RUISBIO ｜ 睿锶生物

ChIP-Seq生信分析报告

1.项目简介

1.1 样本信息

1.2 实验原理及流程

染色质免疫共沉淀高通量测序（Chromatin Immunoprecipitation and high-throughput sequencing，ChIP-Seq）是一种用于分析蛋白质与 DNA 交互作用的研究方法。该技术将染色质免疫沉淀（ChIP）与高通量测序技术结合，可鉴定蛋白在全基因组上的结合位点。
ChIP-Seq 可以鉴定任意具有 DNA 结合活性的蛋白在体内的结合位点，是研究表观遗传学调控和转录调控的最佳方法。
表观遗传学调控的主要方式 DNA 甲基化修饰和组蛋白修饰。组蛋白氨基酸可以发生乙酰化、磷酸化和甲基化等修饰，这些修饰可影响组蛋白与 DNA 的亲和性从而改变染色质的疏松和凝集状态，或改变其他蛋白质与 DNA 结合的亲和性从而间接地影响基因表达。DNA 甲基化修饰可以用 BS-seq 等方法进行研究，而 ChIP-Seq 则是研究组蛋白修饰的黄金方法。
转录调控主要是通过启动子区的顺式作用元件和反式作用的蛋白因子（主要是转录因子）的互作来实现的。ChIP-Seq 技术可以动态的捕捉到生物体内转录因子在基因启动子区的结合情况及强度变化，因此是研究转录调控的有力武器。

1.3 实验流程

ChIP-Seq实验主要步骤有六步:
1)使用甲醛处理组织/细胞，使得DNA与蛋白质交联;
2)破碎细胞，收集细胞核；
3)通过超声打断DNA成200-1000bp的片段;
4)利用抗体将与目的蛋白相结合的DNA片段沉淀下来;
5)解交联释放DNA片段并提取DNA;
6)构建测序文库后进行高通量测序。
ChIP-Seq 文库构建如下图所示：

ChIP-Seq DNA 文库构建流程

进行高通量测序上机前必须要构建与之匹配的测序文库，测序文库的主要构建流程如下：

DNA文库构建

1.4 分析流程

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

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

2. 数据质控

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

将Illumina测序得到的原始图像数据经过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 原始数据质控

ChIP-Seq 测序项目在NGS测序平台上完成，构建双端测序文库（~300bp）进行测序，对获得的测序数据进行质量控制，之后利用生物信息学手段对ChIP-Seq 数据进行分析。我们对所有测序Reads 的每个circle 进行碱基分布和质量波动的统计，可以从宏观上直观地反映出样本的测序质量和文库构建质量。我们对建库产生的片段（fragments）从 5’和 3’端两端分别测序，产生的序列为 paired-end reads，分别为 read 1 和read 2。为了获知测序质量情况，在测序数据下机后必须对测序数据进行质控，质控内容包括碱基质量分布、碱基平衡性分析和重复序列水平等方面。测序数据包含一些带接头、低质量的 reads，为了保证后续信息分析质量，需要对测序数据进行过滤。数据过滤的标准主要包括以下几点：
接头污染去除，去除双端测序建库的接头序列；
reads 过滤，去除含N（N 表示无法确定碱基信息）较多（≥10%）的reads；
质量过滤，对reads按照5 bp滑动窗口质量修剪，如果窗口平均质量低于20则从窗口起始位置剪切；
长度过滤，如果某reads长度小于25 bp，则同时丢弃掉R1和R2。
我们首先使用FastQC（version 0.12.1）(Andrews S, 2010)获取raw data的质量信息，之后使用Fastp(version 0.24.0)(Chen et al., 2018)软件对 raw data 进行数据过滤，之后需要对clean reads 进行新一轮质控，得到clean reads 的数据信息。原始和过滤后质控结果请详见result/1.qc文件夹，raw为原始数据质控结果，clean为过滤后质控结果。

测序数据下机后，获得每个样品的原始测序序列，我们称为 raw data 或raw reads，以 FASTQ（简写为fq）文件格式存储。为了评估测序数据质量情况，我们使用软件FastQC对数据进行评估。重点从碱基质量，碱基分布以及重复序列水平三个方面进行评估测序数据水平。

