本次给学徒讲解的文章是 :The landscape of accessible chromatin in mammalian preimplantation embryos. Nature 2016 Jun 30;534(7609):652-7. PMID: 27309802
查看文章发现数据上传到了GEO，是：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE66581
在SRA数据库可以下载原始测序数据 , 从文章找到数据的ID： https://www.ncbi.nlm.nih.gov/sra?term=SRP055881 把下面的内容保存到文件，命名为 srr.list 就可以使用prefetch这个函数来下载。
linux环境及软件安装
这里首推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后需要设置镜像。
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先进行检索
conda search trim_galore
## 保证所有的软件都是安装在 wes 这个环境下面
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
代码里面提到的软件，都是根据我们对ATAC-seq搜索学习总结的。
值得一提的是自己为了ATAC-seq建立的软件环境可以很方便移植到另外一台电脑!
首先通过activate target_env要分享的环境target_env，然后输入下面的命令会在当前工作目录下生成一个environment.yml文件，
conda env export > environment.yml
小伙伴拿到environment.yml文件后，将该文件放在工作目录下，可以通过以下命令从该文件创建环境
conda env create -f environment.yml
下载作者的数据
前面提到的SRA数据库，该文章配套数据太多，我们节选部分作为练习，文件config.sra 如下：
2-cell-1 SRR2927015
2-cell-2 SRR2927016
2-cell-5 SRR3545580
2-cell-4 SRR2927018
因为conda安装好了sra-toolkit，所以prefetch函数可以直接使用
## 下载数据
# cat srr.list |while read id;do (nohup $prefetch $id -X 100G & );done
## 注意组织好自己的项目
mkdir -p ~/project/atac/
cd ~/project/atac/
mkdir {sra,raw,clean,align,peaks,motif,qc}
cd sra
## vim 或者cat命令创建 srr.list 文件, 里面保存着作为练习使用的4个数据ID
source activate atac
cat srr.list |while read id;do ( nohup prefetch $id & );done
## 默认下载目录：~/ncbi/public/sra/
ls -lh ~/ncbi/public/sra/
## 下载耗时，自行解决，学员使用现成数据：/public/project/epi/atac/sra
## 假如提前下载好了数据。
cd ~/project/atac/
ln -s /public/project/epi/atac/sra sra
总之数据如下：
-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测序数据
通常我们应该是自己的实验数据，自己找公司测序后拿到原始数据，本次讲解使用的是公共数据，所以需要把下载的sra数据转换为fq格式。
## 下面需要用循环
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]}
# 测序数据的sra转fasq
nohup $dump -A $sample -O $analysis_dir --gzip --split-3 sra/$srr.sra &
done
### 如果不只是4个文件，需要使用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数据如下:
-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
第二步，测序数据的质量控制
这个时候选择trim_galore软件进行过滤，双端测序数据的代码如下:
需要自行制作 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
质量控制前后都需要可视化，肯定是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 ./
第三步，比对
比对需要的index，看清楚物种，根据对应的软件来构建，这里直接用bowtie2进行比对和统计比对率, 需要提前下载参考基因组然后使用命令构建索引，或者直接就下载索引文件：下载小鼠参考基因组的索引和注释文件, 这里用常用的mm10
# 索引大小为3.2GB， 不建议自己下载基因组构建，可以直接下载索引文件，代码如下：
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
这些索引文件一个都不能少，而且文件名的前缀很重要，保证一致。
单端测序数据的比对代码如下：
首先可以对测试样本走流程，完善代码:
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功能，避免麻烦。
## https://www.jianshu.com/p/1e6189f641db
samtools markdup -r test.bam test.samtools.rmdup.bam
## 把报错信息在谷歌搜索后，在两个网页上找到了答案。
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要比对到同一条染色体(Proper paired) ，还有高质量的比对结果(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 {}
探索好了整个流程，就可以直接写批处理，代码如下：
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
## 相对目录需要理解
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]}
## 比对过程15分钟一个样本
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比对加入了-X 2000 参数，是最大插入片段，宽泛的插入片段范围(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
上述脚本的步骤都可以拆分运行，比如bam文件构建index或者转为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软件说明书详见：https://www.jianshu.com/p/21e8c51fca23
第5步，计算插入片段长度，FRiP值，IDR计算重复情况
非冗余非线粒体能够比对的fragment、比对率、NRF、PBC1、PBC2、peak数、无核小体区NFR、TSS富集、FRiP 、IDR重复的一致性！
名词解释：https://www.encodeproject.org/data-standards/terms/
参考：https://www.encodeproject.org/atac-seq/
看 sam文件第9列，在R里面统计绘图
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脚本就可以在批量检验bam文件进行出图。
还有NFR:https://github.com/GreenleafLab/NucleoATAC/issues/18
FRiP值的计算：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)
文章其它指标：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()
也可以使用专业软件，IDR 来进行计算出来，同时考虑peaks间的overlap，和富集倍数的一致性 。
source activate atac
# 可以用search先进行检索
conda search idr
source deactivate
## 保证所有的软件都是安装在 wes 这个环境下面
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
结果如下：
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文件转为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附件信号强度：
## 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
### 如果要批处理 ，需要学习好linux命令。
下载 bed文件：https://genome.ucsc.edu/cgi-bin/hgTables 只需要3列坐标格式文件。
查看基因body的信号强度
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注释
统计peak在promoter，exon，intron和intergenic区域的分布
