PCA plot figure 1D

bulkRNAseq/Fig1D.R

+library(RColorBrewer)
+library(DESeq2)
+library(pheatmap)
+library(ggplot2)
+library(gplots)
+library(ggrepel)
+library(dbplyr)
+
+# Metadata definition
+cellcompartment <- c("HSC", "HSC", "MPP", "MPP", "MEP", "MEP", "CMP", "CMP", "GMP", "GMP",
+                   "HSC", "HSC", "MPP", "MPP", "MEP", "MEP", "CMP", "CMP", "GMP", "GMP")
+condition <- c("THAL", "THAL", "THAL", "THAL", "THAL","THAL","THAL","THAL","THAL","THAL",
+               "HD","HD","HD","HD","HD","HD","HD","HD","HD","HD")
+expgroup <- c("HSC_THAL", "HSC_THAL", "MPP_THAL", "MPP_THAL", "MEP_THAL", "MEP_THAL", "CMP_THAL", "CMP_THAL", "GMP_THAL", "GMP_THAL",
+              "HSC_HD", "HSC_HD", "MPP_HD", "MPP_HD", "MEP_HD", "MEP_HD", "CMP_HD", "CMP_HD", "GMP_HD", "GMP_HD")
+
+myCompartment <- data.frame(cellcompartment, condition, expgroup)
+
+myCompartment$expgroup <- as.factor(myCompartment$expgroup)
+
+myDataframe <- read.table(file = "Fig1D_input.txt", header = TRUE, stringsAsFactors = FALSE, row.names = 1)
+myMatrix <- myDataframe[ , 6:ncol(myDataframe)]
+colnames(myMatrix) <- c(paste0("S", 1:20))
+myMatrix <- as.matrix(myMatrix)
+
+#make coldesign
+myGroup <- paste(myCompartment$cellcompartment, myCompartment$condition, sep = "_")
+myCompartment$expgroup <- myGroup
+
+#DESeq2 object
+dds <- DESeqDataSetFromMatrix(countData = myMatrix,
+                              colData = myCompartment,
+                              design = ~ expgroup)
+#pre-filtering 
+dds <- dds[rowSums(counts(dds))>0, ]
+
+#normalization and dispersion estimation
+dds <- DESeq(object = dds, parallel = T)
+
+#r-log transformation for visualization and performing
+rld <- rlogTransformation(dds)
+
+#PCA analysis
+
+pcagenes <- 500
+#dataPCA <- plotPCA(rld, intgroup = c("cellcompartment", "condition"), returnData=TRUE, 
+#                   ntop=pcagenes)
+dataPCA <- plotPCA(rld, intgroup = "expgroup", returnData=TRUE, ntop=pcagenes)
+percentVar <- round(100 * attr(dataPCA, "percentVar"))
+png(filename = "Fig1D.png",width = 6, height = 5, units = "in", res = 300)
+ggplot(dataPCA, aes(PC1, PC2, fill= group, label = name)) + 
+  geom_point(size=2.5, shape = 21, colour= "black") +
+          scale_fill_brewer(palette="Paired") +
+          # scale_colour_brewer(palette="Paired") + 
+          geom_text_repel() +
+          xlab(label = paste0("PC1: ",percentVar[1],"% variance")) +
+          ylab(label = paste0("PC2: ",percentVar[2],"% variance")) +
+          ggtitle(label = "Principal Component Analysis", paste("Top 500 most variable genes - DESeq2")) +
+          theme_bw() +
+          theme(plot.title = element_text(size = 12, face = "bold", hjust = 0.5), 
+                plot.subtitle=element_text(size=8, hjust=0.5, face="italic", color="black"),
+                axis.title = element_text(size=12, hjust=0.5, face="bold", color="black"),
+                legend.text = element_text(size=10, color="black"),
+                legend.title = element_blank(),
+                axis.text = element_text(size=12, hjust=0.5, color="black"))
+dev.off()
+