RUISBIO ｜ 睿锶生物

ATAC-Seq生信分析报告

1.项目简介

1.1 样本信息

1.2 实验原理及流程

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

ATAC-Seq实验原理

1.3 实验流程

ATAC-seq建库主要步骤有四步：
1)收集细胞或组织，制备细胞悬液;
2)加入细胞膜裂解液，获得细胞核;
3)加入Tn5转座酶，对处于开放状态的DNA进行酶切;
4)回收酶切的DNA片段，PCR扩增和文库质控后进行二代高通量测序。

进行高通量测序上机前必须要构建与之匹配的测序文库，测序文库的主要构建流程如下：
1)回收的DNA进行片段化与接头连接；
2)将DNA片段进行PCR扩增，富集带有接头的DNA片段·；
3)利用磁珠纯化回收目的大小片段，获得待测序文库；
4)利用琼脂糖电泳方法对所构建的文库进行质检；
5)利用Qubit4.0对文库进行定量，确定文库浓度是否适合上机;
6)文库质检合格后，把不同文库按照有效浓度及目标下机数据量的需求在测序仪上进行上机测序。

DNA文库构建

1.4 分析流程

获得测序原始数据（raw data）后，首先对原始数据进行过滤，获得高质量的测序数据（clean data），将测序数据（clean data）比对到项目物种的参考基因组上，对比对结果进行鉴定染色质开放区域——峰位点（peak calling），对peak关联基因进行注释以及富集分析， 有生物学重复时进行差异Peak、共识Peak分析以及motif分析。

ATAC-Seq生物信息学分析流程

2. 数据质控

我们交付的原始数据为fastq（简称fq）格式文件的压缩包，文件名后缀通常为 “.fq.gz”。交付数据前我们会计算每个压缩文件的md5值。在您拿到数据之后，请您先校>验每个压缩文件的md5值，Linux下可以在数据目录使用“md5sum -c <*md5.txt>”命令进行校验，Windows下可使用hashmyfiles等校验工具，如发现压缩文件md5值与附在数据文件目录下的md5文档中的不一致则说明文件可能在传输的过程中被损坏。数据文件大小为文件占用磁盘空间的大小，文件的大小通常与磁盘格式、压缩比例等因素有关，与测序数据量（碱基数）的多少无对应关系，因此对应PE测序的 read1和read2两个文件大小也可能不相同。

将高通量测序得到的原始图像数据经过Base Calling 转化为序列数据，即FASTQ格式，得到最原始的测序数据文件。FASTQ 格式文件可记录所测读段（read）的碱基及其质量分数。FASTQ 格式以测序读段为单位进行存储，每条读段占 4 行，第一行是序列标识（read ID）以及相关的描述信息，以“@” 开头；第二行即为碱基序列，长度由测序策略决定；第三行以“+”开头，后面是序列标示符、描述信息，或者什么也不加； 第四行是测序质量值（phred），与第二行一一对应，phred值以ASCII码标记，对应的 ASCII 值减去33，即为第二行对应碱基的测序质量值，示例如下：

@HWI-ST1276:71:C1162ACXX:1:1101:1208:2458 1:N:0:CGATGT
NAAGAACACGTTCGGTCACCTCAGCACACTTGTGAATGTCATGGGATCCAT
+
#55???BBBBB?BA@DEEFFCFFHHFFCFFHHHHHHHFAE0ECFFD/AEHH

测序错误率用e表示， 平台测得数据的碱基质量值用Qphred表示，则有：Qphred=-10log10(e)。软件中碱基识别正确率与Phred分值之间的简明对应关系见下表：

测序Reads的错误率往往会随着测序接近尾声而升高，这是由测序过程中化学试剂的消耗造成共有的特征。

2.1 原始数据质控

