其他
TCGA（转录组）差异分析三大R包及其结果对比
大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~
就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~
这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我！
花花写于2020.1.6
本系列是我的TCGA学习记录，跟着生信技能树B站课程学的，已获得授权（嗯，真的_）。课程链接：https://www.bilibili.com/video/av49363776
目录：
1.准备R包
1if(!require(stringr))install.packages('stringr')
2if(!require(ggplotify))install.packages("ggplotify")
3if(!require(patchwork))install.packages("patchwork")
4if(!require(cowplot))install.packages("cowplot")
5if(!require(DESeq2))BiocManager::install('DESeq2')
6if(!require(edgeR))BiocManager::install('edgeR')
7if(!require(limma))BiocManager::install('limma')
2.准备数据
获取示例数据请在生信技能树公众号聊天窗口回复：”deg“。使用任一counts表达矩阵均可。当然了，你看完B站视频，自己学会下载TCGA-KIRC的miRNA表达矩阵就最好了。
本示例的数据是TCGA-KIRC的miRNA表达矩阵。tcga样本编号14-15位是隐藏分组信息的,详见：
TCGA的样本id里藏着分组信息
这里需要注意的是miRNA也是测序拿到的表达矩阵,所以分析等同于RNA-seq得到表达矩阵，一定要跟芯片数据分析区分开来哦！！！
1rm(list = ls())
2load("tcga_kirc_exp.Rdata") #表达矩阵expr
3dim(expr)
1## [1] 552 593
1group_list <- ifelse(as.numeric(str_sub(colnames(expr),14,15))<10,"tumor","normal")
2group_list <- factor(group_list,levels = c("normal","tumor"))
3table(group_list)
1## group_list
2## normal tumor
3## 71 522
3.三大R包的差异分析
3.1 Deseq2
1library(DESeq2)
2colData <- data.frame(row.names =colnames(expr),
3 condition=group_list)
4dds <- DESeqDataSetFromMatrix(
5 countData = expr,
6 colData = colData,
7 design = ~ condition)
8#参考因子应该是对照组 dds$condition <- relevel(dds$condition, ref = "untrt")
9
10dds <- DESeq(dds)
11# 两两比较
12res <- results(dds, contrast = c("condition",rev(levels(group_list))))
13resOrdered <- res[order(res$pvalue),] # 按照P值排序
14DEG <- as.data.frame(resOrdered)
15head(DEG)
16# 去除NA值
17DEG <- na.omit(DEG)
18
19#添加change列标记基因上调下调
20#logFC_cutoff <- with(DEG,mean(abs(log2FoldChange)) + 2*sd(abs(log2FoldChange)) )
21logFC_cutoff <- 1
22DEG$change = as.factor(
23 ifelse(DEG$pvalue < 0.05 & abs(DEG$log2FoldChange) > logFC_cutoff,
24 ifelse(DEG$log2FoldChange > logFC_cutoff ,'UP','DOWN'),'NOT')
25)
26head(DEG)
27
28DESeq2_DEG <- DEG
3.2 edgeR
1library(edgeR)
2
3dge <- DGEList(counts=expr,group=group_list)
4dge$samples$lib.size <- colSums(dge$counts)
5dge <- calcNormFactors(dge)
6
7design <- model.matrix(~0+group_list)
8rownames(design)<-colnames(dge)
9colnames(design)<-levels(group_list)
10
11dge <- estimateGLMCommonDisp(dge,design)
12dge <- estimateGLMTrendedDisp(dge, design)
13dge <- estimateGLMTagwiseDisp(dge, design)
14
15fit <- glmFit(dge, design)
16fit2 <- glmLRT(fit, contrast=c(-1,1))
17
18DEG=topTags(fit2, n=nrow(expr))
19DEG=as.data.frame(DEG)
20logFC_cutoff <- with(DEG,mean(abs(logFC)) + 2*sd(abs(logFC)) )
21logFC_cutoff <- 1
22DEG$change = as.factor(
23 ifelse(DEG$PValue < 0.05 & abs(DEG$logFC) > logFC_cutoff,
24 ifelse(DEG$logFC > logFC_cutoff ,'UP','DOWN'),'NOT')
25)
26head(DEG)
1## logFC logCPM LR PValue FDR
2## hsa-mir-508 -4.264945 5.3610815 825.7952 1.329948e-181 7.341313e-179
3## hsa-mir-514-3 -4.262325 3.5005425 674.3829 1.112883e-148 3.071556e-146
4## hsa-mir-514-2 -4.258203 3.4771070 658.6855 2.885406e-145 5.309148e-143
5## hsa-mir-506 -5.522829 0.7477531 654.6812 2.143124e-144 2.957511e-142
6## hsa-mir-514-1 -4.271951 3.4852217 642.0128 1.219493e-141 1.346320e-139
7## hsa-mir-514b -5.956182 0.3742949 579.5893 4.606971e-128 4.238413e-126
8## change
9## hsa-mir-508 DOWN
10## hsa-mir-514-3 DOWN
11## hsa-mir-514-2 DOWN
12## hsa-mir-506 DOWN
13## hsa-mir-514-1 DOWN
14## hsa-mir-514b DOWN
1table(DEG$change)
1##
2## DOWN NOT UP
3## 64 368 120
1edgeR_DEG <- DEG
3.limma-voom
1library(limma)
2
3design <- model.matrix(~0+group_list)
4colnames(design)=levels(group_list)
5rownames(design)=colnames(expr)
6
7dge <- DGEList(counts=expr)
8dge <- calcNormFactors(dge)
9logCPM <- cpm(dge, log=TRUE, prior.count=3)
10
11v <- voom(dge,design, normalize="quantile")
12fit <- lmFit(v, design)
13
14constrasts = paste(rev(levels(group_list)),collapse = "-")
15cont.matrix <- makeContrasts(contrasts=constrasts,levels = design)
16fit2=contrasts.fit(fit,cont.matrix)
17fit2=eBayes(fit2)
18
19DEG = topTable(fit2, coef=constrasts, n=Inf)
20DEG = na.omit(DEG)
21#logFC_cutoff <- with(DEG,mean(abs(logFC)) + 2*sd(abs(logFC)) )
22logFC_cutoff <- 1
23DEG$change = as.factor(
24 ifelse(DEG$P.Value < 0.05 & abs(DEG$logFC) > logFC_cutoff,
25 ifelse(DEG$logFC > logFC_cutoff ,'UP','DOWN'),'NOT')
26)
27head(DEG)
1## logFC AveExpr t P.Value adj.P.Val
2## hsa-mir-141 -5.492612 2.990323 -31.22459 1.280624e-127 7.069043e-125
3## hsa-mir-200c -5.333666 5.687063 -30.87441 8.224995e-126 2.270099e-123
4## hsa-mir-3613 2.199074 3.862900 23.32209 5.152888e-86 9.481314e-84
5## hsa-mir-15a 1.335460 7.014647 22.83389 2.008313e-83 2.771472e-81
6## hsa-mir-934 -3.234590 -1.930201 -22.39709 4.148098e-81 4.579500e-79
7## hsa-mir-122 5.554068 3.112250 22.33183 9.192713e-81 8.457296e-79
8## B change
9## hsa-mir-141 280.9396 DOWN
10## hsa-mir-200c 276.7881 DOWN
11## hsa-mir-3613 185.4023 UP
12## hsa-mir-15a 179.4222 UP
13## hsa-mir-934 174.1455 DOWN
14## hsa-mir-122 173.3486 UP
1limma_voom_DEG <- DEG
2#save(DESeq2_DEG,edgeR_DEG,limma_voom_DEG,group_list,file = "DEG.Rdata")
4.差异分析结果的可视化
1rm(list = ls())
2load("tcga_kirc_exp.Rdata")
3load("DEG.Rdata")
4source("3-plotfunction.R")
5logFC_cutoff <- 1
6expr[1:4,1:4]
1## TCGA-A3-3307-01A-01T-0860-13 TCGA-A3-3308-01A-02R-1324-13
2## hsa-let-7a-1 5056 14503
3## hsa-let-7a-2 10323 29238
4## hsa-let-7a-3 5429 14738
5## hsa-let-7b 17908 37062
6## TCGA-A3-3311-01A-02R-1324-13 TCGA-A3-3313-01A-02R-1324-13
7## hsa-let-7a-1 8147 7138
8## hsa-let-7a-2 16325 14356
9## hsa-let-7a-3 8249 7002
10## hsa-let-7b 28984 6909
1dat = log(expr+1)
2pca.plot = draw_pca(dat,group_list)
3
4cg1 = rownames(DESeq2_DEG)[DESeq2_DEG$change !="NOT"]
5cg2 = rownames(edgeR_DEG)[edgeR_DEG$change !="NOT"]
6cg3 = rownames(limma_voom_DEG)[limma_voom_DEG$change !="NOT"]
7
8h1 = draw_heatmap(expr[cg1,],group_list)
9h2 = draw_heatmap(expr[cg2,],group_list)
10h3 = draw_heatmap(expr[cg3,],group_list)
1v1 = draw_volcano(test = DESeq2_DEG[,c(2,5,7)],pkg = 1)
2v2 = draw_volcano(test = edgeR_DEG[,c(1,4,6)],pkg = 2)
3v3 = draw_volcano(test = limma_voom_DEG[,c(1,4,7)],pkg = 3)
4
5library(patchwork)
6(h1 + h2 + h3) / (v1 + v2 + v3) +plot_layout(guides = 'collect')
1#(v1 + v2 + v3) +plot_layout(guides = 'collect')
2ggsave("heat_volcano.png",width = 21,height = 9)
5.三大R包差异基因对比
1# 三大R包差异基因交集
2UP=function(df){
3 rownames(df)[df$change=="UP"]
4}
5DOWN=function(df){
6 rownames(df)[df$change=="DOWN"]
7}
8
9up = intersect(intersect(UP(DESeq2_DEG),UP(edgeR_DEG)),UP(limma_voom_DEG))
10down = intersect(intersect(DOWN(DESeq2_DEG),DOWN(edgeR_DEG)),DOWN(limma_voom_DEG))
11
12hp = draw_heatmap(expr[c(up,down),],group_list)
13
14#上调、下调基因分别画维恩图
15
16up.plot <- venn(UP(DESeq2_DEG),UP(edgeR_DEG),UP(limma_voom_DEG),
17 "UPgene"
18)
1down.plot <- venn(DOWN(DESeq2_DEG),DOWN(edgeR_DEG),DOWN(limma_voom_DEG),
2 "DOWNgene"
3)
4#维恩图拼图，终于搞定
5
6library(cowplot)
7library(ggplotify)
8up.plot = as.ggplot(as_grob(up.plot))
9down.plot = as.ggplot(as_grob(down.plot))
10library(patchwork)
11#up.plot + down.plot
1# 就爱玩拼图
2pca.plot + hp+up.plot +down.plot
1ggsave("deg.png",height = 10,width = 10)
三个R包差异分析结果的交集共有50个上调和51个下调，可以作为最终结果提交。当然，这三个包没有对错之分，你拿其中任意一个包的分析结果都是对的。取交集的方法更可靠，但不是必须的，有些数据取交集会很可怜的。
插个小广告！
再给生信技能树打个call！