第二步：差异分析
前期文件
rm(list=ls())load("gdc.Rdata")dim(expr)dim(clinical)table(group_list)
如前所述，这里使用老师的Rdata文件，存在tumor与normal的分组（但是矩阵与病人信息数量也不相同好像）
gdc.Rdata
预分析--PCA
观察组间差异是否明显
library("FactoMineR")#画主成分分析图需要加载这两个包library("factoextra")expr[1:4,1:4]dat = log(expr+1) #防止为0t(dat)[1:4,1:4] #行名为样本，列名为基因
转置表达矩阵
dat.pca <- PCA(t(dat), graph = FALSE)pca.plot=fviz_pca_ind(dat.pca, geom.ind = "point", col.ind = group_list, addEllipses = TRUE, legend.title = "Groups")pca.plot
如下图，可以看出两组差异还是比较明显的。
主成分分析图
1、三种差异分析方法
(1) Deseq2
library(DESeq2)colData <- data.frame(row.names =colnames(expr), condition=group_list)dds <- DESeqDataSetFromMatrix( countData = expr, colData = colData, design = ~ condition)dds <- DESeq(dds)# 分组比较res <- results(dds, contrast = c("condition",rev(levels(group_list)))) #rev设置下比较方式resOrdered <- res[order(res$pvalue),] # 按照P值排序DEG <- as.data.frame(resOrdered)# 去除NA值DEG <- na.omit(DEG)head(DEG)
head(DEG)
#添加change列标记基因上调下调logFC_cutoff <- with(DEG,mean(abs(log2FoldChange)) + 2*sd(abs(log2FoldChange)) ) #计算上下调标准#logFC_cutoff <- 2DEG$change = as.factor( ifelse(DEG$pvalue < 0.05 & abs(DEG$log2FoldChange) > logFC_cutoff, ifelse(DEG$log2FoldChange > logFC_cutoff ,'UP','DOWN'),'NOT'))head(DEG)table(DEG$change)DESeq2_DEG <- DEG
基因分为下调、无显著差异、上调三类
之前学习RNA-seq时，用的就是此方法。此外还有edgeR，limma两种方法，不太熟悉，仅记录下代码吧。会得到相似的三分类结果，但存在量的差异。
(2) edgeR
library(edgeR)dge <- DGEList(counts=expr,group=group_list)dge$samples$lib.size <- colSums(dge$counts)dge <- calcNormFactors(dge) design <- model.matrix(~0+group_list)rownames(design)<-colnames(dge)colnames(design)<-levels(group_list)dge <- estimateGLMCommonDisp(dge,design)dge <- estimateGLMTrendedDisp(dge, design)dge <- estimateGLMTagwiseDisp(dge, design)fit <- glmFit(dge, design)fit2 <- glmLRT(fit, contrast=c(-1,1)) DEG=topTags(fit2, n=nrow(expr))DEG=as.data.frame(DEG)logFC_cutoff <- with(DEG,mean(abs(logFC)) + 2*sd(abs(logFC)) )#logFC_cutoff <- 2DEG$change = as.factor( ifelse(DEG$PValue < 0.05 & abs(DEG$logFC) > logFC_cutoff, ifelse(DEG$logFC > logFC_cutoff ,'UP','DOWN'),'NOT'))head(DEG)table(DEG$change)edgeR_DEG <- DEG
edgeR三分类
(3) limma
library(limma)design <- model.matrix(~0+group_list)colnames(design)=levels(group_list)rownames(design)=colnames(expr)dge <- DGEList(counts=expr)dge <- calcNormFactors(dge)logCPM <- cpm(dge, log=TRUE, prior.count=3)v <- voom(dge,design, normalize="quantile")fit <- lmFit(v, design)constrasts = paste(rev(levels(group_list)),collapse = "-")cont.matrix <- makeContrasts(contrasts=constrasts,levels = design) fit2=contrasts.fit(fit,cont.matrix)fit2=eBayes(fit2)DEG = topTable(fit2, coef=constrasts, n=Inf)DEG = na.omit(DEG)logFC_cutoff <- with(DEG,mean(abs(logFC)) + 2*sd(abs(logFC)) )#logFC_cutoff <- 2DEG$change = as.factor( ifelse(DEG$P.Value < 0.05 & abs(DEG$logFC) > logFC_cutoff, ifelse(DEG$logFC > logFC_cutoff ,'UP','DOWN'),'NOT'))head(DEG)table(DEG$change)limma_voom_DEG <- DEG
limma三分类
tj = data.frame(deseq2 = as.integer(table(DESeq2_DEG$change)), edgeR = as.integer(table(edgeR_DEG$change)), limma_voom = as.integer(table(limma_voom_DEG$change)), row.names = c("down","not","up") )tjsave(DESeq2_DEG,edgeR_DEG,limma_voom_DEG,group_list,tj,file = "DEG.Rdata")rm(list=ls())
tj
2、差异分析结果可视化
2.1 热图
挑选出显著差异基因
xz_DESeq2 = rownames(DESeq2_DEG)[DESeq2_DEG$change !="NOT"]xz_DESeq2= dat[xz_DESeq2,]xz_edgeR = rownames(edgeR_DEG)[DESeq2_DEG$change !="NOT"]xz_edgeR= dat[xz_edgeR,]xz_limma = rownames(limma_voom_DEG)[DESeq2_DEG$change !="NOT"]xz_limma= dat[xz_limma,]
绘制热图
library(pheatmap)library(ggplot2)n=xz_DESeq2 # n=xz_edgeR # n=xz_limman[1:4,1:4]n = t(scale(t(n)))n[1:4,1:4]n_cutoff=2 #把绝对值对于2的count值改为2或者-2n[n > n_cutoff] = n_cutoffn[n < -n_cutoff] = -n_cutoffn[1:4,1:4]annotation_col = data.frame(group = group_list)rownames(annotation_col) = colnames(n)p.DEseq2 = as.ggplot(pheatmap(n, show_colnames = F, show_rownames = F, annotation_col = annotation_col, legend = T, annotation_legend = T, annotation_names_col = T))
p.DEseq2
2.2 火山图
之前在差异分析是已经增加了down，not、up三分类变化的change列，这里我们直接使用其来绘制
m2d = function(x){ mean(abs(x))+2*sd(abs(x))} dat=DESeq2_DEG# dat=edgeR_DEG# dat=limma_voom_DEGlogFC_cutoff=m2d(dat$log2FoldChange)tile <- paste0("down gene:", nrow(dat[dat$change == "DOWN",]), "\n" , " up gene:" , nrow(dat[dat$change == "UP",]) )ggplot(data = dat, aes(x = log2FoldChange, y = -log10(pvalue))) + geom_point(alpha = 0.4, size = 3.5, aes(color = change)) + scale_color_manual(values = c("blue", "grey", "red")) + geom_vline(xintercept = c(-logFC_cutoff, logFC_cutoff), lty = 4, col = "black", lwd = 0.8) + geom_hline(yintercept = -log10(0.05), lty = 4,col = "black", lwd = 0.8) + theme_bw() + labs(title = tile, x = "logFC", y = "-log10(P.value)") + theme(plot.title = element_text(hjust = 0.5))
DESeq2_DEG火山图
可以利用patch包将上面的六张图拼在一起，类似如下代码
library(patchwork)(h1 + h2 + h3) / (v1 + v2 + v3) +plot_layout(guides = 'collect')
3、三种差异分析方法比较
3.1 韦恩图展示三方法的三分类差异
#定义取类函数UP=function(df){ rownames(df)[df$change=="UP"]}DOWN=function(df){ rownames(df)[df$change=="DOWN"]}#画韦恩图ven.plot=venn.diagram(list(Deseq2 = UP(DESeq2_DEG), edgeR = UP(edgeR_DEG), limma = UP(limma_voom_DEG)), imagetype = "png", filename = NULL, lwd = 1, lty = 1, col = c("#0099CC", "#FF6666", "#FFCC99"), fill = c("#0099CC", "#FF6666","#FFCC99"), cat.col = c("#0099CC", "#FF6666", "#FFCC99"), cat.cex = 1.5, rotation.degree = 0, main = "name", main.cex = 1.5, cex = 1.5, alpha = 0.5, reverse = TRUE)#结果返回是一个list，看不太懂，反正不是图..看了下老师教程，加了如下代码就可以了。#猜测是不同绘图系统导致的吧？library(cowplot)p.up = as.ggplot(as_grob(ven.plot))# 同法也可绘制出 p.down
三种方法韦恩图，316个相同
3.2 用三方法相同差异基因绘制热图
#挑选出共同基因up = intersect(intersect(UP(DESeq2_DEG),UP(edgeR_DEG)),UP(limma_voom_DEG))down = intersect(intersect(DOWN(DESeq2_DEG),DOWN(edgeR_DEG)),DOWN(limma_voom_DEG))#套入绘图流程n=expr[c(up,down),] # n=xz_edgeR # n=xz_limman[1:4,1:4]n = t(scale(t(n)))n[1:4,1:4]n_cutoff=2 #把绝对值对于2的count值改为2或者-2n[n > n_cutoff] = n_cutoffn[n < -n_cutoff] = -n_cutoffn[1:4,1:4]annotation_col = data.frame(group = group_list)rownames(annotation_col) = colnames(n)hp = as.ggplot(pheatmap(n, show_colnames = F, show_rownames = F, annotation_col = annotation_col, legend = T, annotation_legend = T, annotation_names_col = T))#如下图，差异很明显
共有基因热图
3.3 最后拼一张综合图收尾吧
library(patchwork)pca.plot + hp + up.plot + down.plot
PCA+heatmap+veen
