今天来看看ATACseq数据质量评估分析。在前面的笔记：有提到使用 ATACseqQC 这个包来对ATAC-seq数据做质量评估，这个软件在最佳实践中也提到了，今天再来看看~

这个软件提供ATAC-seq特定的质量控制指标和可视化工具，例如：

  • 插入片段大小分布。
  • TSS富集。
  • 聚合转录因子足迹可视化。
  • 进行测序饱和度和文库复杂度分析。
  • 等等。

包的网址：https:///packages/release/bioc/vignettes/ATACseqQC/inst/doc/ATACseqQC.html

导入示例数据

注意，此处使用的示例数据集来自人类。如果要使用来自不同物种或不同组装版本的数据集进行分析，请安装相应的BSgenome、TxDb和phastCons。

rm(list=ls())
## load the library
library(ATACseqQC)

## input the bamFile from the ATACseqQC package 
bamfile <- system.file("extdata", "GL1.bam", package="ATACseqQC", mustWork=TRUE)
bamfile
bamfile.labels <- gsub(".bam", "", basename(bamfile))
bamfile.labels

评估文库复杂度

## 1.评估文库复杂度
#bamQC(bamfile, outPath=NULL)
reads <- readsDupFreq(bamfile)
estimateLibComplexity(readsDupFreq(bamfile))

插入片段大小分布

Fragment size distribution：

  • 首先，应该有大量长度小于100 bp的reads，代表无核小体区域。
  • 其次，片段大小分布应具有明显的周期性，在下面的图中右上角清晰可见，表明核小体的占据情况（通常以核小体长度（约147 bp）的整数倍出现）。
## 2.插入片段大小分布
## generate fragement size distribution
bamfile
fragSize <- fragSizeDist(bamfile, bamfile.labels)

核小体定位

核小体在启动子区域的定位对基因转录起始至关重要。核小体占据的启动子区域会阻碍转录因子和RNA聚合酶的结合，从而抑制基因表达。相反，核小体缺乏区域（Nucleosome-Free Regions, NFRs）则为转录因子和RNA聚合酶提供了结合位点，促进基因的转录。

调整读取起始位点

Tn5转座酶已被证明以二聚体形式结合，并将两个接头插入到间隔为9个碱基对的可接近DNA位置中。

因此，对于下游分析，如 peak-calling 和足迹分析，输入bam文件中的所有读取都需要进行偏移。可以使用 shiftGAlignmentsList 函数来偏移读取。默认情况下，所有与正链比对的读取会向右偏移+4bp，而所有与负链比对的读取会向左偏移-5bp。

调整后的读取将被写入一个新的bam文件中，用于峰值检测或足迹分析。

## 3.碱基偏移
## bamfile tags to be read in
possibleTag <- list("integer"=c("AM", "AS", "CM", "CP", "FI", "H0", "H1", "H2", 
                                "HI", "IH", "MQ", "NH", "NM", "OP", "PQ", "SM",
                                "TC", "UQ"), 
                    "character"=c("BC", "BQ", "BZ", "CB", "CC", "CO", "CQ", "CR",
                                  "CS", "CT", "CY", "E2", "FS", "LB", "MC", "MD",
                                  "MI", "OA", "OC", "OQ", "OX", "PG", "PT", "PU",
                                  "Q2", "QT", "QX", "R2", "RG", "RX", "SA", "TS",
                                  "U2"))
possibleTag

library(Rsamtools)
bamTop100 <- scanBam(BamFile(bamfile, yieldSize = 100),
                     param = ScanBamParam(tag=unlist(possibleTag)))[[1]]$tag
tags <- names(bamTop100)[lengths(bamTop100)>0]
tags

## files will be output into outPath
outPath <- "splited"
dir.create(outPath)

## shift the coordinates of 5'ends of alignments in the bam file
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
seqlev <- "chr1"## subsample data for quick run
seqinformation <- seqinfo(TxDb.Hsapiens.UCSC.hg19.knownGene)
which <- as(seqinformation[seqlev], "GRanges")
gal <- readBamFile(bamfile, tag=tags, which=which, asMates=TRUE, bigFile=TRUE)
shiftedBamfile <- file.path(outPath, "shifted.bam")
gal1 <- shiftGAlignmentsList(gal, outbam=shiftedBamfile)

splited目录会生成调整后的bam：

函数的默认参数：shiftGAlignmentsList(gal, positive = 4L, negative = 5L, outbam)

启动子/转录本体（PT）评分

PT评分是通过启动子区域的覆盖度除以相应转录本体区域的覆盖度计算得出的。PT评分可以显示信号是否在启动子区域富集。

library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txs <- transcripts(TxDb.Hsapiens.UCSC.hg19.knownGene)
pt <- PTscore(gal1, txs)
plot(pt$log2meanCoverage, pt$PT_score, 
     xlab="log2 mean coverage",
     ylab="Promoter vs Transcript")

核小体自由区域（NFR）评分

NFR评分是一种用于评估转录起始位点（TSS）附近核小体缺失区域（Nucleosome Free Regions, NFR）的指标。它通过以下步骤计算：

  1. 将TSS周围的400 bp区域划分为三个子区域：最上游的150 bp（n1）、最下游的150 bp（n2）和中间的100 bp（nf）。
  2. 计算每个TSS中5’端与这三个区域重叠的片段数。
  3. 使用公式NFR-score = log2(nf) – log2((n1 + n2) / 2)计算每个TSS的NFR评分。