ATAC-Seq实验基于第二代测序（NGS）平台完成，采用双端测序文库构建策略。我们需要对原始测序数据进行质量评估与过滤，以确保后续分析的可靠性。首先，我们使用FastQC（version 0.12.1）(Andrews, 2010)对原始测序数据（raw data）进行全局质量分析，包括碱基质量分布（Phred score）、碱基组成平衡性（base content uniformity）、重复序列比例（duplication level）及GC含量偏差等指标，以全面评估测序质量。
我们使用Fastp（version 0.24.0）(Chen et al., 2018)对原始测序数据进行以下过滤操作。
接头序列去除：识别并切除双端reads中的接头序列；
低复杂度序列过滤：剔除含模糊碱基（N碱基占比≥10%）的reads；
动态质量修剪：通过滑动窗口法（5 bp窗口步长）评估局部序列质量，当窗口平均Phred score小于20时，执行3'端截断；
长度筛选：保留长度≥25 bp的paired-end reads，长度不足的reads及其匹配reads（R1/R2）均被排除。
原始和过滤后质控结果请详见result/1.qc文件夹，raw为原始数据质控结果，clean为过滤后质控结果。

Dry_rep1_R1Dry_rep1_R2Dry_rep2_R1Dry_rep2_R2Germinating_rep1_R1Germinating_rep1_R2Germinating_rep2_R1Germinating_rep2_R2

❮

❯

图2.1 各个样本平均测序碱基质量分数，横坐标代表150 bp长度序列中各个位置，纵坐标为该位置平均的碱基质量值Q；盒形图中间的红线表示中位数（median value）；黄色部分代表四分位距（25-75%）；上下分割线代表 90%和 10%的上下临界值；蓝色的线代表碱基质量的平均值。

Dry_rep1_R1Dry_rep1_R2Dry_rep2_R1Dry_rep2_R2Germinating_rep1_R1Germinating_rep1_R2Germinating_rep2_R1Germinating_rep2_R2

❮

❯

图2.2 各个样本碱基平衡性，图中四条线代表A T C G在每个位置平均含量。理论上，A和T应该相等，G和C应该相等，且4种碱基平行且接近分布。正常情况下四种碱基的出现频率应该是接近的，而且没有位置差异。因此好的样本中四条线应该平行且接近。当部分位置碱基的比例出现 bias 时，即四条线波动较大时可能存在测序数据或者文库污染。如果所有位置的碱基比例一致的表现出bias 时，即四条线平行但分开，往往代表文库有 bias (建库过程或本身特点)，或者是测序中的系统误差。测序刚开始由于测序仪状态不稳定，在15bp之前很可能出现波动。

Dry_rep1_R1Dry_rep1_R2Dry_rep2_R1Dry_rep2_R2Germinating_rep1_R1Germinating_rep1_R2Germinating_rep2_R1Germinating_rep2_R2

❮

❯

图2.3 各个样本重复序列水平，测序深度越高，越容易产生一定程度的重复（duplication），这属于正常的现象。但如果duplication 的程度很高，就提示我们可能有 bias 的存在（如建库过程中由于 PCR 扩增引起的duplication）。横坐标为 reads 重复的次数，纵坐标为重复次数对应的 reads 占 unique reads 的比例，以unique reads 的总数作为 100%。这里，我们仅对文件前 2000000 个reads 进行统计：对长度小于75bp 的reads 将其截短为 50bp，用于统计重复。

2.2 过滤后数据质控

这里展示Fastp过滤后的数据质控结果，图片内容与上面raw data类似。

Dry_rep1_R1Dry_rep1_R2Dry_rep2_R1Dry_rep2_R2Germinating_rep1_R1Germinating_rep1_R2Germinating_rep2_R1Germinating_rep2_R2

❮

❯

图2.4 各个样本平均测序碱基质量分数，横坐标代表150 bp长度序列中各个位置，纵坐标为该位置平均的碱基质量值Q；盒形图中间的红线表示中位数（median value）；黄色部分代表四分位距（25-75%）；上下分割线代表 90%和 10%的上下临界值；蓝色的线代表碱基质量的平均值。

Dry_rep1_R1Dry_rep1_R2Dry_rep2_R1Dry_rep2_R2Germinating_rep1_R1Germinating_rep1_R2Germinating_rep2_R1Germinating_rep2_R2

❮

❯

