本次給學(xué)徒講解的文章是 :The landscape of accessible chromatin in mammalian preimplantation embryos. Nature 2016 Jun 30;534(7609):652-7. PMID: 27309802
查看文章發(fā)現(xiàn)數(shù)據(jù)上傳到了GEO，是：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE66581
在SRA數(shù)據(jù)庫可以下載原始測(cè)序數(shù)據(jù) , 從文章找到數(shù)據(jù)的ID： https://www.ncbi.nlm.nih.gov/sra?term=SRP055881 把下面的內(nèi)容保存到文件，命名為 srr.list 就可以使用prefetch這個(gè)函數(shù)來下載店读。
linux環(huán)境及軟件安裝
這里首推conda
# https://mirrors.tuna.tsinghua.edu.cn/help/anaconda/
# https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/
wget https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-latest-Linux-x86_64.sh
bash Miniconda3-latest-Linux-x86_64.sh
source ~/.bashrc
## 安裝好conda后需要設(shè)置鏡像。
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
conda config --set show_channel_urls yes
conda create -n atac -y python=2 bwa
conda info --envs
source activate atac
# 可以用search先進(jìn)行檢索
conda search trim_galore
## 保證所有的軟件都是安裝在 wes 這個(gè)環(huán)境下面
conda install -y sra-tools sambamba
conda install -y trim-galore samtools bedtools
conda install -y deeptools homer meme
conda install -y macs2 bowtie bowtie2
conda install -y multiqc
代碼里面提到的軟件阅悍，都是根據(jù)我們對(duì)ATAC-seq搜索學(xué)習(xí)總結(jié)的。
值得一提的是自己為了ATAC-seq建立的軟件環(huán)境可以很方便移植到另外一臺(tái)電腦!
首先通過activate target_env要分享的環(huán)境target_env昨稼，然后輸入下面的命令會(huì)在當(dāng)前工作目錄下生成一個(gè)environment.yml文件节视，
conda env export > environment.yml
小伙伴拿到environment.yml文件后，將該文件放在工作目錄下假栓，可以通過以下命令從該文件創(chuàng)建環(huán)境
conda env create -f environment.yml
下載作者的數(shù)據(jù)
前面提到的SRA數(shù)據(jù)庫寻行，該文章配套數(shù)據(jù)太多，我們節(jié)選部分作為練習(xí)匾荆，文件config.sra 如下：
2-cell-1 SRR2927015
2-cell-2 SRR2927016
2-cell-5 SRR3545580
2-cell-4 SRR2927018
因?yàn)閏onda安裝好了sra-toolkit拌蜘，所以prefetch函數(shù)可以直接使用
## 下載數(shù)據(jù)
# cat srr.list |while read id;do (nohup $prefetch $id -X 100G & );done
## 注意組織好自己的項(xiàng)目
mkdir -p ~/project/atac/
cd ~/project/atac/
mkdir {sra,raw,clean,align,peaks,motif,qc}
cd sra
## vim 或者cat命令創(chuàng)建 srr.list 文件, 里面保存著作為練習(xí)使用的4個(gè)數(shù)據(jù)ID
source activate atac
cat srr.list |while read id;do ( nohup prefetch $id & );done
## 默認(rèn)下載目錄：~/ncbi/public/sra/
ls -lh ~/ncbi/public/sra/
## 下載耗時(shí)，自行解決牙丽，學(xué)員使用現(xiàn)成數(shù)據(jù)：/public/project/epi/atac/sra
## 假如提前下載好了數(shù)據(jù)简卧。
cd ~/project/atac/
ln -s /public/project/epi/atac/sra sra
總之?dāng)?shù)據(jù)如下：
-rw-r--r-- 1 stu stu 4.2G Aug 25 11:10 SRR2927015.sra
-rw-r--r-- 1 stu stu 5.5G Aug 25 11:13 SRR2927016.sra
-rw-r--r-- 1 stu stu 2.0G Aug 25 11:12 SRR2927018.sra
-rw-r--r-- 1 stu stu 7.0G Aug 25 11:13 SRR3545580.sra
第一步，得到fastq測(cè)序數(shù)據(jù)
通常我們應(yīng)該是自己的實(shí)驗(yàn)數(shù)據(jù)烤芦，自己找公司測(cè)序后拿到原始數(shù)據(jù)举娩，本次講解使用的是公共數(shù)據(jù)，所以需要把下載的sra數(shù)據(jù)轉(zhuǎn)換為fq格式构罗。
## 下面需要用循環(huán)
cd ~/project/atac/
source activate atac
dump=fastq-dump
analysis_dir=raw
mkdir -p $analysis_dir
## 下面用到的 config.sra 文件铜涉，就是上面自行制作的。
# $fastq-dump sra/SRR2927015.sra --gzip --split-3 -A 2-cell-1 -O clean/
cat config.sra |while read id;
do echo $id
arr=($id)
srr=${arr[1]}
sample=${arr[0]}
# 測(cè)序數(shù)據(jù)的sra轉(zhuǎn)fasq
nohup $dump -A $sample -O $analysis_dir --gzip --split-3 sra/$srr.sra &
done
### 如果不只是4個(gè)文件遂唧，需要使用shell腳本批處理芙代。
cut -f 10,13 SRP055881/SraRunTable.txt|\
sed 's/Embryonic stem cell/ESC/'|sed 's/early 2-cell/e2-cell/' |\
perl -alne '{$h{$F[1]}++;print "$_-$h{$F[1]}"}' |tail -n+2|awk '{print $2"\t"$1}'> config.sra
得到的原始fq數(shù)據(jù)如下:
-rw-rw-r-- 1 jmzeng jmzeng 2.6G Aug 24 23:10 2-cell-1_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 2.6G Aug 24 23:10 2-cell-1_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 3.4G Aug 24 23:31 2-cell-2_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 3.7G Aug 24 23:31 2-cell-2_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.2G Aug 24 22:46 2-cell-4_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.2G Aug 24 22:46 2-cell-4_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 4.4G Aug 24 23:52 2-cell-5_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 4.9G Aug 24 23:52 2-cell-5_2.fastq.gz
第二步，測(cè)序數(shù)據(jù)的質(zhì)量控制
這個(gè)時(shí)候選擇trim_galore軟件進(jìn)行過濾蠢箩，雙端測(cè)序數(shù)據(jù)的代碼如下:
需要自行制作 config.raw 文件链蕊， 是3列，第一列占位用谬泌，沒有意義，第二列是fq1的地址逻谦，第3列是fq2的地址掌实。
cd ~/project/atac/
mkdir -p clean
source activate atac
# trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 4 --paired -o clean/ raw/2-cell-1_1.fastq.gz raw/2-cell-1_2.fastq.gz
cat config.raw |while read id;
do echo $id
arr=($id)
fq2=${arr[2]}
fq1=${arr[1]}
sample=${arr[0]}
nohup trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 4 --paired -o clean $fq1 $fq2 &
done
ps -ef |grep trim
得到過濾后的fq文件如下：
-rw-rw-r-- 1 jmzeng jmzeng 2.4G Aug 25 09:35 2-cell-1_1_val_1.fq.gz
-rw-rw-r-- 1 jmzeng jmzeng 2.3G Aug 25 09:35 2-cell-1_2_val_2.fq.gz
-rw-rw-r-- 1 jmzeng jmzeng 3.1G Aug 25 10:10 2-cell-2_1_val_1.fq.gz
-rw-rw-r-- 1 jmzeng jmzeng 3.3G Aug 25 10:10 2-cell-2_2_val_2.fq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.1G Aug 25 08:52 2-cell-4_1_val_1.fq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.1G Aug 25 08:52 2-cell-4_2_val_2.fq.gz
-rw-rw-r-- 1 jmzeng jmzeng 3.7G Aug 25 10:27 2-cell-5_1_val_1.fq.gz
-rw-rw-r-- 1 jmzeng jmzeng 3.9G Aug 25 10:27 2-cell-5_2_val_2.fq.gz
質(zhì)量控制前后都需要可視化，肯定是fastqc+multiqc邦马，代碼如下贱鼻；
cd ~/project/atac/qc
mkdir -p clean
fastqc -t 5 ../clean/2-cell-*gz -o clean
mkdir -p raw
fastqc -t 5 ../raw/2-cell-*gz -o clean
# https://zh.wikipedia.org/wiki/ASCII
## 還有很多其它工具宴卖，比如：
qualimap='/home/jianmingzeng/biosoft/Qualimap/qualimap_v2.2.1/qualimap'
$qualimap bamqc --java-mem-size=20G -bam $id -outdir ./
第三步，比對(duì)
比對(duì)需要的index邻悬，看清楚物種症昏，根據(jù)對(duì)應(yīng)的軟件來構(gòu)建，這里直接用bowtie2進(jìn)行比對(duì)和統(tǒng)計(jì)比對(duì)率, 需要提前下載參考基因組然后使用命令構(gòu)建索引父丰，或者直接就下載索引文件：下載小鼠參考基因組的索引和注釋文件, 這里用常用的mm10
# 索引大小為3.2GB肝谭， 不建議自己下載基因組構(gòu)建，可以直接下載索引文件蛾扇，代碼如下：
mkdir referece && cd reference
wget -4 -q ftp://ftp.ccb.jhu.edu/pub/data/bowtie2_indexes/mm10.zip
unzip mm10.zip
解壓后的索引如下：
848M Jul 5 05:03 /public/reference/index/bowtie/mm10.1.bt2
633M Jul 5 05:00 /public/reference/index/bowtie/mm10.2.bt2
6.0K Jul 5 05:05 /public/reference/index/bowtie/mm10.3.bt2
633M Jul 5 05:05 /public/reference/index/bowtie/mm10.4.bt2
848M Jul 5 04:52 /public/reference/index/bowtie/mm10.rev.1.bt2
633M Jul 5 04:49 /public/reference/index/bowtie/mm10.rev.2.bt2
這些索引文件一個(gè)都不能少攘烛，而且文件名的前綴很重要，保證一致镀首。
單端測(cè)序數(shù)據(jù)的比對(duì)代碼如下：
首先可以對(duì)測(cè)試樣本走流程坟漱，完善代碼:
zcat ../clean/2-cell-1_1_val_1.fq.gz |head -10000 > test1.fq
zcat ../clean/2-cell-1_2_val_2.fq.gz |head -10000 > test2.fq
bowtie2 -x /public/reference/index/bowtie/mm10 -1 test1.fq -2 test2.fq
bowtie2 -x /public/reference/index/bowtie/mm10 -1 test1.fq -2 test2.fq -S test.sam
bowtie2 -x /public/reference/index/bowtie/mm10 -1 test1.fq -2 test2.fq |samtools sort -@ 5 -O bam -o test.bam -
## 建議拋棄 samtools markdup功能，避免麻煩更哄。
## http://www.reibang.com/p/1e6189f641db
samtools markdup -r test.bam test.samtools.rmdup.bam
## 把報(bào)錯(cuò)信息在谷歌搜索后芋齿，在兩個(gè)網(wǎng)頁上找到了答案。
https://github.com/samtools/samtools/issues/765
https://www.biostars.org/p/288496/
## gatk 可以在GitHub下載
/public/biosoft/GATK/gatk-4.0.3.0/gatk MarkDuplicates \
-I test.bam -O test.picard.rmdup.bam --REMOVE_SEQUENCING_DUPLICATES true -M test.log
### picards 被包裝在GATK里面：
### sambamba 文檔： http://lomereiter.github.io/sambamba/docs/sambamba-markdup.html
conda install -y sambamba
sambamba --help
sambamba markdup --help
sambamba markdup -r test.bam test.sambamba.rmdup.bam
samtools flagstat test.sambamba.rmdup.bam
samtools flagstat test.bam
## 接下來只保留兩條reads要比對(duì)到同一條染色體(Proper paired) 成翩，還有高質(zhì)量的比對(duì)結(jié)果(Mapping quality>=30)
## 順便過濾 線粒體reads
samtools view -f 2 -q 30 test.sambamba.rmdup.bam |grep -v chrM|wc
samtools view -f 2 -q 30 test.sambamba.rmdup.bam |wc
samtools view -h -f 2 -q 30 test.sambamba.rmdup.bam |grep -v chrM| samtools sort -O bam -@ 5 -o - > test.last.bam
bedtools bamtobed -i test.last.bam > test.bed
ls *.bam |xargs -i samtools index {}
探索好了整個(gè)流程沟突，就可以直接寫批處理，代碼如下：
ls /home/jmzeng/project/atac/clean/*_1.fq.gz > 1
ls /home/jmzeng/project/atac/clean/*_2.fq.gz > 2
ls /home/jmzeng/project/atac/clean/*_2.fq.gz |cut -d"/" -f 7|cut -d"_" -f 1 > 0
paste 0 1 2 > config.clean ## 供mapping使用的配置文件
cd ~/project/epi/align
## 相對(duì)目錄需要理解
bowtie2_index=/public/reference/index/bowtie/mm10
## 一定要搞清楚自己的bowtie2軟件安裝在哪里捕传，以及自己的索引文件在什么地方；菔谩！庸论！
#source activate atac
cat config.clean |while read id;
do echo $id
arr=($id)
fq2=${arr[2]}
fq1=${arr[1]}
sample=${arr[0]}
## 比對(duì)過程15分鐘一個(gè)樣本
bowtie2 -p 5 --very-sensitive -X 2000 -x $bowtie2_index -1 $fq1 -2 $fq2 |samtools sort -O bam -@ 5 -o - > ${sample}.raw.bam
samtools index ${sample}.raw.bam
bedtools bamtobed -i ${sample}.raw.bam > ${sample}.raw.bed
samtools flagstat ${sample}.raw.bam > ${sample}.raw.stat
# https://github.com/biod/sambamba/issues/177
sambamba markdup --overflow-list-size 600000 --tmpdir='./' -r ${sample}.raw.bam ${sample}.rmdup.bam
samtools index ${sample}.rmdup.bam
## ref:https://www.biostars.org/p/170294/
## Calculate %mtDNA:
mtReads=$(samtools idxstats ${sample}.rmdup.bam | grep 'chrM' | cut -f 3)
totalReads=$(samtools idxstats ${sample}.rmdup.bam | awk '{SUM += $3} END {print SUM}')
echo '==> mtDNA Content:' $(bc <<< "scale=2;100*$mtReads/$totalReads")'%'
samtools flagstat ${sample}.rmdup.bam > ${sample}.rmdup.stat
samtools view -h -f 2 -q 30 ${sample}.rmdup.bam |grep -v chrM |samtools sort -O bam -@ 5 -o - > ${sample}.last.bam
samtools index ${sample}.last.bam
samtools flagstat ${sample}.last.bam > ${sample}.last.stat
bedtools bamtobed -i ${sample}.last.bam > ${sample}.bed
done
其中bowtie2比對(duì)加入了-X 2000 參數(shù)职辅，是最大插入片段，寬泛的插入片段范圍(10-1000bp)
第一步得到的bam文件如下：
-rw-rw-r-- 1 stu stu 3.7G Aug 25 14:17 2-cell-1.bam
-rw-rw-r-- 1 stu stu 4.6G Aug 25 15:32 2-cell-2.bam
-rw-rw-r-- 1 stu stu 1.8G Aug 25 15:47 2-cell-4.bam
-rw-rw-r-- 1 stu stu 5.5G Aug 25 16:49 2-cell-5.bam
過濾后的bam文件是：
3.7G Aug 25 21:08 2-cell-1.bam
490M Aug 25 21:14 2-cell-1.last.bam
776M Aug 25 21:13 2-cell-1.rmdup.bam
4.6G Aug 25 23:51 2-cell-2.bam
678M Aug 25 23:58 2-cell-2.last.bam
1.1G Aug 25 23:57 2-cell-2.rmdup.bam
1.8G Aug 26 00:41 2-cell-4.bam
427M Aug 26 00:43 2-cell-4.last.bam
586M Aug 26 00:43 2-cell-4.rmdup.bam
5.5G Aug 26 03:26 2-cell-5.bam
523M Aug 26 03:33 2-cell-5.last.bam
899M Aug 26 03:32 2-cell-5.rmdup.bam
上述腳本的步驟都可以拆分運(yùn)行聂示，比如bam文件構(gòu)建index或者轉(zhuǎn)為bed的:
ls *.last.bam|xargs -i samtools index {}
ls *.last.bam|while read id;do (bedtools bamtobed -i $id >${id%%.*}.bed) ;done
ls *.raw.bam|while read id;do (nohup bedtools bamtobed -i $id >${id%%.*}.raw.bed & ) ;done
最后得到的bed文件是：
237M Aug 26 08:00 2-cell-1.bed
338M Aug 26 08:01 2-cell-2.bed
203M Aug 26 08:01 2-cell-4.bed
254M Aug 26 08:01 2-cell-5.bed
第4步域携，使用macs2找peaks
# macs2 callpeak -t 2-cell-1.bed -g mm --nomodel --shift -100 --extsize 200 -n 2-cell-1 --outdir ../peaks/
ls *.bed | while read id ;do (macs2 callpeak -t $id -g mm --nomodel --sHit -100 --extsize 200 -n ${id%%.*} --outdir ../peaks/) ;done
## shell 13問
macs2軟件說明書詳見：http://www.reibang.com/p/21e8c51fca23
第5步，計(jì)算插入片段長度鱼喉，F(xiàn)RiP值秀鞭，IDR計(jì)算重復(fù)情況
非冗余非線粒體能夠比對(duì)的fragment、比對(duì)率扛禽、NRF锋边、PBC1、PBC2编曼、peak數(shù)豆巨、無核小體區(qū)NFR、TSS富集掐场、FRiP 往扔、IDR重復(fù)的一致性贩猎！
名詞解釋：https://www.encodeproject.org/data-standards/terms/
參考：https://www.encodeproject.org/atac-seq/
看 sam文件第9列，在R里面統(tǒng)計(jì)繪圖
cmd=commandArgs(trailingOnly=TRUE);
input=cmd[1]; output=cmd[2];
a=abs(as.numeric(read.table(input)[,1]));
png(file=output);
hist(a,
main="Insertion Size distribution",
ylab="Read Count",xlab="Insert Size",
xaxt="n",
breaks=seq(0,max(a),by=10)
);
axis(side=1,
at=seq(0,max(a),by=100),
labels=seq(0,max(a),by=100)
);
dev.off()
有了上面的繪圖R腳本就可以在批量檢驗(yàn)bam文件進(jìn)行出圖萍膛。
還有NFR:https://github.com/GreenleafLab/NucleoATAC/issues/18
FRiP值的計(jì)算：fraction of reads in called peak regions
bedtools intersect -a ../new/2-cell-1.bed -b 2-cell-1_peaks.narrowPeak |wc -l
148928
wc ../new/2-cell-1.bed
5105844
wc ../new/2-cell-1.raw.bed
5105844
### 搞清楚 FRiP值具體定義：
ls *narrowPeak|while read id;
do
echo $id
bed=../new/$(basename $id "_peaks.narrowPeak").raw.bed
#ls -lh $bed
Reads=$(bedtools intersect -a $bed -b $id |wc -l|awk '{print $1}')
totalReads=$(wc -l $bed|awk '{print $1}')
echo $Reads $totalReads
echo '==> FRiP value:' $(bc <<< "scale=2;100*$Reads/$totalReads")'%'
done
Fraction of reads in peaks (FRiP) - Fraction of all mapped reads that fall into the called peak regions, i.e. usable reads in significantly enriched peaks divided by all usable reads. In general, FRiP scores correlate positively with the number of regions. (Landt et al, Genome Research Sept. 2012, 22(9): 1813–1831)
文章其它指標(biāo)：https://www.nature.com/articles/sdata2016109/tables/4
可以使用R包看不同peaks文件的overlap情況吭服。
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")
library('BiocInstaller')
biocLite("ChIPpeakAnno")
biocLite("ChIPseeker")
}
library(ChIPseeker)
library(ChIPpeakAnno)
list.files('project/atac/peaks/',"*.narrowPeak")
tmp=lapply(list.files('project/atac/peaks/',"*.narrowPeak"),function(x){
return(readPeakFile(file.path('project/atac/peaks/', x)))
})
ol <- findOverlapsOfPeaks(tmp[[1]],tmp[[2]])
png('overlapVenn.png')
makeVennDiagram(ol)
dev.off()
也可以使用專業(yè)軟件，IDR 來進(jìn)行計(jì)算出來蝗罗，同時(shí)考慮peaks間的overlap艇棕，和富集倍數(shù)的一致性 。
source activate atac
# 可以用search先進(jìn)行檢索
conda search idr
source deactivate
## 保證所有的軟件都是安裝在 wes 這個(gè)環(huán)境下面
conda create -n py3 -y python=3 idr
conda activate py3
idr -h
idr --samples 2-cell-1_peaks.narrowPeak 2-cell-2_peaks.narrowPeak --plot
結(jié)果如下：
Initial parameter values: [0.10 1.00 0.20 0.50]
Final parameter values: [0.00 1.06 0.64 0.87]
Number of reported peaks - 5893/5893 (100.0%)
Number of peaks passing IDR cutoff of 0.05 - 674/5893 (11.4%)
參考：https://www.biostat.wisc.edu/~kendzior/STAT877/SLIDES/keles3.pdf
第6步绿饵，deeptools的可視化
具體仍然是見：https://mp.weixin.qq.com/s/a4qAcKE1DoukpLVV_ybobA 在ChiP-seq 講解欠肾。
首先把bam文件轉(zhuǎn)為bw文件，詳情：http://www.bio-info-trainee.com/1815.html
cd ~/project/atac/new
source activate atac
#ls *.bam |xargs -i samtools index {}
ls *last.bam |while read id;do
nohup bamCoverage -p 5 --normalizeUsing CPM -b $id -o ${id%%.*}.last.bw &
done
cd dup
ls *.bam |xargs -i samtools index {}
ls *.bam |while read id;do
nohup bamCoverage --normalizeUsing CPM -b $id -o ${id%%.*}.rm.bw &
done
查看TSS附件信號(hào)強(qiáng)度：
## both -R and -S can accept multiple files
mkdir -p ~/project/atac/tss
cd ~/project/atac/tss
source activate atac
computeMatrix reference-point --referencePoint TSS -p 15 \
-b 10000 -a 10000 \
-R /public/annotation/CHIPseq/mm10/ucsc.refseq.bed \
-S ~/project/atac/new/*.bw \
--skipZeros -o matrix1_test_TSS.gz \
--outFileSortedRegions regions1_test_genes.bed
## both plotHeatmap and plotProfile will use the output from computeMatrix
plotHeatmap -m matrix1_test_TSS.gz -out test_Heatmap.png
plotHeatmap -m matrix1_test_TSS.gz -out test_Heatmap.pdf --plotFileFormat pdf --dpi 720
plotProfile -m matrix1_test_TSS.gz -out test_Profile.png
plotProfile -m matrix1_test_TSS.gz -out test_Profile.pdf --plotFileFormat pdf --perGroup --dpi 720
### 如果要批處理 拟赊，需要學(xué)習(xí)好linux命令刺桃。
下載 bed文件：https://genome.ucsc.edu/cgi-bin/hgTables 只需要3列坐標(biāo)格式文件。
查看基因body的信號(hào)強(qiáng)度
source activate atac
computeMatrix scale-regions -p 15 \
-R /public/annotation/CHIPseq/mm10/ucsc.refseq.bed \
-S ~/project/atac/new/*.bw \
-b 10000 -a 10000 \
--skipZeros -o matrix1_test_body.gz
plotHeatmap -m matrix1_test_body.gz -out ExampleHeatmap1.png
plotHeatmap -m matrix1_test_body.gz -out test_body_Heatmap.png
plotProfile -m matrix1_test_body.gz -out test_body_Profile.png
ngsplot也是可以的吸祟。
第7步瑟慈，peaks注釋
統(tǒng)計(jì)peak在promoter，exon屋匕，intron和intergenic區(qū)域的分布
