学徒跟着B站ATAC-seq视频5天完成流程 / 四六文摘

前面我已经学习完了jimmy老师的B站单细胞水平并且提交笔记：处理单细胞? Bioconductor就够用了，希望对大家有帮助哈！

最近刷视频看到了b站jimmy老师又更新了ATAC-seq系列教学指引，赶紧花几天时间follow了一遍！

而且把我自己学习笔记分享给大家，视频的话，文末的阅读原文直达免费学习哈！虽然视频录制是两年前，但是丝毫不影响学习体验！

前提

超大服务器+conda安装软件
fq+trim+bowtie与其他类似，这里省略
比对后需要过滤，需要摸索
mac2找peaks
计算一些指标
deeptools可视化
peaks注释

ATAC-seq

背景知识

真核生物的核DNA并不是裸露的，而是与组蛋白结合形成染色体的基本结构单位核小体，核小体再经逐步的压缩折叠最终形成染色体高级结构（如人的DNA链完整展开约2m长，经过这样的折叠就变成了纳米级至微米级的染色质结构而可以储存在小小的细胞核）。而DNA的复制转录是需要将DNA的紧密结构打开，从而允许一些调控因子结合（转录因子或其他调控因子）。这部分打开的染色质，就叫开放染色质，打开的染色质允许其他调控因子结合的特性称为染色质的可及性（chromatin accessibility）。因此，认为染色质的可及性与转录调控密切相关。ATAC则是原理是通过转座酶Tn5容易结合在开放染色质的特性，然后对Tn5酶捕获到的DNA序列进行测序。来源自简书第1篇：ATAC-seq的背景介绍以及与ChIP-Seq的异同

优势

周期快
样本小
不依赖抗体

实验设计

实战流程

12天入门生物信息学课程
注意一些黑名单（微卫星序列，重复序列)，去除掉不要当做peaks
差异peaks

1.数据下载

GSE与SRA
通过SRP055881 下载原始数据，获得sraruntable与accession list，找到样本对应的信息（例如样本名，分组等)

通过ascp下载，因为样本量太大，这里只做一部分（4个）

制作下载名单

cat >srrlist.txt
/sra/sra-instant/reads/ByRun/sra/SRR/SRR292/SRR2927015/SRR2927015.sra
/sra/sra-instant/reads/ByRun/sra/SRR/SRR292/SRR2927016/SRR2927016.sra
/sra/sra-instant/reads/ByRun/sra/SRR/SRR354/SRR3545580/SRR3545580.sra
/sra/sra-instant/reads/ByRun/sra/SRR/SRR292/SRR2927018/SRR2927018.sra

然后使用下面的命令高速下载（需要注意的是ascp有时候会失效）

nohup /home/xlfang/.aspera/connect/bin/ascp -T -l 200M -P 33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh -T --mode recv --host ftp-private.ncbi.nlm.nih.gov --user anonftp --file-list ./srrlist.txt ./ & ##3803

奇怪的是SRR292开头的数据没有对应的下载链接，SRR354有

先用ascp下载一个

nohup /home/xlfang/.aspera/connect/bin/ascp -T -P 33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh -k 1 -l 200m anonftp@ftp-private.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/sra/SRR/SRR354/SRR3545580/SRR3545580.sra ./ &

还是一样的报错

用ebi源试一下

cat >fq.txt
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR292/005/SRR2927015/SRR2927015_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR292/005/SRR2927015/SRR2927015_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR292/006/SRR2927016/SRR2927016_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR292/006/SRR2927016/SRR2927016_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR292/008/SRR2927018/SRR2927018_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR292/008/SRR2927018/SRR2927018_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR354/000/SRR3545580/SRR3545580_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR354/000/SRR3545580/SRR3545580_2.fastq.gz

cat >fq.command
cat fq.txt|while read id
do ( ascp -QT -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@$id ./ )
done

一样的问题，难道是数据库迁移？

真的是很尴尬啊！

4.10更新

ASCP又可以用了。。重新下载

ls -lh *fastq.gz|cut -d " " -f 5-
2.5G Apr 10 11:15 SRR2927015_1.fastq.gz
2.4G Apr 10 11:20 SRR2927015_2.fastq.gz
3.2G Apr 10 11:26 SRR2927016_1.fastq.gz
3.5G Apr 10 11:32 SRR2927016_2.fastq.gz
1.1G Apr 10 11:34 SRR2927018_1.fastq.gz
1.2G Apr 10 11:36 SRR2927018_2.fastq.gz
4.1G Apr 10 11:44 SRR3545580_1.fastq.gz
4.6G Apr 10 11:52 SRR3545580_2.fastq.gz