图2.5 各个样本碱基平衡性，图中四条线代表A T C G在每个位置平均含量。理论上，A和T应该相等，G和C应该相等，且4种碱基平行且接近分布。正常情况下四种碱基的出现频率应该是接近的，而且没有位置差异。因此好的样本中四条线应该平行且接近。当部分位置碱基的比例出现 bias 时，即四条线在某些位置波动较大时，可能测序数据或者文库存在污染。当所有位置的碱基比例一致的表现出bias 时，即四条线平行但分开，往往代表文库有 bias (建库过程或本身特点)，或者是测序中的系统误差。一般测序的时候，刚开始测序仪状态不稳定，在15bp之前很可能出现波动。

Dry_rep1_R1Dry_rep1_R2Dry_rep2_R1Dry_rep2_R2Germinating_rep1_R1Germinating_rep1_R2Germinating_rep2_R1Germinating_rep2_R2

❮

❯

图2.6 各个样本重复序列水平，测序深度越高，越容易产生一定程度的重复（duplication），这属于正常的现象。但如果duplication 的程度很高，就提示我们可能有 bias 的存在（如建库过程中由于 PCR 扩增引起的duplication）。横坐标为 reads 重复的次数，纵坐标为重复次数对应的 reads 占 unique reads 的比例，以unique reads 的总数作为 100%。这里，我们仅对文件前 2000000 个reads 进行统计：对长度小于75bp 的reads 将其截短为 50bp，用于统计重复。

2.3 数据过滤结果统计

我们对数据过滤结果进行统计，如下表所示：

表 2.1数据过滤结果统计：
Sample：样品名称；
Raw_Total_Reads/Clean_Total_Reads：过滤前后样本总reads数量，单位为百万；
Raw_Total_Bases/Clean_Total_Bases：过滤前后样本总碱基数量，单位为百万；
Raw_Q20_Rate/Clean_Q20_Rate：过滤前后样本Q20碱基比例；
Raw_Q30_Rate/Clean_Q30_Rate：过滤前后样本Q30碱基比例；
Raw_GC_Content/Clean_GC_Content:过滤前后样本GC含量。

3. 比对参考基因组

我们将各样品过滤后的clean data的reads与参考基因组进行比对，获取Reads在参考基因组上的定位信息，这里使用的软件是Bowtie2(version 2.4.5)(Langmead B. et al., 2018)。由于ATAC-Seq建库时可能有较多的细胞器DNA，而在分析时它们会影响call peak结果，因此在比对后去除了比对在细胞器的reads。来自一个DNA片段的多个拷贝，可能会锚定在多个read上，经过测序得到的这些reads就是PCR重复。PCR本身就是为了产生重复序列的。理论上来讲，不同的序列在进行PCR扩增时，扩增的倍数应该是相同的。但是由于聚合酶的偏好性，PCR扩增次数过多的情况下，会导致一些序列持续扩增，而另一些序列扩增到一定程度后便不再进行，也就是我们常说的PCR偏好性。因此，比对完成后我们使用软件Sambamba（version 1.0.1）(Tarasov A. et al., 2015)去除PCR重复，获取unique reads。

3.1 比对参考基因组情况

表 3.1比对结果统计：
Sample：样品名称；
clean_reads：去除细胞器reads后clean reads总数；
PCR_dup：鉴定为PCR重复的reads数；
PCR_dup %：PCR重复reads百分比；
prop_map_reads：完美比对的reads总数，PE两端reads比对到同一条序列，且根据比对结果推断的插入片段大小符合设置的阈值；
prop_map %：完美比对reads百分比；
MAPQ30：MAPQ值大于30的reads数。

3.2 Reads富集情况

我们使用Deeptools(version 3.5.4)(Ramírez F. et al., 2016)软件对reads富集情况进行可视化，绘制信号在基因不同区域（Transcription Start Site,转录起始位点，TSS; Transcription End Site,转录终止位点，TES）的分布。对于可能富集在基因区域或者启动子区域的靶蛋白，其IP信号会富集在基因的TSS上游到TES区域，且显著高于INPUT信号。

Dry_rep1_TSSTESDry_rep2_TSSTESGerminating_rep1_TSSTESGerminating_rep2_TSSTES

❮

❯

图3.1 各样本reads富集情况。横坐标为相对基因位置，纵坐标为按照基因组大小RPGC标准化后reads富集分数。

