CN106033502A - 鉴定病毒的方法和装置 - Google Patents

Abstract

本发明公开了一种鉴定病毒的方法，该方法包括：获取待测样本的RNA测序数据；对测序数据的第一部分进行组装，获得组装序列；将测序数据的第一部分与组装序列进行比对，获得比对结果；基于比对结果，确定组装序列上的突变位点，以及确定各组装序列的(a)‑(c)中的至少之一：(a)平均熵值和中位数熵值至少之一，以及突变位点比例，(b)平均突变率和中位数突变率至少之一，以及突变位点比例，(c)突变位点比例；将组装序列的(a)‑(c)至少之一和与其对应的界限比较，判定落入界限的组装序列来自病毒。本发明还公开一种鉴定病毒的装置。利用本发明的方法和/或装置鉴定病毒，能够不依赖同源序列比对来准确预测未知序列是否为病毒序列。

Description

鉴定病毒的方法和装置

技术领域

本发明涉及生物检测领域，具体的，本发明涉及一种病毒鉴定方法和一种病毒鉴定装置。

背景技术

截止2014年6月30日，国际病毒学分类大会(ICTV)公布的病毒种类有2827种。而2011年是有2484种，2009年是有2285种，每年差不多100种病毒被发现，这个速度远比不上细菌等其他微生物。这是因为病毒本身不能生长，需要寄生在宿主细胞里，很难分离培养，而有研究者预测，只要有细胞就会有病毒，病毒在自然界存在是非常广泛和巨大的，现在已知的病毒，还不及总量的万分之一。随着测序技术的发展，成本不断降低，通量不断提高，使得越来越多的物种被测序，越来越多的样本被发现存在病毒，而病毒种类的数据也在不断增加。测序技术在病毒发现的应用越来越常见(Barzon,L.,E.Lavezzo,etal.(2011)."Applications of next-generation sequencing technologies to diagnostic virology."Int JMol Sci 12(11):7861-7884.)，并且新一代的测序技术在对于样本中一些低丰度的病毒亚型有着独特的先天优势。对于发现新病毒，很典型的例子就是在一个不明原因发热最后死亡的病人身上，传统检测并没有发现致病病原，在病人死后取样，利用新一代的测序技术产生了巨大的数据，之后发现了一株新沙粒病毒才是导致病人发病的原因(Palacios,G.,J.Druce,et al.(2008)."A new arenavirus in a cluster of fatal transplant-associated diseases."NEngl J Med 358(10):991-998.)。同样的技术运用，在非洲发现了一株新的致命病毒Lujovirus(Briese,T.,J.T.Paweska,et al.(2009)."Genetic detection and characterization of Lujovirus,a new hemorrhagic fever-associated arenavirus from southern Africa."PLoS Pathog 5(5):e1000455.)。除了在新病毒的发现外，对于已知的病毒也是有着明显的优势，比如2014年在非洲大爆发的埃博拉疫情。研究人员能在短时间内获得99株最新的埃博拉病毒，并且迅速分析了大量的变异位点，最终确定了传播的可能路径和变异(Gire,S.K.,A.Goba,et al.(2014)."Genomic surveillance elucidates Ebola virus origin and transmission during the 2014outbreak."Science 345(6202):1369-1372.)。

新一代测序技术在病毒发现和分析上，相比传统的技术有了很大通量和阳性率的提高，但是很大程度上依赖于现有病毒序列，要么和已有病毒序列比较相近，发现了有较大的变异、重组，要么和已有的病毒序列相差较大，但是仍然在同一个属水平。目前所用的方法基本都是基于同源序列的比对，对已经构建好的数据库的依赖较大。

同样基于新一代测序技术的基础上，利用小RNA(small RNA)的数据，吴等人通过对小RNA的拼接，从延长的序列里挑出来了病毒的序列(Wu,Q.,Y.Luo,et al.(2010)."Virusdiscovery by deep sequencing and assembly of virus-derived small silencing RNAs."Proc NatlAcad Sci U S A 107(4):1606-1611.)，拼接完后，基于同源比对来鉴定。

前面提到的small RNA在细胞生物里是广泛存在的，包括microRNAs(miRNAs)，小干扰(RNAsmall interfering RNAs(siRNAs)和Piwi-interacting RNAs(piRNAs)，或多或少，直接或者间接的参与了细胞基因表达的调节，保护细胞抵御外来入侵(Ghildiyal,M.and P.D.Zamore(2009)."Small silencing RNAs:an expanding universe."Nat Rev Genet 10(2):94-108.)。在真菌、植物和无脊椎动物中，病毒介导的小RNA(virus-derived small interfering RNAs，vsiRNAs)也是自身免疫系统的重要组成部分(van Mierlo,J.T.,K.W.van Cleef,et al.(2010)."Small Silencing RNAs:Piecing Together a Viral Genome."Cell Host Microbe 7(2):87-89.)。Small RNAs能够直接来源于直接降解的mRNA，或者单链的病毒RNA，和DNA或者mRNA结合后又能影响转录和表达Mlotshwa,S.,G.J.Pruss,et al.(2008)."Small RNAs inviral infection and host defense."Trends Plant Sci 13(7):375-382.)。

发明内容

本发明旨在至少在一定程度上解决上述技术问题之一或至少提供一种商业选择。为此，本发明的目的在于提出鉴定病毒的手段。

依据本发明的一方面，本发明提供一种鉴定病毒的方法，该方法包括：获取待测样本的RNA测序数据，所述测序数据包括多个读段；对所述测序数据的第一部分进行组装，获得组装序列，所述测序数据的第一部分包括不能比对上ncRNA参考序列的读段；将所述测序数据的第一部分与所述组装序列进行比对，获得比对结果；基于所述比对结果，确定所述组装序列上的突变位点，所述突变位点包括SNV，以及确定各条组装序列的(a)-(c)中的至少之一，(a)平均熵值和中位数熵值至少之一，以及突变位点比例，(b)平均突变率和中位数突变率至少之一，以及突变位点比例，(c)突变位点比例，一条组装序列的平均熵值为其上的SNV的熵值的平均值，一条组装序列的中位数熵值为其上的SNV的熵值的中位数，一条组装序列的平均突变率为其上的SNV的突变率的平均值，一条组装序列的中位数突变率为其上的SNV的突变率的中位数，一条组装序列的突变位点比例为其上的SNV的数目所占的比例，SNV的熵值＝-100*∑(Pi*logPi)，Pi为该SNV的各种碱基的深度，SNV的突变率＝支持该SNV的读段数目/比对上该SNV的读段数目；将所述确定的组装序列的(a)-(c)至少之一和与其对应的界限比较，判定落入所述界限的组装序列来自病毒。所称的测序数据通过测序获得，测序方法依据测序平台的不同可选择但不限于Illumina公司的Hisq2000/2500测序平台、Life Technologies公司的Ion Torrent平台和单分子测序平台，测序方式可以选择单端测序，也可以选择双末端测序，获得的下机数据是测读出来的片段，称为读段(reads)。ncRNA(non-coding RNA)为非编码RNA，所称的ncRNA参考序列包括：rRNA参考序列、tRNA参考序列、snRNA参考序列、snoRNA参考序列和microRNA参考序列中的至少之一，所称的参考序列指预先确定的序列，可以是预先获得的待测样本所属或者所包含的生物类别的任意参考模板，例如，若待测样本来源的个体为人类，参考序列可选择NCBI数据库提供的HG19，若待测样本来源的个体宿主为人类，目标核酸是病毒RNA中的ncRNA，参考序列可选择Rfam数据库和/或GeneBank中的人非编码RNA序列，利于将宿主的已知ncRNA和剩余的ncRNA区分开，进一步地，也可以预先配置包含更多参考序列的资源库，例如依据待测样本来源个体的状态、地域等因素选择或是测定组装出更接近的序列作为参考序列。SNV为单核苷酸变异，同SNP，SNV一般由两种碱基组成，是一种二等位基因。所称的深度，即测序深度，指碱基被测序或者读取的平均次数。比对可以利用已知比对软件进行，例如SOAP、BWA和TeraMap等，在比对过程中，一般对比对参数进行设置，设置一条reads最多允许有s个碱基错配(mismatch)，例如设置s≤4，若reads中有超过s个碱基发生错配，则视为该reads无法比对到(比对上)参考序列，而所称的支持某个SNV的读段指比对上该位点的读段的相应位置的碱基与该位点的一致。所称的界限是属于病毒的界限，界限可能是个临界值、数值范围或关系式，但能够界定出病毒。与(a)-(c)中至少之一对应的界限是利用感染病毒的样本的测序数据来确定的，与(a)对应的界限包括平均熵值界限和中位数熵值界限至少之一、以及突变位点比例界限，与(b)对应的界限包括平均突变率界限和中位数突变率界限至少之一、以及突变位点比例界限，与(c)对应的界限包括突变位点比例界限。(a)-(c)共包含组装序列的非重复的五个要素，五个要素对应的界限都可根据已知感染有病毒的样本的核酸测序数据来确定。

