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 outlier
dataPrep2 <- 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 cutoff
xlim = 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遇见更好的你 41篇原创内容 --> 公众号：桓峰基因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.
本文使用 文章同步助手 同步
网友评论