CN111261229B - 一种MeRIP-seq高通量测序数据的生物分析流程 - Google Patents

Abstract

本发明属于生物学技术领域，具体涉及一种MeRIP‑seq高通量测序数据的生物分析流程。该分析流程主要包括数据过滤与质量评估、序列比对分析、单样本Peak分析、组内一致性peak检测与注释、motif分析、组间peak合并以及多样本聚类、组间RNA甲基化修饰有或无的差异及组间RNA甲基化率高低的差异分析。本发明提供的生物信息分析流程分析内容全面详实，深度对测序数据进行了挖掘和有效信息的提取，既能满足标准分析需求，也能满足个性化分析需求，有助于更加快速分析大规模的m6A测序数据，适用于植物、动物RNA甲基化修饰的研究。

Description

一种MeRIP-seq高通量测序数据的生物分析流程

技术领域

本发明属于生物学技术领域，具体涉及一种MeRIP-seq高通量测序数据的生物分析流程。

背景技术

表观遗传学指的是在本身不改变核酸序列的情况下，DNA、RNA和组蛋白通过化学结构修饰，从而影响基因的功能和特性，并通过细胞分裂和增殖遗传给后代。化学修饰是人类基因组丰富的组件，存在于DNA和蛋白质上，其动态变化通过参与基因表达调控，决定细胞的状态，影响细胞的分化和发育。DNA和组蛋白上可逆的化学修饰已经得到了广泛而深入的研究。最新的研究表明，RNA的化学修饰也存在着类似的动态调控机制。

目前已知的RNA转录后修饰形式共有112种，这些修饰广泛分布在古生菌、细菌和真核生物这3大生物中。目前，在所有的RNA修饰中，80％属于或包含甲基化修饰，可见甲基化修饰是一种极为重要的修饰形式。其中，6-甲基腺嘌呤(N6-methyladenosine,m6A)和5-甲基胞嘧啶(C5-methylcytidine,m5C)是最常见的RNA甲基化修饰形式。m6A成为了近年来表观遗传学研究的热点。m6A是发生在碱基A第六位N原子上的甲基化修饰形式，广泛出现在mRNA、tRNA、rRNA以及ncRNA中，其导致了超过80％的RNA碱基甲基化。研究证实，m6A甲基化修饰具有调控基因表达、基因剪切、改变细胞核mRNA输出和mRNA稳定性的功能。随着高通量测序技术的飞速发展，尤其是MeRIP-seq技术的出现，m6A甲基化也得到了广泛的关注和研究。

目前，市面上利用MeRIP-seq测序技术进行m6A甲基化的研究越来越多，测序数据的深入分析和挖掘也是极为重要的。

发明内容

本发明的目的是针对现有技术中普遍存在的缺点，创造性的提供了一种MeRIP-seq高通量测序数据的生物分析流程。本发明提供的生物信息分析流程分析内容全面详实，深度对测序数据进行了挖掘和有效信息的提取，既能满足标准分析需求，也能满足个性化分析需求，有助于更加快速分析大规模的m6A测序数据，适用于植物、动物RNA甲基化修饰的研究。

为了达到上述目的，本发明采用的技术方案为：

一种MeRIP-seq高通量测序数据的生物分析流程，包括如下步骤：

S1、数据过滤与质量评估：对下机得到的原始数据利用fastp进行质控，过滤掉低质量的测序数据，得到clean reads，接着，对得到的clean reads进行碱基的组成和质量分布分析，直观展示数据质量情况；

S2、序列比对分析：将步骤S1过滤得到的clean reads与参考基因组进行比对，确认比对质量合格后，提取比对上唯一位置的序列进行后续信息分析；

S3、单样本Peak分析：利用MACS2分析软件在全基因组范围内进行peak扫描，阈值为Q value＜0.05.，并对peak在基因组上的位置信息、peak区域序列信息等进行分析，筛选出peak相关基因，并进行基因的功能富集分析；

S4、组内一致性peak检测与注释：对组内重复样本间的peak进行过滤，保留overlap＞50％的共有peak，进行后续分析；

S5、将步骤S4保留的overlap＞50％的共有peak进行motif分析；

S6、组间peak合并以及多样本聚类：利用DiffBind软件进行组间peak合并，得到各处理组间的peak并集，并计算几个peak在各样本中的丰度，并进行样本聚类分析；

