# MESSAGE -----------------------------------------------------------------

#

# author: Yulin Lyu

# email: lvyulin@pku.edu.cn

#

# require: R whatever

#

# ---

# 1. Load packages --------------------------------------------------------

setwd("exampleData/RNA")

# grammar

library(tidyverse)

library(magrittr)

library(glue)

library(data.table)

# analysis

library(DESeq2)

library(irlba)

# graphics

library(ggplot2)

library(ggrepel)

library(ggsci)

library(scales)

# 2. Load data ------------------------------------------------------------

vsd <- readRDS("mid/vsd.rds")

vsdMtx <- assay(vsd)

diffData <- fread("mid/ES_vs_Fib.DEG.csv")

colnames(diffData)[1] <- "gene"

# 3. Plot -----------------------------------------------------------------

diffData[, type := "ns"][]

diffData[is.na(padj), padj := 1][]

diffData[log2FoldChange > 3 & padj < 0.05, type := "up"][log2FoldChange < -3 & padj < 0.05, type := "down"][]

usedGene <- diffData[type %in% c("up", "down"), gene]

PCAdata <- prcomp_irlba(t(vsdMtx[usedGene, ]), n = 3, scale. = T)

s <- summary(PCAdata)

s <- s$importance[2, ] %>% round(4)

plotData <- as.data.table(PCAdata$x)

plotData[, id := colnames(vsdMtx)][, type := rep(c("Fib", "CiPS", "ES"), each = 2)][]

usedCol <- pal_npg()(10)[c(1, 4, 3)] %>% set_names(c("Fib", "CiPS", "ES"))

ggplot(plotData, aes(x = PC1, y = PC2)) +

  geom_point(aes(color = type), size = 3, shape = 16, show.legend = F) +

  geom_text_repel(aes(label = id)) +

  scale_color_manual(values = usedCol) +

  labs(

    x = str_c("PC1 (", s[1] * 100, "%)"),

    y = str_c("PC2 (", s[2] * 100, "%)")) +

  coord_fixed() +

  theme(

    panel.background = element_blank(),

    panel.grid = element_blank(),

    panel.border = element_rect(fill = NA)

  )

ggsave("plot/PCA.png", width = 6, height = 2)