ATACSeq 开放染色质分析生信生信

Chip-seq学习记录

2019-11-21  本文已影响0人  小梦游仙境

关于概念问题，很重要

第1篇：ATAC-seq的背景介绍以及与ChIP-Seq的异同https://mp.weixin.qq.com/s/gLcNZkzzQI4AvZs-0DgG5Q

流程

image-20191014150758466

ATAC系列连载：

第1篇：ATAC-seq的背景介绍以及与ChIP-Seq的异同

https://mp.weixin.qq.com/s/gLcNZkzzQI4AvZs-0DgG5Q

【ATAC-Seq简介】

ATAC-seq（Assay for Transposase-Accessible Chromatin with high throughput sequencing） 是2013年由斯坦福大学William J. Greenleaf和Howard Y. Chang实验室开发的用于研究染色质可及性（通常也理解为染色质的开放性）的方法， 原理是通过转座酶Tn5容易结合在开放染色质的特性，然后对Tn5酶捕获到的DNA序列进行测序。
真核生物的核DNA并不是裸露的，而是与组蛋白结合形成染色体的基本结构单位核小体，核小体再经逐步的压缩折叠最终形成染色体高级结构（如人的DNA链完整展开约2m长，经过这样的折叠就变成了纳米级至微米级的染色质结构而可以储存在小小的细胞核）。而DNA的复制转录是需要将DNA的紧密结构打开，从而允许一些调控因子结合（转录因子或其他调控因子）。这部分打开的染色质，就叫开放染色质，打开的染色质允许其他调控因子结合的特性称为染色质的可及性（chromatin accessibility）。因此，认为染色质的可及性与转录调控密切相关。
开放染色质的研究方法有ATAC-seq以及传统的DNase-Seq及FAIRE-seq等，ATAC-Seq由于所需细胞量少，实验简单，可以在全基因组范围内检测染色质的开放状态，目前已经成为研究染色质开放性的首选技术方法。

img

[Nat Methods, 2013. doi: 10.1038/nmeth.2688. Epub 2013 Oct. ]

【ATAC-Seq与CHIP-Seq的异同】

ATAC-Seq与ChIP-Seq的不同的是ATAC-Seq是全基因组范围内检测染色质的开放程度，可以得到全基因组范围内的蛋白质可能结合的位点信息，一般用于不知道特定的转录因子，用此方法与其他方法结合筛查感兴趣的特定调控因子；但是ChIP-Seq是明确知道感兴趣的转录因子是什么，根据感兴趣的转录因子设计抗体去做ChIP实验拉DNA，验证感兴趣的转录因子是否与DNA存在相互作用。ATAC-Seq、ChIP-Seq、Dnase-Seq、MNase-Seq、FAIRE-Seq整体的分析思路一致，找到富集区域，对富集区域进行功能分析。

[An overview of ChIP–seq, DNase-seq, ATAC-seq, MNase-seq and FAIRE–seq experiments]

【CHIP-Seq简介】

ChIP实验（Chromatin immunoprecipitation）即染色质免疫沉淀，根据DNA与蛋白质相互作用的原理，分离富集与感兴趣的蛋白相互作用的DNA。ChIP-Seq即对分离得到的DNA扩增测序，然后通过分析得到DNA的富集区域也称为peaks，同时可以鉴定过表达的序列motif以及进行功能注释分析。

img

[ChIP-Seq Overall]

下面这一部分将会介绍ChIP-Seq数据分析的整个流程，从实验设计到产生原始的测序reads，以及到最后的功能富集分析和motif查找。

img

实验设计和文库构建

文库构建包括以下5步骤：

富集到的DNA片段只有一部分是真实的信号（感兴趣的蛋白结合的DNA区域），这个比例取决于number of active binding sites, the number of starting genomes, and the efficiency of the IP.

ChIP-Seq富集序列存在以下特点：

因此，ChIP-Seq需要有合适的对照组，对照样本需要满足以下其中一个条件：

img img

示例数据介绍

所用到的示例数据是来自于人类胚胎干细胞系（h1-ESC）中Nanog和Pou5f1（Oct4）两个转录因子的结合图谱。这两个转录因子的功能涉及干细胞的多能性，该研究的目标之一是探究这两个转录因子在转录调控中单独和相互的调控作用。

两组重复，每组重复包括3个实验样本信息，共6个样本，数据分析中只用到了12号染色体的信息。

