CN118016159A - 基于drs的生物遗传样本转录组测序分析方法 - Google Patents

Abstract

本发明公开了一种基于DRS的生物遗传样本转录组测序分析方法，涉及基因测序技术领域。本方法包括富集mRNA建立直接测序文库；进行二代测序和三代DRS测序；对测序的数据进行质控和纠错校正；纠错的转录本序列与参考基因比对统计，得到新转录本；对新转录本进行功能注释和编码区预测；优化转录本的结构，进行可变剪切分析、融合转录本分析、lncRNA预测、转录本表达定量分析、差异表达和功能富集分析、poly(A)分析、甲基化分析和假尿苷分析。本方法基于nextflow与docker搭建了一体化分析流程，优化了从新转录本识别到定量差异分析的内容，进一步发挥了RNA直接测序的优势；新增了多种分析过程，并与表达关联，对转录本表达进行系统多角度阐释。

Description

基于DRS的生物遗传样本转录组测序分析方法

技术领域

本发明涉及基因测序技术领域，特别涉及基于DRS的生物遗传样本转录组测序分析方法。

背景技术

对生物样本中的核糖核酸(RNA)进行测序可以获得丰富的信息，包括细菌和病毒的身份、选择性剪接的细微差别或生物体的转录状态。目前大多数的RNA测序技术只能产生相对较短的读取长度，多数方法仍需将RNA反转录成cDNA(互补脱氧核糖核酸)，这不仅会引入错误及偏好，而且效率不高，还会影响RNA自身复杂结构的准确表征。传统的第二代高通量测序平台在分析转录组时存在一些限制，例如无法准确得到或组装出完整的转录本，难以区分和定量同源基因、超基因家族和等位基因表达的转录本，同时也无法直接获得转录本上的甲基化修饰信息和poly(A)尾长信息，从而限制了对生命活动的深层次理解。

