1.分析流程
染色质免疫共沉淀测序 技术是将染色质免疫沉淀技术(Chromatin Immunoprecipitation,ChIP)与新一代测序技术(NGS)相结合,能够高效地在全基因组范围内检测与组蛋白修饰、转录因子等互作的DNA区段。利用该种技术,通过将测序结果精确定位到基因组上,研究人员可获得全基因组范围内与组蛋白、转录因子等互作的DNA区段信息。
生物信息分析流程
2. 分析结果
1. 数据质控
1.1 FastQC(raw_reads)
随着高通量测序技术的不断成熟,通过测序可以获得海量的数据信息。如何从得到的数据中提取出我们所需的部分成为关键。而在数据分析前进行的首项分析是看一下下机数据的质量前比。较流行的数据质量查看软件就是FastQC。采用FastQC对下机后的数据(raw reads)质量进行基本的统计,结果如下:
图 1-1 raw reads of FastQC
上左图 :横轴代表位置,纵轴quality,颜色区分为分位数 ;
上右图 :横轴为位置,纵轴为碱基百分比 ;
中左图 :reads长度的分布 ;
中右图 :统计reads的平均GC含量的分布,红线是实际情况,蓝线是理论分布 ;
下左图 :对所有reads的每个位置,统计N的比率 ;
下右图 :横坐标是 duplication 的次数,纵坐标是 duplicated reads 的数目。
1.2 原始数据 Trimming
在对库检合格样本测序的过程中,随机引物的不完全配对可能会引起reads前几个位置 的核苷酸组成存在一定的偏好,导致测序过程引入错误碱基。同时,由于片段较短导致测序获得的读段中,有一部分会含有接头序列, 需要将接头序列去掉。为了保证信息分析质量,在进行分析以前需要对原始下机数据(raw data)进行筛选,获得clean data,后续分析都基于clean data。
数据处理的步骤如下:
1) 去除低质量reads(质量值sQ小于等于的碱基数占整个read 的50%以上的 reads);
2) 去除N(N 表示无法确定碱基信息)的比例大于15%的reads;
3) 去除有5’ 接头污染的reads;
4) 去除没有3’ 接头序列和插入片段的reads;
5) Trim掉3’ 接头序列;
6) 舍弃修剪后短于18nt的read。
原始数据质控统计
对原始数据进行修剪,切除读段中与 Adpter 序列匹配的序列, 以及 3'端质量值低于 20 的碱基,如果所剩下的序列长度不短于 18nt ,则保留。
Raw reads 指所有的原始读段;
Clean reads 指修剪后被保留的读段,包括 Trimmed 和 Utrimmed 读段;
Clean rate 指 Clean_reads 在 Raw_reads 中的比例。
1.3 FastQC Clean_reads
将经过 Trimming后得到的数据(clean reads)再次进行质量统计,结果如下:
图 1-2 clean reads of FastQC
上左图 :横轴代表位置,纵轴quality,颜色区分为分位数 ;
上右图 :横轴为位置,纵轴为碱基百分比 ;
中左图 :reads长度的分布 ;
中右图 :统计reads的平均GC含量的分布,红线是实际情况,蓝线是理论分布 ;
下左图 :对所有reads的每个位置,统计N的比率 ;
下右图 :横坐标是 duplication 的次数,纵坐标是 duplicated reads 的数目。
2. 序列比对
2.1 基本统计
目前,常用比较软件有Bowtie,BWA,MAQ,TOPhat 等。根据不同的基因组的特征,我们选取相对合适的软件,合适的参数设置,将过滤后的测序序列进行基因组定位分析。由于ChIP-seq打断片段通常偏小,且唯一序列占整体序列数量的百分比是关注的重点,因此,利用BWA(Burrows Wheeler Aligner)可以更准确的把reads比对到参考基因 组上(Li, H. and R. Durbin, 2009)。将读段与参考序列进行比对统计如下表:
2.2 比对质量统计
ChIP-seq分析中,唯一序列占整体序列数量的百分比是关注的重点,通过标记重复读段(Duplicates), 计算比对质量值(MAPQ),选取合适质量值作为唯一比对的阈值。在这里我们选取30作为阈值,即相应读段非唯一比对的概率只有0.01。后续的检峰过程中,对于重复读段一般只保留一条。
图 2- 1 比对质量分布图
图 2-2 读段在基因组上的分布统计
图 2-3 比对质量分布图
2.3 相对基因的位置分布统计
转录因子、组蛋白等对基因的调控机制与蛋白质的结合位置相关,故分析读段相对基因位置分布,有助于我们预测蛋白的功能。将每个基因以及该基因上下游 2K 各自等分成100等分,统计落到每个小格中的读段数,计算每个小格内的reads 数占这些区域总reads数百分比,作为各区段的 reads 密度。
图 4- 1 相对基因位置分布图
横坐标 :基因相对位置;纵坐标:reads 密度。
3. 样本相关性检查
生物学重复是任何生物学实验所必须的,高通量测序技术也不例外(Hansen et al.)。生物学重复主要有两个用途:一个是证明所涉及的生物学实验操作是可以重复的且变异不大,另一个是为了确保后续的差异基因分析得到更可靠的结果。样品间相关性分析是检验实验可靠性和样本选择是否合理的重要指标。相关系数越接近1,表明样品之间表达模式的相似度越高。
图 4-2 样品相关系数热图
样品的相关系数热图,样品间的相关性系数(皮尔逊相关系数), 颜色的深浅表示相关性系数的大小。
4. SCC 统计
我们知道,测得的reads 最终会近似平均地分配到正负链上。这样,我们通过计算正负链间的相关系数(链互相关系数)检测两条链之间的最优距离。通过对 IP 和Input 数据的链互相关性检测不仅可以获得正负链间相关系数,同时也可以监测到IP实验效果。
NSC:相对链互相关系系数,一般不低于 1.05
RSC:reads长度矫正后链相互关系系数, 一般不低于0.8
5. 峰检出
转录因子、组蛋白等蛋白结合位点的注释是探究其调节机制和功能的重要信息。近年来随着新一代测序技术的发展,将染色质免疫沉淀后的DNA直接进行高通量测序,对比基准基因(reference sequence)可直接获得蛋白与DNA的结合位点信息。
利用 MACS2 软件(Yong Zhang,Tao Liu et al., 2008) (阈值为 qvalue=0.05)完成峰检分析(peak calling), 并对峰的个数、宽度、分布等进行统计,筛选出峰的相关基因等。结果如下:
5.1 峰的检出情况
Frip:peak 中 reads 的个数占整个基因组 reads 比率。可以检测 IP 试验效果
5.2 峰的宽度分布
5.3 峰的富集度
计算每个 peak 的显著性(q value)值。 图中横坐标表示峰的-log10 q value ,纵坐标为密度。
横坐标为峰的富集倍数 ,纵坐标为这种宽度的峰密度。
6 模体分析
转录因子、组蛋白等蛋白质与DNA序列的结合并不是随机的,而具有一定的序列偏好性。模体(Motif)分析不仅可以检测到蛋白质特异性结合位点的DNA序列的偏好性, 同时通过模体注释可以获得已知Motif的注释以及蛋白结合位点、Motif序列信息等。采用MEME(Timothy L. Bailey and Charles Elkan,1994)与Dreme(Timothy L. Bailey,2011)软件检测peak序列中显著Motif序列,采用Tomtom(Shobhit Gupta, JA Stamatoyannopolous,2007)软件将得到的 Motif序列与已知Motif数据库进行比对,利用已知Motif对其进行相应注释。
6.1 模体检测
Meme
Dreme
图 6- 1
(注 : 由于结合位点的特异性 ,会出现 Motif 序列只出现在一个区段(<=8 或者>=9),则下图会有一部分没有结果)
6.2 meme 模体注释
采用序列表示图( sequence logo)的方式展示 MEME 检测到的 Motif 与已知 Motif 的比对结果,如下图:
TOMTOM
图 6-2 MEME 模体注释
Meme TOMTOM
图 6-3 MEME 模体注释
6.3 Dreme 模体注释
采用序列表示图( sequence logo)的方式展示 Dreme 检测到的 Motif 与已知 Motif 的比对结果,如下图:
TOMTOM
图 6-4 MEME 模体注释
7. 峰注释
7.1 峰位置注释
通过蛋白结合的特点预测蛋白的调节模式或者功能。每个peak相关基因的TS(transcription start site)位点采用Chip seeker检测,根据peak与TSS位点的距离统计 peak百分比,分析peak与基因TSS位点的距离分布。
图 7- 1 peak 与基因 TSS 位点的距离分布
7.2 GO 富集分析
Peak 重叠基因 GO 富集柱状图,直观的反映出在生物过程(biological process)、 细胞组分(cellular component)和分子功能(molecular function)富集的 GO 项目上 Peak 重叠基因的个数分布情况。
7.3 KEGG 分析
在生物体内,不同基因相互协调行使其生物学功能,通过Pathway显著性富集能确定差异表达基因参与的最主要生化代谢途径和信号转导途径。KEGG(Kyoto Encyclopedia of Genes and Genomes)是系统分析基因功能、基因组信息数据库,它有助于研究者把基因及表达信息作为一个整体网络进行研究。作为Pathway相关的主要公共数据库(Kanehisa,2008),KEGG提供的整合代谢途径 (pathway)查询十分出色,包括碳水化合物、核苷、氨基酸等的代谢及有机物的生物降解,不仅提供了所有可能的代谢途径,而且对催化各步反应的酶进行 了全面的注解,包含有氨基酸序列、PDB 库的链接等等,是进行生物体内代谢分析、代谢网络研究的强有力工具。Pathway 显著性富集分析以 KEGG 数据库中 Pathway 为单位,应用超几何检验,找出与整个基因组背景相比,在差异表达基因中显著性富集的 Pathway。
3. 背景知识
1 实验与建库测序
ChIP-Seq(Chromatin Immunoprecipitation Sequencing ,染色质免疫共沉淀测序)结合了染色质免疫共沉淀(ChIP)与高通量测序,是研究蛋白与 DNA 相互作用的有力手段,广泛用于研究全基因组范围转录因子结合位点以及组蛋白修饰的状态(Landt and Marinov et al., 2012)。 一般的 ChIP-Seq 实验及建库测序流程如下图所示,我们目前提供从建库开始往后的各阶段服务。
使用甲醛固定细胞,将 DNA 与蛋白质交联;裂解细胞、分离染色质,通过超声或酶处理将 染色质随机切割;然后利用抗原抗体的特异性识别反应,将与目的蛋白相结合的 DNA 片段沉淀 下来;解交联释放 DNA 片段并提取 DNA; 构建文库(包括末端修复、加 A 、加接头、长度筛选、 PCR 扩增)并进行高通量测序。对于 ChIP-Seq ,通常留存足够的打断后的染色质,直接用于后续的 DNA 提取,作为对照,称为Input 。也可以使用同型抗体、免疫前血清、或同物种来源 抗体进行模拟的免疫沉淀,所得产物提取 DNA 作为对照,称为 Mock。通常,在建库之前,如果已知一些靶标序列,可以通过 qPCR 等验证 ChIP 的效果。文库构建之后,进行质检,质检过关之后,才进行后续的上机测序。
ChIP-Seq 实验及建库测序流程