H3K27me3_Seedling_rep1_INPUT.R1H3K27me3_Seedling_rep1_INPUT.R2H3K27me3_Seedling_rep1_IP.R1H3K27me3_Seedling_rep1_IP.R2H3K27me3_Seedling_rep2_INPUT.R1H3K27me3_Seedling_rep2_INPUT.R2H3K27me3_Seedling_rep2_IP.R1H3K27me3_Seedling_rep2_IP.R2H3K27me3_Spm_rep1_INPUT.R1H3K27me3_Spm_rep1_INPUT.R2H3K27me3_Spm_rep1_IP.R1H3K27me3_Spm_rep1_IP.R2H3K27me3_Spm_rep2_INPUT.R1H3K27me3_Spm_rep2_INPUT.R2H3K27me3_Spm_rep2_IP.R1H3K27me3_Spm_rep2_IP.R2

❮

❯

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

H3K27me3_Seedling_rep1_INPUT.R1H3K27me3_Seedling_rep1_INPUT.R2H3K27me3_Seedling_rep1_IP.R1H3K27me3_Seedling_rep1_IP.R2H3K27me3_Seedling_rep2_INPUT.R1H3K27me3_Seedling_rep2_INPUT.R2H3K27me3_Seedling_rep2_IP.R1H3K27me3_Seedling_rep2_IP.R2H3K27me3_Spm_rep1_INPUT.R1H3K27me3_Spm_rep1_INPUT.R2H3K27me3_Spm_rep1_IP.R1H3K27me3_Spm_rep1_IP.R2H3K27me3_Spm_rep2_INPUT.R1H3K27me3_Spm_rep2_INPUT.R2H3K27me3_Spm_rep2_IP.R1H3K27me3_Spm_rep2_IP.R2

❮

❯

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

H3K27me3_Seedling_rep1_INPUT.R1H3K27me3_Seedling_rep1_INPUT.R2H3K27me3_Seedling_rep1_IP.R1H3K27me3_Seedling_rep1_IP.R2H3K27me3_Seedling_rep2_INPUT.R1H3K27me3_Seedling_rep2_INPUT.R2H3K27me3_Seedling_rep2_IP.R1H3K27me3_Seedling_rep2_IP.R2H3K27me3_Spm_rep1_INPUT.R1H3K27me3_Spm_rep1_INPUT.R2H3K27me3_Spm_rep1_IP.R1H3K27me3_Spm_rep1_IP.R2H3K27me3_Spm_rep2_INPUT.R1H3K27me3_Spm_rep2_INPUT.R2H3K27me3_Spm_rep2_IP.R1H3K27me3_Spm_rep2_IP.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类似。

H3K27me3_Seedling_rep1_INPUT.R1H3K27me3_Seedling_rep1_INPUT.R2H3K27me3_Seedling_rep1_IP.R1H3K27me3_Seedling_rep1_IP.R2H3K27me3_Seedling_rep2_INPUT.R1H3K27me3_Seedling_rep2_INPUT.R2H3K27me3_Seedling_rep2_IP.R1H3K27me3_Seedling_rep2_IP.R2H3K27me3_Spm_rep1_INPUT.R1H3K27me3_Spm_rep1_INPUT.R2H3K27me3_Spm_rep1_IP.R1H3K27me3_Spm_rep1_IP.R2H3K27me3_Spm_rep2_INPUT.R1H3K27me3_Spm_rep2_INPUT.R2H3K27me3_Spm_rep2_IP.R1H3K27me3_Spm_rep2_IP.R2

❮

❯

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

H3K27me3_Seedling_rep1_INPUT.R1H3K27me3_Seedling_rep1_INPUT.R2H3K27me3_Seedling_rep1_IP.R1H3K27me3_Seedling_rep1_IP.R2H3K27me3_Seedling_rep2_INPUT.R1H3K27me3_Seedling_rep2_INPUT.R2H3K27me3_Seedling_rep2_IP.R1H3K27me3_Seedling_rep2_IP.R2H3K27me3_Spm_rep1_INPUT.R1H3K27me3_Spm_rep1_INPUT.R2H3K27me3_Spm_rep1_IP.R1H3K27me3_Spm_rep1_IP.R2H3K27me3_Spm_rep2_INPUT.R1H3K27me3_Spm_rep2_INPUT.R2H3K27me3_Spm_rep2_IP.R1H3K27me3_Spm_rep2_IP.R2