S7、组间RNA甲基化修饰有或无的差异：包括对组间共有或特有peak分析、组间共有和特有peak相关基因分析及组间共有和特有peak的GO/Pathway富集分析；

S8、组间RNA甲基化率高低的差异分析：对组间peak的甲基化率差异分析、组间差异peak相关基因分析以及组间差异peak相关基因的GO/Pathway富集分析。

优选地，步骤S2的比对过程主要包括核糖体比对、参考基因组比对、Reads在染色体上的分布及Reads比对结果可视化，具体操作过程为：

(1)核糖体比对：使用短reads比对工具bowtie2将步骤S1中的High qualityclean reads比对到该物种的核糖体数据库，去除比对上核糖体RNA的reads；

(2)参考基因组比对：利用比对软件bowtie2将过滤后并去除比对上核糖体RNA的reads比对到参考基因组上，统计各样品比对参考基因组结果，将比对到基因组上唯一位置的reads用于后续的信息分析；

(3)Reads在染色体上的分布：将唯一比对的且去重复后的reads比对到基因组上各个染色体上，并对密度进行统计作图；

(4)Reads比对结果可视化：利用Integrative Genomics Viewer将上述步骤中reads在每条染色体比对结果的bigwig文件进行比对结果可视化展示。

优选地，步骤S3中对peak在基因组上的分析主要包括peak统计、peak深度分布、peak相关reads相对转录本位置分布统计、peak在染色体上的分布、peak宽度分布、peak富集倍数分布和peak显著程度分布。

优选地，步骤S4的具体操作过程包括：

(1)对组内重复样本间的peak进行过滤，保留overlap>50％的共有peak进行后续分析，若组内有2个生物学重复，则保留的peak必须在2个样本中都出现；若组内有3个生物学重复，则保留的peak必须在至少2个样本中都出现；

(2)peak相关基因分析：根据peak在基因组上的区域信息及基因的注释信息，得到peak相关基因，了解发生m6A修饰的基因情况，对每组的peak相关基因进行注释，并且对peak相关基因上的peak数目进行统计，分别统计每组样本peak在5’UTR、start codon、CDS、stop codon、3’UTR等基因功能元件上的个数分布；

(3)进行peak在不同基因功能元件上的富集分析，peak相对基因位置分布分析和peak相关基因GO/Pathway功能富集分析。

优选地，步骤S5具体包括以下过程：组间peak合并以及peak丰度计算、PCA分析、相关性热图分析。首先利用DiffBind软件进行组间peak的合并，得到各个处理组间peak的并集，并计算各个peak在各个样本中的丰度。接着，采用降维处理来进行主成分分析(PCA)以及进行相关性热图分析。

优选地，步骤S6的具体操作如下：

(1)组间共有或特有peak分析：将不存在overlap的peak定义为处理组特有的peak，存在overlap的peak为组间共有的peak。统计出组间共有与特有的peak数目，用韦恩图进行展示；

(2)组间共有和特有peak相关基因分析：对组间共有和特有peak相关基因进行基因注释；

(3)组间共有和特有peak的GO/Pathway富集分析：对组间共有和特有的peak相关基因进行GO/Pathway功能富集分析。

优选地，步骤S8的具体操作过程为：利用MeRIP数据与Input数据，计算每个peak的相对甲基化率，然后用exomePeak软件对比较组的所有peak(两组peak的并集)进行RNA甲基化率差异分析，利用FDR值与log₂FC来筛选差异peak，筛选条件为FDR<0.05且|log₂FC|>1，并对差异peak相关基因进行GO与KEGG功能富集分析。

与现有技术相比，本发明具有如下技术优点：

(1)本发明针对MeRIP-seq得到的测序数据，开发了完整的生物信息分析流程，分析内容全面详实，深度对测序数据进行了挖掘和有效信息的提取，既能满足标准分析需求，也能满足个性化分析需求；

(2)本发明为后续更多的科研人员进行基于MeRIP-seq技术进行RNA的m6A甲基化修饰的研究提供了很好的数据分析研究思路，有助于更加快速分析大规模的m6A测序数据，适用于植物、动物RNA甲基化修饰的研究；

(3)本发明的分析结果以网页版结题报告的形式进行展现，条理清晰，还设置有超链接对各分析项目进行详细的解释说明，有助于用户理解分析报告内容，报告可读性强。

附图说明

图1本发明一种MeRIP-seq高通量测序数据生物信息分析流程图；