本发明的这一方面的病毒鉴定方法，实现非基于同源序列比对来鉴定病毒，且适用于任何新一代高通量测序产出的数据，能有效提高病毒鉴定的检出率和准确率。

根据本发明的一个实施例，该方法中的测序数据的第一部分不包括符合以下(1)-(3)至少之一的读段：(1)包含接头序列，和/或平均单碱基错误率大于0.01，(2)长度小于18nt，和/或长度大于44nt，(3)包含碱基质量值小于10的碱基个数大于2。过滤掉被测序接头污染的和低质量的读段，读段的可靠性提高利于准确鉴定病毒。

根据本发明的一个实施例，该方法中的组装包括：分别对所述测序数据的第一部分中的全部读段和所述测序数据的第一部分中的至少一部分读段进行第一组装和第二组装，获得第一组装结果和第二组装结果；合并所述第一组装结果和第二组装结果。任选的，所述第一组装和/或所述第二组装为混合组装。所称的混合组装指采用多种Kmer分别组装，再合并不同Kmer的组装结果。Kmer指长度为K的一段序列，一条长度为L的reads产生的Kmer数量为L-K+1，组装时，将每个读段分解成其包含的所有长度为K的固定序列，利于快速组装。基于不同数据量和采用不同Kmer分别进行组装，再合并各个组装结果，能够提高组装的精确性。

根据本发明的一个实施例，在获得组装序列之后，去除所述组装序列中的已知序列。例如，通过与公开数据库中收录的物种信息确定的序列进行比对，将比对上的组装序列去除，判定是否比对上可以依据所使用的比对软件或程序的默认设置或者设置软件或程序的参数数值定义。组装序列中的已知序列由于公开数据中已经披露了其来源，可直接鉴定出其来源，无需进行后续步骤来鉴定，去除这些序列减少了数据量，利于利用机器快速运行该病毒鉴定方法。

根据本发明的一个实施例，在确定每条组装序列的(a)-(c)中的至少之一之前，去除该组装序列中深度小于100X的位点。所称的位点为组成该条组装序列的各个位点，包括SNP位点，也包括非SNP位点，将组装序列上的深度小于100X的位点排除在任一(a)-(c)的计算确定过程，有利于准确进行病毒鉴定。

根据本发明的一个实施例，该方法还包括：设计引物，对判定来自病毒的组装序列进行RT-PCR延伸，获得延伸产物；利用延伸产物验证所述组装序列来自病毒。利用判定来自病毒的组装序列进行引物设计，RT-PCR对这条或这些组装序列进行延伸和/或连接，以验证判定的组装序列的确来自病毒。

根据本发明的一些实施例，确定该方法中的界限，包括：获取至少一个已知病毒感染的样本的RNA测序结果，所述测序结果包括多个读段；对所述测序结果的至少一部分进行组装，获得组装片段；将所述组装片段与参考序列进行第一比对，进行组装片段物种注释，获得物种注释结果；基于所述物种注释结果对所述组装片段进行分类，获得第一类组装片段和第二类组装片段，所述第一类组装序列来自病毒，所述第二类组装片段来自宿主；分别将所述测序结果的至少一部分与所述第一类组装片段和所述第二类组装片段进行第二比对，获得第二比对结果；基于所述第二比对结果，确定所述第一类组装片段和所述第二类组装片段上的突变位点，所述突变位点包括SNV，以及确定每条第一类组装片段和每条第二类组装片段的五个因素，所述五个因素为平均熵值、中位数熵值、平均突变率、中位数突变率以及突变位点比例，其中，一条组装片段的平均熵值为其上的SNV的熵值的平均值，一条组装片段的中位数熵值为其上的SNV的熵值的中位数，一条组装片段的平均突变率为其上的SNV的突变率的平均值，一条组装片段的中位数突变率为其上的SNV的突变率的中位数，一条组装片段的突变位点比例为其上的SNV的数目所占的比例，SNV的熵值＝-100*∑(Pi*logPi)，Pi为该SNV的各种碱基的深度，SNV的突变率＝支持该SNV的读段数目/比对上该SNV的读段数目；基于所述第一类组装片段和第二类组装片段、以及每条所述第一类组装片段和每条所述第二类组装片段的五个因素，利用SVM，确定所述五个因素中至少之一的界限。支持向量机(Support Vector Machine，SVM)是比较流行的用来数据分类的方法，它的基本法则可看成，将给定数据定为+1或者-1，然后根据给定数据之外的数据和给定数据的相似程度，判断给定数据之外的数据应归为+1还是-1，例如，根据物种注释结果，将第一类组装片段标记为病毒(+1)，将第二类组装片段标记为宿主(-1)，使用LibSVM软件包(Chih-Chung Chang,C.-J.L.(2011)."LIBSVM:a library for support vectormachines."ACM Transactions on Intelligent Systems and Technology 2:27:1--27:27.)，使用SVC(support vector classification)算法来找出规律即五个要素的界限来区分病毒和宿主序列。在本发明的一个实施例中，要求所述第一类组装片段的数目大于30，利于最后确定出的界限能用于准确分类鉴定。在本发明的一个实施例中，在确定界限后，对确定的界限的分类可靠性进行ROC曲线评估，要求其AUC值不小于0.7。在本发明的一个实施例中，要求确定出的界限对病毒预测的准确率达90％。

根据本发明的一个实施例，所述基于第一类组装片段和第二类组装片段、以及每条第一类组装片段和每条第二类组装片段的五个因素，利用SVM，确定所述五个因素中至少之一的界限，包括：将所述第一类组装片段和第二类组装片段分为多份片段包，每份所述片段包包含多条第一组装片段和多条第二组装片段，利用其中一份片段包中的组装片段的五个因素的至少之一，判断其它各份片段包中的每条组装片段为第一组装片段还是第二组装片段，依据判断的正确率，确定所述五个因素的至少之一的界限。

