scATAC-seq入门

Posted on

Introduction

基本概念：

nucleosomes: 由八个组蛋白构成，包裹147bp的DNA序列。

Promotor：TSS附近区域

Enhancer：promotor上游直到1MB的位置。

转座酶：搬运一段DNA序列到另外一个位置，有转座酶实现插入到DNA序列，需要插入位点染色质是开放的。如果两个相邻的Tn5转座酶切割，长度40bp。文献了解：Rapid, low-input, low-bias construction of shotgun fragment libraries by high-density in vitro transposition

CTCF：由11个锌指蛋白组成，转录结合位点由锌指蛋白，CTCF与自身结合形成同源二聚体，和cohesiin一起作用，导致结合的DNA形成环状。

DNA-loop：It is currently believed that the DNA loops are formed by the “loop extrusion” mechanism, whereby the cohesin ring is actively being translocated along the DNA until it meets CTCF. CTCF has to be in a proper orientation to stop cohesin.

CTCF的结合可以被CpG methylation给破环，另一方面，CTCF结合可能为DNA甲基化的扩散设置边界。

与核小体的关系：CTCF binding sites act as nucleosome positioning anchors so that, when used to align various genomic signals, multiple flanking nucleosomes can be readily identified.所以used as a positive control for assessing if the ATAC-Seq experiment is good quality.Good ATAC-Seq data would have accessible regions both within and outside of TSS, for example, at some CTCF binding sites.

异染色质：是染色质的紧密排列形式，可以沉默基因转录。 异染色质构成端粒、中心周围区域和富含重复序列的区域。 常染色质凝缩较少，含有活性最强的转录基因。active and

基本原理

DNA被核小体（nucleosomes包裹）。当DNA被转录的的时候，DNA会被打开，核小体会松散。很多因子，比如染色质结构，核小体的位置，组蛋白修饰都会对DNA的organization和accessiblity其重要作用。这些因子也会对基因的激活和抑制起到重要作用。ATAC-seq能够检测染色质的开放区域，来作为研究基因调控机制的一种方法。当TF结合enhancer并且连接promoter区域的时候，基因转录或者停止转录。

技术历史

  1. 2013年 Nature Method ATAC-seq第一篇

    Buenrostro, Jason D., Paul G. Giresi, Lisa C. Zaba, Howard Y. Chang, and William J. Greenleaf. “Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin, DNA-binding proteins and nucleosome position.” Nature methods 10, no. 12 (2013): 1213-1218.

  2. 2015年 Nature scATAC-seq第一篇

    Buenrostro, Jason D., Beijing Wu, Ulrike M. Litzenburger, Dave Ruff, Michael L. Gonzales, Michael P. Snyder, Howard Y. Chang, and William J. Greenleaf. “Single-cell chromatin accessibility reveals principles of regulatory variation.” Nature 523, no. 7561 (2015): 486-490.

技术原理和实验过程

转座酶和转座子的区别和共性,转座酶切割的原理？？

技术原理

bulk

single cell: greenleaf，renbin, 10x,

实验过程

一、提取细胞核
二、细胞核完整性观察及计数
三、开放染色质片段化后建库
四、片段分选去除非目的区域的DNA

技术要点

双端测序的必要性：

技术优化

数据分析流程

bulk

1. Input

  • fastq files(fastqsanger.gz),
  • bed files(encodepeak) ,
  • genome: hg38，bed files chrom, start,end, name , score

bed格式；

2. Quality Control

  • Purpose：fastq质控
  • Tool: fastqc
  • Input：R1.fastq，R2.fastq
  • 注意点：
    • 总reads数量
    • per base sequence content: 由于转座酶的strong sequence bias，文献: “Insertion site preference of Mu, Tn5, and Tn7 transposons”
    • sequence duplication levels：PCR duplicates，后续去除
    • overrepresented sequences:adapter序列，需要用cutadapt移除

3. Trimming Reads

  • Purpose：用cutadapt移除adapter序列，用fastqc再次检查序列中adapta的含量
  • Tool: cutadapt,fastqc
  • Input:
    • adapter的序列，
    • R1.fastq，R2.fastq
  • Output:
    • R1_cutAda.fastq,R2_cutAda.fastq及报告