❮

❯

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

H3K27me3_Seedling_rep1_INPUT.R1H3K27me3_Seedling_rep1_INPUT.R2H3K27me3_Seedling_rep1_IP.R1H3K27me3_Seedling_rep1_IP.R2H3K27me3_Seedling_rep2_INPUT.R1H3K27me3_Seedling_rep2_INPUT.R2H3K27me3_Seedling_rep2_IP.R1H3K27me3_Seedling_rep2_IP.R2H3K27me3_Spm_rep1_INPUT.R1H3K27me3_Spm_rep1_INPUT.R2H3K27me3_Spm_rep1_IP.R1H3K27me3_Spm_rep1_IP.R2H3K27me3_Spm_rep2_INPUT.R1H3K27me3_Spm_rep2_INPUT.R2H3K27me3_Spm_rep2_IP.R1H3K27me3_Spm_rep2_IP.R2

❮

❯

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

表 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)。来自一个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：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富集情况进行可视化，绘制信号在基因不同区域（TSS到TES区域以及基因上下游）的分布。对于可能富集在基因区域或者启动子区域的靶蛋白，其IP信号会富集在基因的TSS上游到TES区域，且显著高于INPUT信号。

H3K27me3_Seedling_rep1_TSSTESH3K27me3_Seedling_rep2_TSSTESH3K27me3_Spm_rep1_TSSTESH3K27me3_Spm_rep2_TSSTES

❮

❯

图3.1 各样本reads富集情况。横坐标为相对基因位置，纵坐标为按照基因组大小RPGC标准化后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）。其核心原理是通过对比实验组（IP）与对照组（Input或非特异性抗体样本）的测序信号，构建动态背景噪声模型以区分特异性结合与非随机分布的背景事件。算法首先通过双滑动窗口扫描基因组，利用泊松分布或负二项分布评估局部富集信号的显著性（P-value），并结合片段长度推断（Fragment Size）优化覆盖深度分析；随后通过两步法（粗筛峰与精细调整）定位峰边界，并确定信号峰值中心（Summit）（如果是broad模式则无summit）。本节结果请详见位于report/result/3.peak文件夹中

4.1 Peak信息统计

表 4.1 Peak信息统计：
sample：样品名称；
Peak num：Peak数量；
FRiP（Fraction of Reads in Peaks）值表示映射到峰区的 reads 占总 reads 的比例，反映了 ChIP 实验的富集效果。较高的 FRiP 值表明实验成功地富集了目标区域的 DNA 片段，而较低的 FRiP 值可能表明富集效果差或背景噪声较高。
NCS（Normalized Cross-correlation Score）则衡量富集区域的标准化交叉相关性得分，较高的 NCS 值表示目标区域的富集良好。
RCS（Relative Cross-correlation Score）评估 ChIP-Seq 数据中的伪峰情况，较低的 RCS 值通常意味着背景噪声较少，实验数据质量较高。
Peak reads：映射到峰区域的reads数；Total reads：样本总reads数。

4.2 Call Peak结果

各个样本peak信息结果表部分内容如下，完整信息请查看report/result/3.peak/{样本名称}_peaks.tsv表格。

表 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),narrowPeak存在而broadPeak不存在；
6.pileup，峰顶上的堆积高度(pileup height at peak summit)；
7.-log10(pvalue) for the peak（例如 pvalue =1e-10，那么这个值应该是 10）；
8.fold_enrichment，该峰的富集倍数，与该位置λ的随机泊松分布相对应；
9.-log10(qvalue) for the peak。

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组组内共识峰。

SpmvsSeed_PeakPca2d

❮

❯

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

SpmvsSeed_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 在基因组分布

H3K27me3_Seedling_rep1_FeatureDis_PieH3K27me3_Seedling_rep2_FeatureDis_PieH3K27me3_Spm_rep1_FeatureDis_PieH3K27me3_Spm_rep2_FeatureDis_PieSeed_FeatureDis_PieSpm_FeatureDis_PieSpmvsSeed_down_FeatureDis_PieSpmvsSeed_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_FeatureDistributionsamples_FeatureDistribution

❮

❯

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

cons_Peaks_relative_TSSdiff_Peaks_relative_TSSsamples_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代表组间差异峰注释结果。