依据本发明的另一方面，本发明提供一种鉴定病毒的装置，包括：数据输入单元，用于输入数据；数据输出单元，用于输出数据；存储单元，用于存储数据，其中包括计算机可执行程序；处理器，与所述数据输入单元、数据输出单元和存储单元连接，用于执行所述程序，执行所述程序包括完成上述鉴定病毒的方法。

依据本发明的再一方面，本发明提供一种计算机可读存储介质，用于存储供计算机执行的程序，本领域普通技术人员可以理解，在执行该程序时，通过指令相关硬件可完成上述鉴定病毒的方法的全部或部分步骤。所称存储介质可以包括：只读存储器、随机存储器、磁盘或光盘等。

本发明的方法和/或装置，基于小RNA(small RNAs)，包括病毒介导的小RNA(virus-derived small interfering RNAs，vsiRNAs)，基于病毒和宿主之间的作用机制，利用机器学习和模拟的办法确定未知序列五个要素中至少一个要素的病毒界限，实现非基于同源序列比对来鉴定病毒，且适用于任何新一代高通量测序产出的数据，能有效提高病毒鉴定的检出率和准确率。

附图说明

本发明的上述和/或附加的方面和优点从结合下面附图对实施方式的描述中将变得明显和容易理解，其中：

图1是本发明的一个实施例中的病毒鉴定方法包含的步骤的示意图。

图2是本发明的一个实施例中的病毒鉴定方法的流程图。

图3是本发明的一个实施例中的测序饱和度评估结果示意图。

图4是本发明的一个实施例中的不同组装方式对最大contig影响的评估结果示意图。

图5是本发明的一个实施例中的组装序列的聚类示意图。

图6是本发明的一个实施例中的训练集评估的ROC曲线。

图7是本发明的一个实施例中的病毒序列试验验证示意图。

具体实施方式

下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中，自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，仅用于解释本发明，而不能理解为对本发明的限制。需要说明的，本文中所使用的术语“第一”、“第二”或者“第一部分”等仅为方便描述，不能理解为指示或暗示相对重要性，也不能理解为之间有先后顺序关系。在本发明的描述中，除非另有说明，“多个”的含义是两个或两个以上。在本文中，除非另有明确的规定和限定，术语“相连”、“连接”等术语应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。

如图1所示，依据本发明的一方面，本发明提供一种鉴定病毒的方法，该方法包括：S10获取待测样本的RNA测序数据，所述测序数据包括多个读段；S20对所述测序数据的第一部分进行组装，获得组装序列，所述测序数据的第一部分包括不能比对上ncRNA参考序列的读段；S30将所述测序数据的第一部分与所述组装序列进行比对，获得比对结果；S40基于所述比对结果，确定所述组装序列上的突变位点，所述突变位点包括SNV，以及确定各条组装序列的(a)-(c)中的至少之一：(a)平均熵值和中位数熵值至少之一，以及突变位点比例，(b)平均突变率和中位数突变率至少之一，以及突变位点比例，(c)突变位点比例，一条组装序列的平均熵值为其上的SNV的熵值的平均值，一条组装序列的中位数熵值为其上的SNV的熵值的中位数，一条组装序列的平均突变率为其上的SNV的突变率的平均值，一条组装序列的中位数突变率为其上的SNV的突变率的中位数，一条组装序列的突变位点比例为其上的SNV的数目所占的比例，SNV的熵值＝-100*∑(Pi*logPi)，Pi为该SNV的各种碱基的深度，SNV的突变率＝支持该SNV的读段数目/比对上该SNV的读段数目；S50将所述确定的组装序列的(a)-(c)至少之一和与其对应的界限比较，判定落入所述界限的组装序列来自病毒。其中，所称的测序数据通过测序获得，测序方法依据测序平台的不同可选择但不限于Illumina公司的Hisq2000/2500测序平台、LifeTechnologies公司的Ion Torrent平台和单分子测序平台，测序方式可以选择单端测序，也可以选择双末端测序，获得的下机数据是测读出来的片段，称为读段(reads)。ncRNA(non-coding RNA)为非编码RNA，所称的ncRNA参考序列包括：rRNA参考序列、tRNA参考序列、snRNA参考序列、snoRNA参考序列和microRNA参考序列中的至少之一，所称的参考序列指预先确定的序列，可以是预先获得的待测样本所属或者所包含的生物类别的任意参考模板，例如，若待测样本来源的个体为人类，参考序列可选择NCBI数据库提供的HG19，若待测样本来源的个体宿主为人类，目标核酸是病毒RNA中的ncRNA，参考序列可选择Rfam数据库和/或GeneBank中的人非编码RNA序列，利于将宿主的已知ncRNA和剩余的ncRNA区分开，进一步地，也可以预先配置包含更多参考序列的资源库，例如依据待测样本来源个体的状态、地域等因素选择或是测定组装出更接近的序列作为参考序列。SNV为单核苷酸变异，同SNP，SNV一般由两种碱基组成，是一种二等位基因。所称的深度，即测序深度，指碱基被测序或者读取的平均次数。比对可以利用已知比对软件进行，例如SOAP、BWA和TeraMap等，在比对过程中，一般对比对参数进行设置，设置一条reads最多允许有s个碱基错配(mismatch)，例如设置s≤4，若reads中有超过s个碱基发生错配，则视为该reads无法比对到(比对上)参考序列，而所称的支持某个SNV的读段指比对上该位点的读段的相应位置的碱基与该位点的一致。所称的界限是属于病毒的界限，界限可能是个临界值、数值范围或关系式，但能够界定出病毒。与(a)-(c)中至少之一对应的界限是利用感染病毒的样本的测序数据来确定的，与(a)对应的界限包括平均熵值界限和中位数熵值界限至少之一、以及突变位点比例界限，与(b)对应的界限包括平均突变率界限和中位数突变率界限至少之一、以及突变位点比例界限，与(c)对应的界限包括突变位点比例界限。(a)-(c)共包含组装序列的非重复的五个要素，五个要素至少之一或者其全部或者部分组合能反映出病毒序列自己的特性，虽然非病毒序列可能特性并不单一，但至少和病毒序列是有差别的，序列的五个要素的至少之一或者其全部或部分的组合能够用以判断该序列是否来自病毒，五个要素对应的界限都可根据已知感染有病毒的样本的核酸测序数据来确定。在本发明的一个实施例中，确定了组装序列的(a)-(c)中的至少之二或者全部五个要素。本发明的这一方面的病毒鉴定方法，实现非基于同源序列比对来鉴定病毒，且适用于任何新一代高通量测序产出的数据，能有效提高病毒鉴定的检出率和准确率。

根据本发明的一个实施例，该方法中的测序数据的第一部分不包括符合以下(1)-(3)至少之一的读段：(1)包含接头序列，和/或平均单碱基错误率大于0.01，(2)长度小于18nt，和/或长度大于44nt，(3)包含碱基质量值小于10的碱基个数大于2。包含接头序列指该reads包含测序接头，测读的不是目标区域，为被测序接头污染的读段。碱基错误率和碱基质量值为测序平台赋予读段的值，质量值为-10*lg(p)，这里，p为测错的概率，即单碱基错误率，当一条reads某位置碱基的出错概率为0.1时，其质量值为10。过滤掉被测序接头污染的和低质量的读段，能使测序数据包含的reads整体质量提高，利于减少后续分析鉴定对机器内存的需求，而且读段的可靠性提高利于准确鉴定病毒。

