引言
本系列讲解 CUT&Tag 数据处理和分析教程 教程，持续更新，欢迎关注，转发！
差异分析
- DESeq2：通过 DESeq2 对 RNA-seq 数据进行差异分析，适度估算基因表达的倍数变化和离散度
在高通量测序实验中，分析计数数据的方差与均值之间的关系，并利用负二项分布模型来检测基因表达的差异。
构建样本矩阵
一般来说，差异分析主要是比较同一组蛋白修饰在不同条件下的表达情况。由于本教程使用的数据有限，只能通过比较 H3K27me3 和 H3K4me3 的两组重复样本来进行差异分析的演示。这里以 DESeq2为例。
- 生成一个主峰列表，整合所有样本所检测到的峰。
##=== R command ===##
mPeak = GRanges()
## overlap with bam file to get count
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)
- 统计主峰列表中每个峰的片段数量。
##=== R command ===##
library(DESeq2)
bamDir = paste0(projPath, "/alignment/bam")
countMat = matrix(NA, length(masterPeak), length(histL)*length(repL))
## overlap with bam file to get count
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 = "_")
- 测序深度标准化与差异富集峰的检测
##=== R command ===##
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)
DESeq2 要求输入的矩阵为未经标准化的测序读长计数或估计计数。
DESeq2 模型会自动校正文库大小。
countMatDiff 汇总了差异分析的结果：
- 前 4 列：过滤掉低计数峰区域后的原始读长计数。
- 接下来的 4 列：经过标准化处理的读长计数，消除了文库大小的影响。
- 剩余列：差异检测的结果。
