ATAC-seq分析干货-2

2020-02-22 本文已影响0人生信阿拉丁

作者 | Arno
审稿 | 童蒙
编辑 | amethyst

上一期我们介绍了ATAC-seq相关的背景知识。ATAC-Seq能帮助我们从全基因组范围内推测可能的开放染色质位点，其分析从本质上来说和ChIP-seq区别不大，核心都是peak-calling。如果现在你还不太了解基本分析流程，那快快跟随小编继续学起来吧。

分析步骤

前期我们通过Trimmomatic软件对原始下机数据进行了数据清洗，主要包括：

去除下机数据中的接头污染
去除低质量的reads

从而获得clean reads数据。那接下来就可以将clean数据跟参考基因组进行比对，进行后续的可视化以及peak calling相关分析了。

1.数据比对（Bowtie2,picard-tools）

1.1 比对

#Bowtie2 Mapping
bowtie2 -N 1 -p 4 -X 2000 -q -x ref -1 clean_R1.fq.gz -2 clean_R2.fq.gz 2>alignment.log 
samtools view -bS sanple.sam > sample.bam
#Sort BAM
samtools sort -@ 4 -O BAM -o sample.srt.bam sample.bam
#Stat BAM
samtools flagstat sample.bam >sample.prime.stat

1.2 过滤（去dup，去线粒体，低质量，未比对上的序列）
比对完之后，需要对比对得到的bam文件进行过滤处理。主要考虑过滤以下几个方面：

去除PCR duplicate reads，因为影响后续call peak
线粒体DNA和叶绿体DNA序列。由于线粒体和叶绿体DNA是裸露的，可以被Tn5酶识别切割，因此数据中会有线粒体和叶绿体的污染，比对完之后需要去掉（如果是植物，还需要去掉叶绿体）
低质量以及未比对上的reads

#Mark duplicates（picard-tools）
java -Xmx25g -XX:ParallelGCThreads=4 -jar MarkDuplicates.jar INPUT=sample.srt.bam OUTPUT=sample.srtdup.bam METRICS_FILE=sample.srtdup.qc VALIDATION_STRINGENCY=LENIENT ASSUME_SORTED=true \ REMOVE_DUPLICATES=false
#Remove low MAPQ umap Duplicates chrMT... 
samtools view -h sample.srtdup.bam |grep -v chrM|samtools view -bS -q 30 -F 1804 -f 0x2 - >sample.filt.bam 
## 去除线粒体或者叶绿体，可以自己写脚本或者命令行实现（grep -v chrM）
# Sort BAM 
samtools sort --threads 4 sample.filt.bam  -o sample.final.bam
samtools index sample.final.bam
samtools flagstat sample.final.bam >sample.final.flagstat.qc

1.3 插入片段长度统计（picard-tools）,同时生成call peak的bed文件）
插入片段长度是重要的评估实验好坏的指标，统计出的插入片段长度应该符合实验预期的长度。可以用picard-tools工具的CollectInsertSizeMetrics功能进行统计。

为了更方便的进行peak calling，可以根据插入片段起始终止位置整理成bed文件作为输入（需要注意的是bed文件的染色体起始位置是0起始的，第一个碱基的位置标记为0，而终止位置是从1开始的）。

根据ATAC-seq的实验原理，call peak 时，可以进行Tn5边界调整，将bed文件上下游各移动75bp，从而直接进行calling，而不用再进行reads延伸^[1]，具体call peak操作可以往下看。

#CollectInsertSizeMetrics，使用picard工具统计插入片段长度，获得的插入片段统计文件
java -Xmx25g -XX:ParallelGCThreads=4 -jar CollectInsertSizeMetrics.jar I=sample.final.bam O=sample.InsertSizes.txt H=sampe.InsertSizes.pic METRIC_ACCUMULATION_LEVEL=ALL_READS  
#BAM to BED（可以使用bedtools的bamToBed工具）
bamToBed -i sample.final.bam >sample.temp.bed
#produce insert BED 具体生成bed文件sample.InsertSizes.bed的方式可参考秀丽隐杆菌ATAC-seq分析[1]文献的具体补充细节
# Convert back to a bam, 1bp
samtools view -H sample.final.bam | grep '^@SQ' | awk '{split($2, SN, ":"); split($3, LN, ":"); print SN[2]"\t"LN[2]}' > genome.size
bedToBam -i sample.InsertSizes.bed -g genome.size > sample.temp.bam
samtools sort sample.temp.bam -o sample.Shifted_InsertSizes.bam
samtools index sample.Shifted_InsertSizes.bam
###为后续可直接使用bed文件进行call peak，可以先把bed文件处理为上下游移动75bp sample.Shifted_75bp.bed

