ATAC-seq技术与数据处理流程简述

ATAC-seq技术简介

ATAC-seq全称Assay for Transposase-Accessible Chromatin with high-throughput Sequencing，即利用高活性转座酶Tn5研究染色质可及性的高通量测序技术。

染色质的结构

真核生物的核DNA并不是裸露的，而是有蛋白质即组蛋白与之相结合的。DNA一圈一圈地缠绕在组蛋白上，形成串珠式的结构，而这样的结构还能够进一步折叠、浓缩，并在其他架构蛋白的辅助下形成染色体。

这样做的意义在于，能够将超长的DNA链，折叠成很小很小的结构，从而塞进小小的细胞核里。比如说，人类的DNA链如果完整展开，大概会有2米那么长。而经过这样的折叠，就变成了纳米至微米级染色体/质的结构了。
当基因需要转录时，染色质就需要打开这个基因所在的高级结构，将需要进行复制和转录的DNA裸露，这一过程主要是在组蛋白的修饰（尤其是乙酰化）的作用下完成的。

寻找染色质的开放区域

传统的实验方法主要是借助MNase-seq和DNase I hypersensitivity assay。
这两个实验的主要思路是一致的：染色质变得开放，就意味着DNA和组蛋白的浓聚程度降低，就会有一部分DNA暴露出来。而一旦失去了蛋白质的保护，这部分DNA就可以被DNA酶（MNase或DNase I）所切割。然后，我们再把切割完的DNA拿来测序，和已知的全基因组序列相比较，就能发现被切掉的是哪些地方，没有被切掉的地方又在哪里，从而获知开放的染色质区域。
不过，这两个方法都有明显的缺陷，即耗时费力与重复性差。具体就不展开讲，亲手做过的就知道这两个方法有多么讨厌。
2013年，美国Stanford大学的William Greenleaf教授研发了一种全新的方法，利用DNA转座酶结合高通量测序技术，来研究染色体的，即ATAC-seq。

DNA转座，是一种把DNA序列从染色体的一个区域搬运到另外一个区域的现象，由DNA转座酶来实现。这种转座插入DNA，也是需要插入位点的染色质是开放的，否则就会被一大坨高级结构给卡住。
那么，如上图a，我们只要人为地，将携带已知DNA序列标签的转座复合物（即带着红色蓝色测序标签的转座酶Tn5），加入到细胞核中，再利用已知序列的标签进行PCR后测序，就知道哪些区域是开放染色质了。而这也就是ATAC-seq的原理。
ATAC-seq出来的结果，和传统方法出来的结果具有很强的一致性，同时也和基于组蛋白修饰marker的ChIP-seq有较高的吻合程度。也就是说，ATAC-seq中的peak，往往是启动子、增强子序列，以及一些反式调控因子结合的位点。
相比起来，ATAC-seq的重复性，比MNase-seq和DNase-seq的更强，操作起来也更加简单，而且只需要很少的细胞/组织量，同时出来的信号更加漂亮。目前已经是研究染色质开放性首选的技术方法。

ATAC-seq实验操作

下面，我们简要概述了ATAC-seq实验步骤，该步骤可分为两个主要部分：细胞制备和转座。
细胞制备
第一步需要收获细胞，由于用于ATAC-seq分析的细胞数量对于转座反应和生成的DNA片段的大小分布至关重要，因此对细胞进行技术非常重要。此外，细胞应完好无损，并且是均匀的单细胞悬浮液。收获后，用非离子型去污剂裂解细胞以产生纯细胞核。
转座反应
然后将所得的染色质片段化，并同时用Tn5 转座酶进行测序接头标记，以生成ATAC-seq文库。纯化后，可使用barcode引物通过PCR扩增文库。所得文库随后通过qPCR或NGS进行分析。

ATAC-seq数据处理

质控

fastqc -o 1_fastqc sample_R1.fq.gz sample_R2.fq.gz -t 20
trim_galore --phred33 --fastqc --paired --clip_R1 15 --three_prime_clip_R1 10 --clip_R2 15 --three_prime_clip_R2 10 sample_R1.fq.gz sample_R2.fq.gz -o trim_path