由于小RNA的高通量测序产生的数据庞大，长度都小于50nt，并且混杂了不同类型的小RNA，如：miRNA，siRNA，vsiRNA等，还有属于宿主的各种RNA片段。对于这样的混合数据进行组装，尤其是基于Kmer方式的组装，单一软件都有一定的差别。为提高组装的精确性，根据本发明的一个实施例，该方法中的组装包括：分别对所述测序数据的第一部分中的全部读段和所述测序数据的第一部分中的至少一部分读段进行第一组装和第二组装，获得第一组装结果和第二组装结果；合并所述第一组装结果和第二组装结果。任选的，所述第一组装和/或所述第二组装为混合组装。组装可以利用已知序列组装方法，例如利用soapdenovo、velvet等。所称的混合组装指多个组装混合，如采用多种Kmer分别组装，再合并不同Kmer的组装结果，和/或采用多个数据量分别进行组装。Kmer指长度为K的一段序列，一条长度为L的reads产生的Kmer数量为L-K+1，组装时，将每个读段分解成其包含的所有长度为K的固定序列，利于快速组装。在本发明的一个实施例中，基于所述测序数据的第一部分出4份不同等分数据量分别进行组装，如取25％的测序数据的第一部分数据量，50％的测序数据的第一部分的数据量，75％的测序数据的第一部分的数据量和100％的测序数据得第一部分数据量分别进行组装，以分辨不同等分对组装出的contigs的影响。基于不同数据量和采用不同Kmer分别进行组装，获得不同的组装结果，能够提高组装效果，提高组装的精确性。

根据本发明的一个实施例，在获得组装序列之后，去除所述组装序列中的已知序列。例如，通过与公开数据库中收录的物种信息确定的序列进行比对，将比对上的组装序列去除，判定是否为比对上可以依据所使用的比对软件或程序的默认设置或者设置软件或程序的参数数值定义。组装序列中的已知序列由于公开数据中已经披露了其来源，可直接鉴定出其来源，可以不用进行后续步骤进行鉴定，去除这些序列减少了数据量，利于利用机器快速运行该病毒鉴定方法。在本发明的一个实施例中，使用BLASTn和BLASTx将组装序列与NCBI的nt库和nr库分别进行比对，进行物种注释。与Nt库的比对在同源性85％以上以及序列80％以上能比对上库判定为比对阳性，将比对阴性的序列进行nr库比对，将同源性50％以上以及序列50％以上比上库判定为阳性，这样，比对结果可分成三组：比对上nt库(nt阳性)；没有比对上nt库但比对上nr库(nt阴性但nr阳性)；nt和nr都阴性。nt和nr都阴性的组装序列为完全未知序列，基于目前根据同源比对方法无法确定其来源。

根据本发明的一个实施例，在确定每条组装序列的(a)-(c)中的至少之一之前，去除该组装序列中深度小于100X的位点。所称的位点为组成该条组装序列的各个位点，包括SNP位点，也包括非SNP位点，将组装序列上的深度小于100X的位点排除在任一(a)-(c)的计算确定过程，有利于准确进行病毒鉴定。

根据本发明的一个实施例，该方法还包括：设计引物，对判定来自病毒的组装序列进行RT-PCR延伸，获得延伸产物；利用延伸产物验证所述组装序列来自病毒。利用判定来自病毒的组装序列进行引物设计，RT-PCR对这条或这些组装序列进行延伸，任选的对延伸产物进行顺序确定和连接，以验证判定来自病毒的组装序列来自病毒。

据本发明的一些实施例，确定该方法中的界限，包括：获取至少一个已知病毒感染的样本的RNA测序结果，所述测序结果包括多个读段；对所述测序结果的至少一部分进行组装，获得组装片段；将所述组装片段与参考序列进行第一比对，进行组装片段物种注释，获得物种注释结果；基于所述物种注释结果对所述组装片段进行分类，获得第一类组装片段和第二类组装片段，所述第一类组装序列来自病毒，所述第二类组装片段来自宿主；分别将所述测序结果的至少一部分与所述第一类组装片段和所述第二类组装片段进行第二比对，获得第二比对结果；基于所述第二比对结果，确定所述第一类组装片段和所述第二类组装片段上的突变位点，所述突变位点包括SNV，以及确定每条第一类组装片段和每条第二类组装片段的五个因素，所述五个因素为平均熵值、中位数熵值、平均突变率、中位数突变率以及突变位点比例，其中，一条组装片段的平均熵值为其上的SNV的熵值的平均值，一条组装片段的中位数熵值为其上的SNV的熵值的中位数，一条组装片段的平均突变率为其上的SNV的突变率的平均值，一条组装片段的中位数突变率为其上的SNV的突变率的中位数，一条组装片段的突变位点比例为其上的SNV的数目所占的比例，SNV的熵值＝-100*∑(Pi*logPi)，Pi为该SNV的各种碱基的深度，SNV的突变率＝支持该SNV的读段数目/比对上该SNV的读段数目；基于所述第一类组装片段和第二类组装片段、以及每条所述第一类组装片段和每条所述第二类组装片段的五个因素，利用SVM，确定所述五个因素中至少之一的界限。支持向量机(Support Vector Machine，SVM)是比较流行的用来数据分类的方法，它的基本法则可看成，将给定数据定为+1或者-1，然后根据给定数据之外的数据和给定数据的相似程度，判断给定数据之外的数据应归为+1还是-1。一般的，可以随机选出一些数据作为模型(训练集)，发现其规律，利用其它数据作为测试集来计算检测这个模型的规律与训练集的误差等，从而确定发现的规律是否正确。在本发明的一个实施例，所述基于第一类组装片段和第二类组装片段、以及每条第一类组装片段和每条第二类组装片段的五个因素，利用SVM，确定所述五个因素中至少之一的界限，包括：将所述第一类组装片段和第二类组装片段分为多份片段包，每份所述片段包包含多条第一组装片段和多条第二组装片段，利用其中一份片段包中的组装片段的五个因素的至少之一，判断其它各份片段包中的每条组装片段为第一组装片段还是第二组装片段，依据判断的正确率，调整确定出所述五个因素的至少之一的界限。具体地，根据物种注释结果，将每个样本的contig标记为属于病毒(+1)和属于宿主(-1)，使用LibSVM软件包实现基于SVM的分类和回归(Chih-Chung Chang,C.-J.L.(2011)."LIBSVM:a library for support vector machines."ACMTransactions on Intelligent Systems and Technology 2:27:1--27:27.)，使用SVC(support vectorclassification)算法来区分病毒和宿主序列，为优化分类，在建立训练集的时候，罚分参数C和RBF(Radical Basis Function)内核参数γ的设定将根据5倍交叉验证来确定，即将用来做训练集的数据分5等份，第一次取第一等份的训练集的数据来预测剩下4份的数据，第二次取第二等份的数据预测其他4份数据，依次进行5次循环，确定出分类界限。