4. Mapping

  • Purpose：Mapping Reads to Reference Genome
  • Tool: Bowtie2
  • Input:
    • R1_cutAda.fastq,R2_cutAda.fastq
    • genome.bed
  • Output:BAM 文件
    • 插入片段长度:POS-MPOS+CIGAR
    • POS: leftmost position of where this alignment maps to the reference
    • MPOS: leftmost position of where the next alignment in this group maps to the reference, MPOS or PNEXT.
    • CICAR:string indicating alignment information that allows the storing of clipped
  • Parameters:
    • fragment length:500-1000
    • end_to_end: 不需要剪掉末端
    • mode–very_sensitive
  • 注意点：
    • dovetailing:cutadapt识别至少三个碱基以上的短序列，因此cutadapter之后的序列可能包含 1-2 个接头碱基并超出其配对起始位点。需要加入比对中。
    • 结果参数：unique mapping rate，multi-mapping多的可能原因：1. 使用–very-sensitive的模式，即使第二次命中的质量比第一次低得多，Bowtie2 也会将读取视为多重映射。另一个原因是我们读取了线粒体基因组的映射。线粒体基因组有很多具有相似序列的区域。

5. Filtering Mapped Reads

    1. Filter Uninformative Reads
      • Purpose：线粒体基因组nucleosome-free，并且Tn5可及。过滤线粒体reads，低对比质量reads，not properly paired
      • Input: aligned.bam
      • Tool: FilterBAM
      • Parameters: isProperPair:Yes, reference:!chrM, mapQuality:>=30
      • Output: aligned_filter.bam
      • 注意点：如果想保留比对到重复区域的reads， mapQuality可以降低
      • 线粒体reads的数量：使用Samtools idxstats对aligned.bam统计得到： chromosome name, chromosome length, number of reads mapping to the chromosome, number of unaligned mate whose mate is mapping to the chromosome
    1. Filter Duplicate Reads
      • Purpose：由于Tn5的插入随机，所以两个比对序列完全一样的序列被认为是PCR duplicates
      • Tool: Picard MarkDuplicates
      • Parameter：默认值标记duplicates，需要设置“If true do not write duplicates to the output file instead of writing them with appropriate flags“为yes才是去除。
      • Output：bam文件
      • Duplicates的数量在metric的文件夹里面，UNPAIRED_READS_EXAMINED， READ_PAIR_DUPLICATES
    1. Check Insert Sizes
      • Purpose：Paired-end histogram of insert size frequency，插入片段长度是ATAC-seq一个很重要的指标。
      • Tools: CollectInsertSizeMetrics (picard)
      • 结果解读： The first peak (50bp) corresponds to where the Tn5 transposase inserted into nucleosome-free regions. The second peak (a bit less than 200bp) corresponds to where Tn5 inserted around a single nucleosome. The third one (around 400bp) is where Tn5 inserted around two adjacent nucleosomes and the fourth one (around 600bp) is where Tn5 inserted around three adjacent nucleosomes.

6. Peak calling

  • Call peaks

    • Purpose：将片段信息数字化
    • Input
    • Output
    • Tool：MACS(使用的更加广泛)和Genrich(适用于ATAC，但是more reads，less peaks)
    • Mode： 1.第一个是仅选择片段长度低于 100bp 对应于无核小体区域的配对，并使用峰值调用，就像您对 ChIP-seq 所做的那样，连接配偶之间的信号。 这种方法的缺点是只有在有双端数据时才能使用它，并且会错过只有一个 Tn5 绑定的小开放区域。2.第二个是使用所有读取来更加详尽。 在这种方法中，将每个读取的信号重新集中在 5’ 末端（读取起始位点）非常重要，因为这是 Tn5 切割的地方。 实际上，您希望峰位于核小体周围，而不是直接位于核小体上。
    • Parameters：
    • 注意点：

  • Using MACS2

      1. Convert bam to bed
      • bedtools: Bam to Bed
        Input: markdup.bam
        Output: bed file
      1. Call peaks with MACS2
      • 单端测序需要设置偏移模型，双端测序不需要。在chip_seq等数据中科学家发现在真实的结合位点两侧，正负链的测序深度分布如下图所示，对应峰值的中心距离peak中心有一定的偏移。MACS首先通过一个模型来评估真实的peak中心和测序峰值的偏移距离，给定参数bandwidth和mfold, 采用一个大小为2倍bandwidth的滑动窗口，比较该窗口内真实测序深度的分布与随机测序的差异，如果二者的差异倍数超过了阈值mfold，则认为该窗口是一个peak区域。识别到初始的peak区域之后，随机挑选1000个高可信度的peak区域，分别计算正链和负链的测序深度分布。通过这种方式识别到正负链峰值之间的距离，定义为d。在后续peak calling时，会在初始计算结果的基础上向3’端偏移d/2的距离。

    • -泊松模型，landa可调

      • Input: bed file
        • Parameters:
        • We call peaks with MACS2. In order to get the coverage centered on the 5’ extended 100bp each side we will use –shift -100 and –extend 200:
        • “How many duplicate tags at the exact same location are allowed?”: all

