ATAC-seq(5) -- deeptools可视化及peak

1. deeptools可视化

将bdg文件转为bw文件

ls *treat_pileup.bdg | while read id;
do
sample=${id%%_*}
bedGraphToBigWig ${sample}_treat_pileup.bdg ~/my_data/ref_data/mouse/fasta/mm10.chrom.sizes ../bw/${sample}.bw
done

查看TSS附件信号强度：

cat >tss.sh
ls *.bw | while read id
do
sample=${id%.*}
computeMatrix reference-point  --referencePoint TSS  -p 15  -S ${sample}.bw -b 10000 -a 10000 -R ~/my_data/ref_data/mouse/fasta/mm10.tss.bed --missingDataAsZero --skipZeros  -o ../tss/${sample}_TSS.gz  
done
nohup sh tss.sh &

##     both plotHeatmap and plotProfile will use the output from   computeMatrix
ls *.gz | while read id
do
sample=${id%.*}
plotHeatmap -m ${sample}.gz  --colorMap RdYlBu_r -out ${sample}_Heatmap.png
plotHeatmap -m ${sample}.gz   --colorMap RdYlBu_r -out ${sample}_Heatmap.pdf --plotFileFormat pdf  --dpi 720  
plotProfile -m ${sample}.gz  -out ${sample}_plotProfile.png
plotProfile -m ${sample}.gz  -out ${sample}_plotProfile.pdf --plotFileFormat pdf --perGroup --dpi 720 
done

2-cell-1_TSS_Heatmap.png

查看基因body的信号强度

ls *.bw | while read id
do
sample=${id%.*}
computeMatrix scale-regions  -p 15  \
-S ${sample}.bw -b 10000 -a 10000 -R ~/my_data/ref_data/mouse/fasta/mm10.tss.bed 
 --missingDataAsZero --skipZeros  -o ../tss/${sample}_TSS.gz  
done
#可视化
ls *.gz | while read id
do
sample=${id%.*}
plotHeatmap -m ${sample}.gz  --colorMap RdYlBu_r -out ${sample}_Heatmap.png
plotHeatmap -m ${sample}.gz   --colorMap RdYlBu_r -out ${sample}_Heatmap.pdf --plotFileFormat pdf  --dpi 720  
plotProfile -m ${sample}.gz  -out ${sample}_plotProfile.png
plotProfile -m ${sample}.gz  -out ${sample}_plotProfile.pdf --plotFileFormat pdf --perGroup --dpi 720 
done

2.peaks注释

主要用Y叔的R包ChIPseeker对peaks的位置（如peaks位置落在启动子、UTR、内含子等）以及peaks临近基因的注释。

# 下载老鼠的基因和lincRNA的TxDb对象，#安装所需R包
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("TxDb.Mmusculus.UCSC.mm10.knownGene")

BiocManager::install("org.Mm.eg.db")


#载入各种包
library("ChIPseeker")
library(clusterProfiler)
library("org.Mm.eg.db")
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene

#函数readPeakFile读入summit.bed文件

bedfile <- readPeakFile("2-cell-1_summits.bed")
#注释peaks
peakAnnoList <- annotatePeak(bedfile, tssRegion=c(-5000, 5000), TxDb=txdb, annoDb="org.Mm.eg.db")

promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
tagMatrix <- getTagMatrix(bedfile, 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")

image.png

可视化,展示peak在基因组特征区域的分布以及转录因子在TSS区域的结合。

library(ggplot2)
#饼图
plotAnnoPie(peakAnnoList)

plotDistToTSS(peakAnnoList,title="Distribution of transcription factor-binding loci relative to TSS")
#covplot函数可以直接接受macs2-callpeak的输出bed文件进行画图。
#函数返回的图中，每个竖线表示一个peaks在染色体的位置，线的宽度表示peaks在染色体的分布范围。
covplot(bedfile)

image.png
image.png