NFR评分反映了TSS附近核小体的缺失程度。较高的NFR评分表示中间区域（nf）的信号明显高于两侧区域（n1和n2），这通常提示TSS附近存在核小体缺失区域，这对于基因转录的启动非常重要。通过将NFR评分与400 bp窗口的平均信号绘制图表，可以直观地展示TSS区域的核小体分布情况，类似于基因表达分析中的MA图。

nfr <- NFRscore(gal1, txs)
plot(nfr$log2meanCoverage, nfr$NFR_score, 
     xlab="log2 mean coverage",
     ylab="Nucleosome Free Regions score",
     main="NFRscore for 200bp flanking TSSs",
     xlim=c(-10, 0), ylim=c(-5, 5))

转录起始位点（TSS）富集评分

TSS富集评分是一种用于评估转录起始位点（TSS）附近信号富集程度的指标。它通过以下步骤计算：

ATACseq数据质量如何？来看看 ATACseqQC这个软件~
  1. 收集围绕TSS中心的读取分布，范围为TSS两侧各1000 bp，分为多个100 bp窗口。
  2. 计算每个100 bp窗口的读取深度。
  3. 计算末端侧翼区域（每端100 bp）的平均读取深度。
  4. TSS评分 = 每个100 bp窗口的读取深度 / 末端侧翼区域的平均读取深度。
  5. TSS富集评分 = 所有窗口中TSS评分的最大平均值。

该评分反映了TSS区域的信号富集程度，较高的TSS富集评分表明TSS区域的信号显著高于侧翼区域，这通常提示基因转录的活跃性。TSS富集值的具体数值取决于所使用的参考文件，而高质量数据的截止值可以在ENCODE项目的ATAC-seq数据标准和处理流程页面中找到。

tsse <- TSSEscore(gal1, txs)
tsse$TSSEscore
plot(100*(-9:10-.5), tsse$values, type="b", 
     xlab="distance to TSS",
     ylab="aggregate TSS score")

bam分割

将reads分为四个类别：核小体自由区域、单核小体、双核小体和三核小体。

############### split reads
library(BSgenome.Hsapiens.UCSC.hg19)
library(phastCons100way.UCSC.hg19)
## run program for chromosome 1 only
txs <- txs[seqnames(txs) %in% "chr1"]
genome <- Hsapiens
## split the reads into NucleosomeFree, mononucleosome, 
## dinucleosome and trinucleosome.
## and save the binned alignments into bam files.
objs <- splitGAlignmentsByCut(gal1, txs=txs, genome=genome, outPath = outPath,
                              conservation=phastCons100way.UCSC.hg19)
## list the files generated by splitGAlignmentsByCut.
dir(outPath)


objs <- splitBam(bamfile, tags=tags, outPath=outPath,
                 txs=txs, genome=genome,
                 conservation=phastCons100way.UCSC.hg19)
                 
##  [1] "NucleosomeFree.bam"     "NucleosomeFree.bam.bai" "dinucleosome.bam"      
##  [4] "dinucleosome.bam.bai"   "inter1.bam"             "inter1.bam.bai"        
##  [7] "inter2.bam"             "inter2.bam.bai"         "inter3.bam"            
## [10] "inter3.bam.bai"         "mononucleosome.bam"     "mononucleosome.bam.bai"
## [13] "others.bam"             "others.bam.bai"         "shifted.bam"           
## [16] "shifted.bam.bai"        "trinucleosome.bam"      "trinucleosome.bam.bai"

核小体位置的热图和覆盖曲线

通过对所有活跃的转录起始位点（TSS）的信号进行平均，我们应观察到核小体缺失片段在TSS处富集，而核小体结合片段在活跃TSS的上游和下游均富集，并显示出上游和下游核小体的特征性相位。由于ATAC-seq的读取集中在开放染色质区域，用户应在+1核小体处看到强烈的核小体信号，但在+2、+3和+4核小体处信号会减弱。

library(ChIPpeakAnno)
bamfiles <- file.path(outPath,
                     c("NucleosomeFree.bam",
                     "mononucleosome.bam",
                     "dinucleosome.bam",
                     "trinucleosome.bam"))
## Plot the cumulative percentage of tag allocation in nucleosome-free 
## and mononucleosome bam files.
cumulativePercentage(bamfiles[1:2], as(seqinformation["chr1"], "GRanges"))

计算信号和TSS值：

TSS <- promoters(txs, upstream=0, downstream=1)
TSS <- unique(TSS)
## estimate the library size for normalization
(librarySize <- estLibSize(bamfiles))

## calculate the signals around TSSs.
NTILE <- 101
dws <- ups <- 1010
sigs <- enrichedFragments(gal=objs[c("NucleosomeFree", 
                                     "mononucleosome",
                                     "dinucleosome",
                                     "trinucleosome")], 
                          TSS=TSS,
                          librarySize=librarySize,
                          seqlev=seqlev,
                          TSS.filter=0.5,
                          n.tile = NTILE,
                          upstream = ups,
                          downstream = dws)
## log2 transformed signals
sigs.log2 <- lapply(sigs, function(.ele) log2(.ele+1))
#plot heatmap
featureAlignedHeatmap(sigs.log2, reCenterPeaks(TSS, width=ups+dws),
                      zeroAt=.5, n.tile=NTILE)

结果如下：可以看到核小体缺失片段在TSS处富集

本次分享到这~