一种多重混样直接RNA纳米孔测序方法

一种多重混样直接rna纳米孔测序方法
优先权申请
1.本技术要求2023年3月17日提交的中国发明专利申请cn2023102647091“一种新型直接rna纳米孔测序方法及系统”的优先权，该优先权发明专利申请以引用方式全文并入。
技术领域
2.本发明涉及测序领域，具体涉及一种多重混样直接rna纳米孔测序方法。


背景技术：

3.三代测序技术的发展彻底改变了基因组和转录组的长读测序的能力。它能够在无需碎片化和扩增的情况下，直接对dna和rna及其修饰进行测序。牛津纳米孔技术(ont)的直接rna测序(drs)技术可以通过一系列蛋白质纳米孔在合成膜上直接对rna分子进行测序。
4.当马达酶缓慢地将单个rna分子通过纳米孔时，不同的核苷酸阻塞了纳米孔时会导致电流波动，这可以通过碱基调用算法记录下来以确定rna分子的序列。直接rna测序技术可以产生长读序列(长达21kb)，旨在获得完整的转录本或接近完整的rna基因组。完整的测序使其在测量rnapoly(a)(聚(a))尾长度和分析聚(a)尾对基因表达调控的调节功能方面具有强大的作用。此外，直接rna测序可以直接检测rna分子的修饰，例如m6a、7-甲基鸟苷(m7g)、伪尿嘧啶(ψ)和2'-o-甲基化(nm)，因为这些修饰会产生不同的电流波动。
5.此外，无扩增文库制备方法可以避免扩增或逆转录(rt)偏差，特别适用于识别出外显子内含物(exitrons)，这些外显子内含物很可能源自逆转录造成的伪迹。这些优势使直接rna测序技术在特定亚型的rna结构确定、trna和rrna鉴定以及rna病毒、原核生物和寄生虫等生物体基因组测序方面具有特殊的优势。近期的研究还利用直接rna测序技术在广泛的物种中进行转录组或表观转录组分析，包括人类、动物(线虫、昆虫)、植物、浮游动物和病原体等。直接rna测序技术的另一个应用是通过使用碱基类似物标记新生rna，例如5-乙炔基尿嘧啶和4-硫尿嘧啶，随后进行ont直接rna测序来分析rna代谢动力学。
6.尽管drs具有独特的优点，但它有两个主要的限制，包括缺乏多样本混样测序方法和需要大量rna，通常每次测序需要500ng rna。众所周知，ont提供用于多路复用cdna文库的条形码协议，这些协议依赖于直接连接dna适配器到cdna序列上。在这种情况下，条形码和cdna序列都可以在dna模型下轻松地被识别。然而，在drs的情况下，由于适配器是dna寡核苷酸，而目标片段是rna分子，因此dna适配器不能在rna模型或dna模型下被正确识别。因为dna的过孔速度(450bp/sec)比rna(70nt/sec)快6倍。为了克服这些限制，已经有团队开发了deeplexicon这一方法，将有限的rna样品混入到一张芯片中，生成了多路复用(multiplex)直接rna测序数据集，以此降低单个样品的测序成本。但由于条形码设计和算法的限制，该方法只能在一张芯片中多路复用4个样品。尽管已经有文献报道deeplexicon可以训练新模型，但如果使用不同的序列和多达4个以上的条形码，它将无法准确地解复用读取。因此，提供更多条形码以增加多路复用能力并实现高精度仍然是一个挑战。


技术实现要素：

