轉錄組分析

上一節(jié)，我們簡單介紹了 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 都可以分析

該算法既可以用于多組實驗的統(tǒng)計分析蛔翅，也可以使用廣義線性模型（glm）方法來對多因子實驗數據進行統(tǒng)計分析

不僅可以應用在基因水平敲茄，也可以在外顯子、轉錄本水平進行差異分析山析，我們以基因水平為例

使用 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)