图2本发明实施例中Reads在染色体上的分布图；

图3本发明实施例中peak深度分布图；

图4本发明实施例中peak在不同功能元件上的统计分布图；

图5本发明实施例中peak相对基因位置的分布图；

图6本发明实施例中比较组间共有与特有的peak数目维恩图。

具体实施方式

下面结合具体实施例对本发明作进一步解释，但是应当注意的是，以下实施例仅用以解释本发明，而不能用来限制本发明，所有与本发明相同或相近的技术方案均在本发明的保护范围之内。若未特别指明，实施例中所用的技术手段为本领域技术人员所熟知的常规手段，所用原料为市售商品。

实施例一种MeRIP-seq高通量测序数据生物信息分析流程所述MeRIP-seq高通量测序数据生物信息分析流程，包括如下步骤：

S1、数据过滤与质量评估：首先，对下机得到的原始数据利用fastp进行质控，过滤掉低质量的测序数据，得到clean reads。接着，对得到的clean reads进行碱基的组成和质量分布分析，直观展示数据质量情况。

其中，reads过滤的步骤包括：1)去除含adapter的reads；2)去除含N比例大于10％的reads；3)去除全面都是A碱基的reads；4)去除低质量的reads(质量值Q≤20的碱基数占整条reads的50％以上)。

S2、序列比对分析。具体步骤如下：

(1)核糖体比对。使用短reads比对工具bowtie2将步骤S1中的High qualityclean reads比对到该物种的核糖体数据库，去除比对上核糖体RNA的reads。

(2)参考基因组比对。利用比对软件bowtie2将过滤后并去除比对上核糖体RNA的reads比对到参考基因组上，统计各样品比对参考基因组结果，将比对到基因组上唯一位置的reads用于后续的信息分析。

(3)Reads在染色体上的分布。将唯一比对的且去重复后的reads比对到基因组上各个染色体(分正负链)上，并对密度进行统计作图，如图2所示。

(4)Reads比对结果可视化。利用Integrative Genomics Viewer(IGV)将上述步骤中reads在每条染色体比对结果的bigwig文件进行比对结果可视化展示。

S3、单样本peak分析。peak(富集区域)分析是指从全基因组范围内寻找发生m6A修饰的RNA区域和位点。利用MACS2分析软件在全基因组范围进行peak扫描(peak calling)，阈值为Q value<0.05。并对peak在基因组上的位置信息、peak区域序列信息等进行分析，筛选出peak相关基因，并进行基因的功能富集分析。具体包括：

(1)peak统计：利用MACS2分析软件在全基因组范围进行peak扫描(peakcalling)，阈值为Q value<0.05。将各样本中peak数目、总长度、平均长度、总序列深度、平均序列深度、Peak中reads的个数占所有比对上的reads总数的百分比、Peak基因组比对率(％)等信息进行统计，统计结果见表1；

表1 样本peak数目统计表

(2)peak深度分布：统计出各样本在不同深度下的peak的数目，以peak区域内的测序深度为横坐标、大于某特定测序深度下的peak数目占总peak的比例为纵坐标画图进行peak深度分布展示，如图3所示；

(3)peak相关reads相对转录本位置分布统计：以比对后得到的唯一比对的peak相关reads为分析对象，统计其比对到转录本5’UTR,CDS,3’UTR区的数目，得到reads相对转录本位置的分布结果；

(4)peak在染色体上的分布：对得到的peak在染色体上的分布进行统计，从定位到染色体上的peak个数及分布情况可以反映m6A修饰的分布情况；

(5)peak宽度分布：以peak宽度为指标，统计peak数目，对peak的宽度分布进行统计；

(6)peak富集倍数分布：peak富集倍数即为peak的表达量，将该peak区域内的reads count数对样本总reads数进行归一化后，在IP样本与input样本中的比值。Peak的富集倍数越大，表示富集到该peak中的reads数越多，则该peak的表达量越高；

(7)peak显著程度分布：计算出每个peak的显著性(q value)值，以peak的显著度为指标，统计peak的数目；

S4、组内一致性peak检测与注释。对组内重复样本间的peak进行过滤，保留overlap＞50％的共有peak进行后续分析。具体步骤包括：

(1)peak过滤合并：将组内重复样本间的peak进行过滤，保留overlap＞50％的共有peak进行后续分析。