https://blog.csdn.net/u012110870/article/details/102804191
https://zhuanlan.zhihu.com/p/90180058

7. Visualisation of Coverage

  • 1.Prepare the Datasets

    • Extract CTCF peaks on chr22 in intergenic regions
      ATAC的数据期待会在TSS区域附近富集。同时只有好的ATAC的数据才会测到在intergenic CTCF区域富集。事实上，CTCF 蛋白能够定位核小体并创建一个大约 120bp 的核小体耗尽区。 这比 TSS 周围的 200bp 无核小体区域小，也可能不存在于所有细胞中。

      • 提取chr22上面的CTCF peaks：
      • Filter data on any column using simple expressions
      • bedtools Intersect intervals find overlapping intervals in various ways

    • Convert bedgraph from MACS2 to bigwig

      • 由于MACS2生成的bedgraph文件较大且难以可视化，转化为bigwig的二进制文件。

  • 2.Create heatmap of coverage at TSS with deepTools
    分别做TSS和CTCF附近的富集图

    • Generate computeMatrix
      • Purpose: 确认特殊区域的覆盖度
      • Tools：computeMatrix：evaluate the coverage at each locus we are interested in.
      • Input：bigwig文件
      • Output：compute的matrix
    • Plot with plotHeatmap
      • 热图解读：每一行是一条transcript（非reads），覆盖度no-max，red-blue，TSS放中间，TSS 2kb的被展示。热图上方的是TSS附近的平均信号。一个bigwi件一张热图。coverage并非对称，一般在左边高，因为左边一般是激活基因的promoter的可及。
      • 如果是CTCF的区域，需要intergenic CTCF peaks chr22的reference来进行computeMatrix一步，然后再plotHeatmap。结果更加对称。
  • 3.Visualise Regions with pyGenomeTracks
    • 对于一个特殊区域的可视化，比如几个基因，可以用igv，UCSCC browser或者pyGenomeTracks
    • Tool：pyGenomeTracks
    • 结果；如果peak没有和TSS或者CTCF重合的区域，可能是enhancer，但是需要进一步确认。

single-cell

snaptools是任斌实验室开放的scATAC-seq的工具，包含上游分析和下游分析两块。
网址：https://github.com/r3fang/SnapTools

https://github.com/r3fang/SnapATAC/wiki/FAQs

  1. Index Reference Genome
    1
    2
    3
    4
    5
    6
    7
    8
     $ which bwa
    /opt/biotools/bwa/bin/bwa
    $ snaptools index-genome	\
    	--input-fasta=mm10.fa	\
    	--output-prefix=mm10	\
        --aligner=bwa	\
    	--path-to-aligner=/opt/biotools/bwa/bin/	\
    	--num-threads=5
  1. Alignment
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    $ snaptools align-paired-end	\
    	--input-reference=mm10.fa	\
    	--input-fastq1=demo.R1.fastq.gz	\
    	--input-fastq2=demo.R2.fastq.gz	\
    	--output-bam=demo.bam	\
    	--aligner=bwa	\
    	--path-to-aligner=/opt/biotools/bwa/bin/	\
    	--read-fastq-command=zcat	\
    	--min-cov=0	\
    	--num-threads=5	\
    	--if-sort=True	\
    	--tmp-folder=./	\
    	--overwrite=TRUE
  1. Pre-processing.
    a snap-format (Single-Nucleus Accessibility Profiles) file that contains meta data, cell-by-bin count matrices of a variety of resolutions, cell-by-peak count matrix.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