在不同实施例中，由于构建的训练集的不同，例如选取的已知病毒感染的样本不同、样本数目的不同，使得获得的组装片段、物种注释结果及组装片段的五个要素不同，会使确定出的界限不相同，为确保确定出的界限是病毒界限，能够用于准确鉴定分类出病毒，在本发明的一个实施例中，要求所述第一类组装片段的数目大于30。在本发明的一个实施例中，在确定界限后，对确定的界限的分类效果进行ROC曲线评估，要求其AUC值不小于0.7。ROC曲线(receiver operating characteristic curve，接收者操作特征曲线)，是一种二元分类模型，即输出结果只有两种类别的模型。考虑一个二分问题，即将实例分成正类(positive)或负类(negative)，对一个二分问题来说，会出现四种情况：如果一个实例是正类并且也被预测成正类，即为真正类(True positive，TP)，如果实例是负类被预测成正类，称之为假正类(False positive，FP)，相应地，如果实例是负类被预测成负类，称之为真负类(True negative，TN)，正类被预测成负类则为假负类(false negative，FN)。TP：正确肯定的数目；FN：漏报，没有正确找到的匹配的数目；FP：误报，给出的匹配是不正确的；TN：正确拒绝的非匹配对数。在一个二分类模型中，对于所得到的连续结果，这边的连续结果指五个要素的至少之一或者其全部或部分的组合对组装序列来源于病毒还是非病毒的分类结果，假设已确定(c)中的突变位点比例的界限，比如为阈值0.5，大于这个值的组装序列划归为病毒(正类)，小于这个值则划到宿主(负类)。如果减小阀值，减到0.3，固然能识别出更多的来源于病毒的序列，也就是提高了识别出的正类占所有正类的比例，即TPR(true positive rate，真正类率)，但同时也将更多的负类当作了正类，即提高了FPR(false positive rate，假正类率)。引入ROC能够形象化这一变化，ROC曲线可以用于评价一个分类器，即评价这一界限。AUC(Area Under roc Curve)为ROC曲线下方的面积，AUC值介于0.5到1.0之间，AUC越大，分类器分类效果越好。

依据本发明的另一方面，本发明提供一种鉴定病毒的装置，包括：数据输入单元，用于输入数据；数据输出单元，用于输出数据；存储单元，用于存储数据，其中包括计算机可执行程序；处理器，与所述数据输入单元、数据输出单元和存储单元连接，用于执行所述程序，执行所述程序包括完成上述鉴定病毒的方法。

依据本发明的再一方面，本发明提供一种计算机可读存储介质，用于存储供计算机执行的程序，本领域普通技术人员可以理解，在执行该程序时，通过指令相关硬件可完成上述鉴定病毒的方法的全部或部分步骤。所称存储介质可以包括：只读存储器、随机存储器、磁盘或光盘等。

本发明的方法和/或装置，基于小RNA(small RNAs)，包括病毒介导的小RNA(virus-derived small interfering RNAs，vsiRNAs)，基于病毒和宿主之间的作用机制，实现非基于同源序列比对来鉴定病毒，利用机器学习和模拟的办法确定未知序列五个要素中至少一个要素的病毒界限，适用于任何新一代高通量测序产出的数据，能有效提高病毒鉴定的检出率和准确率。

以下结合具体实施例对本发明病毒鉴定方法和/或装置进行详细的描述。除另有交待，以下实施例中涉及的未特别交待的试剂、序列(接头、标签和引物)、软件及仪器，都是常规市售产品或者开源的，例如购买Illumina的转录组文库构建试剂盒。

实施例一

图2示意病毒鉴定过程以及界限的确定过程，主要包括：

1.样本选取

该实施例选取了16个植物样本，样品名称为：Cooks_footf、Grass100f、Poplar100f、TCV_add、TCV、TCV-TYMV_add、TCV-TYMV、TGM-CK、TYMV-2、TYMV、Willow100f、GSM548932、GSM548933、peach_flower、peach_fruit、peach_leaf。其中，TGM-CK是实验室纯培养的无菌幼苗，TCV、TCV-TYMV、TYMV是分别在纯培养的幼苗基础上人工主动感染病毒，TCV(Turnip crinkle virus)为芜菁皱缩病毒，TYMV(Turnip yellow mosaic virus)为芜菁黄花叶病毒，样本TCV感染TCV病毒，样本TCV-TYMV感染TCV和TYMV两种病毒，样本TYMV感染TYMV病毒，样本TYMV-2为实验重复，样本TCV_add和TCV-TYMV_add为测序数据重复(技术重复)，其他样本为野外采集样本，其中样本Grass100f为多种草的混合物。样本选择包括了实验室的和野外的，有病毒感染的和纯净的，涉及到了实验重复和技术重复。

2.样本处理

将上述样本组织分别加液氮、然后碾磨，后用Trizol(Invitrogen)提取总RNA。然后将总RNA用15％的琼脂糖凝胶电泳，在maker为50nt左右位置，切约1g琼脂糖凝胶，做纯化回收。接着在5’和3’加测序接头(adaptor)，最后RT-PCR反转成cDNA，HiSeq 2000测序(Liang,C.,X.Zhang,et al.(2010)."Identification of miRNA from Porphyra yezoensis byhigh-throughput sequencing and bioinformatics analysis."PLoS One 5(5):e10698.)。

3.数据预处理