(2)peak相关基因分析：根据peak在基因组上的区域信息及基因的注释信息，得到peak相关基因，了解发生m6A修饰的基因情况，对每组的peak相关基因进行注释。

(3)基因上的peak数目统计：计算peak相关基因上peak的数目，可以了解样本基因发生m6A修饰的程度，并且可以比较不同处理组样本的基因发生m6A修饰的变化。统计peak相关基因上peak的数目，并画柱状图展示。

(4)peak在不同基因功能元件上的数目统计：分别统计每组样本peak在5’UTR、start codon、CDS、stop codon、3’UTR等基因功能元件上的个数分布，统计结果以统计表格的形式进行展示并可绘制饼图，如图4。

(5)peak在不同基因功能元件上的富集分析：首先计算peak在不同基因功能元件上的enrichment score,enrichment score＝n/(N*p),其中，n:每个基因功能元件上的peak数目；N:样本peak的总数目；p:每种分类占基因组长度的比例。接着，进行卡方检验。最后，将计算出的peak在不同功能元件上的富集结果绘制出柱状图进行结果展示。

(6)peak相对基因位置分布分析：计算peak在编码基因5’UTR,CDS,3’UTR上的数目及百分比，以组为单位画图展示，分析结果如图5。

(7)peak相关基因GO/Pathway功能富集分析。

S5、motif分析：m6A主要富集在mRNA的转录起始位点、终□密码子区，并且有特定的结合序列，通过结合到特定的位置上，从而在基因表达调控中发挥作用。我们利用MEMESuite中的MEME-chip识别peak中显著的motif序列。MEME-chip整合了MEME(检测8～30bp)与Dreme(检测≤8bp)功能，可同时检测长motif和短motif。

S6、组间peak合并以及多样本聚类。具体操作过程如下：

(1)组间peak合并以及peak丰度计算：利用DiffBind软件进行组间peak合并，得到各处理组间的peak并集，并计算几个peak在各样本中的丰度。

(2)PCA分析：

(3)相关性热图分析：用唯一比对的reads进行两两样本间相关性的计算，具体地是将基因组分成10kb的窗口，统计每个窗口内的reads数目，计算所有窗口上reads数目的皮尔逊相关系数。相关系数越大，说明样品之间RNA甲基化模式的相似度越高。

S7、组间RNA甲基化修饰有或无的差异。具体包括：首先，对组间共有或特有peak分析：在组间，将不存在overlap的peak定义为处理组特有peak，存在overlap的peak定义为组间共有的peak。统计组间共有与特有的peak数目，用韦恩图来展示，不同比较组间共有和特有的peak,如图6。接着，组间共有和特有peak相关基因分析：对组间共有和特有的peak相关基因进行注释。最后，进行组间共有和特有peak的GO/Pathway富集分析。

S8、组间RNA甲基化率高低的差异分析。具体步骤如下：

(1)组间peak的甲基化率差异分析：利用MeRIP数据与Input数据，计算每个peak的相对甲基化率，然后用exomePeak软件对比较组的所有peak(两组peak的并集)进行RNA甲基化率差异分析，利用FDR值与log₂FC来筛选差异peak，筛选条件为FDR<0.05且|log₂FC|>1。得到组间差异peak,并对数目进行统计，结果如表2所述。

表2组间差异peak数目统计表

Name	Up	Down	All
				IP-vs UV-IP	0	0	0
IP-vs 1-IP	0	0	0
				IP-vs 2-IP	0	0	0
IP-vs 3-IP	130	0	130

(2)组间差异peak相关基因分析：对上述步骤S8.1中筛选出的差异peak相关基因进行注释。

(3)组间差异peak相关基因的GO/Pathway富集分析。

以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

Claims (5)

1.一种MeRIP-seq高通量测序数据的生物分析方法，其特征在于，包括如下步骤：

S1、数据过滤与质量评估：对下机得到的原始数据利用fastp进行质控，过滤掉低质量的测序数据，得到clean reads，接着，对得到的clean reads进行碱基的组成和质量分布分析，直观展示数据质量情况；

S2、序列比对分析：将步骤S1过滤得到的clean reads与参考基因组进行比对，确认比对质量合格后，提取比对上唯一位置的序列进行后续信息分析；

S3、单样本Peak分析：利用MACS2分析软件在全基因组范围内进行peak扫描，阈值为Qvalue＜0.05，并对peak在基因组上的位置信息、peak区域序列信息进行分析，筛选出peak相关基因，并进行基因的功能富集分析；