$ snaptools snap-pre  \
	--input-file=demo.bam  \
	--output-snap=demo.snap  \
	--genome-name=mm10  \
	--genome-size=mm10.chrom.size  \
	--min-mapq=30  \
	--min-flen=0  \
	--max-flen=1000  \
	--keep-chrm=TRUE  \
	--keep-single=TRUE  \
	--keep-secondary=False  \
	--overwrite=True  \
	--max-num=1000000  \
	--min-cov=100  \
	--verbose=True
  1. Cell-by-Bin Matrix
    1
    2
    3
    4
    $ snaptools snap-add-bmat  \
    	--snap-file=demo.snap  \
    	--bin-size-list 5000 10000  \
    	--verbose=True

技术应用

技术缺陷和展望

与其他组学的关系

H3K4me1，H3K4me2，H3K4me3和H3K27ac与染色质开放的关系：与H3K27ac相比，ATAC-Seq峰与H3K4me1的重叠具有更紧密的相关性，这表明增强子处H3K4me1的得失会影响可及染色质的形成。此外，ATAC-Seq能够识别造血过程中新的血细胞谱系调控因子。

比较：

bulk ATAC-seq的缺陷

    1. Tn5插入DNA将测序接头连接到剪断的DNA片段末端。因为两端接头连接随机， 不同接头的DNA片段才可以用于富集扩增及测序，相同接头的DNA片段无法利用。
    1. 长片段无法PCR扩增
    1. Tn5的切割活性
    1. 植物细胞难做
    1. 对细胞数量很敏感，对于转座反应和生成的DNA片段的大小分布至关重要。

Chip-seq, FAIRE-Seq and DNase-Seq的基本原理和ATAC-seq的比较

Meyer, Clifford A., and X. Shirley Liu. “Identifying and mitigating bias in next-generation sequencing methods for chromatin biology.” Nature Reviews Genetics 15, no. 11 (2014): 709-721.

传统使用的的实验方法主要是有MNase-seq和DNase-seq ，这两种实验方法的主要思路是：染色质变得开放，就意味着DNA和组蛋白的聚集程度降低，就会有一部分DNA暴露出来。而一旦失去了蛋白质的保护，这部分DNA就可以被DNA酶（MNase或DNase I）所切割。然后，我们再把切割完的DNA拿来测序，和已知的全基因组序列相比较，就能发现被切割的是哪些序列，没有被切掉的基因序列又在哪里，就知道开放的染色质区域在哪里了。不过，这两个方法有明显的缺陷，即耗时费力与重复性差。虽然FAIRE-seq 不依赖酶和抗体，但其检测背景较高，测序信噪比低，甲醛交联时间不好把握等缺陷，限制其使用范围。ATAC-seq出来的结果，和传统方法出来的结果具有很强的一致性，同时也和ChIP-seq有较高的吻合程度。而相比较而言，ATAC-seq的重复性，比MNase-seq和DNase-seq的更强，操作起来也更加简便，而且只需要很少的细胞/组织量，同时测序信号更加好。

生信分析基本流程

生物信息学数据分析由我们内部的生物信息学专家团队进行。标准的生物信息学分析包括：

    1. 序列分析：将测序读段比对到基因组，并删除重复的reads。
      峰发现：
    1. 使用MACS 2.1.0，将双端测序中的两个reads用于peak calling。
    1. 片段密度的确定： 为了鉴定基因组的转座事件的密度，将基因组分为32 bp的条带，并确定每个条带中的片段数。为此，将reads扩展到200 bp以使数据平滑。
    1. 通过随机采样对所有样本的比对reads数进行归一化，以使每一个样本包含所有样本中比对reads数最少的样本相同数目的reads。
    1. 活动区域分析： 为了比较两个或多个样本之间的峰，将重叠的峰分组为“活动区域”，这是由最上游峰的起始坐标和最下游峰的终止坐标定义的专有度量。

参考

国外教程
https://informatics.fas.harvard.edu/atac-seq-guidelines.html#design
https://training.galaxyproject.org/training-material/topics/epigenetics/tutorials/atac-seq/tutorial.html#introduction

国内入门：

https://www.activemotif.com.cn/blog-atac-seq
https://www.plob.org/article/13950.html
https://www.jianshu.com/p/87bc2002e82c 个人博客，有关于scATAC-seq的分析很多