這一章主要介紹了ATAC-seq的基本分析方法
ATAC-seq
這種測序主要用于檢測全基因組的染色體開放程度的一種方法
它的基本原理是用Tn5 transposase隨機(jī)插入開放的染色體區(qū)域（未開放的區(qū)域Tn5插入不進(jìn)去）
那么PCR擴(kuò)增出來的是兩個(gè)Tn5之間的序列，并且但部分片段為100bp
類似于Chip-seq，當(dāng)我們測序的reads重新mapping到參考基因組以后，開放的區(qū)域所測得reads會(huì)形成一個(gè)peak，那么這個(gè)peak所在區(qū)域即為染色質(zhì)開放的區(qū)域
所以整個(gè)上游分析：
1.mapping：bowtie
2.peak calling：MACS
3.peak annotation
QC
ATACseqQC是針對Paired-End測序開發(fā)的
并且對bam文件進(jìn)行QC
## 加載包
library(ATACseqQC)
## input the bamFile from the ATACseqQC package
#讀取file
bamfile <- system.file("extdata", "GL1.bam",
package="ATACseqQC", mustWork=TRUE)
bamfile.labels <- gsub(".bam", "", basename(bamfile))
## generate fragement size distribution
fragSize <- fragSizeDist(bamfile, bamfile.labels)
這幅圖主要是看reads的分布
一般read長度集中于100bp左右
那么我們想看染色質(zhì)開放區(qū)域是否在TSS區(qū)域，這樣TF就可以很輕易的結(jié)合上去，從而發(fā)生生物學(xué)過程
但是Tn5 transposase會(huì)帶來正鏈4堿基負(fù)鏈5堿基的adaptor偏移，所以我們要消除這樣的誤差
#這是bam文件里面的標(biāo)簽
tags <- c("AS", "XN", "XM", "XO", "XG", "NM", "MD", "YS", "YT")
#輸出路徑
outPath <- "splited"
dir.create(outPath)
#矯正5’末端
library(BSgenome.Hsapiens.UCSC.hg19)
seqlev <- "chr1" ##以chr1為例
which <- as(seqinfo(Hsapiens)[seqlev], "GRanges")
gal <- readBamFile(bamfile, tag=tags, which=which, asMates=TRUE)
gal1 <- shiftGAlignmentsList(gal)
shiftedBamfile <- file.path(outPath, "shifted.bam")
export(gal1, shiftedBamfile)
矯正完成以后，我們可以計(jì)算Promoter/Transcript body (PT) score來查看信號是否在promoter區(qū)域富集；利用Nucleosome Free Regions (NFR) score來查看TSS附近打開的情況。
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")
經(jīng)過比較，我們可以看到在整個(gè)基因組上，在TSS附近富集的情況了
同樣，我們也可以在基因組范圍查看一下整體的開放程度
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))
查看NFR（Nucleosome Free Regions）上下游分布情況
#分類
txs <- txs[seqnames(txs) %in% "chr1"]
genome <- Hsapiens
## split the reads into NucleosomeFree, mononucleosome,
## dinucleosome and trinucleosome.
objs <- splitGAlignmentsByCut(gal1, txs=txs, genome=genome)
## Save the binned alignments into bam files.
null <- writeListOfGAlignments(objs, outPath)
這里做分類是為了依據(jù)片段的大小將reads分類為nucleosome free, mononucleosome, dinucleosome, 和trinucleosome，方便后期查看NFR上下游分布情況以及nucleosome的分布狀況
接下來就可以用熱圖可視化
library(ChIPpeakAnno)
bamfiles <- file.path(outPath,
c("NucleosomeFree.bam",
"mononucleosome.bam",
"dinucleosome.bam",
"trinucleosome.bam"))
##這里是四種類型的nucleosome的bam文件
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)
接下來我們可以畫nucleosome-free and nucleosome-bound regions
#獲取nucleosome-free and nucleosome-bound regions的信號
out <- featureAlignedDistribution(sigs, reCenterPeaks(TSS, width=ups+dws), zeroAt=.5, n.tile=NTILE, type="l", ylab="Averaged coverage")
#標(biāo)準(zhǔn)化
range01 <- function(x){(x-min(x))/(max(x)-min(x))}
out <- apply(out, 2, range01)
matplot(out, type="l", xaxt="n",
xlab="Position (bp)",
ylab="Fraction of signal")
axis(1, at=seq(0, 100, by=10)+1,
labels=c("-1K", seq(-800, 800, by=200), "1K"), las=2)
abline(v=seq(0, 100, by=10)+1, lty=2, col="gray")
reCenterPeaks(TSS, width=ups+dws),這幅圖說明了在TSS上下1k的局域內(nèi)peak的情況，是否開放實(shí)在TSS上富集，或是其他區(qū)域
footprint
footprints可以幫助我們查看轉(zhuǎn)錄因子在全基因組上結(jié)合的狀態(tài)。因?yàn)閠ranscription factor結(jié)合在打NFR區(qū)域，所有Tn5 transposase并沒有辦法在TF結(jié)合的區(qū)域插入，所以客觀上形成一個(gè)被保護(hù)的區(qū)域，一個(gè)低coverage區(qū)域
library(MotifDb)
CTCF <- query(MotifDb, c("CTCF"))
CTCF <- as.list(CTCF)
print(CTCF[[1]], digits=2)
sigs <- factorFootprints(shiftedBamfile, pfm=CTCF[[1]],
genome=genome,
min.score="90%", seqlev=seqlev,
upstream=100, downstream=100)
結(jié)合motif，我們可以看到由于Tn5 transposase沒有辦法插入到TF結(jié)合的位置，相當(dāng)于該位置是關(guān)閉著的 ，所以還可以去做一些motif方面的分析
或用V-plot表示：
vp <- vPlot(shiftedBamfile, pfm=CTCF[[1]],
genome=genome, min.score="90%", seqlev=seqlev,
upstream=200, downstream=200,
ylim=c(30, 250), bandwidth=c(2, 1))
distanceDyad(vp, pch=20, cex=.5)
其中CTCF為CTCF蛋白結(jié)合位點(diǎn)
該圖中間位置即為結(jié)合的蛋白
注釋
同Chip-seq
