这里是佳奥！新的一年，ATAC-Seq的学习也进入了尾声。
让我们开始吧！
1 peaks注释
统计peak在promoter，exon，intron和intergenic区域的分布
QQ截图20230101155625.png
if(F){
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
source("http://bioconductor.org/biocLite.R")
BiocManager::install('TxDb.Mmusculus.UCSC.mm10.knownGene')
BiocManager::install('org.Mm.eg.db')
}
bedPeaksFile = '2-ceLL-1_peaks.narrowPeak';
bedPeaksFile
## loading packages
require(ChIPseeker)
require(TxDb.Mmusculus.UCSC.mm10.knownGene)
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
require(clusterProfiler)
peak <- readPeakFile( bedPeaksFile )
##去除含_的染色体
keepChr= !grepl('_',seqlevels(peak))
seqlevels(peak, pruning.mode="coarse") <- seqlevels(peak)[keepChr]
peakAnno <- annotatePeak(peak, tssRegion=c(-3000, 3000),
TxDb=txdb, annoDb="org.Mm.eg.db")
peakAnno_df <- as.data.frame(peakAnno)
promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
tagMatrix <- getTagMatrix(peak, windows=promoter)
# 然后查看这些peaks在所有基因的启动子附近的分布情况，热图模式
tagHeatmap(tagMatrix, xlim=c(-3000, 3000), color="red")
# 然后查看这些peaks在所有基因的启动子附近的分布情况，信号强度曲线图
plotAvgProf(tagMatrix, xlim=c(-3000, 3000),
xlab="Genomic Region (5'->3')", ylab = "Read Count Frequency")
plotAnnoPie(peakAnno)
QQ截图20230101160301.png
QQ截图20230101160325.png
QQ截图20230101160332.png
可以载入IGV看看效果，检测软件找到的peaks是否真的合理，还可以配合rmarkdown来出自动化报告。
https://ke.qq.com/course/274681
我们可以看到Tcea1基因转录起始位置有peaks富集
QQ截图20230101160746.png
也可以使用其它代码进行下游分析；
https://github.com/jmzeng1314/NGS-pipeline/tree/master/CHIPseq
Homer 可以做，但是需要下载数据库
# perl ~/miniconda3/envs/atac/share/homer-4.9.1-5/configureHomer.pl -install mm10
# ln -s /home/jmzeng/miniconda3/envs/chipseq/share/homer-4.9.1-5/data/genomes/ genomes
# cp /home/jmzeng/miniconda3/envs/chipseq/share/homer-4.9.1-5/config.txt /home/stu/miniconda3/envs/atac/share/homer-4.9.1-5/config.txt
## 保证数据库下载是OK
ls -lh ~/miniconda3/envs/atac/share/homer-4.9.1-5/data/genomes
mkdir -p ~/project/atac/peaks
source activate atac
cd ~/project/atac/peaks
ls *.narrowPeak |while read id;
do
echo $id
awk '{print $4"\t"$1"\t"$2"\t"$3"\t+"}' $id >{id%%.*}.homer_peaks.tmp
annotatePeaks.pl {id%%.*}.homer_peaks.tmp mm10 1>${id%%.*}.peakAnn.xls
2>${id%%.*}.annLog.txt
done
Bedtools也可以做
https://bedtools.readthedocs.io/en/latest/content/tools/annotate.html
2 motif寻找及注释
Homer可以做
ls -lh ~/miniconda3/envs/atac/share/homer-4.9.1-5/data/genomes
mkdir -p ~/project/atac/motif
cd ~/project/atac/motif
source activate atac
ls ../peaks/*.narrowPeak |while read id;
do
file=$(basename $id )
sample=${file%%.*}
echo $sample
awk '{print $4"\t"$1"\t"$2"\t"$3"\t+"}' $id > ${sample}.homer_peaks.tmp
nohup findMotifsGenome.pl ${sample}.homer_peaks.tmp mm10 ${sample}_motifDir -len 8,10,12 &
done
meme 也可以做 ，首先利用.bed获取.fa序列：
https://github.com/jmzeng1314/NGS-pipeline/blob/master/CHIPseq/step7-peaks2sequence.R
##usage: Rscript peakView.R peaks.bed IP.sorted.bam input.sorted.bam 10
#options(echo=TRUE) # if you want see commands in output file
args <- commandArgs(trailingOnly = TRUE)
if(length(args) != 1 ){
print(" usage: Rscript peakAnno.R peaks.bed ")
}
bedPeaksFile = args[1] ;
##自这开始，.bed文件要和R Project文件在同一目录下
bedFiles=list.files(pattern = '*.bed')
> bedFiles
[1] "2-ce11-2_summits.bed" "2-ce11-4_summits.bed" "2-ce11-5_summits.bed" "2-ceLL-1_summits.bed"
BiocManager::install("BSgenome.Mmusculus.UCSC.mm10")
library(BSgenome.Mmusculus.UCSC.mm10)
library(ChIPpeakAnno)
##生成.fa文件
bedPeaksFile=bedFiles[2]##第二个文件即2-ce11-4_summits.bed，要下一个就[3]
sampleName=strsplit(bedPeaksFile,'\\.')[[1]][1]
peak <- toGRanges(bedPeaksFile, format="BED")
keepChr= !grepl('_',seqlevels(peak))
#seqlevels(peak, force=TRUE) <- seqlevels(peak)[keepChr]
seq <- getAllPeakSequence(peak, upstream=20, downstream=20, genome=Mmusculus)
write2FASTA(seq, paste0(sampleName,'.fa'))
使用网页端注释
https://meme-suite.org/meme/
QQ截图20230101170341.png
QQ截图20230101170412.png
R包，比如
motifmatchr包 也可以做。
https://bioconductor.org/packages/release/bioc/html/motifmatchr.html
3 多组学整合分析
RNS-Seq、ChIP-Seq、ATAC-Seq
以及一些整合的R包：esATAC
ATAC-Seq的实战流程学习至此结束。
但个性化的分析还有很多要钻研的地方，尤其是官方文档。
新年快乐！我们下一个篇章再见！
