TCGA 数据分析实战 —— 差异基因
2021-07-16 本文已影响0人
名本无名
转录组分析
上一节，我们简单介绍了 CNV 数据的处理以及突变数据可视化。下面我们简单介绍一下转录组数据分析中必不可少的差异基因分析，以及通路富集分析
1. 数据准备
我们来分析 LGG 和 GBM 之间的转录组差异，首先从 GDC 中获取原位癌 read count 数据
get_count <- function(cancer) {
query <- GDCquery(
project = cancer,
data.category = "Gene expression",
data.type = "Gene expression quantification",
platform = "Illumina HiSeq",
file.type = "results",
sample.type = c("Primary Tumor"),
legacy = TRUE
)
# 选择 20 个样本
query$results[[1]] <- query$results[[1]][1:20,]
GDCdownload(query)
# 获取 read count
exp.count <- GDCprepare(
query,
summarizedExperiment = TRUE,
)
return(exp.count)
}
gbm.exp <- get_count("TCGA-GBM")
lgg.exp <- get_count("TCGA-LGG")
dataPrep_GBM <- TCGAanalyze_Preprocessing(
object = gbm.exp,
cor.cut = 0.6,
datatype = "raw_count"
)
dataPrep_LGG <- TCGAanalyze_Preprocessing(
object = lgg.exp,
cor.cut = 0.6,
datatype = "raw_count"
)
# 合并数据并使用 gcContent 方法进行标准化
dataNorm <- TCGAanalyze_Normalization(
tabDF = cbind(dataPrep_LGG, dataPrep_GBM),
geneInfo = TCGAbiolinks::geneInfo,
method = "gcContent"
)
# 分位数过滤
dataFilt <- TCGAanalyze_Filtering(
tabDF = dataNorm,
method = "quantile",
qnt.cut = 0.25
)
# 将数据拆分
dataLGG <- subset(dataFilt, select = gbm.exp$barcode)
dataGBM <- subset(dataFilt, select = lgg.exp$barcode)
2. edgeR
edgeR 可以对 RNA-seq、SAGE 或 ChIP-Seq 等数据进行差异表达分析，任何从基因组特征上产生的 read count 都可以分析
该算法既可以用于多组实验的统计分析，也可以使用广义线性模型（glm）方法来对多因子实验数据进行统计分析
不仅可以应用在基因水平，也可以在外显子、转录本水平进行差异分析，我们以基因水平为例
使用 TCGAbiolinks 提供的差异表达分析方法，可以很容易地获取差异基因列表
DEGs.edgeR <- TCGAanalyze_DEA(
mat1 = dataLGG,
mat2 = dataGBM,
Cond1type = "LGG",
Cond2type = "GBM",
fdr.cut = 0.01 ,
logFC.cut = 1,
method = "glmLRT"
)
edgeR 包含许多分析的变体，其中 glm 方法较经典的方法更加灵活，而 glm 方法又包含两种检验方法：似然率检验和准似然率 F 检验，其中准似然率的方法适用于 RNA-seq，似然率检验方法适用于单细胞 RNA-seq 和无生物学或技术重复的数据。
其中 method 参数支持两种方法：
-
glmLRT：似然率检验 -
exactTest：经典方法，两组间配对均值差的精确检验
也可以使用 edgeR 包提供的函数来识别差异基因
library(edgeR)
# 合并 read counts 数据
x <- cbind(dataPrep_LGG, dataPrep_GBM)
group <- ifelse(colnames(x) %in% gbm.exp$barcode, 1, 2)
y <- DGEList(counts = x, group = group)
# 过滤低表达基因
keep <- filterByExpr(y)
y <- y[keep,,keep.lib.sizes=FALSE]
# RNA-seq 数据，推荐使用 TMM 标准化
y <- calcNormFactors(y)
# 构建设计矩阵，由于设置分组
design <- model.matrix(~group)
# 离散评估
y <- estimateDisp(y,design)
# quasi-likelihood F-tests
fit <- glmQLFit(y,design)
qlf <- glmQLFTest(fit,coef=2)
准似然率检验的 top 基因
> topTags(qlf)
Coefficient: group
logFC logCPM F PValue FDR
FTHL3 6.199229 3.276914 368.7956 3.511379e-22 5.467568e-18
NACAP1 5.696704 1.979549 282.2026 4.742793e-20 3.692501e-16
LOC728643 5.803125 1.175056 172.4759 2.527741e-16 1.311982e-12
H3F3C 1.966811 4.100946 160.0586 8.745783e-16 3.404515e-12
RPL17 2.141402 6.987182 146.1642 3.856715e-15 1.201058e-11
TMOD2 1.926996 7.926813 135.9716 1.231274e-14 3.195360e-11
RPL13AP3 5.331669 -1.467946 128.9576 2.848142e-14 6.335488e-11
VKORC1 -1.464038 6.009896 125.9088 4.145942e-14 8.069558e-11
LYPLA1 -1.560197 5.911095 124.9255 4.686569e-14 8.108285e-11
CPEB3 2.157688 3.651139 120.9170 7.784031e-14 1.212051e-10
# likelihood ratio tests
fit <- glmFit(y,design)
lrt <- glmLRT(fit,coef=2)
似然率检验的 top 基因
> topTags(lrt)
Coefficient: group
logFC logCPM LR PValue FDR
FTHL3 6.198600 3.276914 372.0878 6.570886e-83 1.023153e-78
NACAP1 5.696150 1.979549 345.9657 3.203870e-77 2.494373e-73
LOC728643 5.803098 1.175056 273.6598 1.808411e-61 9.386255e-58
RPL13AP3 5.330410 -1.467946 151.4426 8.387839e-35 3.265176e-31
H3F3C 1.966764 4.100946 143.7204 4.090057e-33 1.273726e-29
RPL17 2.141400 6.987182 140.4016 2.174685e-32 5.513899e-29
SLC6A10P 4.275208 -1.038397 140.1416 2.478794e-32 5.513899e-29
TMOD2 1.926998 7.926813 128.9986 6.786503e-30 1.320908e-26
GCSH 3.854903 1.669420 121.2087 3.439727e-28 5.951109e-25
LYPLA1 -1.560212 5.911095 114.7787 8.798646e-27 1.370037e-23
3. limma
limma 也是一个可以对芯片和 RNA-seq 数据进行差异表达分析的包，它的线性模型和差异表达函数可以应用于任何基因表达定量技术，也包括定量 PCR，还可以处理单通道和双通道的数据。
limma 包集成了很多功能，包括数据的读取、预处理（如背景矫正、组内或组件标准化等）和差异表达分析
使用 TCGAanalyze_DEA 函数，只需指定 pipeline = "limma"，便会使用 limma 包中的函数来识别差异表达基因，例如
DEGs.limma <- TCGAanalyze_DEA(
mat1 = dataLGG,
mat2 = dataGBM,
pipeline = "limma",
Cond1type = "LGG",
Cond2type = "GBM",
fdr.cut = 0.01 ,
logFC.cut = 1
)
运行上面的代码会返回一张图
如果使用 limma 包来分析，可以先用 edgeR 的几个函数来过滤低表达的基因，然后进行 TMM 标准化，例如
counts <- cbind(dataPrep_LGG, dataPrep_GBM)
dge <- DGEList(counts=counts)
keep <- filterByExpr(dge, design)
dge <- dge[keep,,keep.lib.sizes=FALSE]
dge <- calcNormFactors(dge)
如果样本之间的测序深度高度一致，则最简单最鲁棒的方法是使用 limma-trend 来进行差异表达分析。该方法需要使用 edgeR 的 cpm 函数先将 count 数转换为 logCPM
logCPM <- cpm(dge, log=TRUE, prior.count=3)
prior.count 用于控制低 count 的对数方差，获取的 logCPM 值可以应用于任何标准的 limma 流程，例如
fit <- lmFit(logCPM, design)
fit <- eBayes(fit, trend=TRUE)
topTable(fit, coef=ncol(design))
> topTable(fit, coef=ncol(design))
logFC AveExpr t P.Value adj.P.Val B
NACAP1|83955 5.019785 0.2889119 26.09854 1.784215e-27 2.795507e-23 50.54761
FTHL3|2498 5.674651 1.1845469 25.41154 5.094355e-27 3.990918e-23 49.64564
LOC728643|728643 4.654882 -0.4465236 22.68152 4.279208e-25 2.234888e-21 45.76199
GCSH|2653 3.754083 0.4104193 16.75698 3.830152e-20 1.500271e-16 35.31615
SLC6A10P|386757 2.962269 -1.6739395 16.15459 1.446361e-19 4.532316e-16 34.06532
RPL13AP3|645683 2.695707 -2.2357252 14.03958 2.075289e-17 5.419270e-14 29.34302
PPIAL4C|653598 2.681801 -1.3227442 13.72591 4.523101e-17 1.012399e-13 28.59619
TMOD2|29767 2.139542 7.4007515 12.49070 1.090567e-15 2.020331e-12 25.53092
HNRNPA3P1|10151 2.239508 0.7594537 12.46728 1.160517e-15 2.020331e-12 25.47083
RPL17|6139 2.173207 6.5201863 12.32842 1.680040e-15 2.632287e-12 25.11311
可以在对基因排序时，给 FC 添加更高的权重
fit <- lmFit(logCPM, design)
fit <- treat(fit, lfc=log2(1.2), trend=TRUE)
> topTreat(fit, coef=ncol(design))
logFC AveExpr t P.Value adj.P.Val
NACAP1|83955 5.019785 0.2889119 24.73099 7.502004e-27 1.175414e-22
FTHL3|2498 5.674651 1.1845469 24.23365 1.676295e-26 1.313210e-22
LOC728643|728643 4.654882 -0.4465236 21.39985 2.035471e-24 1.063058e-20
GCSH|2653 3.754083 0.4104193 15.58288 2.656256e-19 1.040455e-15
SLC6A10P|386757 2.962269 -1.6739395 14.72014 1.993135e-18 6.245687e-15
RPL13AP3|645683 2.695707 -2.2357252 12.66966 3.402740e-16 8.885689e-13
PPIAL4C|653598 2.681801 -1.3227442 12.37965 7.334959e-16 1.641773e-12
RPL23P8|222901 3.021467 0.5510455 11.01340 3.164181e-14 5.657453e-11
HNRNPA3P1|10151 2.239508 0.7594537 11.00297 3.249749e-14 5.657453e-11
TMOD2|29767 2.139542 7.4007515 10.95510 3.723658e-14 5.834227e-11
如果样本间的文库大小差别较大，则 voom 方法理论上相较于 limma-trend 的性能更好。例如，传递标准化和过滤后的对象
v <- voom(dge, design, plot=TRUE)
也可以直接传入未表转化的 count 数据
v <- voom(counts, design, plot=TRUE)
如果背景噪音较大，可以使用芯片间的标准化方法来矫正
v <- voom(counts, design, plot=TRUE, normalize="quantile")
会生成这样一个均值方差图
最后，识别差异表达基因
fit <- lmFit(v, design)
fit <- eBayes(fit)
> topTable(fit, coef=ncol(design))
logFC AveExpr t P.Value adj.P.Val B
NACAP1|83955 5.487935 0.03746277 20.76468 1.329296e-23 2.082741e-19 38.49041
LOC728643|728643 5.591719 -0.93701291 17.67713 5.582431e-21 4.373276e-17 33.18987
GCSH|2653 3.978260 0.25749396 15.95254 2.350618e-19 1.227650e-15 31.48577
FTHL3|2498 5.944472 1.02921281 15.30391 1.036540e-18 4.060127e-15 30.54694
SLC6A10P|386757 4.002446 -2.30901427 14.71063 4.188955e-18 1.312651e-14 27.61085
RPL13AP3|645683 4.460776 -3.31934440 14.30732 1.106812e-17 2.890256e-14 26.63276
HNRNPA3P1|10151 2.350619 0.68521834 12.89222 3.887209e-16 8.700684e-13 25.56282
RPL17|6139 2.189649 6.51891837 12.26764 2.020113e-15 3.165113e-12 24.99887
TMOD2|29767 2.122613 7.40007617 12.26784 2.019028e-15 3.165113e-12 24.99688
PPIAL4C|653598 3.288928 -1.72427506 12.79122 5.057666e-16 9.905439e-13 24.02883
4. DESeq2
DEseq2 是 DEseq 的升级版，能够对 RNA-seq 流程产生的 count 数据进行差异表达分析，也可以对其他芯片类型的数据进行分析，如 ChIP-Seq、HiC、shRNA 等。
该算法的核心是使用负二项广义线性模型来检验基因表达的差异。
简单示例
counts <- cbind(dataLGG, dataGBM)
coldata <- data.frame(
row.names = colnames(counts),
group = factor(ifelse(colnames(counts) %in% gbm.exp$barcode, "gbm", "lgg"))
)
# 构造 DESeqDataSet 数据结构
dds <- DESeqDataSetFromMatrix(
countData = counts,
colData = coldata,
design = ~ group
)
# 差异分析
dds <- DESeq(dds)
#
resultsNames(dds)
# [1] "Intercept" "group_lgg_vs_gbm"
# 获取结果
res <- results(dds, name="group_lgg_vs_gbm")
# 筛选差异基因
DEGs.DESeq2 <- res[res$padj < 0.01 & abs(res$log2FoldChange) >= 1, ]
> DEGs.DESeq2
log2 fold change (MLE): group lgg vs gbm
Wald test p-value: group lgg vs gbm
DataFrame with 3338 rows and 6 columns
baseMean log2FoldChange lfcSE stat pvalue padj
<numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
100134869 36.0292 1.27684 0.346283 3.68726 2.26678e-04 6.52739e-04
A1BG 244.3233 -1.09141 0.315552 -3.45874 5.42710e-04 1.41740e-03
A2BP1 740.5992 2.92719 0.511364 5.72427 1.03877e-08 8.67522e-08
A2LD1 171.9329 -1.69698 0.320965 -5.28710 1.24268e-07 7.90876e-07
AATK 2436.6236 1.90075 0.355487 5.34689 8.94790e-08 5.90934e-07
... ... ... ... ... ... ...
ZSCAN23 193.559 1.37127 0.250557 5.47287 4.42808e-08 3.18100e-07
ZWILCH 632.742 -1.29730 0.194161 -6.68157 2.36401e-11 3.88757e-10
ZWINT 716.560 -1.23933 0.275902 -4.49192 7.05831e-06 2.90102e-05
ZYX 11344.017 -1.28341 0.225298 -5.69648 1.22304e-08 1.00619e-07
psiTPTE22 720.162 2.75603 0.717245 3.84252 1.21780e-04 3.72260e-04
绘制 MA 图
plotMA(res)
5. 差异基因交集
我们使用 UpSetR 包来对这三种方法所识别出的差异基因进行比较
library(UpSetR)
library(RColorBrewer)
g.edgeR <- rownames(DEGs.edgeR)
g.limma <- rownames(DEGs.limma)
g.DESeq2 <- rownames(DEGs.DESeq2)
set_list <- list(
edgeR = g.edgeR,
limma = g.limma,
DESeq2 = g.DESeq2
)
upset(
fromList(set_list),
order.by = "freq",
sets.bar.color = brewer.pal(7, "Set2")[1:3],
matrix.color = brewer.pal(4, "Set1")[2],
main.bar.color = brewer.pal(7, "Set2"),
)
或者使用 VennDiagram 绘制韦恩图
library(VennDiagram)
grid.newpage()
venn_ploy <- venn.diagram(
x = list(
edgeR = g.edgeR,
limma = g.limma,
DESeq2 = g.DESeq2
),
filename = NULL,
fill = brewer.pal(3, "Set2")
)
grid.draw(venn_ploy)
