本发明属于生物学
技术领域：
，具体涉及一种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进行质控，过滤掉低质量的测序数据，得到cleanreads，接着，对得到的cleanreads进行碱基的组成和质量分布分析，直观展示数据质量情况；s2、序列比对分析：将步骤s1过滤得到的cleanreads与参考基因组进行比对，确认比对质量合格后，提取比对上唯一位置的序列进行后续信息分析；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富集分析。优选地，步骤s2的比对过程主要包括核糖体比对、参考基因组比对、reads在染色体上的分布及reads比对结果可视化，具体操作过程为：(1)核糖体比对：使用短reads比对工具bowtie2将步骤s1中的highqualitycleanreads比对到该物种的核糖体数据库，去除比对上核糖体rna的reads；(2)参考基因组比对：利用比对软件bowtie2将过滤后并去除比对上核糖体rna的reads比对到参考基因组上，统计各样品比对参考基因组结果，将比对到基因组上唯一位置的reads用于后续的信息分析；(3)reads在染色体上的分布：将唯一比对的且去重复后的reads比对到基因组上各个染色体上，并对密度进行统计作图；(4)reads比对结果可视化：利用integrativegenomicsviewer将上述步骤中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、startcodon、cds、stopcodon、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值与log2fc来筛选差异peak，筛选条件为fdr<0.05且|log2fc|>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进行质控，过滤掉低质量的测序数据，得到cleanreads。接着，对得到的cleanreads进行碱基的组成和质量分布分析，直观展示数据质量情况。其中，reads过滤的步骤包括：1)去除含adapter的reads；2)去除含n比例大于10％的reads；3)去除全面都是a碱基的reads；4)去除低质量的reads(质量值q≤20的碱基数占整条reads的50％以上)。s2、序列比对分析。具体步骤如下：(1)核糖体比对。使用短reads比对工具bowtie2将步骤s1中的highqualitycleanreads比对到该物种的核糖体数据库，去除比对上核糖体rna的reads。(2)参考基因组比对。利用比对软件bowtie2将过滤后并去除比对上核糖体rna的reads比对到参考基因组上，统计各样品比对参考基因组结果，将比对到基因组上唯一位置的reads用于后续的信息分析。(3)reads在染色体上的分布。将唯一比对的且去重复后的reads比对到基因组上各个染色体(分正负链)上，并对密度进行统计作图，如图2所示。(4)reads比对结果可视化。利用integrativegenomicsviewer(igv)将上述步骤中reads在每条染色体比对结果的bigwig文件进行比对结果可视化展示。s3、单样本peak分析。peak(富集区域)分析是指从全基因组范围内寻找发生m6a修饰的rna区域和位点。利用macs2分析软件在全基因组范围进行peak扫描(peakcalling)，阈值为qvalue<0.05。并对peak在基因组上的位置信息、peak区域序列信息等进行分析，筛选出peak相关基因，并进行基因的功能富集分析。具体包括：(1)peak统计：利用macs2分析软件在全基因组范围进行peak扫描(peakcalling)，阈值为qvalue<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区域内的readscount数对样本总reads数进行归一化后，在ip样本与input样本中的比值。peak的富集倍数越大，表示富集到该peak中的reads数越多，则该peak的表达量越高；(7)peak显著程度分布：计算出每个peak的显著性(qvalue)值，以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、startcodon、cds、stopcodon、3’utr等基因功能元件上的个数分布，统计结果以统计表格的形式进行展示并可绘制饼图，如图4。(5)peak在不同基因功能元件上的富集分析：首先计算peak在不同基因功能元件上的enrichmentscore,enrichmentscore＝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值与log2fc来筛选差异peak，筛选条件为fdr<0.05且|log2fc|>1。得到组间差异peak,并对数目进行统计，结果如表2所述。表2组间差异peak数目统计表nameupdownallip-vsuv-ip000ip-vs1-ip000ip-vs2-ip000ip-vs3-ip1300130(2)组间差异peak相关基因分析：对上述步骤s8.1中筛选出的差异peak相关基因进行注释。(3)组间差异peak相关基因的go/pathway富集分析。以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。当前第1页1 2 3
技术特征：1.一种merip-seq高通量测序数据的生物分析流程，其特征在于，包括如下步骤：
s1、数据过滤与质量评估：对下机得到的原始数据利用fastp进行质控，过滤掉低质量的测序数据，得到cleanreads，接着，对得到的cleanreads进行碱基的组成和质量分布分析，直观展示数据质量情况；
s2、序列比对分析：将步骤s1过滤得到的cleanreads与参考基因组进行比对，确认比对质量合格后，提取比对上唯一位置的序列进行后续信息分析；
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富集分析。
2.如权利要求1所述merip-seq高通量测序数据的生物分析流程，其特征在于，步骤s2的比对过程主要包括核糖体比对、参考基因组比对、reads在染色体上的分布及reads比对结果可视化，具体操作过程为：
(1)核糖体比对：使用短reads比对工具bowtie2将步骤s1中的highqualitycleanreads比对到该物种的核糖体数据库，去除比对上核糖体rna的reads；
(2)参考基因组比对：利用比对软件bowtie2将过滤后并去除比对上核糖体rna的reads比对到参考基因组上，统计各样品比对参考基因组结果，将比对到基因组上唯一位置的reads用于后续的信息分析；
(3)reads在染色体上的分布：将唯一比对的且去重复后的reads比对到基因组上各个染色体上，并对密度进行统计作图；
(4)reads比对结果可视化：利用integrativegenomicsviewer将上述步骤中reads在每条染色体比对结果的bigwig文件进行比对结果可视化展示。
3.如权利要求1所述merip-seq高通量测序数据的生物分析流程，其特征在于，步骤s3中对peak在基因组上的分析主要包括peak统计、peak深度分布、peak相关reads相对转录本位置分布统计、peak在染色体上的分布、peak宽度分布、peak富集倍数分布和peak显著程度分布。
4.如权利要求1所述merip-seq高通量测序数据的生物分析流程，其特征在于，步骤s4的具体操作过程包括：
(1)对组内重复样本间的peak进行过滤，保留overlap>50％的共有peak进行后续分析，若组内有2个生物学重复，则保留的peak必须在2个样本中都出现；若组内有3个生物学重复，则保留的peak必须在至少2个样本中都出现；
(2)peak相关基因分析：根据peak在基因组上的区域信息及基因的注释信息，得到peak相关基因，了解发生m6a修饰的基因情况，对每组的peak相关基因进行注释，并且对peak相关基因上的peak数目进行统计，分别统计每组样本peak在5’utr、startcodon、cds、stopcodon、3’utr等基因功能元件上的个数分布；(3)进行peak在不同基因功能元件上的富集分析，peak相对基因位置分布分析和peak相关基因go/pathway功能富集分析。
5.如权利要求1所述merip-seq高通量测序数据的生物分析流程，其特征在于，步骤s5具体包括以下过程：组间peak合并以及peak丰度计算、pca分析、相关性热图分析。首先利用diffbind软件进行组间peak的合并，得到各个处理组间peak的并集，并计算各个peak在各个样本中的丰度。接着，采用降维处理来进行主成分分析(pca)以及进行相关性热图分析。
6.如权利要求1所述merip-seq高通量测序数据的生物分析流程，其特征在于，步骤s6的具体操作如下：
(1)组间共有或特有peak分析：将不存在overlap的peak定义为处理组特有的peak，存在overlap的peak为组间共有的peak。统计出组间共有与特有的peak数目，用韦恩图进行展示；
(2)组间共有和特有peak相关基因分析：对组间共有和特有peak相关基因进行基因注释；
(3)组间共有和特有peak的go/pathway富集分析：对组间共有和特有的peak相关基因进行go/pathway功能富集分析。
7.如权利要求1所述merip-seq高通量测序数据的生物分析流程，其特征在于，步骤s8的具体操作过程为：利用merip数据与input数据，计算每个peak的相对甲基化率，然后用exomepeak软件对比较组的所有peak(两组peak的并集)进行rna甲基化率差异分析，利用fdr值与log2fc来筛选差异peak，筛选条件为fdr<0.05且|log2fc|>1，并对差异peak相关基因进行go与kegg功能富集分析。
技术总结本发明属于生物学技术领域，具体涉及一种MeRIP‑seq高通量测序数据的生物分析流程。该分析流程主要包括数据过滤与质量评估、序列比对分析、单样本Peak分析、组内一致性peak检测与注释、motif分析、组间peak合并以及多样本聚类、组间RNA甲基化修饰有或无的差异及组间RNA甲基化率高低的差异分析。本发明提供的生物信息分析流程分析内容全面详实，深度对测序数据进行了挖掘和有效信息的提取，既能满足标准分析需求，也能满足个性化分析需求，有助于更加快速分析大规模的m6A测序数据，适用于植物、动物RNA甲基化修饰的研究。
技术研发人员：夏昊强;周煌凯;高川;陶勇;艾鹏
受保护的技术使用者：广州基迪奥生物科技有限公司
技术研发日：2020.01.17
技术公布日：2020.06.09