为了克服这些限制，Direct RNA Sequencing(DRS)技术应运而生。DRS是一种基于牛津纳米孔公司(Oxford Nanopore Technologies，ONT)的第三代测序平台的直接RNA测序方法。与传统方法不同，DRS不需要对RNA进行反转录或扩增，可以避免引入假阳性转录本。DRS可以检测单碱基的甲基化修饰(m5C、m6A、假尿苷等)及其修饰率、Poly(A)尾长、可变剪切、融合基因等信息，反映RNA的最原始信息，表征最全面的RNA特征。DRS一次测序，原始数据可终身使用，随着分析算法的更新，后续可用原始数据分析其他类型的甲基化修饰(m1A、2'-O-甲基化等)。

现有的基于DRS的转录组三代测序方法，由于读长的限制，不能有效地实现较大RNA分子的全长结构分析，如信使RNA和核糖体RNA。需要PCR扩增也意味着通常不能捕获RNA分子上的原生修饰，不能在单分子水平上分析RNA结构。

发明内容

针对现有技术中存在的上述问题，本发明提供了一种基于DRS的生物遗传样本转录组测序分析方法，利用Nextflow与docker软件串行了一套完整、可移植的DRS分析技术，能够更准确、更全面地分析RNA特征，深入挖掘生命过程中的转录信息，从而为生命科学研究、医学诊断和药物开发等领域提供更多有价值的数据和应用。具体通过以下技术实现。

本发明提供的基于DRS的生物遗传样本转录组测序分析方法，步骤包括：富集样本mRNA，建立直接测序文库；

基于Illumina进行二代测序，基于ONT进行三代DRS测序，分别得到二代原始测序数据和三代原始测序数据；

对所述二代原始测序数据和三代原始测序数据分别进行质控，基于所述二代原始测序数据对所述三代原始测序数据进行纠错校正；将纠错后的转录本序列与参考基因组进行比对，统计得到一致性序列；

过滤所述旧转录本中的全长reads，然后将一致性序列与参考基因组的已知转录本进行比对，得到新转录本和新基因；

识别所述新转录本的潜在编码区序列，预测所述新转录本和新基因的编码区序列，对所述新转录本和新基因进行功能注释；

对转录本的结构进行优化，进行可变剪切分析、融合转录本分析、lncRNA预测、转录本表达定量分析、差异表达和功能富集分析、poly(A)分析、甲基化分析和假尿苷分析。

进一步地，对所述二代原始测序数据的质控方式为：

去除N碱基含量＞5％的reads；

并且，去除质量值≤5，碱基数目达到50％的reads；

并且，去除有adapter污染的reads；

并且，去除PCR扩增造成的重复序列；

对所述三代原始测序数据的质控方式为：过滤Q＜9且长度＜100bp的reads。进一步地，基于所述二代原始测序数据，使用LRECE软件中的run_correction_tools.sh脚本对所述三代原始测序数据进行纠错校正；

使用minimap2软件将纠错校正后的转录本序列与参考转录本进行比对，使用samtools软件统计比对结果，得到一致性序列。

进一步地，使用LAFITE软件在所述一致性序列中引入poly(A)鉴定标签过滤去除全长reads；使用gffcompare软件将所述一致性序列与基因组的已知转录本进行比对，发现新转录本和新基因。

进一步地，使用TransDecoder软件对所述新转录本和新基因进行潜在编码区序列识别，使用orfpy软件对所述新转录本和新基因的编码区进行预测。

进一步地，对所述新转录本和新基因进行Nr、Pfam、Uniprot、KEGG、GO、KOG/COG和PATHWAY七个数据库的转录本功能注释。

更进一步地，功能注释的方法包括：

使用diamond blastp软件将新转录本中转录本编码的蛋白序列与现有蛋白质数据库uniprot和Nr进行比对，获得序列的功能信息，以及蛋白可能参与的代谢通路信息；基于数据库之间的关联，得到KOG/COG注释结果，进行KOG/COG的分类统计及绘图；

使用Pfam数据库和hmmscan软件进行所述新转录本的结构域预测；使用KOfam数据库和kofam_scan软件进行所述新转录本的同源性搜索。

进一步地，使用gffcompare软件，将一致性序列与基因组已知转录本进行比较；如果存在转录本在原有转录本边界之外的区域，则将转录本的非翻译区向上下游延伸，修正转录本的边界；通过比较结果，对基因的5’端或3’端进行延长，完成对转录本的结构优化。

进一步地，使用LASER软件完成转录本的可变剪切分析；

使用CTAT-LR-Fusion软件比对、寻找融合转录本；

使用CNCI软件、CPC2软件和PLEK软件对所述新转录本进行lncRNA预测；

采用TPM作为衡量表达水平的指标，使用bamboo软件进行所述新转录本表达定量分析；

使用DESeq2软件进行差异表达分析，筛选阈值为padj<0.05，且|log2FoldChange|>1；若显著差异转录本数目＜10时，筛选阈值为pvalue<0.05，且|log2FoldChange|>1；基于Gene Ontology和KEGG Pathway方法，使用clusterProfiler软件进行功能富集分析；

使用NanoPolish软件对有效reads的poly(A)进行计算，使用minimap2软件将reads比对到参考基因组序列后，从bam文件中提取比对到参考基因组的终点位置，再使用Quantifypoly(A)软件进行poly(A)位点的鉴定、聚类与注释；

使用Tombo软件预测出RNA分子序列中m5C修饰位点，使用m6Anet软件预测出RNA分子序列中m6A修饰位点，使用R语言计算并绘制甲基化位点的位置、分布、5bp的motif图，使用methylkit软件进行差异甲基化位点分析；若每组有若干样本则采用Logistic回归检验；

使用Nanopsu软件获得RNA分子序列中的假尿苷修饰的P值，基于P>0.9进行P值筛选。

与现有技术相比，本发明的有益之处在于：

1、在转录组测序分析过程中，本发明在原有的三代转录组分析基础上，基于nextflow与docker搭建了一套一体化的分析流程，使得DRS分析具有更优异的可移植、可重复、便携的特点；使用LAFITE、bamboo、LASER等准确性更高的软件优化了从新转录本识别到定量分析、可变剪切、差异分析等分析内容，使其更充分利用了DRS的特点，更好的发挥RNA直接测序的特长；自主研发了poly(A)与表达关联与可视化、甲基化分布与motif的可视化及差异分析、假尿苷的鉴定分布与差异分析与可视化等脚本，系统化、多角度地对转录本的表达进行阐释。

2、本发明基于DRS进行测序分析，整个分析过程无需进行反转录和扩增，无测序偏好性，有着储存并同时检测RNA上的甲基化修饰位点信息，准确分析可变剪接、融合基因和鉴定新异构体的特点；还可对poly(A)尾长进行相对准确的估算，还原真实RNA特征，并且实现转录本的表达水平准确定量；

3、本发明在样本mRNA富集和反转录引物连接环节添加了RNA酶抑制剂，有效的降低了RNA在建库过程中发生降解的情况；选用T4 DNA ligase提高了单链RNA和双链接头的连接效率，使得测序所得数据更接近样本转录本的真实表达情况。

附图说明

图1为三代Direct RNA文库构建及测序流程图；

图2为二代建库流程图；

图3为CDS的长度分布图；

图4为主要的转录本可变剪接类型示意图；

图5为鉴定的差异表达转录本数目；

图6为假尿苷单样本位点venn图。

具体实施方式

下面将对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动条件下所获得的所有其它实施例，都属于本发明保护的范围。

以下具体上述方式提供的转录组测序分析方法，具体步骤如下所示。

一、实验前准备和样本准备

本具体实施方式采用超低温冻存常规动植物RNA核酸样本。具体地，样本为小鼠(Mus musculus)的肝脏RNA样本和水稻(Oryza sativa)的叶片RNA样本，通过商业购买原材料，由公司实验人员自行提取RNA。

1、清理擦拭超净工作台，紫外消杀、打冰、擦拭桌面与金属浴；

2、取出富集试剂盒，将富集磁珠平衡至室温；

3、将样本取出，置于冰盒上；

4、待样本融化后按建库起始量要求取样，加入1μl RiboLock RNase Inhibitor(RNA酶抑制剂)，置于冰上备用。

二、富集磁珠准备

1、取磁珠平衡至室温，取磁珠于1.5μl离心管中；

2、将离心管置于磁力架上处理，澄清后弃去上清；加入RNA Binding Buffer，取下离心管，重悬磁珠；重复本步骤一次；

3、将离心管至于磁力架上，弃去上清，加入RNA Binding Buffer再次重悬磁珠，备用。

三、mRNA富集

1、将备用的样本与磁珠混合，移液器吹打混匀；65℃处理6min后冰浴迅速降温至4℃；恢复至室温后重悬磁珠，静置20min，每5min混匀磁珠一次；

2、将离心管置于磁力架上，弃去上清；使用wash buffer清洗磁珠两次(需重悬磁珠)；弃去上清后取下离心管，瞬时离心后重新上架，吸取弃去多余废液；

3、加入tris buffer，移液器吹打混匀；80℃处理2min，然后室温静置10min；

4、加入RNA Binding Buffer，移液器吹打混匀；静置20min，每5min混匀磁珠一次；

5、将离心管置于磁力架上，弃去上清；使用wash buffer清洗磁珠一次(需重悬磁珠)；弃去上清后取下离心管，瞬时离心后重新上架，吸取弃去多余废液；

6、加入规定体积tris buffer，移液器吹打混匀；80℃处理2min，然后室温静置10min；迅速上架吸出上清至新的离心管；

7、吸取1μl做qubit检测(荧光染料法)，其余样本置于冰上备用。

四、建库

1、在第三步富集的样本中加入下表1所示的试剂，24℃反应30min，每15min轻弹混匀一次；备用。

表1

NEBNext Quick Ligation Reaction Buffer	3μL
		RNA CS(RCS)，110nM	0.5μL
T4 DNA Ligase(T4 DNA连接酶)	1.5μL
		RT Adapter(RTA)	1μL
Total volume	15μL
		RiboLock RNase Inhibitor	1μL

2、取新的1.5ml离心管，按下表2所示的配方配置Mix。

表2

Nuclease-free water	9μL
		10mM dNTPs	2μL
5×first-strand buffer	8μL
		0.1M DTT	4μL

3、在第1步的样本液中加入第2步所得Mix，轻弹混匀后加入2μl SuperScript IIIreverse transcriptase，混匀，离心，置于金属浴上；

4、50℃反应60min，70℃反应10min，转4℃；

5、取出磁珠恢复至室温；向上述反应体系中加入磁珠，纯化一次；

6、向所得样本中加入下表3所示的试剂，24℃反应10min，轻弹混匀一次。

表3

7、从金属浴取出样本加入磁珠，纯化一次；将所得产物置于冰上，准备上机loading芯片操作。

五、DRS测序(三代试验测序)

使用ONT测序芯片进行三代DRS测序，得到三代原始测序数据。三代Direct RNA文库构建及测序流程如图1所示。测序序列数据量统计结果如下表4所示。

表4三代原始测序数据量统计

六、二代试验测序

如图2所示，RNA样品检测合格后，使用Oligo dT磁珠富集mRNA，然后进行片段化，之后反转录cDNA第一链以及第二链，再经末端修复、加A尾、加测序接头、纯化、PCR扩增等步骤完成整个文库制备工作。

文库检测合格后，按照有效浓度及目标下机数据量的需求将不同文库pooling至Flow cell(即合并到流动池)，cBOT成簇后使用Illumina高通量测序平台IlluminaNovaSeq进行测序，得到二代原始测序数据。

测序序列数据量统计结果如下表5所示。

表5二代原始测序数据量统计

七、测序数据的质量控制

Nanopore测序的下机数据的原始数据格式为包含所有原始测序信号的fast5格式，通过GUPPY(版本：5.0.16)软件进行base calling后将fast5格式数据转换为fastq格式。

由于原始测序数据可能包含低质量序列、接头序列等，为了保证信息分析结果的可靠性，原始测序数据需要经过过滤，后得到高质量、长度长的有效reads，仍然以FASTQ格式存储(该序列文件可直接用于发表，公共数据库提交等，后续分析都基于此数据)，过滤标准如下：

1、按照以下标准过滤二代原始测序数据：

(1)去除N碱基含量超过5％的reads；

(2)去除低质量(质量值小于等于5)碱基数目达到50％的reads；

(3)去除有adapter污染的reads；

(4)去除PCR扩增造成的重复序列。

2、按照以下标准过滤三代原始测序数据：过滤Q小于9，且长度＜100bp的reads。

八、二代数据纠错及比对参考基因组

基于二代Illumina测序数据使用LRECE中的run_correction_tools.sh脚本(参数：default)对三代有效数据数据进行纠错校正，提升数据质量。

根据minimap2(版本：2.17-r941；参数：-ax splice-uf-k14)将纠错后的转录本序列与参考基因比对，使用samtools(版本：1.11；参数：flagstat)进行比对结果统计，得到一致性序列。如下表6所示。

表6纠错后比对结果

九、一致性转录本

使用LAFITE(版本：1.0.1；参数：-t 20)软件引入了poly(A)鉴定标签来过滤全长reads。基于表6的比对结果得到一致性转录本序列(新转录本)，样品的一致性序列信息统计如下表7所示。

表7一致性转录本序列统计

使用LAFITE软件引入了poly(A)鉴定标签来过滤全长reads，结合二代有效数据比对结果使junction更精确，使得低丰度全长转录本的得以实现。

利用gffcompare(版本：0.12.1；参数：-R-C-K-M)将一致性序列与基因组的已知转录本进行比较，发现新的转录本和新基因。新转录本数量和结果如下表8所示。

表8新转录本数量和结果如下表

使用TransDecoder(版本：5.5.0；参数：-m 50，-single_best_only)软件对新转录本和新基因进行潜在编码区序列(Coding Sequence，CDS)识别。

TransDecoder基于开放阅读框(Open Reading Frame，ORF)长度、对数似然函数值(Log-likelihood Score)等信息，能够从新转录本序列中识别可靠的潜在编码区序列。

使用orfpy(版本：0.0.4；参数：--strand f--start ATG)软件对新转录本的编码区进行预测，取正链上最长的ORF作为新转录本的编码区进行后续分析和预测。CDS的长度分布图如图3所示；图中，横坐标表示reads长度；纵坐标表示该长度范围内的reads数目，其中红色的虚线表示N50的长度。

为获得新转录本和新基因的全面的功能信息，对新转录本进行以下七个数据库的新转录本功能注释，包括：Nr、Pfam、Uniprot、KEGG、GO(Gene Ontology)、KOG/COG和PATHWAY。新转录本的功能注释结果统计如下表9所示。

表9新转录本的功能注释结果统计

功能注释主要有如下2种方法：

1、序列相似性搜索

将全长转录本中转录本编码的蛋白序列与现有蛋白质数据库Uniprot、Nr进行diamond blastp(版本：2.0.6.144；参数：-evalue 1e-5)比对，获得序列的功能信息，以及蛋白可能参与的代谢通路信息。

Uniprot数据库记录了每个蛋白质家族与Gene Ontology中的功能节点的对应关系，通过此系统预测转录本编码的蛋白序列所执行的生物学功能。

基于数据库之间的关联，得到KOG/COG注释结果，进行KOG/COG的分类统计及绘图。

2、结构域相似性搜索：

蛋白质中，一般由一个或多个功能区构成，这些区通常被称为域。结构域的不同组合方式产生了多种多样的蛋白质，因此蛋白结构域的鉴别对分析蛋白质的功能来说尤其重要。

Pfam数据库是一个蛋白质家族大集合，依赖于多序列比对和隐马尔可夫模型。使用Pfam数据库和hmmscan(版本：3.3.2；参数：e-value 0.01)软件进行所述新转录本和新基因的结构域预测，获得蛋白质的保守序列、motif和结构域等。

KOfam数据库是以隐马尔可夫模型构建的具有预先计算的自适应得分阈值的KEGGOrthologs数据库，使用KOfam数据库和kofam_scan(版本：1.3.0；参数：-E 1e-5)软件进行所述新转录本和新基因的同源性搜索，可以准确地将KEGG Orthologs(KOs)分配给蛋白质序列。

十、新基因和新转录本结构分析

1、转录本结构优化

由于使用的软件或数据本身的局限性，导致所选参考基因组的注释往往不够精确，因此有必要对原有注释的转录本结构进行优化。

使用gffcompare(版本：0.12.1；参数：-R-C-K-M)软件，将新转录本与基因组的已知转录本进行比较，获得新转录本和新基因，对现有注释进行补充，得到一致性序列在调整前与调整后的区间；如果存在转录本在原有转录本边界之外的区域，将转录本的非翻译区(Untranslated Region，UTR)向上下游延伸，从而修正转录本的边界。通过比较结果，对基因的5’端或3’端进行延长。

例如，基因GENE1，既新鉴定了个转录本Gene1.t1(属于上述第九点)，也是对老转录本Trans1进行了结构的调整，区间从参考的(100-1000)变成了(100-1500)。

结构优化展示如下表10所示。

表10结构优化展示结果

2、可变剪切分析

基因转录生成的前体mRNA(pre-mRNA)，有多种剪接方式，选择不同的外显子，产生不同的成熟mRNA，从而翻译为不同的蛋白质，构成生物性状的多样性。这种转录后的mRNA加工过程称为可变剪接或选择性剪接(Alternative splicing)。

使用LASER(https://github.com/hilgers-lab/LASER；版本：v0.1.0；参数：默认)软件获取每个样品存在的可变剪接类型。主要的转录本可变剪接类型如下图4所示。

对转录本发生可变剪接事件情况进行统计，预测的可变剪接事件数量统计见下表11所示。

表11预测的可变剪接事件数量统计结果

3、融合转录本分析

融合转录本是指两个或多个转录本的编码区首尾相连，置于同一套调控序列(包括启动子、增强子、核糖体结合序列、终止子等)控制之下构成的嵌合体转录本。

目前研究表明，许多疾病、癌症与融合转录本的发生相关。使用CTAT-LR-Fusion(https://github.com/NCIP/Trinity_CTAT，版本：1.0.0；参数：default)软件比对及寻找融合转录本。

检测融合转录本的分析原理为：使用ctat-minimap2进行嵌合只读搜索，基于读取和基因组比对位置来筛选嵌合读取比对，以定义候选者。对于每个候选，构建有序和定向融合对的模型。将候选嵌合读数重新排列到这些融合重叠群的数据库中，并根据读数排列断点(也称为融合转录物断点)对每个融合对的读数支持进行评分。

4、lncRNA预测

lncRNA(长链非编码RNA)是一类转录本长度超过200nt，不编码蛋白质的RNA分子。使用CNCI(版本：2.0；参数：default)、CPC2(版本：standalone_python3v1.0.1)软件和PLEK软件对新鉴定的转录本进行编码潜能预测。

CPC2编码潜能计算软件将转录本与已知蛋白数据库做BLAST比对，CPC2依据转录本各个编码框的生物学序列特征，通过支持向量机的分类器来评估转录本的编码潜能。

CNCI根据相邻三核苷酸的频谱有效地区分编码和非编码的序列，CNCI可以有效地对不完整的转录本和反义转录本对进行编码潜能预测。

PLEK软件通过序列的kmer构成来区分编码和非编码转录本。

十一、转录本表达定量

为了消除假阳性及批次效应，使得不同样品的不同基因/转录本表达水平具有可比性，采用TPM(Transcripts Per Kilobase Million)作为衡量表达水平的指标。

在转录组测序中，常见的几种reads count值的标准化方法基本都是基于测序深度和转录本长度进行标准化。在TPM的计算过程中，对于每个转录本，将其reads count值除以其转录本的长度(以千碱基为单位)，获得转录本的每千碱基reads覆盖数(reads permillion，RPK)。计算样本中所有RPK值，然后除以1000000，获得百万缩放因子。将RPK值除以百万缩放因子，获得TPM值。

使用TPM时，每个样本中所有TPM的总和是相同的，这样可以更轻松地比较每个样本中映射到转录本的reads比例。

使用bamboo(版本：3.2.5)软件对转录本进行表达定量。

十二、差异表达转录本分析

利用表达定量得到的转录本在各个样品中的表达量reads count数据，进行差异表达分析。差异表达分析的软件为DESeq2(版本：1.26.0)，筛选阈值为padj<0.05，且|log2FoldChange|>1；如若显著差异转录本数目＜10时，筛选阈值为pvalue<0.05，且|log2FoldChange|>1。鉴定的差异表达转录本数目如图5所示；图中，横坐标为各比较组，纵坐标为上下调数目。

富集分析，是指利用统计学方法，研究差异表达的转录本功能是否集中在某些特有的功能分类上。常用的富集分析包括Gene Ontology富集分析和KEGG Pathway富集分析。富集分析统计学方法一般使用超几何检验，找出差异转录本相对于所有有注释的转录本显著富集的GO条目或KEGG通路。

差异表达转录本富集分析软件为clusterProfiler(版本：3.14.3)，qvalue是做过多重假设检验校正之后的pvalue，qvalue的取值范围为[0,1]，越接近于零，表示富集越显著。

十三、poly(A)分析

大多数真核mRNA的3’非翻译区域下游含有一个非模板化的poly(A)尾，这个尾有助于mRNA的稳定和向细胞质内的运输。poly(A)分析一般分为APA分析和poly(A)尾长分析。APA是指一个基因有多个poly(A)尾的插入位点。poly(A)尾长度和插入位点会影响mRNA的翻译效率。因此，poly(A)分析可解析mRNA降解和翻译的动态调控过程，也能够发现一些RNA裂解和加尾的未知特点。

使用NanoPolish(版本：0.13.2；参数：poly(A))软件对clean_reads的poly(A)进行计算。

DRS测序得到的每条read对应一个转录本，且转录本可以被测到完整的3’端，因此比对参考基因组的终点就是poly(A)插入位点。

使用minimap2软件将reads比对到参考基因组序列后，从bam文件中提与比对到的终点位置，再使用Quantifypoly(A)软件进行poly(A)位点的鉴定、聚类与注释。

本实施方式的样品中poly(A)长度统计结果如表12所示。

表12样品中poly(A)长度统计结果

十四、甲基化分析(甲基化鉴定、分布、差异与关联)

甲基化在核酸和蛋白质中是一种非常重要的修饰，调节基因的表达和关闭，与癌症、衰老、老年痴呆等许多疾病密切相关，是表观遗传学的重要研究内容之一。甲基化是指从活性甲基化合物(如S-腺苷基甲硫氨酸)上将甲基催化转移到其他化合物的过程。可形成各种甲基化合物，或是对某些蛋白质或核酸等进行化学修饰形成甲基化产物。在生物系统内，甲基化是经酶催化的，这种甲基化涉及基因修饰、基因表达的调控、蛋白质功能的调节以及核糖核酸(RNA)加工。Direct RNA测序能在读取RNA碱基序列的同时获得碱基修饰信息。

RNA甲基化是在甲基转移酶的催化作用下，在RNA分子的某一个原子上添加一个甲基基团(CH3)，主要包括m5C RNA甲基化、m6A RNA甲基化、m7G RNA甲基化等。目前最热门的是m6A RNA甲基化，即RNA分子腺嘌呤第6位氮原子上的甲基化修饰(N6-methyladenosine，m6A)，是真核生物mRNA最常见的一种转录后修饰。碱基修饰如m6A，可调节RNA分子的活性和稳定性，并与多种人类疾病和抗菌素耐药性相关。与传统技术不同，纳米孔技术可以对天然RNA分子进行测序，无需扩增或逆转录，可直接鉴定核苷酸序列上的碱基修饰。

Tombo(版本：1.5)(Stoiber et al.,2016)是nanopore官网推出用于从原始纳米孔测序数据中鉴定核苷酸修饰的一套工具。使用Tombo软件预测出RNA序列中m5C修饰位点。使用m6Anet(Christopher Hendra,et al，参数：默认)软件预测出RNA序列中m6A修饰位点。使用R语言、ggplot2、ggseqlogo用来计算并绘制(统计和可视化)甲基化位点的位置、分布、5bp的motif图。

使用methylkit软件进行差异甲基化位点(Differential methylation loci，DML)分析。如果每组有多个样本，采用Logistic回归检验。

对每个转录本的甲基化进行差异检验。当每组仅有一个样本时，使用Fisher软件精确测试进行差异甲基化的检测。默认情况下，差异水平设置为10，q值设定为0.01，使用SLIM(sliding linear model)方法矫正P值。

十五、假尿苷分析

假尿苷修饰(Ψ)是目前已知丰度最高的RNA修饰，广泛存在于多个物种的多类RNA中。作为尿苷的5位异构体，假尿苷在真核生物中的形成机制主要有两种：依赖于假尿苷合酶，或是依赖于H/ACA核糖核蛋白复合体。假尿苷修饰在多个生物学过程中发挥作用，同时不同位点的假尿苷修饰具有不同的功能。而人为地在mRNA终止密码子中引入假尿苷修饰则可以使其具有编码能力，这改变了传统的中心法则。纳米孔RNA测序是一种有前景的方法，可用于区分和鉴定天然RNA中不同的RNA修饰。现有研究发现假尿苷会导致纳米孔中特征性的碱基识别“错误”，因此可以通过捕捉错误的碱基识别用来鉴定假尿苷修饰。

位点的鉴定使用了Nanopsu软件(默认参数)，软件使用十二种错配相关的特征值作为输入，建立了AUC高达0.94的随机森林模型。结果会输出假尿苷修饰的P值。对P值进行筛选(>0.9)，即可得到结果，如图6所示。

上述具体实施方式中，实验端的样本富集和反转录引物连接环节添加了RNA酶抑制剂，有效的降低了RNA在建库过程中发生降解的情况；并将连接反应中的quick T4 DNALigase更换为T4 DNA ligase，一定程度上提高了单链RNA和双链接头的连接效率，最终使得测序所得数据更接近于样本转录本的真实表达情况

转录组测序和结构分析方面，在原有三代转录组分析基础上，基于nextflow与docker搭建了一套一体化的分析流程，使得DRS的分析具有可移植、可重复、便携的特点。与此同时更新了从新转录本识别到定量差异分析的分析内容，使其充分利用了DRS的特点，更好的发挥RNA直接测序的特长。新增了poly(A)、甲基化、假尿苷的鉴定、分布与差异分析，并且将其与表达进行关联，系统化、多角度地对转录本的表达进行阐释。

以上具体实施方式详细描述了本发明的实施，但是，本发明并不限于上述实施方式中的具体细节。在本发明的权利要求书和技术构思范围内，可以对本发明的技术方案进行多种简单改型和改变，这些简单变型均属于本发明的保护范围。

Claims (9)

1.基于DRS的生物遗传样本转录组测序分析方法，其特征在于，步骤包括：

富集样本mRNA，建立直接测序文库；

基于Illumina进行二代测序，基于ONT进行三代DRS测序，分别得到二代原始测序数据和三代原始测序数据；

对所述二代原始测序数据和三代原始测序数据分别进行质控，基于所述二代原始测序数据对所述三代原始测序数据进行纠错校正；将纠错后的转录本序列与参考基因组进行比对，统计得到一致性序列；

过滤所述旧转录本中的全长reads，然后将一致性序列与参考基因组的已知转录本进行比对，得到新转录本和新基因；

识别所述新转录本的潜在编码区序列，预测所述新转录本和新基因的编码区序列，对所述新转录本和新基因进行功能注释；

对转录本的结构进行优化，进行可变剪切分析、融合转录本分析、lncRNA预测、转录本表达定量分析、差异表达和功能富集分析、poly(A)分析、甲基化分析和假尿苷分析。

2.根据权利要求1所述的基于DRS的生物遗传样本转录组测序分析方法，其特征在于，对所述二代原始测序数据的质控方式为：

去除N碱基含量＞5％的reads；

并且，去除质量值≤5，碱基数目达到50％的reads；

并且，去除有adapter污染的reads；

并且，去除PCR扩增造成的重复序列；

对所述三代原始测序数据的质控方式为：过滤Q＜9且长度＜100bp的reads。

3.根据权利要求1所述的基于DRS的生物遗传样本转录组测序分析方法，其特征在于，基于所述二代原始测序数据，使用LRECE软件中的run_correction_tools.sh脚本对所述三代原始测序数据进行纠错校正；

使用minimap2软件将纠错校正后的转录本序列与参考转录本进行比对，使用samtools软件统计比对结果，得到一致性序列。

4.根据权利要求1所述的基于DRS的生物遗传样本转录组测序分析方法，其特征在于，使用LAFITE软件在所述一致性序列中引入poly(A)鉴定标签过滤去除全长reads；使用gffcompare软件将所述一致性序列与基因组的已知转录本进行比对，发现新转录本和新基因。

5.根据权利要求1所述的基于DRS的生物遗传样本转录组测序分析方法，其特征在于，使用TransDecoder软件对所述新转录本和新基因进行潜在编码区序列识别，使用orfpy软件对所述新转录本和新基因的编码区进行预测。

6.根据权利要求1所述的基于DRS的生物遗传样本转录组测序分析方法，其特征在于，对所述新转录本和新基因进行Nr、Pfam、Uniprot、KEGG、GO、KOG/COG和PATHWAY七个数据库的转录本功能注释。

7.根据权利要求6所述的基于DRS的生物遗传样本转录组测序分析方法，其特征在于，功能注释的方法包括：

使用diamond blastp软件将新转录本中转录本编码的蛋白序列与现有蛋白质数据库Uniprot和Nr进行比对，获得序列的功能信息，以及蛋白可能参与的代谢通路信息；基于数据库之间的关联，得到KOG/COG注释结果，进行KOG/COG的分类统计及绘图；

使用Pfam数据库和hmmscan软件进行所述新转录本的结构域预测；使用KOfam数据库和kofam_scan软件进行所述新转录本的同源性搜索。

8.根据权利要求1所述的基于DRS的生物遗传样本转录组测序分析方法，其特征在于，使用gffcompare软件，将一致性序列与基因组已知转录本进行比较；如果存在转录本在原有转录本边界之外的区域，则将转录本的非翻译区向上下游延伸，修正转录本的边界；通过比较结果，对基因的5’端或3’端进行延长，完成对转录本的结构优化。

9.根据权利要求1所述的基于DRS的生物遗传样本转录组测序分析方法，其特征在于，使用LASER软件完成转录本的可变剪切分析；

使用CTAT-LR-Fusion软件比对、寻找融合转录本；

使用CNCI软件、CPC2软件和PLEK软件对所述新转录本进行lncRNA预测；

采用TPM作为衡量表达水平的指标，使用bamboo软件进行所述新转录本表达定量分析；

使用DESeq2软件进行差异表达分析，筛选阈值为padj<0.05，且|log2FoldChange|>1；若显著差异转录本数目＜10时，筛选阈值为pvalue<0.05，且|log2FoldChange|>1；基于Gene Ontology和KEGG Pathway方法，使用clusterProfiler软件进行功能富集分析；

使用NanoPolish软件对有效reads的poly(A)进行计算，使用minimap2软件将reads比对到参考基因组序列后，从bam文件中提取比对到参考基因组的终点位置，再使用Quantifypoly(A)软件进行poly(A)位点的鉴定、聚类与注释；

使用Tombo软件预测出RNA分子序列中m5C修饰位点，使用m6Anet软件预测出RNA分子序列中m6A修饰位点，使用R语言计算并绘制甲基化位点的位置、分布、5bp的motif图，使用methylkit软件进行差异甲基化位点分析；若每组有若干样本则采用Logistic回归检验；

使用Nanopsu软件获得RNA分子序列中的假尿苷修饰的P值，基于P>0.9进行P值筛选。