Dry_rep1_heatmapDry_rep2_heatmapGerminating_rep1_heatmapGerminating_rep2_heatmap

❮

❯

图3.2 各样本reads富集热图。下方热图代表基因上下游Reads富集情况，每一行代表一个基因上下游区域reads富集程度。

3.3 比对可视化

软件比对所得结果为bam格式文件(位于report/result/2.map文件夹中)，bam文件是压缩的⼆进制⽂件，无法直接作为文本打开查看。由于bam文件数据较大，我们将其转为较小的bw格式文件。客户可以结合物种参考基因组和注释文件使用IGV (Integrative Genomics Viewer) 浏览器对bam、bw、bed等文件进行可视化浏览。IGV浏览器使用方法可参考我们提供的使用说明文档IGV快速上手

4. 峰鉴定

我们使用MACS3（version 3.0.0）(Zhang Y. et al. 2008)进行peak鉴定，即找到染色质开放性高的基因组区域。MACS(Model-based Analysis of ChIP-Seq)是一种基于统计学模型的算法，专用于从表观组学测序数据中精准识别蛋白富集区域或染色质高开放性区域（Peaks）。其核心原理是通过构建动态背景噪声模型以区分特异性结合与非随机分布的背景事件。算法首先通过双滑动窗口扫描基因组，利用泊松分布或负二项分布评估局部富集信号的显著性（P-value），并结合片段长度推断（Fragment Size）优化覆盖深度分析；随后通过两步法（粗筛峰与精细调整）定位峰边界，并确定信号峰值中心（Summit）。使用Phantompeakqualtools(version 1.2.2)(Landt SG. et al. 2012)进行基于phantom peak的数据质量评估，计算NSC和RSC。本节结果请详见位于report/result/3.peak文件夹中

4.1 Peak信息统计

表 4.1 Peak信息统计：
sample：样品名称；
Peak num：Peak数量；
FRiP（Fraction of Reads in Peaks）值表示映射到峰区的 reads 占总 reads 的比例，反映了ATAC-Seq实验的Tn5酶富集效果。较高的 FRiP 值表明实验成功地富集了目标区域的 DNA 片段，而较低的 FRiP 值可能表明富集效果差或背景噪声较高。
NSC（Normalized strand cross-correlation coefficent）衡量富集区域的标准化交叉相关性得分。NSC值越大表明富集效果越好，NSC值低于1.1表明较弱的富集，小于等于1表示无富集。
RSC（Relative strand cross-correlation coefficient）评估ATAC-Seq数据中的伪峰情况，RSC是片段长度相关值减去背景相关值除以phantom-peak相关值减去背景相关值。RSC的最小值可能是0，表示无信号；富集好的实验RSC值大于1；低于0.8表示质量偏低。
Peak reads：映射到峰区域的reads数；Total reads：样本总reads数。

4.2 Call Peak结果

各个样本peak信息结果表部分内容如下，完整信息请查看report/result/3.peak/{样本名称}_peaks.tsv表格。“.broadPeak”或“.narrowPeak”文件为MACS3的输出文件，用于描述峰区域信息，可在IGV浏览器中打开。

表 4.2 Call Peak结果。样本_peaks.tsv是一个表格文件，其中包含有关被调用峰的信息。您可以在excel/WPS中打开它并使用函数进行排序/过滤。各列信息为：
1.chr，染色体名称；
2.start，peak起始位置；
3.end，peak的结束位置；
4.length，peak长度；
5.abs_summit,峰顶的位置(absolute peak summit position)；
6.pileup，峰顶上的堆积高度(pileup height at peak summit)；
7.-log10(pvalue) for the peak（例如 pvalue =1e-10，那么这个值应该是 10）；
8.fold_enrichment，该峰的富集倍数，与该位置λ的随机泊松分布相对应，peak文件中signalValue列等于该列；
9.-log10(qvalue) for the peak,peak文件中score列是该列数值x10。

5. 差异Peak分析

在存在多个分组且组内有生物学重复的情况下，可以对组间进行差异Peak(differential peak)分析，以确定哪些Peak在组间存在显著差异，同时获取组内共识峰（consensus peak）。如果没有差异分析则本节内容为空。

