TCGA 数据分析实战 —— 差异基因
转录组分析
上一节，我们简单介绍了 CNV 数据的处理以及突变数据可视化。下面我们简单介绍一下转录组数据分析中必不可少的差异基因分析，以及通路富集分析
1. 数据准备
我们来分析 LGG 和 GBM 之间的转录组差异，首先从 GDC 中获取原位癌 read count 数据
get_count <- function(cancer) {
query <- GDCquery(
project = cancer,
data.category = "Transcriptome Profiling",
data.type = "Gene Expression Quantification",
workflow.type = "STAR - Counts",
sample.type = c("Primary Tumor"),
)
# 选择 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 = "unstranded"
)
dataPrep_LGG <- TCGAanalyze_Preprocessing(
object = lgg.exp,
cor.cut = 0.6,
datatype = "unstranded"
)
# 合并数据并使用 gcContent 方法进行标准化
dataNorm <- TCGAanalyze_Normalization(
tabDF = cbind(dataPrep_LGG, dataPrep_GBM),
geneInfo = TCGAbiolinks::geneInfoHT,
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
# ENSG00000243498.2 -3.846834 -1.7460275 365.4453 2.779593e-29 5.976958e-25
# ENSG00000188219.14 -5.786513 -1.0088548 340.2569 3.747256e-28 4.028863e-24
# ENSG00000241511.1 -3.895737 -0.9897558 420.2361 1.873504e-27 1.342866e-23
# ENSG00000226540.2 -4.435408 -1.6694688 281.4515 5.171097e-27 2.779852e-23
# ENSG00000232801.1 -4.296052 -2.0512485 253.6942 1.489829e-25 6.407158e-22
# ENSG00000222038.4 -4.005681 -1.9434298 248.9198 1.955087e-25 7.006706e-22
# ENSG00000249210.1 -4.319108 -0.5711147 332.5780 6.650366e-25 2.042897e-21
# ENSG00000204434.6 -4.158661 -0.9833792 277.2299 4.878320e-24 1.311231e-20
# ENSG00000268034.1 -4.041803 -1.0947910 275.5744 5.967367e-24 1.425737e-20
# ENSG00000244485.1 -4.076731 -2.0762422 211.3509 1.614016e-23 3.470619e-20
使用似然率检验
# likelihood ratio tests
fit <- glmFit(y,design)
lrt <- glmLRT(fit,coef=2)
似然率检验的 top 基因
topTags(lrt)
# Coefficient: group
# logFC logCPM LR PValue FDR
# ENSG00000230897.1 -4.798073 1.5590102 312.5829 5.978120e-70 1.285475e-65
# ENSG00000248180.1 -4.663868 0.2823858 305.5012 2.085817e-68 2.242566e-64
# ENSG00000188219.14 -5.791169 -1.0088548 298.3306 7.611857e-67 5.455926e-63
# ENSG00000218582.2 -3.866119 0.7310565 295.3294 3.430518e-66 1.844161e-62
# ENSG00000234664.1 -4.589833 2.2933366 284.0840 9.675023e-64 4.160840e-60
# ENSG00000235587.2 -3.886623 1.9078041 276.8738 3.604657e-62 1.291849e-58
# ENSG00000249210.1 -4.322112 -0.5711147 272.5557 3.147089e-61 9.667407e-58
# ENSG00000248626.1 -3.982143 0.2286603 269.7605 1.279710e-60 3.439699e-57
# ENSG00000242299.1 -3.255345 3.5054320 268.1920 2.811643e-60 6.717639e-57
# ENSG00000220472.1 -4.529480 0.8686615 267.8744 3.297548e-60 7.090717e-57
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
# ENSG00000248180.1 -4.017196 -0.8532763 -21.33534 3.512245e-24 4.271878e-20 44.81005
# ENSG00000230897.1 -4.461000 0.1622144 -21.26095 4.018806e-24 4.271878e-20 44.68001
# ENSG00000188219.14 -3.595875 -1.9644400 -20.94124 7.203313e-24 4.271878e-20 44.11627
# ENSG00000231848.1 -2.661537 -2.6473631 -20.88783 7.946571e-24 4.271878e-20 44.02133
# ENSG00000220472.1 -4.071414 -0.3489775 -19.97565 4.393957e-23 1.639845e-19 42.36440
# ENSG00000241511.1 -2.937408 -1.6392051 -19.94236 4.682510e-23 1.639845e-19 42.30265
# ENSG00000234664.1 -4.469849 0.8903101 -19.87390 5.338284e-23 1.639845e-19 42.17537
# ENSG00000226540.2 -2.835009 -2.2250009 -19.74411 6.851226e-23 1.680581e-19 41.93296
# ENSG00000249210.1 -3.384945 -1.4213241 -19.72748 7.074458e-23 1.680581e-19 41.90180
# ENSG00000248626.1 -3.505510 -0.6927554 -19.65368 8.158554e-23 1.680581e-19 41.76322
可以在对基因排序时，给 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
# ENSG00000230897.1 -4.461000 0.1622144 -20.00734 2.089884e-23 2.552641e-19
# ENSG00000248180.1 -4.017196 -0.8532763 -19.93836 2.374219e-23 2.552641e-19
# ENSG00000188219.14 -3.595875 -1.9644400 -19.40941 6.583007e-23 4.718480e-19
# ENSG00000231848.1 -2.661537 -2.6473631 -18.82353 2.091660e-22 9.965311e-19
# ENSG00000234664.1 -4.469849 0.8903101 -18.70440 2.685454e-22 9.965311e-19
# ENSG00000220472.1 -4.071414 -0.3489775 -18.68512 2.780629e-22 9.965311e-19
# ENSG00000249210.1 -3.384945 -1.4213241 -18.19452 7.514892e-22 1.618909e-18
# ENSG00000248626.1 -3.505510 -0.6927554 -18.17898 7.763922e-22 1.618909e-18
# ENSG00000218582.2 -3.559710 -0.2117952 -18.17638 7.807892e-22 1.618909e-18
# ENSG00000241511.1 -2.937408 -1.6392051 -18.15659 8.112117e-22 1.618909e-18
如果样本间的文库大小差别较大，则 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
# ENSG00000235587.2 -3.851652 0.793595296 -20.22831 4.045452e-22 2.453971e-17 40.02019
# ENSG00000248180.1 -4.539996 -1.199882809 -19.67450 1.074928e-21 3.260257e-17 38.81154
# ENSG00000218582.2 -3.695739 -0.361025227 -18.65475 6.905945e-21 1.217588e-16 37.18162
# ENSG00000220472.1 -4.343932 -0.544061553 -18.50139 9.200248e-21 1.217588e-16 36.88809
# ENSG00000248626.1 -3.735037 -0.907388005 -18.45509 1.003617e-20 1.217588e-16 36.75868
# ENSG00000242299.1 -3.237173 2.670786286 -18.10944 1.931877e-20 1.819198e-16 36.47328
# ENSG00000234664.1 -4.595866 0.806614940 -18.06594 2.099305e-20 1.819198e-16 36.24159
# ENSG00000237350.1 -3.408635 -0.377114901 -17.61462 5.020317e-20 3.806655e-16 35.29454
# ENSG00000231848.1 -4.524592 -4.000702094 -17.29721 9.367443e-20 5.479215e-16 34.79879
# ENSG00000230897.1 -4.638526 -0.008763245 -17.32766 8.819939e-20 5.479215e-16 34.79152
4. DESeq2
DEseq2 是 DEseq 的升级版，能够对 RNA-seq 流程产生的 count 数据进行差异表达分析，也可以对其他芯片类型的数据进行分析，如 ChIP-Seq、HiC、shRNA 等。
该算法的核心是使用负二项广义线性模型来检验基因表达的差异。
简单示例
library(DESeq2)
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 <- dplyr::filter(as.data.frame(res),
!is.na(padj) & padj < 0.01 & abs(log2FoldChange) >= 1)
差异基因
DEGs.DESeq2[1:3,]
# baseMean log2FoldChange lfcSE stat pvalue padj
# ENSG00000000971 16034.7189 -1.282929 0.3916925 -3.275347 1.055324e-03 4.164254e-03
# ENSG00000001036 1671.8820 -1.273786 0.2598955 -4.901147 9.527894e-07 9.330880e-06
# ENSG00000001617 109.4029 -1.594253 0.3100178 -5.142458 2.711670e-07 3.067991e-06
绘制 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)
扫一扫
2839
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