2.质检 trim

下载完成后对数据进行质检 trim 构建小鼠索引

###conda activate rna
###fastqc
nohup fastqc -t 2 -o ./ ./*.fastq.gz &
###multiqc
multiqc ./*zip -o ./
###trim
for i in `ls *_1.fastq.gz`
do
i=${i/_1.fastq.gz/}
echo "trim_galore --phred33 -q 25 --length 35 --stringency 4 -e 0.1 --fastqc --paired -o ../2.clean ${i}_1.fastq.gz ${i}_2.fastq.gz"
done > trim_galore.command
cat trim_galore.command

3.bowtie2比对

构建索引

质检完成后就需要用bowtie2去比对，第一步需要构建索引。官网有构建好的常见物种的索引（人，小鼠，斑马鱼等）就可以直接下载用，如果是非模式生物就需要自己构建（不推荐）

###下载小鼠的mm10基因索引
nohup wget -4 -q ftp://ftp.ccb.jhu.edu/pub/data/bowtie2_indexes/mm10.zip &
ls -lh|cut -d " " -f 5-
#3.2G Apr 11 03:14 mm10.zip 大小为3.2G
unzip mm10.zip ##解压

解压完成后有以下7个文件

比对

决定拆成3部分来做

第一部分 bowtie2进行比对

cd 2.clean
## 相对目录需要理解
index=/zju/phf5a/data/bwindex/mm10
## 一定要搞清楚自己的bowtie2软件安装在哪里，以及自己的索引文件在什么地方！！！
#conda activate rna
for i in `ls *_1_val_1.fq.gz`
do
i=${i/_1_val_1.fq.gz/}
echo "bowtie2 -p 2 --very-sensitive -X 2000 -x ${index} -1 ${i}_1_val_1.fq.gz -2 ${i}_2_val_2.fq.gz |samtools sort -O bam -T ${i} -@ 2 -o - > ${i}.raw.bam"
done > mapping.command

这里一直报错。。。(ERR): bowtie2-align exited with value 141

搞得我满头雾水？？应该不会有错的呀。

后来才发现原来是bowtie2版本问题，升级一下版本，就不会报错了（https://www.biostars.org/p/229653/）；

另外samtools排序需要 -T 参数加前缀，不然也会报错

8个样本居然比对了24小时，看来确实是very sensitive了。。（其实主要是sam文件转为bam文件的时候非常慢）排序好的bam文件如下

ls -lh *raw.bam|cut -d " " -f 5-
3.7G Apr 11 22:55 SRR2927015.raw.bam
4.6G Apr 12 07:19 SRR2927016.raw.bam
1.7G Apr 12 09:50 SRR2927018.raw.bam
5.4G Apr 12 19:23 SRR3545580.raw.bam

第二部分 sambamba去重

对ATAC-seq数据来说，大部分教程都推荐去除PCR重复，所以我们加上这个步骤：

conda activate rna
for i in `ls *.raw.bam`
do
i=${i/.raw.bam/}
echo " samtools index ${i}.raw.bam | bedtools bamtobed -i ${i}.raw.bam > ${i}.raw.bed|samtools flagstat ${i}.raw.bam > ${i}.raw.stat|sambamba markdup --overflow-list-size 600000 --tmpdir='./' -r ${i}.raw.bam ${i}.rmdup.bam |samtools index ${i}.rmdup.bam"
done > rmdup.command

###
## ref:https://www.biostars.org/p/170294/
## Calculate %mtDNA:
mtReads=$(samtools idxstats  ${i}.rmdup.bam | grep 'chrM' | cut -f 3)
totalReads=$(samtools idxstats  ${i}.rmdup.bam | awk '{SUM += $3} END {print SUM}')
echo '==> mtDNA Content:' $(bc <<< "scale=2;100*$mtReads/$totalReads")'%'
echo "samtools flagstat  ${i}.rmdup.bam > ${i}.rmdup.stat|samtools view  -h  -f 2 -q 30 ${i}.rmdup.bam |grep -v chrM |samtools sort  -O bam  -@ 2 -o - > ${i}.last.bam|samtools index  ${i}.last.bam|samtools flagstat  ${i}.last.bam > ${i}.last.stat|bedtools bamtobed -i ${i}.last.bam  > ${i}.bed"

8个样本大概运行了2个小时，结果如下，可以看到去重后的bam文件比原bam文件小得多

打开一个flagstat文件看看统计

cat SRR3545580.raw.stat
# 105342880 + 0 in total (QC-passed reads + QC-failed reads)
# 0 + 0 secondary
# 0 + 0 supplementary
# 0 + 0 duplicates
# 103697403 + 0 mapped (98.44%:-nan%)
# 105342880 + 0 paired in sequencing
# 52671440 + 0 read1
# 52671440 + 0 read2
# 100212070 + 0 properly paired (95.13%:-nan%)
# 103230528 + 0 with itself and mate mapped
# 466875 + 0 singletons (0.44%:-nan%)
# 477956 + 0 with mate mapped to a different chr
# 18927 + 0 with mate mapped to a different chr (mapQ>=5)

QC pass的reads的数量为105342880，未通过QC的reads数量为0，意味着一共有105342880条reads；
0条reads比对到第二个地方
0条reads只比对上部分
重复reads的数量，QC pass和failed
比对到参考基因组上的reads数量；
paired reads数据数量；
read1的数量；
read2 的数量；
正确地匹配到参考序列的reads数量；
一对reads都比对到了参考序列上的数量，但是并不一定比对到同一条染色体上；
一对reads中只有一条与参考序列相匹配的数量；
一对reads比对到不同染色体的数量；
一对reads比对到不同染色体的且比对质量值大于5的数量。

见[https://www.biostars.org/p/149883/]

第三部分获取只比对到常染色体上的高质量序列

其实就是去除在chrM染色体的reads，代码如下：

for i in `ls *.raw.bam`
do
i=${i/.raw.bam/}
mtReads=$(samtools idxstats  ${i}.rmdup.bam | grep 'chrM' | cut -f 3)
totalReads=$(samtools idxstats  ${i}.rmdup.bam | awk '{SUM += $3} END {print SUM}')
echo '==> mtDNA Content:' $(bc <<< "scale=2;100*$mtReads/$totalReads")'%'
echo " samtools flagstat  ${i}.rmdup.bam > ${i}.rmdup.stat|samtools view  -h  -f 2 -q 30 ${i}.rmdup.bam   |grep -v chrM |samtools sort -O bam  -T ${i} -@ 2 -o - > ${i}.last.bam|samtools index ${i}.last.bam|samtools flagstat ${i}.last.bam > ${i}.last.stat"
done > quality.command

可以看到线粒体DNA含量确实很高，差不多占%20。需要去除掉，不然对结果造成线粒体污染！

前面漏了一步，转成为bed格式。。

for i in `ls *.raw.bam`
do
i=${i/.raw.bam/}
echo "bedtools bamtobed -i ${i}.last.bam> ${i}.bed "
done > bed.command

最后得到的bed文件如下

ls -lh *.bed|cut -d " " -f 5-
247M Apr 14 16:40 SRR2927015.bed
352M Apr 14 16:41 SRR2927016.bed
211M Apr 14 16:41 SRR2927018.bed
264M Apr 14 16:41 SRR3545580.bed

跟Jimmy老师的结果比起来，平均要大10M。。不知道为什么

4. MACS2找peak

安装 MACS2

同样是conda搞定一切。

这里注意MACS2只支持python3.6以上版本，注意版本环境！！

conda create -n actc python=3 ###重新创建爬虫3的环境
conda activate actc
conda install -c bioconda macs2
macs2 -h##能调用就说明没有问题

运行

同样是批量

conda activate atac
cat >macs2.command
ls *.bed | while read id ;do (macs2 callpeak -t $id -g mm --nomodel --shift -100 --extsize 200 -n ${id%%.*} --outdir ../3.macs2/) ;done
###参数含义

结果如下

ls -lh|cut -d " " -f 5-
# 644K Apr 14 17:06 SRR2927015_peaks.narrowPeak
# 704K Apr 14 17:06 SRR2927015_peaks.xls
# 441K Apr 14 17:06 SRR2927015_summits.bed
# 1.2M Apr 14 17:07 SRR2927016_peaks.narrowPeak
# 1.3M Apr 14 17:07 SRR2927016_peaks.xls
# 835K Apr 14 17:07 SRR2927016_summits.bed
# 374K Apr 14 17:07 SRR2927018_peaks.narrowPeak
# 409K Apr 14 17:07 SRR2927018_peaks.xls
# 256K Apr 14 17:07 SRR2927018_summits.bed
# 1.1M Apr 14 17:07 SRR3545580_peaks.narrowPeak
# 1.2M Apr 14 17:07 SRR3545580_peaks.xls
# 714K Apr 14 17:07 SRR3545580_summits.bed

其中的excel文件就是结果了，可以看见前几行是运行命令及信息，

后面从左到右分别是染色体，peak的起始终止位置，长度，峰值位置，高度，pvalue值记忆富集程度

后面就是对结果进行分析

5.统计indel插入长度的分布

在上一部生成的bam文件的第9列，就是插入indel的长度啦，需要将4个bam的第9列取出来

conda activate rna
ls  *.last.bam >1
ls  *.last.stat >2
paste 1 2 >config_bam
cat config_bam
#SRR2927015.last.bam SRR2927015.last.stat
#SRR2927016.last.bam SRR2927016.last.stat
#SRR2927018.last.bam SRR2927018.last.stat
#SRR3545580.last.bam SRR3545580.last.stat
###bash脚本
cat >length.sh
cat config_bam |while read id;
do
arr=($id)
sample=${arr[0]}
sample_name=${arr[1]}
samtools view $sample |awk '{print $9}'  > ${sample_name}_length.txt
done

结果如下

ls -lh *last.stat_length.txt
-rw-rw-r-- 1 xlfang xlfang 22M Apr 16 10:48 SRR2927015.last.stat_length.txt
-rw-rw-r-- 1 xlfang xlfang 32M Apr 16 10:49 SRR2927016.last.stat_length.txt
-rw-rw-r-- 1 xlfang xlfang 19M Apr 16 10:49 SRR2927018.last.stat_length.txt
-rw-rw-r-- 1 xlfang xlfang 24M Apr 16 10:49 SRR3545580.last.stat_length.txt

然后在R里面新建一个R脚本画长度直方图

indel_length_distribution.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()

制作一个config

ls *last.stat_length.txt >3
ls *last.stat_length.txt|cut -d "." -f 1-3>4
paste 3 4 >config.indel_length_distribution
cat config.indel_length_distribution
# SRR2927015.last.stat_length.txt SRR2927015.last.stat_length
# SRR2927016.last.stat_length.txt SRR2927016.last.stat_length
# SRR2927018.last.stat_length.txt SRR2927018.last.stat_length
# SRR3545580.last.stat_length.txt SRR3545580.last.stat_length

bash

cat >indel_length_distribution.sh
cat config.indel_length_distribution |while read id;
do
arr=($id)
input=${arr[0]}
output=${arr[1]}
Rscript indel_length_distribution.R $input $output
done

但是我还是用不惯linux的R，还是用windows吧。

a=read.delim("SRR2927015.last.stat_length.txt")
a=abs(as.numeric(a[,1]))
png(file=output)
hist(a,
  main="SRR2927015_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()

其中一个结果如下面

可以很直观的看到，插入片段在150-300之间是counts数最大，随着长度增大，counts越不断缩减

6.FRiP值的计算

FRiP 是 The fraction of reads in called peak regions的缩写，指的是mapping到peak里的reads 除以总mapped reads

conda activate rna
for i in `ls *.narrowPeak`
do
i=${i/_peaks.narrowPeak/}
bed=../2.clean/${i}.bed
Reads=$(bedtools intersect -a ${bed} -b ${i}_peaks.narrowPeak |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 > frip.command

结果如下

cat frip.command
148214 5105856
==> FRiP value: 2.90%
320405 7292136
==> FRiP value: 4.39%
90854 4399724
==> FRiP value: 2.06%
258981 5466470
==> FRiP value: 4.73%

7.样本间peaks值交集

比较粗狂的办法，这里我默认了4个样本都是重复样本（其实不是。。）

library(ChIPseeker)
library(ChIPpeakAnno)
setwd("E://desktop/sept/ATAC-seq_practice/find_peaks_overlaping/")
list.files('./',"*.narrowPeak")
tmp = lapply(list.files('./',"*.narrowPeak"),function(x){
return(readPeakFile(file.path('./', x)))
})
tmp
ol <- findOverlapsOfPeaks(tmp[[1]],tmp[[2]],tmp[[3]],tmp[[4]])
png('overlapVenn.png')
makeVennDiagram(ol)
dev.off()

但这种方法比较粗犷，还是使用IDR吧

##安装IDR
conda activate atac
conda isntall -y idr
ls -lh *narrowPeak|cut -d " " -f 9|tr -s "\n" " "
#SRR2927015_peaks.narrowPeak SRR2927016_peaks.narrowPeak SRR2927018_peaks.narrowPeak SRR3545580_peaks.narrowPeak
###2个样本的overlap peaks IDR只能对于2个重复样本进行比较，所以分为2组
idr --samples SRR2927015_peaks.narrowPeak SRR2927016_peaks.narrowPeak \
--input-file-type narrowPeak \
--rank p.value \
--output-file group1-idr \
--plot \
--log-output-file group1.idr.log

idr --samples SRR2927018_peaks.narrowPeak SRR3545580_peaks.narrowPeak \
--input-file-type narrowPeak \
--rank p.value \
--output-file group2-idr \
--plot \
--log-output-file group2.idr.log

输出结果如下

ls -lh group*|cut -d " " -f 5-
# 740K Apr 19 22:00 group1-idr
# 205 Apr 19 22:00 group1.idr.log
# 341K Apr 19 22:00 group1-idr.png
# 315K Apr 19 22:00 group2-idr
# 202 Apr 19 22:00 group2.idr.log
# 218K Apr 19 22:00 group2-idr.png

cat group1.idr.log
# Initial parameter values: [0.10 1.00 0.20 0.50]
# Final parameter values: [0.59 1.05 0.62 0.79]
# Number of reported peaks - 5869/5869 (100.0%)
# Number of peaks passing IDR cutoff of 0.05 - 1571/5869 (26.8%)

cat group2.idr.log
# Initial parameter values: [0.10 1.00 0.20 0.50]
# Final parameter values: [0.35 1.06 0.47 0.50]
# Number of reported peaks - 2505/2505 (100.0%)
# Number of peaks passing IDR cutoff of 0.05 - 38/2505 (1.5%)

这里主要是idr.log与idr.png
idr.log指通过IDR的0.05的共有peaks数，可以看到比之前用narrowPeak结果直接取交集少了很多peaks

而给的图标则表示：

左上：Rep1 peak ranks vs Rep2 peak ranks, 没有通过特定IDR阈值的peaks显示为红色。
右上：Rep1 log10 peak scores vs Rep2 log10 peak scores，没有通过特定IDR阈值的peaks显示为红色。
下面两个图：Peak rank vs IDR scores，箱线图展示了IDR值的分布，默认情况下，IDR值的阈值为-1E-6。

8.deeptools可视化

介绍主要看Jimmy老师的推文如何使用deeptools处理BAM数据

###安装deeptools
conda activate atac ###激活小环境
conda install -y deeptools ###安装
cd ../../zju/phf5a/atac/
mkdir deeptools
cd deeptools/
ln -s ../2.clean/*.last.bam ./
###构建索引
ls *.bam |xargs -i samtools index {}
###将最终的bam转换成bigwig文件
ls *last.bam |while read id;do
nohup bamCoverage -b ${id} -p 5 --normalizeUsing RPKM -o ${id%%.*}.last.bw &
done
###将去重的bam转换成bigwig文件
ls *rmdup.bam |while read id;do
nohup bamCoverage --normalizeUsing CPM -b $id -o ${id%%.*}.rm.bw &
done

下载mm10的Refgene文件

curl 'http://genome.ucsc.edu/cgi-bin/hgTables?hgsid=646311755_P0RcOBvAQnWZSzQz2fQfBiPPSBen&boolshad.hgta_printCustomTrackHeaders=0&hgta_ctName=tb_ncbiRefSeq&hgta_ctDesc=table+browser+query+on+ncbiRefSeq&hgta_ctVis=pack&hgta_ctUrl=&fbQual=whole&fbUpBases=200&fbExonBases=0&fbIntronBases=0&fbDownBases=200&hgta_doGetBed=get+BED' >
###16m大小下载了18min。。

查看TSS附件信号强度

cut -f 1-3 mm10.refseq.bed > ref.bed
computeMatrix reference-point  --referencePoint TSS  -p 15  \
-b 10000 -a 10000 \
-R ./ref.bed  \
-S ./*.bw  \
--skipZeros  -o matrix1_test_TSS.gz  \
--outFileSortedRegions regions1_test_genes.bed

遇到报错。。应该是python版本问题吧。