S4、组内一致性peak检测与注释：对组内重复样本间的peak进行过滤，保留overlap＞50％的共有peak，进行后续分析；

S5、将步骤S4保留的overlap＞50％的共有peak进行motif分析；

S6、组间peak合并以及多样本聚类：利用DiffBind软件进行组间peak合并，得到各处理组间的peak并集，并计算几个peak在各样本中的丰度，并进行样本聚类分析；

S7、组间RNA甲基化修饰有或无的差异：包括对组间共有或特有peak分析、组间共有和特有peak相关基因分析及组间共有和特有peak的GO/Pathway富集分析；

S8、组间RNA甲基化率高低的差异分析：对组间peak的甲基化率差异分析、组间差异peak相关基因分析以及组间差异peak相关基因的GO/Pathway富集分析；步骤S4的具体操作过程包括：

(1)对组内重复样本间的peak进行过滤，保留overlap>50％的共有peak进行后续分析，若组内有2个生物学重复，则保留的peak必须在2个样本中都出现；若组内有3个生物学重复，则保留的peak必须在至少2个样本中都出现；

(2)peak相关基因分析：根据peak在基因组上的区域信息及基因的注释信息，得到peak相关基因，了解发生m6A修饰的基因情况，对每组的peak相关基因进行注释，并且对peak相关基因上的peak数目进行统计，分别统计每组样本peak在5’UTR、start codon、CDS、stopcodon、3’UTR基因功能元件上的个数分布；(3)进行peak在不同基因功能元件上的富集分析，peak相对基因位置分布分析和peak相关基因GO/Pathway功能富集分析；步骤S8的具体操作过程为：利用MeRIP数据与Input数据，计算每个peak的相对甲基化率，然后用exomePeak软件对比较组的所有peak，进行RNA甲基化率差异分析，利用FDR值与log₂FC来筛选差异peak，筛选条件为FDR<0.05且|log₂FC|>1，并对差异peak相关基因进行GO与KEGG功能富集分析。

2.如权利要求1所述MeRIP-seq高通量测序数据的生物分析方法，其特征在于，步骤S2的比对过程主要包括核糖体比对、参考基因组比对、Reads在染色体上的分布及Reads比对结果可视化，具体操作过程为：

(1)核糖体比对：使用短reads比对工具bowtie2将步骤S1中的High qualitycleanreads比对到与原始数据相对应物种的核糖体数据库，去除比对上核糖体RNA的reads；

(2)参考基因组比对：利用比对软件bowtie2将过滤后并去除比对上核糖体RNA的reads比对到参考基因组上，统计各样品比对参考基因组结果，将比对到基因组上唯一位置的reads用于后续的信息分析；

(3)Reads在染色体上的分布：将唯一比对的且去重复后的reads比对到基因组上各个染色体上，并对密度进行统计作图；

(4)Reads比对结果可视化：利用Integrative Genomics Viewer将上述步骤中reads在每条染色体比对结果的bigwig文件进行比对结果可视化展示。

3.如权利要求1所述MeRIP-seq高通量测序数据的生物分析方法，其特征在于，步骤S3中对peak在基因组上的分析主要包括peak统计、peak深度分布、peak相关reads相对转录本位置分布统计、peak在染色体上的分布、peak宽度分布、peak富集倍数分布和peak显著程度分布。

4.如权利要求1所述MeRIP-seq高通量测序数据的生物分析方法，其特征在于，步骤S5具体包括以下过程：组间peak合并以及peak丰度计算、PCA分析、相关性热图分析；首先利用DiffBind软件进行组间peak的合并，得到各个处理组间peak的并集，并计算各个peak在各个样本中的丰度；接着，采用降维处理来进行PCA分析以及进行相关性热图分析。

5.如权利要求1所述MeRIP-seq高通量测序数据的生物分析方法，其特征在于，步骤S6的具体操作如下：

(1)组间共有或特有peak分析：将不存在overlap的peak定义为处理组特有的peak，存在overlap的peak为组间共有的peak，统计出组间共有与特有的peak数目，用韦恩图进行展示；

(2)组间共有和特有peak相关基因分析：对组间共有和特有peak相关基因进行基因注释；

(3)组间共有和特有peak的GO/Pathway富集分析：对组间共有和特有的peak相关基因进行GO/Pathway功能富集分析。

Images