测序产生的序列(读段，reads)可能会有接头污染以及质量值低的情况，因此，去掉有接头污染的序列，去掉测序接头自连的序列以及确保每条序列的平均单碱基错误率在0.01以下。最后，丢弃那些长度在18nt以下以及44nt以上的序列，同时，丢弃每条序列有大于两个碱基的质量值在10以下的序列。过滤完数据之后，所有序列通过BLAST 2.2.23来和非编码RNA数据库比对，包括Rfam(rRNA，tRNA，snRNA，snoRNA，http://www.sanger.ac.uk/software/Rfam)和Genbank的非编码RNA比对(http://www.ncbi.nlm.nih.gov/)，比对条件过滤：允许两个碱基错配，然后，所有序列拿来和microRNA数据库比对，参考数据版本是miRBase release 18，比对条件过滤：无错配。

4.读段组装和比对

将过滤完低质量和非编码RNA的序列数据分成4份，分别是25％的数据量，50％的数据量，75％的数据量和100％的数据量，依次使用软件SOAPdenovo-Trans v1.0进行基于Kmer15和Kmer17的混合组装(http://soap.genomics.org.cn/SOAPdenovo-Trans.html)(Li,R.,H.Zhu,et al.(2010)."De novo assembly of human genomes with massively parallel short readsequencing."Genome Res 20(2):265-272.)。将4份的组装结果混合，使用软件PHRAP进行下一步组装，原理及软件参数含义参考(de la Bastide,M.and W.R.McCombie(2007)."Assembling genomic DNA sequences with PHRAP."Curr Protoc Bioinformatics Chapter 11:Unit1114.)，使用的参数为50-100overlap match,vector bound 30,max gap 5。将组装好的长序列(contigs)进行两步操作，一步是使用BLASTn和BLASTx比对到NCBI的nt库和nr库，进行物种注释。Nt库的比对在同源性85％以上以及序列80％以上能比上库判定为比对阳性，将比对阴性的序列进行nr库比对，将同源性50％以上以及序列50％以上比上库判定为阳性，将比对结果分成三组：A组，nt阳性；B组，nt阴性但nr阳性；C组，nt和nr都阴性。另一步是将过滤完的原始短序列即读段用软件BOWTIE v0.12.7(Langmead,B.,C.Trapnell,et al.(2009)."Ultrafast and memory-efficient alignment of short DNA sequences to thehuman genome."Genome Biol 10(3):R25.)比对到组装好的contigs序列上，允许2个错配。

5.突变位点计算

将原始短序列和contig比对结果使用软件Samtools(Li,H.,B.Handsaker,et al.(2009)."The Sequence Alignment/Map format and SAMtools."Bioinformatics 25(16):2078-2079.)进行排序，生成一个pipeup文件，基于这个文件计算每条contig的SNV，并且计算每个SNV的熵值S：S＝-100*Sum(Pi*logPi)，Pi为每个SNV位点的各碱基深度数值(http://www.fludb.org/brcDocs/documents/IRD_FluPolymorphism.pdf)。每个突变位点的覆盖深度小于100则不参与熵值和突变率的计算。该位点错配碱基的数量除以总的覆盖深度，取百分数后则为该位点的突变率。突变位点的比例计算为：该条contig的序列的突变位点数除以该条contig所有不小于100覆盖深度的位点，取百分数。最后，每条序列都会生成5个属性值，即前述的5个要素：平均熵值，中位数熵值，平均突变率，中位数突变率，突变位点比例。

6.机器学习方法构建

SVM(Support Vector Machine)方法是比较流行的用来数据分类的方法，他的基本法则是，给定数据制定为+1或者-1，然后根据给定数据之外的数据和给定数据的相似度，判断归为+1还是-1。根据第4步的比对注释结果，将每个样本的contigs标记为病毒(+1)和宿主(-1)，使用LibSVM软件包(Chih-Chung Chang,C.-J.L.(2011)."LIBSVM:a library forsupport vector machines."ACM Transactions on Intelligent Systems and Technology2:27:1--27:27.)，使用SVC(support vector classification)算法来区分病毒和宿主序列；为优化分类，在建立训练集的时候，罚分参数C和RBF(Radical Basis Function)内核参数γ的设定将根据5倍交叉来确定，即将用来做训练集的数据分5等分，第一次取第一等份的训练集的数据来预测剩下4份的数据，第二次取第二等份的数据预测其他4份数据，依次进行5次循环，确定出分类界限。

7.模型评价

为了评估模型的可靠性，进行了如下几个参数的评价：i)灵敏度(Sensitivity)：正确预测为病毒序列的百分比；ii)特异性(Specificity)：正确预测为宿主序列的数量占预测为宿主序列的总数的比例；iii)准确率(Accuracy)：正确预测病毒和宿主序列的比例；iv)二值相关系数(MCC，Matthews correlation coefficient)是评估预测质量的值(Matthews,B.W.(1975)."Comparison of the predicted and observed secondary structure of T4phage lysozyme."Biochim Biophys Acta 405(2):442-451.)。MCC系数越靠近1，预测结果越完美，越靠近0，预测结果越随机。如果一个实例是正类并且也被预测成正类，即为真正类(True positive，TP)，如果实例是负类被预测成正类，称之为假正类(False positive，FP)，相应地，如果实例是负类被预测成负类，称之为真负类(True negative，TN)，正类被预测成负类则为假负类(false negative，FN)，这4个评估参数可以通过如下公式算得：

$\mathrm{Sensitivity}\left(\mathrm{Sn}\right)=\frac{\mathrm{TP}}{(\mathrm{TP}+\mathrm{FN})}\×100$

$\mathrm{Specificity}\left(\mathrm{Sp}\right)=\frac{\mathrm{TN}}{(\mathrm{TN}+\mathrm{FP})}\×100$

$\mathrm{Accuracy}\left(\mathrm{Acc}\right)=\frac{(\mathrm{TP}+\mathrm{TN})}{(\mathrm{TP}+\mathrm{TN}+\mathrm{FP}+\mathrm{FN})}\×100$

$\mathrm{MCC}=\frac{(\mathrm{TP}\×\mathrm{TN})-(\mathrm{FP}\×\mathrm{FN})}{\sqrt{(\mathrm{TP}+\mathrm{FP})(\mathrm{TP}+\mathrm{FN})(\mathrm{TN}+\mathrm{FP})(\mathrm{TN}+\mathrm{FN})}}$

同时，为了评估预测模型可靠性，ROC(Receiver Operating Characteristic)曲线和AUC(Area Under Curve)曲线将用来评估模型结果。我们使用软件包LibSVM里的plotroc.py脚本进行评估。较佳的，要求预测模型(即确定出的界限)的预测准确率达到90％，在保证准确度的前体下，尽量提高灵敏度；和/或要求AUC值不小于0.7。

参照图2，总结以上，界限的确定和病毒的鉴定包括：1)用已知病毒感染实验室培育的无菌幼苗，待出现病毒感染特征后，取有感染特征的叶子组织少许，用液氮冷却碾磨，2)提取总RNA。用2％凝胶电泳分离50nt以下的片段，回收凝胶，3)得到small RNA样本，4)通过高通量测序，生成大量的小RNA序列，将小RNA序列通过软件组装，5)得到50nt以上的长序列(contigs)。6)将长序列和已知数据库比对，得到每条长序列属于宿主和病毒的信息。7)将测序出来的短序列重新比对到组装好的长序列上，8)在比对结果中得到5个因素的信息，形成每条长序列对应的5个因素的信息。9)用每条长序列的这5个因素的信息和属于宿主还是病毒的信息，进行训练集的构建，训练集完成自身评估后备用。10)同样方法处理和测序出来的未知病毒感染样本或者野生样本，得到50nt以上的长序列，同样也生成5因素的信息，11)再用训练集预测该样本，得到每条序列是否是病毒或者宿主，完成对未知样本的预测。

实施例二

病毒鉴定相关因素考量与确定，包括测序数据量、组装方式、五个要素判断标准、预测模型的建议和评估。

1.测序数据量确定

为尽可能能找到样本中的病毒，我们期望测更多的数据量，但是测序量增加成本也随之上升，况且，如果测序量已经能覆盖样品中的全部序列，测更多的数据也是一种浪费，因此，需要大致确定small RNA的测序数据量。

为了评估这个量，我们测了100M的序列，分别取不同数据量中没有冗余的序列的数量比上非冗余序列数量的比值做为评估测序数据饱和度的参考。

图3显示测序饱和度评估图，图中a线表示非冗余序列条数比上测序总条数的比值(uniq/total)，如左纵坐标轴所示，b线表示没有冗余的序列的条数比上非冗余序列的总数的比值(Expr1/Uniq)，如右纵坐标轴所示。随着测序数据的增加，uniq/total的比值在减小，说明越来越多的序列被重复测到。但同时，从1M的序列到16M序列，没有冗余的序列比上非冗余序列的比值在增加，说明测序量越大，不断有新的序列被测到，并且新的序列被测到的趋势要大于已测到的序列被重复测到的概率。从图3的结果来看，每个样品的测序数据大于8M较佳。

2.组装方式

由于small RNA的高通量测序产生的数据庞大，但是都长度都小于50nt，并且混杂了不同类型的small RNA，如：miRNA，siRNA，vsiRNA等，还有属于宿主的各种RNA片段。因此，对于这样的混合数据进行组装，尤其是基于Kmer方式的组装，单一软件都有一定的差别。我们在此基础上，为了提高组装的精确性，我们尝试了不同数据量和不同数据类型的组装效果。

