CUT&Tag 分析教程 | 完结撒花

引言

本系列讲解 CUT&Tag 数据处理和分析教程 教程，持续更新，欢迎关注，转发！

CUT&Tag 数据处理和分析教程（9）

CUT&Tag 数据处理和分析教程（8）

CUT&Tag 数据处理和分析教程（7）

CUT&Tag 数据处理和分析教程（6）

CUT&Tag 数据处理和分析教程（5）

CUT&Tag 数据处理和分析教程（4）

CUT&Tag 数据处理和分析教程（3）

CUT&Tag 数据处理和分析教程（2）

CUT&Tag 数据处理和分析教程（1）

差异分析

• 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 列：经过标准化处理的读长计数，消除了文库大小的影响。
  • 剩余列：差异检测的结果。

本文由mdnice多平台发布