7.第一方面，本发明提供了一种直接rna纳米孔测序方法(简称为decoder)，其特征在于，所述方法用于对n种测试样本rna多核苷酸进行测序，其中2≤n≤24且n为整数，所述方法包括：
8.s1针对每一种测试样本rna多核苷酸，提供一种双链dna多核苷酸来捕获一种待测rna序列并形成第一杂交结构；所述双链dna多核苷酸为逆转录接头的前向链和后向链组装成的双链结构，所述前向链从5’至3’方向包括条形码序列和第一连接序列，所述后向链从5’至3’方向包括第二连接序列、条形码序列和poly(t)序列，所述待测rna序列包括所述测试样本rna多核苷酸和poly(a)序列，所述双链dna多核苷酸通过所述poly(t)序列捕获所述待测rna序列；所述条形码序列对所述每一种测试样本rna多核苷酸具有特异性，所述条形码序列选自seq id nos:1-24中的任意一种；
9.s2将所述第一杂交结构逆转录，形成第二杂交结构；
10.s3将测序接头与所述第二杂交结构连接，形成第三杂交结构；
11.s4对所述第三杂交结构进行纳米孔测序，获得原始测序数据；
12.s5提取所述原始测序数据对应的电流信号数据；
13.s6将所述电流信号数据中所述条形码序列的电流信号数据进行分割，并调用基于所述条形码序列训练得到的分类模型进行分析，获得所述条形码序列的序列信息，进而获得所述n种测试样本rna多核苷酸的分类信息。
14.在一些实施例中，所述双链dna多核苷酸经预退火生成。
15.在一些实施例中，所述预退火的方法包括将每个所述逆转录接头的前向链和后向链混合在缓冲液中，按照以下程序进行：95℃变性5分钟，65℃、50℃、37℃、22℃分别退火30分钟，4℃保持。
16.在一些实施例中，所述第一连接序列包括tagtaggttc。
17.在一些实施例中，所述第二连接序列包括gaggcgagcggtcaatttt。
18.在一些实施例中，所述poly(t)序列包括10-30个重复t序列。
19.在一些实施例中，所述poly(a)序列包括10-30个重复a序列。
20.在一些实施例中，所述poly(t)序列包括10个重复t序列。
21.在一些实施例中，所述poly(a)序列包括10个重复a序列。
22.在一些实施例中，所述分类模型的训练方法包括将训练样本数据中条形码序列的电流信号数据进行分割，利用机器学习算法对分割后的条形码序列的电流信号数据进行训练学习，从而确定所述分类模型，其中所述训练样本数据包括一种或多种训练样本rna多核苷酸和对应的条形码序列的电流信号数据，所述条形码序列选自seq id nos:1-24中的任意一种。
23.在一些实施例中，所述分类模型为基于选自seq id nos:1-24中的任意2种、3种、4种、5种、6种、7种、8种、9种、10种、11种、12种、13种、14种、15种、16种、17种、18种、19种、20种、21种、22种、23种或24种的所述条形码序列训练得到的分类模型。
24.在一些实施例中，所述训练样本rna多核苷酸包括体外转录rna。
25.在一些实施例中，所述分割为将所述条形码序列的电流信号数据分割成70-100列矩阵。
26.在一些实施例中，所述分割为将所述条形码序列的电流信号数据分割成100列矩阵。
27.在一些实施例中，所述机器学习算法包括k最近邻、神经网络、朴素贝叶斯分类、回归树、adaboost和随机森林中的一种或多种。
28.在一些实施例中，所述机器学习算法包括随机森林。
29.在一些实施例中，所述测试样本rna多核苷酸的来源包括病原体、植物、浮游动物、昆虫、哺乳动物和肿瘤中的一种或多种。
30.在一些实施例中，所述病原体包括rna病毒、细菌、真菌和寄生虫中的一种或多种。
31.在一些实施例中，所述体外转录rna包括表2所示的基因中的一种或多种。
32.在一些实施例中，所述分类的分类可能性的阈值》0.3。
33.在一些实施例中，所述s5进一步包括提取所述原始测序数据对应的原始电流信号数据，并对所述原始电流信号数据进行预处理，获得预处理后的电流信号数据，所述预处理包括降噪处理、抛光处理和归一化处理中的一种或多种。
34.在一些实施例中，所述方法进一步包括将所述一种或多种训练样本rna多核苷酸的电流信号数据与对应的参考序列比对，以比对质量大于60、比对至唯一参考序列为条件将不满足条件的读取序列过滤。
35.第二方面，本发明提供了一种用于直接rna测序的试剂盒，其特征在于，所述试剂盒包括如seq id nos:1-24所示的逆转录接头的条形码序列。
36.第三方面，本发明提供了一种直接rna纳米孔测序系统，其特征在于，包括：
37.测序数据接收模块，用于接收和存储测序数据，所述测序数据包括双链多核苷酸的测序数据，所述双链多核苷酸包括第一链和第二链，所述第一链从5’至3’方向包括rna多核苷酸、poly(a)序列、条形码序列、第一连接序列和测序接头，所述第二链从5’至3’方向包括测序接头、第二连接序列、条形码序列、poly(t)序列和所述rna多核苷酸的逆转录序列；所述测序数据包括测试样本数据和训练样本数据；所述测试样本数据包括n种测试样本rna多核苷酸对应的测序数据，其中2≤n≤24且n为整数，所述条形码序列对每一种测试样本rna多核苷酸具有特异性，所述条形码序列选自seq id nos:1-24中的任意一种；
38.信号处理模块，用于将所述测序数据处理成电流信号数据；
39.分类模块，用于调用基于所述测试样本数据中的所述条形码序列构建的分类模型，对所述测试样本数据的电流信号数据进行分类并输出所述n种测试样本rna多核苷酸的分类结果。
40.在一些实施例中，所述系统进一步包括模型训练模块，所述模型训练模块利用机器学习算法对所述训练样本数据的电流信号数据进行训练学习，从而确定所述分类模型。
41.在一些实施例中，所述信号处理模块对所述测序数据进行预处理和分割。
42.在一些实施例中，所述预处理包括降噪处理、抛光处理和归一化处理中的一种或多种。
43.在一些实施例中，所述分割为将所述条形码序列的电流信号数据分割成70-100列矩阵。
44.在一些实施例中，所述分割为将所述条形码序列的电流信号数据分割成100列矩阵。
45.在一些实施例中，所述机器学习算法包括k最近邻、神经网络、朴素贝叶斯分类和回归树、adaboost和随机森林中的一种或多种。
46.在一些实施例中，所述机器学习算法包括随机森林。
47.在一些实施例中，所述系统进一步包括过滤模块，用于基于设置的分类可能性的阈值，将低于所述阈值的读取过滤掉。
48.在一些实施例中，所述训练样本数据包括体外转录rna的测序数据。
49.在一些实施例中，所述测试样本rna多核苷酸的来源包括病原体、植物、浮游动物、昆虫、哺乳动物和肿瘤中的一种或多种。
50.在一些实施例中，所述病原体包括rna病毒、细菌、真菌和寄生虫中的一种或多种。
51.在一些实施例中，所述阈值》0.3。
52.在一些实施例中，所述体外转录rna包括表2所示的基因中的一种或多种。
53.在一些实施例中，所述训练样本数据中用作训练集和测试集的比例为7:3。
54.在一些实施例中，所述模型训练模块用于将所述训练样本数据的电流信号数据与对应的参考序列比对，以比对质量大于60、比对至唯一参考序列为条件将不满足条件的读取序列过滤。
55.有益效果
56.现有技术中纳米孔直接rna混样测序分类器deeplexicon能够将有限的rna样本混入到一张芯片中，生成多路复用(multiplex)直接rna测序数据集，以降低单个样本的测序成本。然而，deeplexicon只能在同一芯片上拆分4个样本，具有混样样本通量低、拆分准确率低、实现准确拆分需要舍弃大量的测序数据等诸多限制，造成测序成本大大提高。
57.本发明提供的方法(简称decoder)提供了24种条形码以在一张芯片中标记多种rna样本，并通过对直接rna-seq数据解复用(demultiplex)和机器学习分类，实现在一张芯片中高精度地多路复用更多样本且检测限低，克服了rna样本量低和测序成本高的限制。本发明采用的将条形码电流分割成70-100列矩阵，能够降低电流之间的干扰，降低了噪音。经验证，decoder的恢复率和准确率可以达到100％和92.2％。将decoder应用于多种病原体(rna病毒、细菌、真菌和疟原虫)和胶质瘤样本，均得到了良好的验证结果。此外，本发明通过decoder还揭示了不同胶质瘤样本的转录组和突变的一些有意义的发现，这表明本发明提供的方法具有较大的实际应用潜力。本发明提供的方法在纳米孔直接rna测序中可以实现在一张芯片中同时对多个生物样本测序，降低了生物样本起始投入量，降低了多样本测序成本，能够实现微量生物样本的测序，已在肿瘤患者样本中得到应用，并具有在其他临床样本中应用的潜力。
58.综上所述，本发明提供了一种完整、成熟、高通量、低成本的直接rna混样测序方法，提供更多条形码以增加多路复用能力并实现对纳米孔drs中的多种生物样本的高精度的测序和分析，有望在临床样本中实际应用。
59.如本文所使用，“核酸”、“核酸分子”、“多核苷酸”和“寡核苷酸”可互换地使用，是指共价连接的核苷酸序列(即用于rna的核糖核苷酸和用于dna的脱氧核糖核苷酸)，其中一个核苷酸的戊糖的3’位置通过磷酸二酯基团连接到下一个核苷酸的戊糖的5’位置。“多核苷酸”包括单链和双链多核苷酸，例如各种dna、rna分子、或其任何杂交体或其衍生物或组合或片段。
60.如本文所使用，“测试样本rna多核苷酸”是指期望被测序的rna。所述rna可以是编码和/或非编码rna，例如mrna、sirna、rrna、mirna、trna、lncrna、snorna、snrna、exrna、pirna。
61.如本文所使用，“互补”是指两条核酸序列能够根据碱基配对原则(waston-crick原则)在彼此之间形成氢键，并由此形成双链结构。在本发明中，“互补”包括实质上互补和完全互补；其中完全互补是指一条核酸序列中的每一个碱基都能够与另一条核酸序列中的碱基配对，而不存在错配或缺口；实质上互补是指一条核酸序列中的大部分碱基都能够与另一条核酸序列中的碱基配对，其允许错配或缺口存在。通常，在允许核酸杂交、退火或扩增的条件下，互补(例如实质上互补或完全互补)的两条核酸序列将选择性地/特异性地发生杂交或退火，并形成双链结构。相应地，“不互补”(或“不匹配”)是指两条核酸序列在允许核酸杂交、退火或扩增的条件下不能发生杂交或退火。
62.如本文所使用，“杂交”和“退火”可互换地使用，是指互补的单链核酸分子形成双链核酸的过程。通常，完全互补或实质上互补的两条核酸序列可发生杂交或退火。
63.如本文所使用，“后端”(或“下游”)用于描述两条核酸序列的相对位置关系，并且具有本领域技术人员通常理解的含义。例如，“一条核酸序列在另一条核酸序列的后端”是指，当以5’端至3’端方向排列时，与后者相比，前者位于更靠后的位置(即更接近3’端的位置)。相应地，“前端”(或“上游”)具有与“后端”相反的含义。
64.如本文所使用，“条形码”是一种特定的核苷酸序列，用于标记期望被测序的rna。在本发明中，条形码的序列是经设计的。换句话说，条形码的序列是已知的。在一些实施例中，所述条形码的序列包括seq id nos:1-24中的一种或多种。在一次测序过程中，每一种期望被测序的rna对应的条形码序列是唯一且特定的(也可以理解为具有特异性)。此外，条形码序列在本发明的测序结构的位置是固定的(例如在poly(a)的下游)，因此可以通过poly(a)的连续a信号来定位读取序列中的条形码序列。本发明可以通过机器学习算法对基于不同条形码组合(例如seq id nos:1-24中的中的任意1种、2种、3种、4种、5种、6种、7种、8种、9种、10种、11种、12种、13种、14种、15种、16种、17种、18种、19种、20种、21种、22种、23种或24种)的训练样本数据进行训练学习，以构建分类模型。进而，在一次测序过程中(例如一张芯片中)，基于期望被测序的rna相连的条形码序列，可以选择相应条形码的组合训练得到的分类模型来分析测序数据，进而实现对期望被测序的rna的分类。在一些实施例中，在一次测序过程中，本发明提供的方法可以实现对1种、2种、3种、4种、5种、6种、7种、8种、9种、10种、11种、12种、13种、14种、15种、16种、17种、18种、19种、20种、21种、22种、23种或24种期望被测序的rna的测序。
65.如本文所使用，“读取”是指通过测序方法获得的核酸分子的任何部分或全部的测序信息。读取可以存储在存储介质中并被适当地处理以确定其是否与参考序列匹配或符合其他标准。可以从测序装置直接获得读取，或者从存储的关于样本的序列信息间接获得。
66.如本文所使用，“参考序列”是指读取与其比对的任何已知序列。在一些实施例中，参考序列可以对应于生物体的基因组或转录组的全部或仅一部分。
67.如本文所使用，“比对”是指将读取与参考序列进行比较的过程，由此确定参考序列是否含有读取序列。如果参考序列含有读取，则读取可以被定位到参考序列，或者被定位到参考序列中的特定位置。比对中读取的匹配可以是100％序列匹配或小于100％(非完全
匹配)。
68.如本文所使用，“混样”是指含有不同基因组的核酸分子混合物的样本。
69.如本文所使用，“库”或“文库”是指源自一种或多种核酸样本的核酸分子的集合。在一些实施例中，所述核酸分子包括可鉴别的序列标记(例如条形码)。在一些实施例中，可以将两个或多个文库合并以创建文库池。
70.如本文所使用，“复用”或“多路复用”是指一种或多种核酸样本或其所衍生的库分子在一个池、试管或反应中(例如在一张芯片中)的集合。
71.如本文所使用，“样本”是指包含生物分子的样本，典型地源自生物学流体、细胞、组织、器官或生物体。“生物样本”是指包含生物分子的任何生物来源样本。在本发明中，所述生物分子包括核酸或核酸混合物，优选为rna。在一些实施例中，“样本”或“生物样本”可以是液体样本(例如体液，包括血液、血浆、血清、尿液、阴道液、水囊肿液、阴道冲洗液、胸膜液、腹水液、脑脊髓液、唾液、汗液、泪液、痰液、支气管肺泡灌洗液、乳头排出液、来自身体不同部位(例如，甲状腺、乳房)的抽吸液)或固体样本(例如组织样本、活检样本或从中衍生的组织培养物或细胞及其子代)。在一些实施例中，“样本”或“生物样本”可以包括富集特定类型分子的样本，例如核酸。在一些实施例中，“样本”或“生物样本”还可以包括临床样本，例如通过手术切除获得的组织、通过活检获得的组织、培养中的细胞、细胞上清液、细胞裂解物、组织样本、器官、骨髓、血液、血浆、血清等。在一些实施例中，“样本”或“生物样本”还可以包括肿瘤样本，例如来自癌症患者的细胞(例如癌细胞)的样本或从癌症患者的癌细胞(例如，细胞裂解物、细胞提取物)中获得的核酸样本。在一些实施例中，“样本”或“生物样本”中的核酸可以是游离的。在一些实施例中，“样本”或“生物样本”的来源可以是病原体，例如rna病毒、细菌、真菌、寄生虫。在一些实施例中，“样本”或“生物样本”的来源可以是植物、浮游动物、昆虫、哺乳动物(例如人类、狗、猫、马、山羊、绵羊、牛、猪、大鼠、小鼠等)。在一些实施例中，“样本”或“生物样本”的来源还可以是肿瘤。“样本”或“生物样本”可以在从来源获得后直接使用，或者在对其进行预处理(例如，固定、嵌入介质、切片、清洗、富集)后使用。在一些实施例中，还需提取“样本”或“生物样本”中的核酸(例如rna多核苷酸)。
72.如本文所使用，“测试样本”是指包含期望被测序的核酸的样本或生物样本。在一些实施例中，所述核酸优选为rna多核苷酸。在一些实施例中，所述测试样本来源于生物。在一些实施例中，所述测试样本还可以来源于环境，例如土壤、湖、河等。
附图说明
73.为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍。在所有附图中，类似的元件或部分一般由类似的附图标记标识。附图中，各元件或部分并不一定按照实际的比例绘制。显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其它的附图。
74.图1示出了decoder的示意图、示例性流程和算法比较；
75.图2示出了本发明的decoder与现有技术的分类器的比较；
76.图3示出了不同条形码的比较；
77.图4示出了本发明的decoder在真实样本中的性能效果；
78.图5示出了本发明的decoder检测来自病毒样本rna；
79.图6示出了本发明的decoder对肿瘤样本的突变和m6a分析；
80.图7示出了本发明提供的一种直接rna纳米孔测序方法的流程图；
81.图8示出了本发明提供的一种多条形码直接rna纳米孔测序分类系统的模块示意图。
具体实施方式
82.为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
83.本文中，使用用于表示元件的诸如“模块”、“部件”或“单元”的后缀仅为了有利于本发明的说明，其本身没有特定的意义。因此，“模块”、“部件”或“单元”可以混合地使用。
84.本文中，术语“上”、“下”、“内”、“外”“前”、“后”、“一端”、“另一端”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性。
85.本文中“和/或”包括任何和所有一个或多个列出的相关项的组合。
86.本文中“多个”意指两个或两个以上，即其包含两个、三个、四个、五个等。
87.需要说明的是，在本文中，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者装置不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者装置所固有的要素。在没有更多限制的情况下，由语句“包括一个
……”
限定的要素，并不排除在包括该要素的过程、方法、物品或者装置中还存在另外的相同要素。
88.如在本说明书中使用的，术语“大约”，典型地表示为所述值的+/-5％，更典型的是所述值的+/-4％，更典型的是所述值的+/-3％，更典型的是所述值的+/-2％，甚至更典型的是所述值的+/-1％，甚至更典型的是所述值的+/-0.5％。
89.在本说明书中，某些实施方式可能以一种处于某个范围的格式公开。应该理解，这种“处于某个范围”的描述仅仅是为了方便和简洁，且不应该被解释为对所公开范围的僵化限制。因此，范围的描述应该被认为是已经具体地公开了所有可能的子范围以及在此范围内的独立数字值。例如，范围1～6的描述应该被看作已经具体地公开了子范围如从1到3，从1到4，从1到5，从2到4，从2到6，从3到6等，以及此范围内的单独数字，例如1，2，3，4，5和6。无论该范围的广度如何，均适用以上规则。
90.实施例一
91.1.1条形码设计及合成
92.如表1所示，本实施例纳入本发明设计的40个新的逆转录接头(reverse transcript adapter，rta)的条形码序列(rta-01至rta-40，长度在20-28bp之间)，以及来自ont官方的rta的条形码序列(rta-00)、deeplexicon的3个rta的条形码序列(rta-41、
rta-42和rta-43)和来自poreplex的4个rta的条形码序列(rta-44、rta-45、rta-46和rta-47)，共计48个rta的条形码序列。
93.本发明所基于的测序结构如图1a所示。每一对rta由前向链(forward strands)和后向链(reverse strands)组成。按照5’端至3’端的方向，前向链序列包括条形码序列(x1)以及接在条形码后端的第一连接序列(x2，也可以称为第一不匹配多核苷酸序列)，后向链包括第二连接序列(x2，也可以称为第二不匹配多核苷酸序列)、条形码序列(x1)以及接在条形码后端的poly(t)序列(例如10-30个重复t序列)。将rta的前向链和后向链序列进行杂交并组装成双链结构。根据碱基互补配对原则，rta包括条形码区域(x1)、不匹配多核苷酸区域(x2)以及突出末端的poly(t)序列。接着，将组装好的双链结构捕获并杂交带有poly(a)序列的待测rna序列，经逆转录形成rna/dna杂化双链结构。然后，再将rna/dna杂化双链结构与ont测序接头(x3)连接，最终形成本发明基于的完整的测序结构。在一些实施例中，所述第一不匹配多核苷酸序列包括tagtaggttc(seq id no:49)，第二不匹配多核苷酸序列包括gaggcgagcggtcaatttt(seq id no:50)。
94.接下来，本实施例使用oligoanalyzer
tm
工具评估了本发明所形成的dna寡核苷酸的gc含量、熔点和二级结构。本发明设计的条形码的长度被控制20-28nt，避免了超长的条形码可能会干扰磁珠纯化步骤中去除过量的rta的问题。
95.表1自定义rta序列
96.1.2体外转录(ivt)rna的制备
97.为了测试条形码并进行模型训练，本实施例生成了52个体外转录(in-vitro transcribed，ivt)rna，其序列来自27个小鼠基因、22个人类基因和3个sars-cov-2基因(表2)。针对这些基因序列，使用primer3web工具(version4.1.0)设计了前向引物和后向引物，并在前向引物中添加了t7 rna聚合酶(t7 rnapolymerase)启动子序列，在反向引物中添加了15个重复t序列(oligo(dt)15)序列。
98.使用逆转录酶，将人类和小鼠总rna以及sars-cov-2基因组rna逆转录为第一链cdna。使用上述设计的引物，对不同的dna片段进行pcr扩增，使用琼脂糖凝胶电泳确定pcr产物的质量。使用1μg纯化的pcr产物作为模板，通过t7 rna聚合酶生成ivtrna，并使用rna clean&concentrator-5按照制造商的建议纯化。使用qsep100确定ivtrna的长度，并使用nanodrop 2000测量每个ivtrna产物的浓度和纯度。
99.表2体外转录rna序列信息
100.1.3混样直接rna测序建库流程
101.自定义rta的组装。将每个rta的前向链(1.54μm)和后向链(1.4μm)分别混合在缓冲液中(10mm tris-hcl，ph 7.5，50mm kcl)(表3)。本实施例采用过量的前向链的原因在于，过量的前向链可以确保后向链的完全消耗，避免了后向链的存在可能会降低rna-rta的连接产物收率的问题。自定义rta的组装遵循如下pcr程序进行：95℃变性5分钟，65℃、50℃、37℃、22℃分别退火30分钟，4℃保持(表4)。最后，通过4％琼脂糖凝胶电泳确认rta是否完成组装。该组装产物能作为带有不同条形码的rta，用于本发明的混样直接rna测序方案。
102.表3自定义rta的组装
103.表4 pcr程序设计
104.混样直接rna测序文库制备。使用nebnext快速连接模块(neb，e6056)，将每个带有poly(a)序列的rna样本单独与预退火的自定义rta进行杂交和连接，并且记录其一一对应的关系。并使用superscriptiv逆转录酶(thermofisher，18090200)将其逆转录，形成rna/dna稳定双链并破坏rna的二级结构。使用1.8xrnaclean xpbeads(beckman，a63987)对产物进行纯化。然后，将每个反应的逆转录产物按相同的摩尔质量进行混合，并在室温下将rna适配体混合物(rmx)连接到rna/dna杂交体的3'端。连接后，将反应混合物使用1xrnaclean xpbeads进行纯化(beckman，a63987)，用缓冲液(wsb)洗涤两次，并用洗脱缓冲液(eb)洗脱产物。接下来，将遵从制造商的说明书(ont，sqk-rna002)构建上机混合液，根据minion或promethion测序仪适配的文库体积进行调整。接下来，将上机混合液加载至质检合格的ont测序芯片中，设置minknow软件中的测序参数，测序持续48h以上以获取足够的测序数据。
105.1.4decoder机器学习模型的训练
106.首先，使用了48种rta连接ivtrna的各种组合进行了5次独立的实验，获得了5次drs原始电信号数据集，在碱基调用(base-calling)和与参考序列比对后，生成了总共6,276,168个有效reads(范围从623k到2.30m)。质量控制后，由于生成的读取(reads)数量较少，24个rta条形码被过滤。
107.接下来，将ivtrna下机的原始信号数据(fast5)使用guppy(v3.1.5)在高精度模型下进行碱基调用。使用minimap2(v2.17-r941)将碱基调用reads(fastq)与ivtrna参考序列比对。接下来，以比对质量(mapq)大于60，比对至唯一参考序列为条件过滤比对reads，并根据它们比对的连接序列将其分组。使用r包rhdf5(https://github.com/grimbough/rhdf5)从fast5文件中提取原始信号。还可以使用r包smoother(https://cran.r-project.org/package＝smoother)使得原始电流信号更平滑，以减少原始数据的噪声。
108.鉴于poly(a)电流信号是连续的碱基a信号，本实施例首先抛光(polish)和归一化(normalize)电流信号，并从rna 3'端的信号的均值和方差中确定poly(a)的位置。于是，获得了ont测序接头(也可以理解为“测序接头”、“适配体”、“适配器”、“rna适配体混合物”等)、rta条形码、poly(a)和体外转录rna序列产生的一系列电流信号(图1b)。由于每个读取的测序接头序列是相同的，因此在分析时没有删除ont测序接头的信号。然后，使用changepoint包(https://cran.r-project.org/package＝changepoint)中的cpt.meanvar函数，根据均值和方差对poly(a)前面的rta条形码的电流信号进行分割，每个条形码信号被分成70-100段，优选为100段(图1c)。
109.以分割成100列矩阵为例，将所有条形码生成的100列矩阵最终用于模型训练。r包caret(https://cran.r-project.org/package＝caret)的实施是为了简化基于随机森林分类模型生成预测模型的过程。为避免过拟合，本实施例采用重复k折交叉验证(10折，重复
10次)的方法进行重采样。为了提高模型训练的计算效率，采用了r包中doparallel包支持的并行处理框架。由于每个条形码的序列读取数不同，因此对于每个分类器，本实施例首先使用caret包中的downsample函数随机抽取一个数据集，以便所有分类组与少数类具有相同的频率。然后，将这些reads的70％用于训练集，其余所有reads用作测试集。本实施例训练了2、4、6、8、10、12、14、16、18、20、24种条形码组合的模型。应当理解，还可以训练其他种类条形码组合(例如3种、5种、11种、19种、22种等)来产生更多的模型。
110.综上所述，如图7所示，本发明提供的多重混样直接rna纳米孔测序方法包括：首先，针对每一种测试样本rna多核苷酸，将rta的前向链和后向链序列进行杂交并组装成双链结构(即双链dna多核苷酸)。将组装好的双链结构捕获并杂交带有poly(a)序列的待测rna序列，形成第一杂交结构(s1)。第一杂交结构经逆转录形成rna/dna杂化双链结构(即第二杂交结构)(s2)。再将rna/dna杂化双链结构与ont测序接头连接，形成第三杂交结构(s3)。将第三杂交结构从带有马达蛋白的3’端序列开始测序，获得原始测序数据(s4)，进而获得包括ont测序接头、条形码、poly(a)和测试样本rna多核苷酸产生的一系列原始电流信号并可选地对其进行处理(例如提取、降噪、归一化)(s5)。然后，对条形码的电流信号数据进行分割，选择基于特定数量条形码的训练好的分类模型对电流信号数据进行预测，可选地选择恰当的过滤阈值(》0.3)，最终获得每条电信号的所属条形码和分类概率(s6)。
111.基于上述方法，本发明还提供了一种多条形码直接rna纳米孔测序系统100。
112.下面说明本发明的测序系统的示例性结构，在一些实施例中，如图8所示，所述测序系统100可以包括：
113.测序数据接收模块102，用于接收和存储测序数据，所述测序数据包括双链多核苷酸的测序数据，所述双链多核苷酸包括第一链和第二链，所述第一链从5’至3’方向包括rna多核苷酸、poly(a)序列、条形码序列、第一连接序列和测序接头，所述第二链从5’至3’方向包括测序接头、第二连接序列、条形码序列、poly(t)序列和所述rna多核苷酸的逆转录序列；所述测序数据包括测试样本数据和训练样本数据。
114.在一些实施例中，所述测试样本数据包括n种测试样本rna多核苷酸对应的测序数据，其中2≤n≤24且n为整数，所述条形码序列对每一种测试样本rna多核苷酸具有特异性，所述条形码序列选自seq id nos:1-24中的任意一种。
115.在一些实施例中，所述测试样本rna多核苷酸的来源包括病原体、植物、浮游动物、昆虫、哺乳动物和肿瘤中的一种或多种。在一些实施例中，所述病原体包括rna病毒、细菌、真菌和寄生虫中的一种或多种。
116.信号处理模块104，用于将所述测序数据处理成电流信号数据。
117.在一些实施例中，所述信号处理模块104对所述测序数据进行预处理和分割。
118.在一些实施例中，所述预处理包括降噪处理、抛光处理和归一化处理中的一种或多种。
119.在一些实施例中，所述分割为将所述条形码序列的电流信号数据分割成70-100列矩阵。
120.在一些实施例中，所述分割为将所述条形码序列的电流信号数据分割成100列矩阵。
121.分类模块106，用于调用基于所述测试样本数据中的所述条形码序列构建的分类
模型，对所述测试样本数据的电流信号数据进行分类并输出所述n种测试样本rna多核苷酸的分类结果。
122.在一些实施例中，所述测序系统100进一步包括模型训练模块108，所述模型训练模块108利用机器学习算法对所述训练样本数据的电流信号数据进行训练学习，从而确定所述分类模型。
123.在一些实施例中，所述机器学习算法包括k最近邻、神经网络、朴素贝叶斯分类、回归树、adaboost和随机森林中的一种或多种。
124.在一些实施例中，所述机器学习算法包括随机森林。
125.在一些实施例中，所述训练样本数据包括体外转录rna的测序数据。
126.在一些实施例中，所述训练样本数据中用作训练集和测试集的比例为7:3。
127.在一些实施例中，所述体外转录rna包括表2所示的基因中的一种或多种。
128.在一些实施例中，所述模型训练模块108用于将所述训练样本数据的电流信号数据与对应的参考序列比对，以比对质量大于60、比对至唯一参考序列为条件将不满足条件的读取序列过滤。
129.在一些实施例中，所述测序系统100进一步包括过滤模块110，用于基于设置的分类可能性的阈值，将低于所述阈值的读取过滤掉。
130.在一些实施例中，所述阈值》0.3。
131.在一些实施例中，所述过滤模块110和所述分类模块106可以整合为一个分析模块112。
132.1.5decoder的条形码识别分类流程
133.在一些实施例中，decoder的具体工作流程如下：
134.根据建库实验的条形码组合类型，选择不同的经过训练的模型，输入至model参数；
135.选择minknow软件保存并过滤后的原始电信号数据，一般在fast5_pass路径下，输入至fast5参数；
136.选择0-1范围内值作为阈值，对分类的标准进行过滤，输入至cutoff参数；
137.选择运行的cpu数目，不超过运算机器最大核心数，输入至nt参数。
138.在一些实施例中，decoder的默认使用方法如下：decoder(fast5,model,nt＝1,cutoff＝0,include.lowest＝false)，输出值将记录每条readid以及其所属的条形码信息和分类的可能性(probability，0-1)，如表5所示。
139.表5decoder输出结果示例
140.根据分类可能性过滤掉分类不太可靠的readid，并以readid为条件分类原始
fast5信号和fastq序列。
141.实施例二
142.2.1decoder算法与其他算法的比较
143.使用验证集评估了6种不同分类算法的准确性，包括k最近邻(k-nearest neighbour，knn)、神经网络(neural network，nnet)、朴素贝叶斯(bayes，nb)分类、回归树(classificationandregression tree，cart)、adaboost和随机森林(random forest，rf)(表6)。尽管6种方法都实现了接收操作特征曲线下的面积(area under the receiver operating characteristic curve，auroc)高于0.92(图1d)，并且相关的精确召回曲线下的面积(area under the precision-recall curve，auprc)大于0.83(图1e)，但随机森林算法优于其他方法，在6种方法中实现了最高的auroc(0.9961)和最高的auprc(0.9906)(图1d、图1e)。此外，随机森林算法的精确度(accuracy)、灵敏度(sensitivity)、特异性(specificity)、准确性(precision)、召回率(recall)和f1分数(f1 score)均显著高于其他五种方法(图1f，表6)。因此，本发明选择基于随机森林算法建立测序方法(简称为decoder)(图1g)，用于进一步的训练、测试和应用，分别针对rna病毒、细菌、真菌、寄生虫和肿瘤样本的rna进行了性能评估和实施例分析(图1h)。
144.表6多种算法的比较
145.2.2decoder与其他条形码解复用方法的比较
146.目前，已经有两种直接rna测序(drs)条形码解复用(demultiplex)计算方法，分别为基于残差神经网络分类器算法的图像识别的deeplexicon(https://github.com/psy-fer/deeplexicon)和尚未正式发表的poreplex(https://github.com/hyeshik/poreplex)。
147.为了比较decoder与其他方法的性能，本实施例评估了每种方法的拆分精度和运行时间。首先，分别对deeplexicon和poreplex的8种条形码连接ivtrna进行了5次独立的实验，获得了5次drs原始电信号数据集(表7)。接着，使用decoder对其中3次数据(训练集)进行训练模型，并使用该模型预测另外2次数据(验证集)。类似地，使用deeplexicon和poreplex对相同的验证集进行条形码拆分的预测。接收操作特征(receiver operating characteristic，roc)分析，曲线下面积(area under curve，auc)指标和精度召回(recall)曲线由r包multiroc(https://cran.r-project.org/package＝multiroc)产生。cpu时间由gnutime命令评估，并将用户时间和系统时间相加。
148.表7解复用方法比较数据集
149.与deeplexicon、poreplex两种方法相比，decoder在解复中使用poreplex预设的4种条形码中的分类准确度范围为95.82至97.83％(96.98
±
0.84)，而poreplex的分类准确度范围为83.4到93.18％(89.4
±
4.21)(图2a，unknown表示未鉴定)，而且准确性、召回率和f1-分数均高于poreplex方法。同时，decoder在解复用deeplexicon条形码中的分类准确率
范围为90.57至94.43％(92.29
±
1.60)，相比之下，deeplexicon的分类准确率范围为51.33至85.22％(68.39
±
15.16)(图2b，unknown表示未鉴定)。此外，与四个条形码的多路分解相比，decoder实现了显著更高的灵敏度、召回率和f1得分(图2b，表8)。decoder获得了比deeplexicon的0.946更高的准确度auc值0.997，与0.888相比更高的精度召回auc值0.993，以及显著更高的分类精确度(accuracy)和恢复率(图2c、图2d和图2e，表8)。最后，为了评估解复用的效率，还计算了decoder和deeplexicon在处理相同数量(1-150k)的drs读取时的cpu时间，发现decoder的运行速度比deeplexicon快10倍(图2f)。
150.表8decdoer与其他解复用方法比较
151.2.3decoder的性能效果
152.为了评估decoder在解复用(demultiplex)直接rna测序方面的性能，使用ivt-rna测试了它在解复用2、4、6、8、10和多达24个条形码时灵敏度、特异性、召回率等指标，解复用的全局精确度从98.4％下降到92.2％(图3c)，而特异性的所有auc值都大于0.99(图3a、图3b，表9)。精确度和恢复率虽然随着条形码数量从2个增加到24个有所下降，但始终保持在较高的水平(图3c、图3d)。
153.表9不同条形码下decoder的性能表现
154.鉴于decoder在多路分解ivt-rna的drs数据方面的理想性能，接下来将decoder应用于从病原体分离的真实的rna样本，包括2种rna病毒：senecavalley病毒(svv)和猪繁殖与呼吸综合征病毒(prrsv)；2种细菌：大肠杆菌(e.coli)和肠炎链球菌(s.enteritidis)；1种真菌：酿酒酵母(s.cerevisiae)；1种寄生虫：伯氏疟原虫(p.berghei)。对不同的rna样本进行了三次decoder建库实验和直接rna测序，并生成了三批数据。roc分析显示，所有三次
重复实验具有高auc值的特异性(0.999、0.998、0.976)，灵敏度分别为99.2％、93.3％、90.6％，特异性分别为99.4％、98.4％、98.4％和最大精确度分别为97.6％、95.9％、92.5％(图4a，表10)。decoder在所有三批数据中也表现出高auc值的精确召回率(0.998、0.994、0.937)(图4b)。在阈值(cutoff)和精确度(accuracy)方面，我们发现cutoff值越高，精确度越高，而召回率越低(图4c)。样本不平衡问题导致这些模型的精度值和f1分数不同，严重的样本不平衡问题与低精度和f1分数相关(表10)。
155.表10真实样本下decoder的性能表现
156.实施例三
157.3.1decoder构建rna病毒参考基因组
158.为了评估decoder是否可以促进病毒基因组的构建，对三次decoder建库实验的rna病毒数据进行了参考基因组比对和比对reads长度的统计，发现decoder测得的reads具有几乎完全覆盖svv和prrsv的参考基因组的能力(图5a、图5b)，这使得rna病毒基因组的组装变得快速和容易。具体而言，724条reads(9.68％，实验1)和111条reads(1.42％，实验2)接近完整的svv基因组(》7,000nt)，42条reads(2.26％)接近完整的prrsv基因组(》15,000nt)(图5b)。接下来分析了svv基因组中的单碱基突变或单核苷酸多态性(snp)，并在实验1和实验2中分别发现了414和423个突变，其中411个是2次实验中的共有突变(图5e)，并且这些突变的突变比率具有高度可重复性(图5c)。举例来说，本实施例检测到了svv病毒中5700位点c到t的突变，在两次实验中确实存在，且突变频率分别为0.98和0.99(图5d)。decoder也具有检测高频缺失(deletion)的能力，svv的缺失(例如4022-4024位点的3nt缺失)在两次实验中均被准确识别(图5d)。总体而言，decoder可用于以低成本效益的方式检测病原体rna，并有可能应用于鉴定新型或突变型的rna病毒。
159.3.2decoder鉴定的肿瘤样本的突变和m6a修饰
160.为了举例说明decoder在癌症中的潜在应用，对弥漫性中线胶质瘤-h3k27m突变体(dmg)和胶质母细胞瘤(gbm)共4个样本进行了decoder混样建库实验和rna直接测序，并获得了全转录组的突变、isoform鉴定和rnam6a修饰的概况(图6a)。在突变检测分析中，共有17,587个单碱基突变在dmg和gbm不同样本中被decoder检测到，并且其中97-98％的突变被直接rna测序(drs)得到验证(图6b)。与二代测序(ngs)和sanger测序相比，decoder在检测肿瘤突变方面表现出相似的准确性。在通过sanger测序验证的突变位点中，decoder和ngs之间检测到的共有突变的突变频率存在高度相关性(r＝0.93，pearson相关性，图6c)。举例来说，hla-a基因上的chr6:29913430:c》t突变位点在三个dmg样本中表现出较高的突变频率，在dmg1、dmg2和dmg3中分别为0.87、0.46和0.84，而gbm样本中的突变频率为0，同时在ngs和sanger测序中也检测到突变位点，突变频率与decoder相似(图6d)。相反的是，atp6v0c基因上的chr16:2519933:c》a突变位点在gbm中具有较高的突变频率，而在dmg三个样本中几乎不突变，同样得到了ngs和sanger测序的验证(图6e)。decoder还鉴定了不同肿瘤类型中的保守缺失区域，例如snord137基因上chr9:12972603-12972607位点具有5nt的
缺失，在所有肿瘤样本中保守出现(图6f)。
161.在m6a修饰分析中，decoder在dmg和gbm不同样本中共检测到2,158个高度置信的m6a位点，其中39％得到drs和m6aaltas公共数据库的支持，60％仅通过drs验证，只有1％未验证(图6g)。decoder在基因表达和m6a检测方面也显示出与drs的高度一致性，相关性分别为0.91和0.92(图6h)。总体而言，decoder可用于肿瘤的疾病样本的突变和m6arna修饰的检测，而且检出结果真实可靠，被一代测序、二代测序验证，并且与直接rna测序一致性较高。
162.通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到上述实施例方法可借助软件加必需的通用硬件平台的方式来实现，当然也可以通过硬件，但很多情况下前者是更佳的实施方式。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质(如rom/ram、磁碟、光盘)中，包括若干指令用以使得一台计算机终端(可以是手机，计算机，服务器，或者网络设备等)执行本发明各个实施例所述的方法。
163.上面结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。

技术特征：
1.一种直接rna纳米孔测序方法，其特征在于，所述方法用于对n种测试样本rna多核苷酸进行测序，其中2≤n≤24且n为整数，所述方法包括：s1针对每一种测试样本rna多核苷酸，提供一种双链dna多核苷酸来捕获一种待测rna序列并形成第一杂交结构；所述双链dna多核苷酸为逆转录接头的前向链和后向链组装成的双链结构，所述前向链从5’至3’方向包括条形码序列和第一连接序列，所述后向链从5’至3’方向包括第二连接序列、条形码序列和poly(t)序列，所述待测rna序列包括所述测试样本rna多核苷酸和poly(a)序列，所述双链dna多核苷酸通过所述poly(t)序列捕获所述待测rna序列；所述条形码序列对所述每一种测试样本rna多核苷酸具有特异性，所述条形码序列选自seq id nos:1-24中的任意一种；s2将所述第一杂交结构逆转录，形成第二杂交结构；s3将测序接头与所述第二杂交结构连接，形成第三杂交结构；s4对所述第三杂交结构进行纳米孔测序，获得原始测序数据；s5提取所述原始测序数据对应的电流信号数据；s6将所述电流信号数据中所述条形码序列的电流信号数据进行分割，并调用基于所述条形码序列训练得到的分类模型进行分析，获得所述条形码序列的序列信息，进而获得所述n种测试样本rna多核苷酸的分类信息。2.如权利要求1所述的方法，其特征在于，所述分类模型的训练方法包括将训练样本数据中条形码序列的电流信号数据进行分割，利用机器学习算法对分割后的条形码序列的电流信号数据进行训练学习，从而确定所述分类模型，其中所述训练样本数据包括一种或多种训练样本rna多核苷酸和对应的条形码序列的电流信号数据，所述条形码序列选自seq id nos:1-24中的任意一种。3.如权利要求1或2所述的方法，其特征在于，所述分割为将所述条形码序列的电流信号数据分割成70-100列矩阵。4.如权利要求2所述的方法，其特征在于，所述训练样本rna多核苷酸包括体外转录rna。5.如权利要求2所述的方法，其特征在于，所述机器学习算法包括k最近邻、神经网络、朴素贝叶斯分类、回归树、adaboost和随机森林中的一种或多种。6.如权利要求1所述的方法，其特征在于，所述测试样本rna多核苷酸的来源包括病原体、植物、浮游动物、昆虫、哺乳动物和肿瘤中的一种或多种。7.如权利要求1所述的方法，其特征在于，所述分类的分类可能性的阈值>0.3。8.如权利要求1或2所述的方法，其特征在于，所述分割为将所述条形码序列的电流信号数据分割成100列矩阵。9.如权利要求1所述的方法，其特征在于，所述s5进一步包括提取所述原始测序数据对应的原始电流信号数据，并对所述原始电流信号数据进行预处理，获得预处理后的电流信号数据，所述预处理包括降噪处理、抛光处理和归一化处理中的一种或多种。10.如权利要求2所述的方法，其特征在于，所述方法进一步包括将所述一种或多种训练样本rna多核苷酸的电流信号数据与对应的参考序列比对，以比对质量大于60、比对至唯一参考序列为条件将不满足条件的读取序列过滤。

技术总结
本发明涉及测序领域，具体涉及一种多重混样直接RNA纳米孔测序方法。本发明提供的方法用于对n种测试样本RNA多核苷酸进行测序(2≤n≤24且n为整数)，包括：对所述n种测试样本RNA多核苷酸对应的测序结构进行纳米孔测序并提取对应的电流信号数据；将其中的条形码序列的电流信号数据进行分割，并调用基于所述条形码序列训练得到的分类模型进行分析，获得所述条形码序列的序列信息，进而获得所述n种测试样本RNA多核苷酸的分类信息；其中所述条形码序列对每一种测试样本RNA多核苷酸具有特异性，所述条形码序列选自SEQIDNOs:1-24中的任意一种。本发明提供的方法可以实现在一张芯片中同时对多种生物样本测序，降低了生物样本起始投入量和多样本测序成本。入量和多样本测序成本。入量和多样本测序成本。


技术研发人员：陈路 林静雯 宋俊伟 唐超 陈钏 耿佳
受保护的技术使用者：四川大学
技术研发日：2023.05.30
技术公布日：2023/8/24