以下表1显示基于不同整体分割的序列(读段)的组装结果。其中，A:Expr1，表示过滤完数据后，所有冗余数为1的序列。B:Filter或者Filtered，表示过滤完数据后，去除已知非编码RNA的数据。C:Full，为过滤完原始数据后剩下的质量好的全部数据。D:No.expr1过滤完数据后去除冗余数为1的短序列。Cut：表示将数据累计等分的数量，如4＝1/4+2/4+3/4+4/4的数据混合组装。

表1

从表1可看出，不同组合的组装结果中，取不同的数据量对组装的序列条数影响不大，但是不同类型的数据的组装结果序列条数是有差别的。类型D在取不同数据量的情况下，差异最小，虽然在数据量少的时候最长序列长度是最长的，但是数据增加基本没带来影响。由于知道的病毒序列基因组在7k多，因此选择最大长度在7k左右的。太长的，可能有引入组装错误，太短的组装不完整。

图4为不同方式组装在最大contig长度水平的评估结果，其中，Expr1表示过滤完数据后，所有冗余数为1的序列。Filtered或Filter表示过滤完数据后，去除已知非编码RNA的数据。Full，表示过滤完原始数据后剩下的质量好的全部数据。No.expr1表示过滤完数据后去除冗余数为1的短序列。M-reads：将数据累计等分的数量，如4＝1/4+2/4+3/4+4/4的数据混合组装。

从图4得知，高通量测序数据的分类组装法能有效提高针对small RNA数据的组装效果。full样品组装出的最长contig增长趋势。no.expr1样品在缺乏expr1序列的支持下，其组装出的最长contig落后于full样品在。expr1的数据对于支持contig连接成更长的序列起至关重要的作用。对于filtered样品，即过滤非编码RNA的样品来说，分成4份以后最长contig基本上没有增加。说明：初级组装的序列有可能通过一些非编码RNA连接。基于以上结果，较佳的，可选取分切4份、全部数据组装作为后续组装的标准。

3.病毒的5因素判断标准

用原始的短序列(读段，reads)map(比对)到组装好的contig上，每条序列map最多允许2个错配。最后每条contig都有5因素属性。针对单个对每条序列进行注释，并根据5因素进行距离聚类，对聚类的方式不作限定，图5显示的聚类图，是将5因素越相近的contigs归越近的结果。

图5显示出，根据每条contig的5因素特征，病毒和非病毒序列能聚在一起。黑色框内是病毒序列(V开头标记)，其他序列为宿主(非病毒，H开头标记)序列。

由图5可以看到，V开头的标记序列在一个分支上，H开头的序列在其他分支上，其他分支根据具体情况会有不同分支，这也证明了病毒序列有自己特性，而非病毒序列可能特性并不单一，但至少和病毒序列是有差别的。

4.预测模型的建立和评估

从实施例一中提到的16个样本的测序数据中选取了11个样本(Cooks_footf、Grass100f、Poplar100f、TCV_add、TCV、TCV-TYMV_add、TCV-TYMV、TGM-CK、TYMV-2、TYMV、Willow100f)构建了一个训练集(命名为training11)，为了评估训练集的准确性，对训练集进行了交叉验证。画了ROC曲线(Receiver Operating Characteristic)和AUC曲线(Area UnderCurve)，如图6。对于训练集来说，AUC值到了0.8956，说明准确性比较高。图6是训练集评估的ROC曲线(Receiver Operating Characteristic)图。图6中的小表格是对训练集进行的分割验证。将总的训练集数据分成5份，分别独立进行验证，表格里列出的是每次的准确率和准确的条数。

同时，为了评估训练集的可靠性，一些指标如Sensitivity(灵敏度，Sn)、Specificity(特异性，Sp)、Accuracy(准确率，Acc)、Matthews correlation coefficient(二值相关系数，Mcc)和Precision(精确度)用来计算评估，其中，上面未提及的Precision＝TP/(TP+FP)，反映了被分类器判定的正例中真正的正例样本的比重。训练集重头分(等分)或者随机分成100份，如重头分割为：1％、2％、…、100％，然后分别预测自己本身，一百份的每一份都预测后计算那些指标值。Sensitivity是预测为病毒的比例，也可以理解为对于病毒的召回能力，但是这个值的增加，往往会带来Precision的降低。从训练集的稳定性的评估结果来看，包括数据等分的评估和随机数据量的评估结果，训练集training11预测本身，重头分割(等分)的所有指标的指都稳定在80％-100％，而随机分割的稳定范围变大，在30％-100％。此外，precision的值在重头分割的部分要比在随机分割的部分高，而sensitivity值却相反。但是，不管这些值如何，在接近100％的地方这些值都是较高的，因此较佳的，可用100％的数据做训练集。

实施例三

对未知序列是否来自病毒的预测和验证。

1)低质量数据过滤。将测序得到的fq文件进行去接头和低质量，生成fasta格式文件，并统计冗余数，即，每条序列的ID包含该条序列的重复数，如：t00001200。其中t00001是ID，200是该条序列的重复数。

2)已知ncRNA过滤。将初步过滤完成的fasta文件和已知的ncRNA数据库比对，包括：http://www.sanger.ac.uk/software/Rfam，GenBank的非编码RNA数据库(noncoding RNAdatabase)(http://www.ncbi.nlm.nih.gov/)，使用BLAST 2.2.23软件，要求e值<0.01。再和miRBase(release 18)数据库比对，要求完全一样比对上。得到进一步过滤的fasta格式文件。

3)数据组装。如果该物种已经有基因组序列，那么将生成好的fasta文件使用BOWTIEv0.12.7软件与基因组进行比对，要求完全匹配。滤除能比上基因组的序列，将剩下的序列分成4份，分别为总序列条数的25％，50％，75％和100％，4份数据分别进行组装。组装软件使用SOAPdenovo-Trans v1.0(http://soap.genomics.org.cn/SOAPdenovo-Trans.html)，使用Kmer 15和17，将两个Kmer组装结果混合后去除冗余，定为组装1。4份数据分别生成组装1，组装2，组装3，组装4，将4份混合再去冗余，去除50bp以下长度的序列，最后使用软件PHRAP或者minimo进行最后的组装，选择50–100个碱基长度作为overlap，生成最后的组装结果，contig的fasta文件。

4)Contig和sRNA比对。将第二步生成的fasta文件和第三步生成的最终Contig文件使用BOWTIE v0.12.7软件比对，要求2个mismatch。比对完成后，根据比对结果，统计每条contig的每个位点的碱基组成，计算每个位点的熵值(参照实施例一描述)，突变率(单个位点除了最大频率碱基的比率外的比率)，整条contig最终得到突变位点率(突变位点的个数除以总的位点个数的百分比)、平均熵值、中位数熵值、平均突变率、中位数突变率。