#首先利用fastqc是对原始数据进行质量检查，数据质量结果以网页形式输出，具体解读参考https://editor.csdn.net/md/?articleId=127923307
#如果fastqc结果显示数据质量不好，需要用trim_galore进行3端和5端的剪切和去接头操作，具体剪切多少碱基根据fastqc结果视情况而定

比对

bowtie2-build -f hg38.fa --threads 24 index_path/hg38.fa
bowtie2 -p 8 --very-sensitive --end-to-end --no-mixed --no-discordant --phred33 -I 10 -X 700  -x index_path/hg38.fa -1 sample_R1.fq.gz -2 sample_R2.fq.gz |  samtools sort -O bam -@ 8 -o -> bowtie2_result_path/sample.bam
samtools index bowtie2_result_path/sample.bam
samtools flagstat bowtie2_result_path/sample.bam > bowtie2_result_path/sample.stat

#bowtie2建立索引进行比对
#samtools对bam文件建立索引和比对质量文件

bam文件过滤

去除PCR重复
gatk MarkDuplicates -MAX_FILE_HANDLES 1000 -I bowtie2_result_path/sample.bam -O rmDuplicate_path/sample.rmdup.bam --REMOVE_DUPLICATES true -M rmDuplicate_path/sample.rdup.bam.metrics.txt
samtools index  rmDuplicate_path/sample.rmdup.bam
samtools flagstat  rmDuplicate_path/sample.rmdup.bam >  rmDuplicate_path/sample.rmdup.stat
	
去除线粒体DNA和比对质量不好的reads
samtools view -h -f 2 -q 30 rmDuplicate_path/sample.rmdup.bam  |grep -v chrM | samtools sort -O bam -@ 8 -o ->clean_path/sample.clean.bam
samtools index clean_path/sample.clean.bam
samtools flagstat clean_path/sample.clean.bam > clean_path/sample.clean.stat

peak calling

macs2 callpeak -t clean_path/sample.clean.bam -g hs --nomodel --shift -100 --extsize 200 -n sample -q 0.01 --outdir peakCalling_path

IGV可视化

genomesize=faCount.py hg38.fa
bamCoverage --numberOfProcessors 8 --binSize 10 --normalizeUsing RPKM --effectiveGenomeSize ${genomesize} --bam clean_path/sample.clean.bam -o bw_path/sample.clean.bw

faCount.py:

import sys
aList=[]
fa_file = sys.argv[1]
with open(fa_file,'r') as f:
    for line in f:
        line = line.strip()
        line = line.upper()
        if not line.startswith(">"):
            baseA = line.count("A")
            baseT = line.count("T")
            baseC = line.count("C")
            baseG = line.count("G")
            aList.extend([baseA, baseT, baseC, baseG])
            # print(aList)
    print(sum(aList))
#bw文件即可在IGV软件中可视化

重复样本分析

样本相关性分析

multiBigwigSummary bins -b sample_replicate1.clean.bw sample_replicate2.clean.bw -o correlation_path/sample_number.of.bins.npz
plotCorrelation -in correlation_path/sample_number.of.bins.npz \
		--corMethod pearson \
		--skipZeros \
		--plotTitle "Pearson Correlation of Average Scores Per Transcript of sample" \
		--plotFileFormat png \
		--whatToPlot scatterplot \
		--log1p \
		-o correlation_path/sample_scatterplot_PearsonCorr_bigwigScores.png \
		--outFileCorMatrix correlation_path/sample_PearsonCorr_bigwigScores.tab 

peak重复性分析

for s in sample_replicate1 sample_replicate2
do
	sort -k8,8nr peakCalling_path/${s}_peaks.narrowPeak > peakCalling_path/${s}_peaks.sorted.narrowPeak	
done
idr --samples peakCalling_path/sample_replicate1_peaks.sorted.narrowPeak peakCalling_path/sample_replicate2_peaks.sorted.narrowPeak --input-file-type narrowPeak --rank p.value --output-file idr_path/sample_idr --plot --log-output-file idr_path/sample.idr.log