RNA 4. SCI 文章中基于TCGA 差异表达之 edgeR

01. 前言

在不久的将来，预计新兴的数字基因表达(DGE) 技术在许多功能基因组学应用方面将超越微阵列技术。其中的一个基础数据分析任务，尤其是基因表达研究，涉及确定是否有证据表明对于转录本或外显子在实验中显著不同。

edgeR 是用于检查的 Bioconductor 软件包重复计数数据的差异表达。过度分散的泊松模型用于解释生物学和技术变化性。经验贝叶斯方法用于调节度数跨转录本的过度分散，提高可靠性推理。该方法甚至可以使用最小的复制水平，提供至少一种表型或实验条件被复制。该软件可能有其他应用程序超越测序数据，例如蛋白质组肽计数数据。

本篇继续展示 R 包 edgeR 的差异基因分析流程。类似 limma, DESeq2，edgeR 作为被广泛使用的 R 包，文献中经常能看到它的身影。基于转录组测序获得的定量表达值，识别差异表达变化的基因或其它非编码 RNA 分子，实际上方法还是非常多的。

02. 数据文件读取

数据的读取我们仍然使用的是 TCGA-COAD 的数据集，表达数据的读取以及临床信息分组的获得我们上期已经提过，不过这里还是附上源代码，方便大家对整体分析有个全面的理解，如下：

library(TCGAbiolinks)
query <- GDCquery(project ="TCGA-COAD",data.category = "Transcriptome Profiling",data.type ="Gene Expression Quantification" ,workflow.type="HTSeq - Counts")samplesDown <- getResults(query,cols=c("cases"))  dataSmTP <- TCGAquery_SampleTypes(barcode = samplesDown,typesample = "TP")dataSmNT <- TCGAquery_SampleTypes(barcode = samplesDown,typesample = "NT")group<-as.data.frame(list(Sample=c(dataSmTP,dataSmNT),Group=c(rep("TP",length(dataSmTP)),rep("NT",length(dataSmNT)))))###设置barcodes参数
queryDown <- GDCquery(project = "TCGA-COAD",data.category = "Transcriptome Profiling",data.type = "Gene Expression Quantification", workflow.type = "HTSeq - Counts", barcode = c(dataSmTP, dataSmNT))# 首次下载数据需要执行，这里已经下载过了，不在执行，默认存放位置为当前工作目录下的GDCdata文件夹中。
#GDCdownload(queryDown,method = "api",
#            directory = "GDCdata",#           files.per.chunk = 10)# GDCprepare():Prepare GDC data,准备GDC数据，使其可用于R语言中进行分析
dataPrep1 <- GDCprepare(query = queryDown, save = TRUE, save.filename = "COAD_case.rda")# 函数功能描述：Array Array Intensity correlation (AAIC) and correlation boxplot to define outlierdataPrep2 <- TCGAanalyze_Preprocessing(object = dataPrep1,cor.cut = 0.6,datatype = "HTSeq - Counts")
dim(dataPrep2)

03. 数据预处理

在数据处理这部分，我们使用 edgeR 的数据过滤模式，该文章提到的过滤 low count 数据，例如 CPM 标准化（推荐），然后标准化，以 TMM 标准化为例，如下：

library(edgeR)
#（1）构建表格
dgelist <- DGEList(counts = dataPrep2, group = group$Group)#（2）过滤
keep <- rowSums(cpm(dgelist) > 1 ) >= 2
dgelist <- dgelist[keep, , keep.lib.sizes = FALSE]#（3）标准化
dgelist_norm <- calcNormFactors(dgelist, method = 'TMM')

04. 差异表达分析

首先根据分组信息构建试验设计矩阵，分组信息中一定要是对照组在前，处理组在后，然后）估算基因表达值的离散度，模型拟合，edgeR 提供了2种拟合算法包括：

负二项广义对数线性模型；
拟似然负二项广义对数线性模型。

