本发明属于生物学技术领域，具体涉及一种MeRIP‑seq高通量测序数据的生物分析流程。该分析流程主要包括数据过滤与质量评估、序列比对分析、单样本Peak分析、组内一致性peak检测与注释、motif分析、组间peak合并以及多样本聚类、组间RNA甲基化修饰有或无的差异及组间RNA甲基化率高低的差异分析。本发明提供的生物信息分析流程分析内容全面详实，深度对测序数据进行了挖掘和有效信息的提取，既能满足标准分析需求，也能满足个性化分析需求，有助于更加快速分析大规模的m6A测序数据，适用于植物、动物RNA甲基化修饰的研究。
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＞
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值与log2FC来筛选差异peak，筛选条件为FDR<0.05且|log2FC|>1，并对差异peak相关基因进行GO与KEGG功能富集分析。
2.如权利要求1所述MeRIP-seq高通量测序数据的生物分析方法，其特征在于，步骤S2的比对过程主要包括核糖体比对、参考基因组比对、Reads在染色体上的分布及Reads比对结果可视化，具体操作过程为：(1)核糖体比对：使用短reads比对工具bowtie2将步骤S1中的High quality cleanreads比对到与原始数据相对应物种的核糖体数据库，去除比对上核糖体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功能富集分析。
一种MeRIP-seq高通量测序数据的生物分析流程
技术领域
[0001] 本发明属于生物学技术领域，具体涉及一种MeRIP-seq高通量测序数据的生物分析流程。
背景技术
[0002] 表观遗传学指的是在本身不改变核酸序列的情况下，DNA、RNA和组蛋白通过化学结构修饰，从而影响基因的功能和特性，并通过细胞分裂和增殖遗传给后代。化学修饰是人类基因组丰富的组件，存在于DNA和蛋白质上，其动态变化通过参与基因表达调控，决定细胞的状态，影响细胞的分化和发育。DNA和组蛋白上可逆的化学修饰已经得到了广泛而深入的研究。最新的研究表明，RNA的化学修饰也存在着类似的动态调控机制。
[0003] 目前已知的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甲基化也得到了广泛的关注和研究。
[0004] 目前，市面上利用MeRIP-seq测序技术进行m6A甲基化的研究越来越多，测序数据的深入分析和挖掘也是极为重要的。
发明内容
[0005] 本发明的目的是针对现有技术中普遍存在的缺点，创造性的提供了一种MeRIP-seq高通量测序数据的生物分析流程。本发明提供的生物信息分析流程分析内容全面详实，深度对测序数据进行了挖掘和有效信息的提取，既能满足标准分析需求，也能满足个性化分析需求，有助于更加快速分析大规模的m6A测序数据，适用于植物、动物RNA甲基化修饰的研究。
[0006] 为了达到上述目的，本发明采用的技术方案为：
[0007] 一种MeRIP-seq高通量测序数据的生物分析流程，包括如下步骤：
[0008] S1、数据过滤与质量评估：对下机得到的原始数据利用fastp进行质控，过滤掉低质量的测序数据，得到clean reads，接着，对得到的clean reads进行碱基的组成和质量分布分析，直观展示数据质量情况；
[0009] S2、序列比对分析：将步骤S1过滤得到的clean reads与参考基因组进行比对，确认比对质量合格后，提取比对上唯一位置的序列进行后续信息分析；
[0010] S3、单样本Peak分析：利用MACS2分析软件在全基因组范围内进行peak扫描，阈值为Q value＜0.05.，并对peak在基因组上的位置信息、peak区域序列信息等进行分析，筛选出peak相关基因，并进行基因的功能富集分析；
[0011] S4、组内一致性peak检测与注释：对组内重复样本间的peak进行过滤，保留overlap＞50％的共有peak，进行后续分析；
[0012] S5、将步骤S4保留的overlap＞50％的共有peak进行motif分析；
[0013] S6、组间peak合并以及多样本聚类：利用DiffBind软件进行组间peak合并，得到各处理组间的peak并集，并计算几个peak在各样本中的丰度，并进行样本聚类分析；
[0014] S7、组间RNA甲基化修饰有或无的差异：包括对组间共有或特有peak分析、组间共有和特有peak相关基因分析及组间共有和特有peak的GO/Pathway富集分析；
[0015] S8、组间RNA甲基化率高低的差异分析：对组间peak的甲基化率差异分析、组间差异peak相关基因分析以及组间差异peak相关基因的GO/Pathway富集分析。
[0016] 优选地，步骤S2的比对过程主要包括核糖体比对、参考基因组比对、Reads在染色体上的分布及Reads比对结果可视化，具体操作过程为：
[0017] (1)核糖体比对：使用短reads比对工具bowtie2将步骤S1中的High quality clean reads比对到该物种的核糖体数据库，去除比对上核糖体RNA的reads；
[0018] (2)参考基因组比对：利用比对软件bowtie2将过滤后并去除比对上核糖体RNA的reads比对到参考基因组上，统计各样品比对参考基因组结果，将比对到基因组上唯一位置的reads用于后续的信息分析；
[0019] (3)Reads在染色体上的分布：将唯一比对的且去重复后的reads比对到基因组上各个染色体上，并对密度进行统计作图；
[0020] (4)Reads比对结果可视化：利用Integrative Genomics Viewer将上述步骤中reads在每条染色体比对结果的bigwig文件进行比对结果可视化展示。
[0021] 优选地，步骤S3中对peak在基因组上的分析主要包括peak统计、peak深度分布、peak相关reads相对转录本位置分布统计、peak在染色体上的分布、peak宽度分布、peak富集倍数分布和peak显著程度分布。
[0022] 优选地，步骤S4的具体操作过程包括：
[0023] (1)对组内重复样本间的peak进行过滤，保留overlap>50％的共有peak进行后续分析，若组内有2个生物学重复，则保留的peak必须在2个样本中都出现；若组内有3个生物学重复，则保留的peak必须在至少2个样本中都出现；
[0024] (2)peak相关基因分析：根据peak在基因组上的区域信息及基因的注释信息，得到peak相关基因，了解发生m6A修饰的基因情况，对每组的peak相关基因进行注释，并且对peak相关基因上的peak数目进行统计，分别统计每组样本peak在5’UTR、start codon、CDS、stop codon、3’UTR等基因功能元件上的个数分布；
[0025] (3)进行peak在不同基因功能元件上的富集分析，peak相对基因位置分布分析和peak相关基因GO/Pathway功能富集分析。
[0026] 优选地，步骤S5具体包括以下过程：组间peak合并以及peak丰度计算、PCA分析、相关性热图分析。首先利用DiffBind软件进行组间peak的合并，得到各个处理组间peak的并集，并计算各个peak在各个样本中的丰度。接着，采用降维处理来进行主成分分析(PCA)以及进行相关性热图分析。
[0027] 优选地，步骤S6的具体操作如下：
[0028] (1)组间共有或特有peak分析：将不存在overlap的peak定义为处理组特有的peak，存在overlap的peak为组间共有的peak。统计出组间共有与特有的peak数目，用韦恩图进行展示；
[0029] (2)组间共有和特有peak相关基因分析：对组间共有和特有peak相关基因进行基因注释；
[0030] (3)组间共有和特有peak的GO/Pathway富集分析：对组间共有和特有的peak相关基因进行GO/Pathway功能富集分析。
[0031] 优选地，步骤S8的具体操作过程为：利用MeRIP数据与Input数据，计算每个peak的相对甲基化率，然后用exomePeak软件对比较组的所有peak(两组peak的并集)进行RNA甲基化率差异分析，利用FDR值与log2FC来筛选差异peak，筛选条件为FDR<0.05且|log2FC|>1，并对差异peak相关基因进行GO与KEGG功能富集分析。
[0032] 与现有技术相比，本发明具有如下技术优点：
[0033] (1)本发明针对MeRIP-seq得到的测序数据，开发了完整的生物信息分析流程，分析内容全面详实，深度对测序数据进行了挖掘和有效信息的提取，既能满足标准分析需求，也能满足个性化分析需求；
[0034] (2)本发明为后续更多的科研人员进行基于MeRIP-seq技术进行RNA的m6A甲基化修饰的研究提供了很好的数据分析研究思路，有助于更加快速分析大规模的m6A测序数据，适用于植物、动物RNA甲基化修饰的研究；
[0035] (3)本发明的分析结果以网页版结题报告的形式进行展现，条理清晰，还设置有超链接对各分析项目进行详细的解释说明，有助于用户理解分析报告内容，报告可读性强。
附图说明
[0036] 图1本发明一种MeRIP-seq高通量测序数据生物信息分析流程图；
[0037] 图2本发明实施例中Reads在染色体上的分布图；
[0038] 图3本发明实施例中peak深度分布图；
[0039] 图4本发明实施例中peak在不同功能元件上的统计分布图；
[0040] 图5本发明实施例中peak相对基因位置的分布图；
[0041] 图6本发明实施例中比较组间共有与特有的peak数目维恩图。
具体实施方式
[0042] 下面结合具体实施例对本发明作进一步解释，但是应当注意的是，以下实施例仅用以解释本发明，而不能用来限制本发明，所有与本发明相同或相近的技术方案均在本发明的保护范围之内。若未特别指明，实施例中所用的技术手段为本领域技术人员所熟知的常规手段，所用原料为市售商品。
[0043] 实施例一种MeRIP-seq高通量测序数据生物信息分析流程所述MeRIP-seq高通量测序数据生物信息分析流程，包括如下步骤：
[0044] S1、数据过滤与质量评估：首先，对下机得到的原始数据利用fastp进行质控，过滤掉低质量的测序数据，得到clean reads。接着，对得到的clean reads进行碱基的组成和质量分布分析，直观展示数据质量情况。
[0045] 其中，reads过滤的步骤包括：1)去除含adapter的reads；2)去除含N比例大于10％的reads；3)去除全面都是A碱基的reads；4)去除低质量的reads(质量值Q≤20的碱基数占整条reads的50％以上)。
[0046] S2、序列比对分析。具体步骤如下：
[0047] (1)核糖体比对。使用短reads比对工具bowtie2将步骤S1中的High quality clean reads比对到该物种的核糖体数据库，去除比对上核糖体RNA的reads。
[0048] (2)参考基因组比对。利用比对软件bowtie2将过滤后并去除比对上核糖体RNA的reads比对到参考基因组上，统计各样品比对参考基因组结果，将比对到基因组上唯一位置的reads用于后续的信息分析。
[0049] (3)Reads在染色体上的分布。将唯一比对的且去重复后的reads比对到基因组上各个染色体(分正负链)上，并对密度进行统计作图，如图2所示。
[0050] (4)Reads比对结果可视化。利用Integrative Genomics Viewer(IGV)将上述步骤中reads在每条染色体比对结果的bigwig文件进行比对结果可视化展示。
[0051] S3、单样本peak分析。peak(富集区域)分析是指从全基因组范围内寻找发生m6A修饰的RNA区域和位点。利用MACS2分析软件在全基因组范围进行peak扫描(peak calling)，阈值为Q value<0.05。并对peak在基因组上的位置信息、peak区域序列信息等进行分析，筛选出peak相关基因，并进行基因的功能富集分析。具体包括：
[0052] (1)peak统计：利用MACS2分析软件在全基因组范围进行peak扫描(peak calling)，阈值为Q value<0.05。将各样本中peak数目、总长度、平均长度、总序列深度、平均序列深度、Peak中reads的个数占所有比对上的reads总数的百分比、Peak基因组比对率(％)等信息进行统计，统计结果见表1；
[0053] 表1 样本peak数目统计表
[0054]
[0055] (2)peak深度分布：统计出各样本在不同深度下的peak的数目，以peak区域内的测序深度为横坐标、大于某特定测序深度下的peak数目占总peak的比例为纵坐标画图进行peak深度分布展示，如图3所示；
[0056] (3)peak相关reads相对转录本位置分布统计：以比对后得到的唯一比对的peak相关reads为分析对象，统计其比对到转录本5’UTR,CDS,3’UTR区的数目，得到reads相对转录本位置的分布结果；
[0057] (4)peak在染色体上的分布：对得到的peak在染色体上的分布进行统计，从定位到染色体上的peak个数及分布情况可以反映m6A修饰的分布情况；
[0058] (5)peak宽度分布：以peak宽度为指标，统计peak数目，对peak的宽度分布进行统计；
[0059] (6)peak富集倍数分布：peak富集倍数即为peak的表达量，将该peak区域内的reads count数对样本总reads数进行归一化后，在IP样本与input样本中的比值。Peak的富集倍数越大，表示富集到该peak中的reads数越多，则该peak的表达量越高；
[0060] (7)peak显著程度分布：计算出每个peak的显著性(q value)值，以peak的显著度为指标，统计peak的数目；
[0061] S4、组内一致性peak检测与注释。对组内重复样本间的peak进行过滤，保留overlap＞50％的共有peak进行后续分析。具体步骤包括：
[0062] (1)peak过滤合并：将组内重复样本间的peak进行过滤，保留overlap＞50％的共有peak进行后续分析。
[0063] (2)peak相关基因分析：根据peak在基因组上的区域信息及基因的注释信息，得到peak相关基因，了解发生m6A修饰的基因情况，对每组的peak相关基因进行注释。
[0064] (3)基因上的peak数目统计：计算peak相关基因上peak的数目，可以了解样本基因发生m6A修饰的程度，并且可以比较不同处理组样本的基因发生m6A修饰的变化。统计peak相关基因上peak的数目，并画柱状图展示。
[0065] (4)peak在不同基因功能元件上的数目统计：分别统计每组样本peak在5’UTR、start codon、CDS、stop codon、3’UTR等基因功能元件上的个数分布，统计结果以统计表格的形式进行展示并可绘制饼图，如图4。
[0066] (5)peak在不同基因功能元件上的富集分析：首先计算peak在不同基因功能元件上的enrichment score,enrichment score＝n/(N*p),其中，n:每个基因功能元件上的peak数目；N:样本peak的总数目；p:每种分类占基因组长度的比例。接着，进行卡方检验。最后，将计算出的peak在不同功能元件上的富集结果绘制出柱状图进行结果展示。
[0067] (6)peak相对基因位置分布分析：计算peak在编码基因5’UTR,CDS,3’UTR上的数目及百分比，以组为单位画图展示，分析结果如图5。
[0068] (7)peak相关基因GO/Pathway功能富集分析。
[0069] S5、motif分析：m6A主要富集在mRNA的转录起始位点、终□密码子区，并且有特定的结合序列，通过结合到特定的位置上，从而在基因表达调控中发挥作用。我们利用MEME Suite中的MEME-chip识别peak中显著的motif序列。MEME-chip整合了MEME(检测8～30bp)与Dreme(检测≤8bp)功能，可同时检测长motif和短motif。
[0070] S6、组间peak合并以及多样本聚类。具体操作过程如下：
[0071] (1)组间peak合并以及peak丰度计算：利用DiffBind软件进行组间peak合并，得到各处理组间的peak并集，并计算几个peak在各样本中的丰度。
[0072] (2)PCA分析：
[0073] (3)相关性热图分析：用唯一比对的reads进行两两样本间相关性的计算，具体地是将基因组分成10kb的窗口，统计每个窗口内的reads数目，计算所有窗口上reads数目的皮尔逊相关系数。相关系数越大，说明样品之间RNA甲基化模式的相似度越高。
[0074] S7、组间RNA甲基化修饰有或无的差异。具体包括：首先，对组间共有或特有peak分析：在组间，将不存在overlap的peak定义为处理组特有peak，存在overlap的peak定义为组间共有的peak。统计组间共有与特有的peak数目，用韦恩图来展示，不同比较组间共有和特有的peak,如图6。接着，组间共有和特有peak相关基因分析：对组间共有和特有的peak相关基因进行注释。最后，进行组间共有和特有peak的GO/Pathway富集分析。
[0075] S8、组间RNA甲基化率高低的差异分析。具体步骤如下：
[0076] (1)组间peak的甲基化率差异分析：利用MeRIP数据与Input数据，计算每个peak的相对甲基化率，然后用exomePeak软件对比较组的所有peak(两组peak的并集)进行RNA甲基化率差异分析，利用FDR值与log2FC来筛选差异peak，筛选条件为FDR<0.05且|log2FC|>1。得到组间差异peak,并对数目进行统计，结果如表2所述。
[0077] 表2组间差异peak数目统计表
[0078]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
[0079] (2)组间差异peak相关基因分析：对上述步骤S8.1中筛选出的差异peak相关基因进行注释。
[0080] (3)组间差异peak相关基因的GO/Pathway富集分析。
[0081] 以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。
