CUT&Tag 数据处理与分析教程 七：差异分析

本章节将介绍 CUT&Tag 数据的差异分析。

• DESeq2: 2014-Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2

估计高通量测序分析计数数据的方差 - 均值依赖性，并基于使用负二项分布的模型对差异表达进行检验。

创建 peak x 样本矩阵

通常，差异分析比较相同组蛋白修饰的两种或两种以上的情况。在本教程中，由于演示数据的限制，我们将通过比较 H3K27me3 的两个重复和 H3K4me3 的两个重复来说明差异分析。我们将使用 DESeq2(完整教程 2020-Analyzing RNA-seq data with DESeq2^[1]) 作为说明。

创建一个最终 Peak 列表，合并每个示例 call 的所有 Peak 为一个

##=== R 命令 ===## 
mPeak = GRanges()
## 合并所有 Peak
for(hist in histL){
  for(rep in repL){
    peakRes = read.table(paste0(projPath, "/peakCalling/SEACR/", hist, "_", rep, "_seacr_control.peaks.stringent.bed"), header = FALSE, fill = TRUE)
    mPeak = GRanges(seqnames = peakRes$V1, IRanges(start = peakRes$V2, end = peakRes$V3), strand = "*") %>% append(mPeak, .)
  }
}
masterPeak = reduce(mPeak)

获取最终 Peak 文件中每一个 Peak 中片段数

##=== R 命令 ===## 
library(DESeq2)
bamDir = paste0(projPath, "/alignment/bam")
countMat = matrix(NA, length(masterPeak), length(histL)*length(repL))
## overlap with bam file to get count
## 计算 Peak 中的片段数
i = 1
for(hist in histL){
  for(rep in repL){
    
    bamFile = paste0(bamDir, "/", hist, "_", rep, "_bowtie2.mapped.bam")
    fragment_counts <- getCounts(bamFile, masterPeak, paired = TRUE, by_rg = FALSE, format = "bam")
    countMat[, i] = counts(fragment_counts)[,1]
    i = i + 1
  }
}
colnames(countMat) = paste(rep(histL, 2), rep(repL, each = 2), sep = "_")

测序深度标准化和差异富集 Peak 分析

##=== R 命令 ===## 
selectR = which(rowSums(countMat) > 5) ## remove low count genes
dataS = countMat[selectR,]
condition = factor(rep(histL, each = length(repL)))
dds = DESeqDataSetFromMatrix(countData = dataS,
                              colData = DataFrame(condition),
                              design = ~ condition)
DDS = DESeq(dds)
normDDS = counts(DDS, normalized = TRUE) ## normalization with respect to the sequencing depth
colnames(normDDS) = paste0(colnames(normDDS), "_norm")
res = results(DDS, independentFiltering = FALSE, altHypothesis = "greaterAbs")
countMatDiff = cbind(dataS, normDDS, res)
head(countMatDiff)

DataFrame with 6 rows and 14 columns
  K27me3_rep1 K4me3_rep1 K27me3_rep2 K4me3_rep2 K27me3_rep1_norm
    <numeric>  <numeric>   <numeric>  <numeric>        <numeric>
1           6          2           1          6         1.408657
2           1          0         242        182         0.234776
3           0          0         176         88         0.000000
4           0          0         274        194         0.000000
5           3          4           0          1         0.704328
6           0          1         109         59         0.000000
  K4me3_rep1_norm K27me3_rep2_norm K4me3_rep2_norm  baseMean log2FoldChange
        <numeric>        <numeric>       <numeric> <numeric>      <numeric>
1        0.620403            4.170         18.2724   6.11787       3.496854
2        0.000000         1009.141        554.2634 390.90978      12.510325
3        0.000000          733.921        267.9955 250.47905      13.297304
4        0.000000         1142.581        590.8082 433.34733      14.089840
5        1.240806            0.000          3.0454   1.24763       0.846266
6        0.310202          454.530        179.6788 158.62986      11.189689
      lfcSE      stat      pvalue        padj
  <numeric> <numeric>   <numeric>   <numeric>
1   1.19893  2.916635 3.53829e-03 4.22134e-02
2   1.50039  8.338074 7.55102e-17 2.18197e-15
3   1.58547  8.386969 4.98837e-17 1.50548e-15
4   1.55196  9.078730 1.09850e-19 6.73560e-18
5   2.18326  0.387617 6.98300e-01 9.72755e-01
6   1.53046  7.311313 2.64545e-13 4.33546e-12

• DESeq2 要求输入矩阵应为非标准化计数或序列 reads 的估计计数
• DESeq2 模型内部修正了文库的大小
• countMatDiff 总结了差异分析结果:
• 前 4 列：在过滤低计数的 Peak 区域后的原始 reads 计数
• 第二个 4 列：标准化后 reads 计数消除文库大小差异
• 其余列：差异分析结果

其他差异分析相关包

• limma: 2015-limma powers differential expression analyses for RNA-sequencing and microarray studies^[2]

limma 是一个用于分析 基因表达芯片数据的 R 软件包，特别是用于分析设计的实验和评估差异表达的线性模型。 limma 提供了在 任意复杂的设计实验中同时分析比较多个 RNA 目标的能力。经验贝叶斯方法被用来提供稳定的结果，即使是数组的数目很小。 Limma 可以扩展到 Peak 区域差异片段富集分析的研究。值得注意的是， limma 可以同时处理 固定效应模型和随机效应模型。

• edgeR: 2012-Differential Expression Analysis of Multifactor RNA-Seq Experiments With Respect to Biological Variation^[3]

利用生物重复 RNA-seq 表达谱的差异表达分析。实施一系列基于 负二项分布的统计方法，包括经验贝叶斯估计，精确检验，广义线性模型和拟似然检验。以及 RNA-seq，它被应用于产生 reads 计数的其他类型的基因组数据的差分信号分析，包括 ChIP-seq, ATAC-seq, Bisulfit-seq, SAGE 和 CAGE。能处理 多因素问题。

参考资料

[1]

2020-Analyzing RNA-seq data with DESeq2: http://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#why-un-normalized-counts

[2]

2015-limma powers differential expression analyses for RNA-sequencing and microarray studies: https://academic.oup.com/nar/article/43/7/e47/2414268

[3]

2012-Differential Expression Analysis of Multifactor RNA-Seq Experiments With Respect to Biological Variation: https://academic.oup.com/nar/article/40/10/4288/2411520