2.比对结果可视化(deeptools,samtools)

得到比对结果之后，我们就可以图形化的展示peak的富集情况，peak reads在染色体的分布情况等等。另外，通常bam文件是比较大的，所以为了更加方便的使用IGV可视化的展示reads情况，通常会把bam文件转为bigwig文件形式。
2.1 bam转bigwig文件（deeptools,bamCoverage）
bam的格式转换可以使用deeptools的bamCoverage工具，具体使用代码参考如下：

## 使用deeptools的bamCoverage工具，方便的将bam转为bw格式。
bamCoverage --bam sample.final.bam --outFileName sample.bw --binSize 10 --numberOfProcessors 4 --normalizeUsingRPKM --extendReads
### -bs/--binSize：可以对基因组进行等宽分箱（划bin）,统计每个bin的reads，进行统计
### --region/-r CHR:START:END: 也可以使用此参数选取某个区域统计，默认为全基因组范围内
### --normalizeUsingRPKM：对每条染色体进行标准化的方式，如果不需要标准化使用--ignoreForNormalization
### -e/--extendReads：使用此参数可以拓宽原来的reads长度

2.2 覆盖度(samtools）
使用samtools的mpileup工具统计base\reads在各个染色体的覆盖度情况，并进行可视化展示。

#单碱基以及染色体的覆盖度统计
samtools mpileup sample.final.bam # 对得到的结果进行简单处理可以直接用R进行碱基或者染色体覆盖情况的可视化
#染色体reads分布统计
samtools view -F 0x0010 sample.final.bam |awk '{print $3"\t"int($4/1000)}'  ## 0x0010代表read reverse strand
samtools view -f 0x0010 sample.final.bam |awk '{print $3"\t"int($4/1000)}' ## 可进行划窗统计
##整理完以上统计文件后，便可以直接使用R进行简单的绘图展示了

3.TSS、TES区域可视化

计算TSS或者TES附近的信号强度，可以使用deeptools^[2]的computeMatrix工具进行。computeMatrix具有两个模式:scale-region和reference-point。

scale-region用来计算某一个区域内信号分布，比如如下代码设置的TSS和TES区域以及上下游各2K的区域内的信号分布计算；
reference-point可以计算相对于某一个点的信号分布情况。

computeMatrix工具计算信号分布情况结果，可以使用plotProfile以折线图的方式对覆盖区域信号分布进行可视化，也可以使用plotHeatmap以热图的方式对覆盖区域信号分布进行可视化。

##BigWig to matrix
computeMatrix scale-regions --regionsFileName Gene.bed --scoreFileName sample.bw --outFileNameMatrix outFileNameMatrix --regionBodyLength 6000 --beforeRegionStartLength 3000 --afterRegionStartLength 3000 --startLabel TSS --endLabel TES --skipZeros --numberOfProcessors 4 --outFileName plotMatrix.gz
##plot line plotProfile
plotProfile -m plotMatrix.gz -out sample.lineProfile.pdf --samplesLabel sample_name --plotType lines --startLabel TSS --endLabel TES --plotFileFormat pdf --plotHeight 10 --plotWidth 14 --yAxisLabel "Reads Density" --perGroups
##plot heatmap plotHeatmap
plotHeatmap -m plotMatrix.gz -out sample.heatmapProfile.pdf --samplesLabel sample_name --plotType heatmap --startLabel TSS --endLabel TES --plotFileFormat pdf

TSS\TES区域可视化结果示例图

4.链互相关性SCC分析(phantompeakqualtools)

SCC（Strand Cross Correlation analysis）是可用于评估ATAC/Chip实验质量好坏的一个重要指标（转录因子富集），其原理基于peak的reads富集分布规律^[3]。

第一：peak位点附近的正负链上reads分布相同；

第二：reads分布的中心点和peak的中心点存在偏移，如果将reads的位置移动一定的距离之后，正负链的中心重合，上下成对称分布。

用泊松相关系数来分析正负链测序深度的相关性，当正负链的中心点重合时，相关系数最高，可以有效衡量偏移过程。由此，可以得到偏移距离和相关系数之间的对应关系。因此，一方面，计算SCC有助于了解正负链间的相关性系数，另一方面可以检测Reads富集的程度来富集实验是否成功。

Chip-seq数据检测蛋白结合位点

典型的链交叉相关图会产生两个峰：

一个富集峰与主要的插入片段长度相关（高相关性），为peak峰。
另一个与测序read长度相关，被称为"phantom "peak。

一个高质量的转录因子富集数据，peak对应的峰最高，phantom peak对应的峰较矮，如果两种峰都能够观测到，而phantom peak最高，则说明抗体还是富集到了部分序列，但是背景噪声太高了，不利于后续分析；而如果观测不到 peak峰，则说明富集实验是失败的。

Examples of Cross-Correlation Plots

为了更精准的进行判断，研究者提出两个量化指标：

Normalized strand cross-correlation coefficent (NSC)是最大交叉相关值除以背景交叉相关的比率，NSC值越大表明富集效果越好，NSC值低于1.1，表明较弱的富集，小于1表示无富集。NSC值稍微低于1.05，表明有较低的信噪比或很少的峰；
Relative strand cross-correlation coefficient (RSC)，片段长度相关值减去背景相关值除以phantom-peak相关值减去背景相关值。RSC的最小值可能是0，表示无信号；富集好的实验RSC值大于1；低于1表示质量低。


##将bam文件转为tagAlign文件 $2>16的即为负链，小于等于16的为正义链
samtools view -F 0x0204 -o - sample.final.bam|awk 'BEGIN{OFS="\t"}{if (and($2,16) > 0) {print $3,($4-1),($4-1+length($10)),"N","1000","-"} else {print $3,($4-1),($4-1+length($10)),"N","1000","+"} }' | gzip -c >sample.tagAlign.gz
### samtools view Reads.bam | gawk '(and(16, $2))' > forwardStrandReads.sam
##Strand cross-correlation analysis
Rscript ./SPP_1.10.1/run_spp.R -c=sample.tagAlign.gz -s=0:5:1000 -p=4 -tmpdir=./tmp -odir=./sample -savp=./sample.SCC.pdf -out sample.SCC_report.txt -ccs=sample.SCC_data.txt -samname=sample_name -rf

5.Peak Calling

Model-based Analysis of ChIP-Seq (MACS2) 是用来检测DNA片段富集的软件，尽管当时是针对Chip-seq开发的，但是其非常好的适用于ATAC-seq。现在ATAC-seq主流的call peak软件依然是MACS2。

蛋白或转录因子结合位点附近会有reads富集，MACS2根据一定的统计模型，扫描全基因组，构建双峰模型，结合对照（background）来检测富集片段中真正的结合位点。MACS2可以鉴定Chip-seq中两种主要类型的富集模式：

组蛋白修饰的broad peaks或者broad domains
转录因子结合的narrow peaks。

而对于ATAC-seq而言，捕获的是染色质开放区域，可类似于Chip-seq中转录因子的检测模式。

MACS2 Peak Calling原理

由于ATAC-seq检测固定的染色质开放区域，所以MACS2检测ATAC peak时，可以使用--nomode1参数，不需要MACS去构建模型，同时需要设置--extsize参数和--shift参数。

--extsize参数和--shift参数是非常重要可以明显左右检测结果的参数，可以根据实验方法以及目的，设置相应的值。具体的脚本可以参考C. elegans的ATAC-seq分析^[1]或者文档ATACSeq Pipeline，相关的参数和结果文件说明可以参考MACS2作者的github，MACS2。

#MACS2 Peak Calling
macs2 callpeak -t sample.Shifted_75bp.bed -f BED -n sample_name -g $genome_size -q 0.05 --nomodel --extsize 150 --keep-dup all --call-summits --outdir ./
## --nomodel：此参数表明不需要MACS去构建模型，因此需要设置--extsize参数
## --extsize：MACS使用这个参数将read以5'-> 3'拓展，如染色体开放区域为150bp，可以设置--extsize 150
## --shift：read绝对偏移量，会先于--extsize前对read进行延伸，由于前期计算插入片段长度时，已进行reads延伸，所以此处可不再进行处理
## --keep-dup：dup reads的处理，默认为保留1条，由于前期已过滤掉dup reads，所以此处可直接全部保留

考虑到ATAC-seq是通过Tn5酶结合的染色质开放区域，传统的使用MACS2软件（使用--extsize,--shift参数）检测ATAC peak的方式可能并不是最优的，比对上的成对的reads可能会被丢弃，从而忽略很多有效的序列信息。基于此，Harvard John M. Gaspar团队^[5]推荐使用Genrich(https://github.com/jsh58/Genrich)软件，Genrich有特定的针对ATAC的分析模式，以Tn5转座酶切割位点为中心进行拓展，而不是通过拓展reads来检测peak。Genrich对线粒体reads、PCR dup reads、多重比对reads、生物学重复都有较好的处理，感兴趣的可以尝试使用此软件分析。

6.FRiP的统计

检峰之后，可以对检测到的peak数进行基本的个数、长度等统计。FRiP(Fraction of reads in peaks)^[6]，peaks中的reads与总reads的比例，即文库中结合位点片段占背景reads的比例，可理解为“信噪比”，也是样本富集效果的评价指标，可一定程度上反应富集效果。

通常，一个典型质量好的TF富集FRiP值约5%或者更高,大部分(787 of 1052)ENCODE数据集中FRiP值在1%以上，但低于此阈值的并不能说明富集实验是失败的，如ZNF274以及人类RNA聚合酶III有非常少的结合位点，所以具体需要综合实验以及研究的对象结合FRiP进行评估。

#Peak Stat
## peak峰的个数
Peak_number=`wc -l sample_peaks.narrowPeak|cut -f1`
## 峰顶的个数
Summit_number=`cut -f1-3 sample_peaks.narrowPeak|sort -u |wc -l`
## 峰平均长度
cut -f1-3 sample_peaks.narrowPeak|awk -F '\t' 'BEGIN{sum=0}{sum+=$3-$2}END{print("%.2f\n",sum/"'$Peak_number'")}'

#Calculate FRiP
Total_reads=`samtools view -c  -F 0x0100 sample.final.bam` &&\
Peak_reads=`samtools view -c -L sample_peaks.narrowPeak -F 0x0100 sample.final.bam` &&\
echo "scale=4; 100*$Peak_reads/$Total_reads" | bc | awk '{print "FRiP(%)""\t"$1}'

经过以上分析后，ATAC-seq分析的主要结果就已经有了，后续可以基于分析得到的peak结合自己的实验课题进行各种可视化展示、peak注释、差异分析、功能富集以及Motif分析等等，筛选出跟课题相关的结合位点。现在，感兴趣的小伙伴就可以上手啦~

参考文献：

Daugherty A C, Yeo R W, Buenrostro J D, et al. Chromatin accessibility dynamics reveal novel functional enhancers in C. elegans[J]. Genome research, 2017, 27(12): 2096-2107.
https://deeptools.readthedocs.io/en/develop/index.html
Kharchenko P V, Tolstorukov M Y, Park P J. Design and analysis of ChIP-seq experiments for DNA-binding proteins[J]. Nature biotechnology, 2008, 26(12): 1351.
https://github.com/hbctraining/In-depth-NGS-Data-Analysis-Course/blob/master/sessionV/lessons
https://informatics.fas.harvard.edu/atac-seq-guidelines.html(2017).
Landt SG, Marinov GK, Kundaje A, et al. ChIP-seq guidelines and practices of the ENCODE and modENCODE consortia. Genome Res. 2012;22(9):1813–1831. doi:10.1101/gr.136184.111

关注阿拉丁，为科研实用而生