design <- model.matrix(~group$Group)#install.packages("statmod")
library(statmod)
dge <- estimateDisp(dgelist_norm, design, robust = TRUE)#（2）模型拟合，edgeR 提供了多种拟合算法
#负二项广义对数线性模型
fit1 <- glmFit(dge, design, robust = TRUE)
lrt1 <- topTags(glmLRT(fit1), n = nrow(dgelist$counts))
lrt1<-as.data.frame(lrt1)
head(lrt1)##                     logFC      logCPM        LR        PValue           FDR
## ENSG00000142959 -5.924830  3.20085238 1457.1137 8.159939e-319 2.592984e-314
## ENSG00000163815 -4.223787  2.59792097 1145.7948 3.679815e-251 5.846674e-247
## ENSG00000107611 -5.131288  1.71477711 1122.0643 5.289173e-246 5.602468e-242
## ENSG00000162461 -4.101967  1.51480635 1085.9423 3.752089e-238 2.980753e-234
## ENSG00000163959 -4.295806  3.39558390 1080.7407 5.067873e-237 3.220836e-233
## ENSG00000144410 -6.284258 -0.02616284  916.3497 2.739233e-201 1.450744e-197#拟似然负二项广义对数线性模型
fit2 <- glmQLFit(dge, design, robust = TRUE)
lrt2 <- topTags(glmQLFTest(fit2), n = nrow(dgelist$counts))

05. 绘制火山图

火山图的绘制仍然使用 EnhancedVolcano 包，因为非常好用，只需要修改几个参数，比如输入矩阵的变量，x， y 变量所使用的列名称即可，输出火山图，这里我们使用的是负二项广义对数线性模型获得的差异基因结果，如下：

require(EnhancedVolcano)
EnhancedVolcano(lrt1,lab = rownames(lrt1),labSize = 2,x = "logFC",y = "PValue",#selectLab = rownames(lrt)[1:4],xlab = bquote(~Log[2]~ "fold change"),ylab = bquote(~-Log[10]~italic(P)),pCutoff = 0.001,## pvalue闃堝€?FCcutoff = 1,## FC cutoffxlim = c(-5,5),ylim = c(0,5),colAlpha = 0.6,legendLabels =c("NS","Log2 FC"," p-value"," p-value & Log2 FC"),legendPosition = "top",legendLabSize = 10,legendIconSize = 3.0,pointSize = 1.5,title ="edgeR results",subtitle = 'Differential Expression Genes',
)

06. 筛选差异基因

筛选差异基因，首先对表格排个序，按 FDR 值升序排序，相同 FDR 值下继续按 log2FC 降序排序，然后标记三种情况的基因：

上调的基因：log2FC≥1 & FDR<0.01 标识 up；
下调的基因：log2FC≤-1 & FDR<0.01 标识 down；
非差异的基因：其余标识 none。

lrt1 <- lrt1[order(lrt1$FDR, lrt1$logFC, decreasing = c(FALSE, TRUE)), ]lrt1[which(lrt1$logFC >= 1 & lrt1$FDR < 0.01),'sig'] <- 'up'
lrt1[which(lrt1$logFC <= -1 & lrt1$FDR < 0.01),'sig'] <- 'down'
lrt1[which(abs(lrt1$logFC) <= 1 | lrt1$FDR >= 0.01),'sig'] <- 'none'
#输出选择的差异基因总表
res <- subset(lrt1, sig %in% c('up', 'down'))
head(res)##                     logFC      logCPM        LR        PValue           FDR  sig
## ENSG00000142959 -5.924830  3.20085238 1457.1137 8.159939e-319 2.592984e-314 down
## ENSG00000163815 -4.223787  2.59792097 1145.7948 3.679815e-251 5.846674e-247 down
## ENSG00000107611 -5.131288  1.71477711 1122.0643 5.289173e-246 5.602468e-242 down
## ENSG00000162461 -4.101967  1.51480635 1085.9423 3.752089e-238 2.980753e-234 down
## ENSG00000163959 -4.295806  3.39558390 1080.7407 5.067873e-237 3.220836e-233 down
## ENSG00000144410 -6.284258 -0.02616284  916.3497 2.739233e-201 1.450744e-197 down

差异基因个数，还是蛮多的，估计是正常样本的数据量挺少的，导致偏向性非常高。

table(res$sig)
##
## down   up
## 3308 8164

基于基因表达的三个软件包 limma, DESeq2，edgeR 都已经介绍了，筛选自己的数据用起来吧，后面在讲文章中差异表达的数据在 SCI 文章中怎样展示给出大家，敬请期待！

关注公众号，每日有更新！

桓峰基因

生物信息分析，SCI文章撰写及生物信息基础知识学习：R语言学习，perl基础编程，linux系统命令，Python遇见更好的你

40篇原创内容

公众号

Reference:

Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139-140.
Agarwal A, Koppstein D, Rozowsky J, et al. Comparison and calibration of transcriptome data from RNA-Seq and tiling arrays. BMC Genomics. 2010;11:383.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.
Maurya NS, Kushwaha S, Chawade A, Mani A. Transcriptome profiling by combined machine learning and statistical R analysis identifies TMEM236 as a potential novel diagnostic biomarker for colorectal cancer. Sci Rep. 2021;11(1):14304.