5)病毒和宿主预测。得到contig的5因素信息后，使用R语言的LibSVM软件包里的SVC(support vector classification)工具，以我们构建好的训练集为模版，对contig进行分类。结果中标记为1的是病毒，标记为-1的为非病毒。具体的，对样本中组装得到的未知的序列，即，nt和nr比对都是阴性的contig，进行5因素预测处理。根据训练集的5因素标准，判断未知contig的5因素综合结果是否符合病毒特征，若符合，判断为病毒，反之，判断为非病毒。若未知序列的突变位点率，或者突变位点率、以及平均熵值和中位数熵值至少之一，或者突变位点率、以及平均突变率和中位数突变率的至少之一，或者突变位点率、平均熵值和中位数熵值至少之一、以及平均突变率和中位数突变率的至少之一，或者全部5个因素都落在其对应的标准因素的±30％内，则判定为符合标准，判定该序列来自病毒。表2显示利用建立好的训练集对未知序列的非同源预测结果。

表2

为了验证未知序列(contig)是否真的是病毒的，根据样本Grass100f被预测的50条序列设计特异的引物，验证并延伸得到的contig。根据预测的结果，样本Grass100f中Contig9、Contig10、Contig11三条contig被预测为病毒，并且都是在nt和nr库里比对阴性的。分别对三条序列设计引物，RT-PCR验证。图7显示验证结果，如图7中的A部分，显示congtig9-11都有明显条带。同时，对congtig9-10、congtig10-11、congtig9-11分别做扩展延伸验证，最后验证contig10-11和contig9-11为阳性，如图中的A和B部分所示，验证三条序列的关系为congtig9-11-10，如图7中的C部分所示，黑色加深的为3730验证区域。将三条序列验证、连接后，能鉴定为Bell pepper endornavirus(甜椒内源RNA病毒)，最初组装得到的同源性很低，不到50％。因此鉴定了病毒在样本中的存在，验证了预测出来的无法通过同源比对鉴定的病毒。

以上所述仅为本发明的较佳实施例，应当理解，这些实施例仅用以解释本发明，并不用于限定本发明。对于本领域的一般技术人员，依据本发明的思想，可以对上述具体实施方式进行变化。

Claims (10)

1.一种病毒鉴定方法，其特征在于，包括，

获取待测样本的RNA测序数据，所述测序数据包括多个读段；

对所述测序数据的第一部分进行组装，获得组装序列，所述测序数据的第一部分包括不能比对上ncRNA参考序列的读段；

将所述测序数据的第一部分与所述组装序列进行比对，获得比对结果；

基于所述比对结果，确定所述组装序列上的突变位点，所述突变位点包括SNV，以及确定各条组装序列的(a)-(c)中的至少之一，

(a)平均熵值和中位数熵值至少之一，以及突变位点比例，

(b)平均突变率和中位数突变率至少之一，以及突变位点比例，

(c)突变位点比例，

一条组装序列的平均熵值为其上的SNV的熵值的平均值，

一条组装序列的中位数熵值为其上的SNV的熵值的中位数，

一条组装序列的平均突变率为其上的SNV的突变率的平均值，

一条组装序列的中位数突变率为其上的SNV的突变率的中位数，

一条组装序列的突变位点比例为其上的SNV的数目所占的比例，

SNV的熵值＝-100*∑(Pi*logPi)，Pi为该SNV的各种碱基的深度，

SNV的突变率＝支持该SNV的读段数目/比对上该SNV的读段数目，

将所述确定的组装序列的(a)-(c)至少之一和与其对应的界限比较，判定落入所述界限的组装序列来自病毒。

2.权利要求1的方法，其特征在于，所述测序数据的第一部分不包括符合(1)-(3)至少之一的读段，

(1)包含接头序列，和/或平均单碱基错误率大于0.01，

(2)长度小于18nt，和/或长度大于44nt，

(3)包含碱基质量值小于10的碱基个数大于2。

3.权利要求1的方法，其特征在于，所述ncRNA参考序列包括，rRNA参考序列、tRNA参考序列、snRNA参考序列、snoRNA参考序列和microRNA参考序列。

4.权利要求1的方法，其特征在于，所述组装包括，

分别对所述测序数据的第一部分中的全部读段和所述测序数据的第一部分中的至少一部分读段进行第一组装和第二组装，获得第一组装结果和第二组装结果，

合并所述第一组装结果和第二组装结果；

任选的，

所述第一组装和/或所述第二组装为混合组装。

5.权利要求1的方法，其特征在于，在获得组装序列之后，去除所述组装序列中的已知序列。

6.权利要求1的方法，其特征在于，在确定每条组装序列的(a)-(c)中的至少之一之前，去除该组装序列中深度小于100X的位点。

7.权利要求1的方法，其特征在于，还包括，

设计引物，对判定来自病毒的组装序列进行RT-PCR延伸，获得延伸产物；

利用延伸产物验证所述组装序列来自病毒。

8.权利要求1-7任一方法，其特征在于，所述界限的确定，包括，

获取至少一个已知病毒感染的样本的RNA测序结果，所述测序结果包括多个读段；

对所述测序结果的至少一部分进行组装，获得组装片段；

将所述组装片段与参考序列进行第一比对，进行组装片段物种注释，获得物种注释结果；

基于所述物种注释结果对所述组装片段进行分类，获得第一类组装片段和第二类组装片段，所述第一类组装序列来自病毒，所述第二类组装片段来自宿主；

分别将所述测序结果的至少一部分与所述第一类组装片段和所述第二类组装片段进行第二比对，获得第二比对结果；

基于所述第二比对结果，确定所述第一类组装片段和所述第二类组装片段上的突变位点，所述突变位点包括SNV，以及确定每条第一类组装片段和每条第二类组装片段的五个因素，所述五个因素为平均熵值、中位数熵值、平均突变率、中位数突变率以及突变位点比例，其中，

一条组装片段的平均熵值为其上的SNV的熵值的平均值，

一条组装片段的中位数熵值为其上的SNV的熵值的中位数，

一条组装片段的平均突变率为其上的SNV的突变率的平均值，

一条组装片段的中位数突变率为其上的SNV的突变率的中位数，

一条组装片段的突变位点比例为其上的SNV的数目所占的比例，

SNV的熵值＝-100*∑(Pi*logPi)，Pi为该SNV的各种碱基的深度，

SNV的突变率＝支持该SNV的读段数目/比对上该SNV的读段数目；

基于所述第一类组装片段和第二类组装片段、以及每条所述第一类组装片段和每条所述第二类组装片段的五个因素，利用SVM，确定所述五个因素中至少之一的界限；

任选的，

所述第一类组装片段的数目大于30；

任选的，

所述界限的ROC曲线评估的AUC值不小于0.7；

任选的，

所述界限的病毒预测的准确率达90％。

9.权利要求8的方法，其特征在于，所述基于第一类组装片段和第二类组装片段、以及每条第一类组装片段和每条第二类组装片段的五个因素，利用SVM，确定所述五个因素中至少之一的界限，包括，

将所述第一类组装片段和第二类组装片段分为多份片段包，每份所述片段包包含多条第一类组装片段和多条第二类组装片段，

利用其中一份片段包中的组装片段的五个因素的至少之一，判断其它各份片段包中的每条组装片段为第一组装片段还是第二组装片段，依据判断的正确率，调整确定出所述五个因素的至少之一的界限。

10.一种鉴定病毒的装置，其特征在于，包括，

数据输入单元，用于输入数据；

数据输出单元，用于输出数据；

存储单元，用于存储数据，其中包括计算机可执行程序；

处理器，与所述数据输入单元、数据输出单元和存储单元连接，用于执行所述程序，执行所述程序包括完成权利要求1-9任一方法。