5.1 差异Peak分析结果

存在组内生物学重复时，我们使用软件DiffBind(version 3.10)(Stark,R., & Brown,G.,2012)对样本peaks进行分析。结果详见report/result/4.peak。后文中提到的“diff”代表组间差异，“cons”代表组内交集。
DiffPeak:“sampAvssampB_res.csv”为各组样本差异分析结果，sampA代表实验组，sampB代表对照组； DiffPeak_sampAvssampB_up.bed为sampA和sampB比较，结合强度上调的peak; DiffPeak_sampA vs sampB_down.bed为sampA和sampB比较，结合强度下调的peak; sampA_consensus_peaks.bed为A组组内共识峰。

GerminatingvsDry_PeakPca2d

❮

❯

图5.1 比较组PCA图。主成分分析是将原来较多维度的指标 (peak 的分布特征)，降维到较低的维度（二维），来研究样品间的主成分关系。二维PCA分析结果中，会展示主成分1(PC1) 和主成分2(PC2)分别作为 X 轴和 Y 轴的散点图，每个点代表 1 个样本。坐标轴上百分比代表主成分的贡献率，贡献率越大，说明该主成分对样本差异的解释能力越强。如果两个样本距离越远，则说明样本 peaks 分布的差异越大。 反之，则说明相应样本peaks整体分布模式越接近。所以，PCA 分析常用于评估样本重复性的好坏。理想情况下，生物学重复的样本应该聚类在一起，而不同组间应该可以区分开。

GerminatingvsDry_vol

❮

❯

图5.2 差异Peak火山图。横坐标为log2(Fold Change)，纵坐标为-log10(FDR)，蓝色为显著性下调的峰，红色为显著性上调的峰，灰色为非显著性差异的峰。

6. 基因组注释

为了进一步探讨peak结合位点特征，理解染色质开放区域对基因调控的机制， 使用R包ChIPseeker(version 1.36)(Wang et al., 2022)对Peak区域进行注释，我们统计Peak在各基因功能元件分布情况，并将各个peak与基因关联。本节结果请详见位于report/result/5.anno文件夹。

6.1 Peak 在基因组分布

Dry_FeatureDis_PieDry_rep1_FeatureDis_PieDry_rep2_FeatureDis_PieGerminating_FeatureDis_PieGerminating_rep1_FeatureDis_PieGerminating_rep2_FeatureDis_PieGerminatingvsDry_down_FeatureDis_PieGerminatingvsDry_up_FeatureDis_Pie

❮

❯

图6.1 Peak在基因功能元件上分布饼图。
一般来说，peaks最多的区域是位于转录起始点(TSS)上游1kb的启动子区域“promoter（<=1kb）”，它与基因的表达调控密切相关；“promoter（1~2kb）”代表TSS上游1~2kb的启动子区域，“promoter（2~3kb）”代表TSS上游2~3kb的启动子区域。
5'非翻译区（5' UTR）和外显子区域（Exon）与mRNA的稳定性或基因表达的调控有关。
3'非翻译区（3' UTR）、内含子（Intron）、远端基因间区（Distal Intergenic）以及TSS下游区（Downstream），这些区域的调控活动可能涉及长距离的基因调控或影响基因的后续处理和表达。

cons_FeatureDistributiondiff_FeatureDistributionsample_FeatureDistribution

❮

❯

图6.2 各样本Peak在基因功能元件上分布比例堆叠条状图，samples代表单个样本，cons代表组内共识峰，diff代表组间差异峰。

cons_Peaks_relative_TSSdiff_Peaks_relative_TSSsample_Peaks_relative_TSS

❮

❯

图6.3 各样本Peak在TSS(转录起始位点)侧翼分布比例堆叠条状图，samples代表单个样本，cons代表组内共识峰，diff代表组间差异峰，各元件内容含义见图 6.1。

6.2 Peak关联基因注释

各个样本Peak关联基因注释结果表部分内容如下，完整信息请查看report/result/5.anno/{样本名称}_PeakAnno.csv表格。{组名}_PeakAnno.csv代表组内共识峰注释结果，{比较组}_{up/down}_PeakAnno.csv代表组间差异峰注释结果。