img

分析环境配置

这个课程提供了示例数据和分析代码，可以参考这里连接他们的服务器，我没有连接成功，不知道是不是打开方式不对，大家可以尝试下，如果连接成功，这一部分就是配置服务器的环境，准备数据；如果也连接不上可以用自己的数据或者下载公共数据。

img

第2篇：原始数据的质控、比对和过滤

https://mp.weixin.qq.com/s/WmMBMYFCsXkl4QZ7y8AoEg

这部分内容包括对原始测序数据质控，然后比对过滤，这是所有NGS数据处理的上游分析。

学习目标

测序reads 的质控流程示意图

img

img

FASTQC

首先对拿到的原始测序数据（fastq或fastq.gz格式）进行质控检测，直接用fastqc软件，再加上multiqc将多个检测结果一起展示。
如:

fastqc -o out_dir raw_data/*gz
multiqc *fastqc.zip --ignore *.html

Trimmomatic

Trimmomatic 可以用于去除接头，过滤低质量数据。相同功能的软件还有很多，如trim_galore、cutadapt等，个人比较喜欢trim_galore可以自动识别接头类型。

# 课程中给出的Trimmomatic 的用法（单端测序）
$ java -jar /opt/Trimmomatic-0.33/trimmomatic-0.33.jar SE \
-threads 2 \
-phred33 \
H1hesc_Input_Rep1_chr12.fastq \
../results/trimmed/H1hesc_Input_Rep1_chr12.qualtrim20.minlen36.fq \
LEADING:20 \
TRAILING:20 \
MINLEN:36

Trimmomatic参数含义：可以参考NGS 数据过滤之 Trimmomatic 详细说明
trim_galore使用示例

trim_galore -q 20 --phred33 --stringency 3 --length 20 -e 0.1 --paired fq1 fq2  --gzip -o input_data_dir
# 重新用fastqc检测进行过滤后的reads质量
fastqc -o out_dir *fq.gz
multiqc *fastqc.zip --ignore *.html

比对

Bowtie2是一个快速精确的比对工具，基于Burrows-Wheeler Transform 构建基因组的FM 索引，比对过程所耗内存少。Bowtie2支持局部、双端、缺口比对模式，对大于50bp的reads比对效果更好（小于50bp的reads用Bowtie1）。

创建Bowtie2索引

bowtie2-build <path_to_reference_genome.fa> <prefix_to_name_indexes>
# Can find indexes for the entire genome on Orchestra using following path: /groups/shared_databases/igenome/Homo_sapiens/UCSC/hg19/Sequence/Bowtie2Index/

Bowtie2 比对

## 课程中给出的代码是单端比对
bowtie2 -p 2 -q \
-x ~/ngs_course/chipseq/reference_data/chr12\
-U ~/ngs_course/chipseq/results/trimmed/H1hesc_Input_Rep1_chr12.qualtrim20.minlen36.fq \
-S ~/ngs_course/chipseq/results/bowtie2/H1hesc_Input_Rep1_chr12_aln_unsorted.sam

NOTE: 如果fastq文件是没有经过trim的，可以用局部比对执行soft-clipping，加上参数--local

过滤reads

首先将sam文件转为bam格式，再对bam文件进行排序，接着过滤唯一比对的reads，去除线粒体reads。
转化为bam格式
使用samtools转换格式

samtools view -h -S -b \
-o H1hesc_Input_Rep1_chr12_aln_unsorted.bam \
H1hesc_Input_Rep1_chr12_aln_unsorted.sam

对bam文件排序
对bam文件按照基因组坐标排序，可以直接使用samtools，也可以使用Sambamba。sambamba快速处理bam和sam文件。

sambamba sort -t 2 \
-o H1hesc_Input_Rep1_chr12_aln_sorted.bam \
H1hesc_Input_Rep1_chr12_aln_unsorted.bam

过滤唯一比对的reads

sambamba view -h -t 2 -f bam \
-F "[XS] == null and not unmapped " \
H1hesc_Input_Rep1_chr12_aln_sorted.bam > H1hesc_Input_Rep1_chr12_aln.bam

去除PCR重复
PCR扩增和一些重复序列（如微卫星、着丝粒）会产生重复，干扰真实的富集信号，所以在call peaks前需要先去除重复，这里先用picard去除PCR重复。

java -jar picard-tools-1.119/MarkDuplicates.jar REMOVE_DUPLICATES=true I=H1hesc_Input_Rep1_chr12_aln.bam O=H1hesc_Input_Rep1_chr12_aln.dedup.bam M=H1hesc.duplicates.log

过滤线粒体reads

samtools index H1hesc_Input_Rep1_chr12_aln.dedup.bam
samtools idxstats H1hesc_Input_Rep1_chr12_aln.dedup.bam > H1hesc_Input_Rep1_chr12_aln.dedup.mitochondrial.stats
samtools view -h H1hesc_Input_Rep1_chr12_aln.dedup.bam | grep -v 'chrM' | samtools view -bS -o H1hesc.final.bam

上面给出的仅是示例代码，和参考课程不一样，实际运行需要修改相应文件。
此时就得到了唯一比对且已经去除过线粒体的比对文件，可以用于接下来的peaks calling。

第3篇：用MACS2软件call peaks

https://mp.weixin.qq.com/s/IjbRisuSnfUWsYs1dSn5_g

学习目标

1

Peak Calling

Peak calling即利用计算的方法找出ChIP-seq或ATAC-seq中reads富集的基因组区域。

img

如下图所示，比对结果的文件中reads在正负链不均匀分布，但在结合位点聚集。正负链5‘末端的reads各形成一组合，通过统计学的方法评估这些组合的分布并和对照组比较，确定这些结合位点是否是显著的。

img

NOTE：ChIP-seq的分析方法可以鉴定两种类型的富集模式：broad domains和narrow peaks。broad domains，如组蛋白修饰在整个基因body区域的分布；narrow peak，如转录因子的结合。narrow peak相对于broad 或者分散的marks更易被检测到。也有一些混合的结合图谱，如PolII包括narrow和broad信号。

2

MACS2

peaks calling 有不同的方法，MACS2是最常用的call peaks工具。 MACS全称Model-based Analysis of ChIP-Seq，最初的设计是用来鉴定转录因子的结合位点，但是它也可以用于其他类型的富集方式测序。

MACS通过整合序列标签位置信息和方向信息提高结合位点的空间分辨率。MACS的工作流如下所示：

img

MACS2的用法，call peaks的参数及输出文件的解读可以参考MACS2文档学习。

了解相关参数：

输入文件参数：

输出文件参数：

peak calling 参数

q值与峰宽有一定的联系。理想情况下，如果放宽阈值，您将简单地获得更多的峰值，但是使用MACS2放松阈值也会导致更宽的峰值。

Shift 模型参数：

  1. 想找富集剪切位点，如DNAse-seq，所有5'端的序列reads应该从两个方向延伸，如果想设置移动的窗口是200bp，参数设置如下：
    --nomodel --shift -100 --extsize 200

  2. 对nucleosome-seq数据，用核小体大小的一半进行extsize,所以参数设置如下：
    --nomodel --shift 37 --extsize 73

    --call-summits

MACS利用此参数重新分析信号谱，解析每个peak中包含的subpeak。对相似的结合图谱，推荐使用此参数，当使用此参数时，输出的subpeak会有相同的peak边界，不同的绩点和peak summit poisitions.

3 ATAC-Seq call peaks示例

ATAC-seq关心的是在哪切断，断点才是peak的中心，所以使用shift模型，--shift -75或-100

对人细胞系ATAC-seq 数据call peak的参数设置如下：

macs2 callpeak -t H1hesc.final.bam -n sample --shift -100 --extsize 200 --nomodel -B --SPMR -g hs --outdir Macs2_out 2> sample.macs2.log

MACS2输出文件解读

**Coordinates in XLS is 1-based which is different with BED format**
XLS里的坐标和bed格式的坐标还不一样，起始坐标需要减1才与narrowPeak的起始坐标一样。

![img](https://img.haomeiwen.com/i11316862/3fde17133be280fe.jpg?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)

Remove the beginning track line if you want to analyze it by other tools.???

summits.bed, narrowPeak, bdg, xls四种输出文件的比较

img

第4篇：对ATAC-Seq/ChIP-seq的质量评估（一）——phantompeakqualtools

https://mp.weixin.qq.com/s/b-GbjioGQH2JXwcvtRt51w

学习目标

ChIP-Seq质量评估

在下游分析前，最好是先对peak calling 后的ChIP-Seq数据进行质量评估。

img

链交叉相关（Strand cross-correlation）

链交叉相关是一个有效的评估ChIP-Seq质量的方法，它不依赖于peak calling，而是基于ChIP-Seq实验。如果ChIP-Seq实验成功，DNA富集序列标签（蛋白质相互作用的序列）会在reads的双峰富集中产生显著的聚集。
产生reads的双峰富集的原因如下：
在ChIP-Seq实验中，DNA被片段化，蛋白质结合的片段会被免疫沉淀，所以产生了有蛋白质结合的DNA片段（fragments ）。
DNA的正链从5'端开始被测序（如下图红色reads），DNA负链也从5’末端被测序产生如下图所示的蓝色reads。

Nat Biotechnol. 2008 Dec; 26(https://tva1.sinaimg.cn/large/006y8mN6gy1g7xxlfc3ufj307l08t0sq.jpg): 1351–1359

Nat Biotechnol. 2008 Dec; 26(12): 1351–1359

由于从DNA片段的5′末端测序，使+链reads的富集（下图中的蓝色部分）与负链reads的富集（下图红色部分）有少量的相互抵消区域。我们我们可以用计算产生最大相关的位移。

img

交叉相关性度量值
交叉相关度量是在Watson移动k个碱基后，计算Crick链与Watson链之间的Pearson线性相关。如下面的示意图：
首先在链位移为0时，两个向量之间的Pearson相关值为0.539。

img

在链位移5bp，两个向量之间的Pearson相关值为0.931。

img

继续移动这些向量，对于每个链位移计算一个相关值。

img

最后，我们将有一个每个碱基对移位与皮尔森相关值的对应表。这是针对每个染色体的每一个峰计算的，然后该值乘以一个缩放因子，再对所有染色体的值相加，就可以绘交叉相关值（y轴）相对于移位值（x轴）生成的交叉相关图。典型的交叉相关图会产生两个峰****片段长度（predominant fragment length）****read 长度（read length）

交叉相关性质量评估度量值

交叉相关谱图可以计算评估ChIP_Seq实验信噪比的度量值，并且实验设计确保fragment length准确。低信噪比和不准确的fragment length 表明ChIP-Seq实验可能有问题。
Normalized strand cross-correlation coefficent (NSC)：
NSC是最大交叉相关值除以背景交叉相关的比率(所有可能的链转移的最小交叉相关值)。NSC值越大表明富集效果越好，NSC值低于1.1 表明较弱的富集，小于1表示无富集。NSC值稍微低于1.05，有较低的信噪比或很少的峰，这肯能是生物学真实现象，比如有的因子在特定组织类型中只有很少的结合位点；也可能确实是数据质量差。
Relative strand cross-correlation coefficient (RSC):
RSC是片段长度相关值减去背景相关值除以phantom-peak相关值减去背景相关值。RSC的最小值可能是0，表示无信号；富集好的实验RSC值大于1；低于1表示质量低。

phantompeakqualtools

phantompeakqualtools 是一个用于计算ChIP-Seq数据富集和质量度量值的一个工具包。我们将使用该包来计算基于链交叉相关峰的主要插入大小(fragment length)和基于相对phantom peak的数据质量度量值。phantompeakqualtools是一个R包，依赖samtools。
下载phantompeakqualtools

wget https://storage.googleapis.com/google-code-archive-downloads/v2/code.google.com/phantompeakqualtools/ccQualityControl.v.1.1.tar.gz
tar -xzf ccQualityControl.v.1.1.tar.gz
cd phantompeakqualtools
# 查看README
less README.txt

Linux下安装

R 
install.packages("caTools", lib="~/R/library")

运行phantompeakqualtools

mkdir -p logs qual

for bam in bam_dir/sample1.final.bam bam_dir/sample2.final.bam
do 
bam2=`basename $bam .final.bam`
Rscript run_spp_nodups.R -c=$bam -savp -out=qual/${bam2}.qual > logs/${bam2}.Rout
done

参数含义：

输出文件解读
输出文件会产生一个tab分割的名为qual的文件,包含的信息如下：

ATAC-seq数据特有的fragment size分布

在第一篇ATAC-seq文章里面，如图：

img

可以使用atacQC 对 ATAC文库出 Fragment size distribution ，示例图如下：

img

最后友情宣传生信技能树

上一篇下一篇

猜你喜欢

热点阅读