鑒定病毒的方法和裝置的制造方法
【專利摘要】本發(fā)明公開了一種鑒定病毒的方法，該方法包括：獲取待測(cè)樣本的RNA測(cè)序數(shù)據(jù)；對(duì)測(cè)序數(shù)據(jù)的第一部分進(jìn)行組裝，獲得組裝序列；將測(cè)序數(shù)據(jù)的第一部分與組裝序列進(jìn)行比對(duì)，獲得比對(duì)結(jié)果；基于比對(duì)結(jié)果，確定組裝序列上的突變位點(diǎn)，以及確定各組裝序列的(a)?(c)中的至少之一：(a)平均熵值和中位數(shù)熵值至少之一，以及突變位點(diǎn)比例，(b)平均突變率和中位數(shù)突變率至少之一，以及突變位點(diǎn)比例，(c)突變位點(diǎn)比例；將組裝序列的(a)?(c)至少之一和與其對(duì)應(yīng)的界限比較，判定落入界限的組裝序列來自病毒。本發(fā)明還公開一種鑒定病毒的裝置。利用本發(fā)明的方法和/或裝置鑒定病毒，能夠不依賴同源序列比對(duì)來準(zhǔn)確預(yù)測(cè)未知序列是否為病毒序列。
【專利說明】
鑒定病毒的方法和裝置
技術(shù)領(lǐng)域
[0001] 本發(fā)明涉及生物檢測(cè)領(lǐng)域，具體的，本發(fā)明涉及一種病毒鑒定方法和一種病毒鑒 定裝置。
【背景技術(shù)】
[0002] 截止2014年6月30日，國(guó)際病毒學(xué)分類大會(huì)（ICTV)公布的病毒種類有2827 種。而2011年是有2484種，2009年是有2285種，每年差不多100種病毒被發(fā)現(xiàn)，這個(gè)速 度遠(yuǎn)比不上細(xì)菌等其他微生物。這是因?yàn)椴《颈旧聿荒苌L(zhǎng)，需要寄生在宿主細(xì)胞里，很 難分離培養(yǎng)，而有研究者預(yù)測(cè)，只要有細(xì)胞就會(huì)有病毒，病毒在自然界存在是非常廣泛和巨 大的，現(xiàn)在已知的病毒，還不及總量的萬分之一。隨著測(cè)序技術(shù)的發(fā)展，成本不斷降低，通 量不斷提高，使得越來越多的物種被測(cè)序，越來越多的樣本被發(fā)現(xiàn)存在病毒，而病毒種類的 數(shù)據(jù)也在不斷增加。測(cè)序技術(shù)在病毒發(fā)現(xiàn)的應(yīng)用越來越常見（Barzon，L.，E. Lavezzo，et al. (2011) · ''Applications of next-generation sequencing technologies to diagnostic virology. 〃Int J Mol Sci 12(11):7861_7884·)，并且新一代的測(cè)序技術(shù)在對(duì) 于樣本中一些低豐度的病毒亞型有著獨(dú)特的先天優(yōu)勢(shì)。對(duì)于發(fā)現(xiàn)新病毒，很典型的例子就 是在一個(gè)不明原因發(fā)熱最后死亡的病人身上，傳統(tǒng)檢測(cè)并沒有發(fā)現(xiàn)致病病原，在病人死后 取樣，利用新一代的測(cè)序技術(shù)產(chǎn)生了巨大的數(shù)據(jù)，之后發(fā)現(xiàn)了一株新沙粒病毒才是導(dǎo)致病 人發(fā)病的原因（Palacios，G.，J.Druce，et al. (2008). "A new arenavirus in a cluster of fatal transplant-associated diseases. "N Engl J Med 358(10):991_998·)。同 樣的技術(shù)運(yùn)用，在非洲發(fā)現(xiàn)了一株新的致命病毒Lujo virus (Briese, T.，J. T. Paweska, et al. (2009) ·''Genetic detection and characterization of Lujo virus, a new hemorrhagic fever-associated arenavirus from southern Africa. ^PLoS Pathog 5(5) :e100 0455.)。除了在新病毒的發(fā)現(xiàn)外，對(duì)于已知的病毒也是有著明顯的優(yōu)勢(shì)，比如 2014年在非洲大爆發(fā)的埃博拉疫情。研究人員能在短時(shí)間內(nèi)獲得99株最新的埃博拉病 毒，并且迅速分析了大量的變異位點(diǎn)，最終確定了傳播的可能路徑和變異（Gire，S. K.，A. Goba,et al· (2014) · "Genomic surveillance elucidates Ebola virus origin and transmission during the 2014outbreak.''Science 345 (6202) : 1369-1372.) 〇
[0003] 新一代測(cè)序技術(shù)在病毒發(fā)現(xiàn)和分析上，相比傳統(tǒng)的技術(shù)有了很大通量和陽(yáng)性率的 提高，但是很大程度上依賴于現(xiàn)有病毒序列，要么和已有病毒序列比較相近，發(fā)現(xiàn)了有較大 的變異、重組，要么和已有的病毒序列相差較大，但是仍然在同一個(gè)屬水平。目前所用的方 法基本都是基于同源序列的比對(duì)，對(duì)已經(jīng)構(gòu)建好的數(shù)據(jù)庫(kù)的依賴較大。
[0004] 同樣基于新一代測(cè)序技術(shù)的基礎(chǔ)上，利用小RNA (smal 1 RNA)的數(shù)據(jù)，吳 等人通過對(duì)小RNA的拼接，從延長(zhǎng)的序列里挑出來了病毒的序列（Wu，Q.，Y.Luo，et al· (2010) ·''Virus discovery by deep sequencing and assembly of virus-derived small silencing RNAs."Proc Natl Acad Sci U S A 107(4): 1606-1611·)，拼接完后，基 于同源比對(duì)來鑒定。
[0005] 前面提到的small RNA在細(xì)胞生物里是廣泛存在的，包括microRNAs (miRNAs)，小 干擾（RNAsmall interfering RNAs(siRNAs)和 Piwi-interacting RNAs(piRNAs)，或多或 少，直接或者間接的參與了細(xì)胞基因表達(dá)的調(diào)節(jié)，保護(hù)細(xì)胞抵御外來入侵（Ghildiyal，Μ. and P. D. Zamore(2009) · "Small silencing RNAs: an expanding universe. ^Nat Rev Genet 10(2) :94-108.)。在真菌、植物和無脊椎動(dòng)物中，病毒介導(dǎo)的小RNA (virus-derived small interfering RNAs，vsiRNAs)也是自身免疫系統(tǒng)的重要組成部分（￥&11]\^61'1〇，]\ T. , K. ff. 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結(jié)合后又能影響轉(zhuǎn)錄和表達(dá)Mlotshwa, S.，G. J. Pruss, et al. (2008) ·''Small RNAs in viral infection and host defense.''Trends Plant Sci 13(7):375-382.)〇

【發(fā)明內(nèi)容】

[0006] 本發(fā)明旨在至少在一定程度上解決上述技術(shù)問題之一或至少提供一種商業(yè)選擇。 為此，本發(fā)明的目的在于提出鑒定病毒的手段。
[0007] 依據(jù)本發(fā)明的一方面，本發(fā)明提供一種鑒定病毒的方法，該方法包括：獲取待測(cè)樣 本的RNA測(cè)序數(shù)據(jù)，所述測(cè)序數(shù)據(jù)包括多個(gè)讀段；對(duì)所述測(cè)序數(shù)據(jù)的第一部分進(jìn)行組裝，獲 得組裝序列，所述測(cè)序數(shù)據(jù)的第一部分包括不能比對(duì)上ncRNA參考序列的讀段；將所述測(cè) 序數(shù)據(jù)的第一部分與所述組裝序列進(jìn)行比對(duì)，獲得比對(duì)結(jié)果；基于所述比對(duì)結(jié)果，確定所述 組裝序列上的突變位點(diǎn)，所述突變位點(diǎn)包括SNV，以及確定各條組裝序列的（a)-(c)中的至 少之一，（a)平均熵值和中位數(shù)熵值至少之一，以及突變位點(diǎn)比例，（b)平均突變率和中位 數(shù)突變率至少之一，以及突變位點(diǎn)比例，（c)突變位點(diǎn)比例，一條組裝序列的平均熵值為其 上的SNV的熵值的平均值，一條組裝序列的中位數(shù)熵值為其上的SNV的熵值的中位數(shù)，一條 組裝序列的平均突變率為其上的SNV的突變率的平均值，一條組裝序列的中位數(shù)突變率為 其上的SNV的突變率的中位數(shù)，一條組裝序列的突變位點(diǎn)比例為其上的SNV的數(shù)目所占的 比例，SNV的熵值=-100* Σ (Pi*logPi)，Pi為該SNV的各種堿基的深度，SNV的突變率=支 持該SNV的讀段數(shù)目/比對(duì)上該SNV的讀段數(shù)目；將所述確定的組裝序列的（a) - (c)至少之 一和與其對(duì)應(yīng)的界限比較，判定落入所述界限的組裝序列來自病毒。所稱的測(cè)序數(shù)據(jù)通過 測(cè)序獲得，測(cè)序方法依據(jù)測(cè)序平臺(tái)的不同可選擇但不限于Illumina公司的Hisq2000/2500 測(cè)序平臺(tái)、Life Technologies公司的Ion Torrent平臺(tái)和單分子測(cè)序平臺(tái)，測(cè)序方式可 以選擇單端測(cè)序，也可以選擇雙末端測(cè)序，獲得的下機(jī)數(shù)據(jù)是測(cè)讀出來的片段，稱為讀段 (reads)。ncRNA (non-coding RNA)為非編碼RNA，所稱的ncRNA參考序列包括：rRNA參考序 列、tRNA參考序列、snRNA參考序列、snoRNA參考序列和microRNA參考序列中的至少之一， 所稱的參考序列指預(yù)先確定的序列，可以是預(yù)先獲得的待測(cè)樣本所屬或者所包含的生物類 別的任意參考模板，例如，若待測(cè)樣本來源的個(gè)體為人類，參考序列可選擇NCBI數(shù)據(jù)庫(kù)提 供的HG19,若待測(cè)樣本來源的個(gè)體宿主為人類，目標(biāo)核酸是病毒RNA中的ncRNA，參考序列 可選擇Rfam數(shù)據(jù)庫(kù)和/或GeneBank中的人非編碼RNA序列，利于將宿主的已知ncRNA和 剩余的ncRNA區(qū)分開，進(jìn)一步地，也可以預(yù)先配置包含更多參考序列的資源庫(kù)，例如依據(jù)待 測(cè)樣本來源個(gè)體的狀態(tài)、地域等因素選擇或是測(cè)定組裝出更接近的序列作為參考序列。SNV 為單核苷酸變異，同SNP，SNV-般由兩種堿基組成，是一種二等位基因。所稱的深度，即測(cè) 序深度，指堿基被測(cè)序或者讀取的平均次數(shù)。比對(duì)可以利用已知比對(duì)軟件進(jìn)行，例如SOAP、 BWA和TeraMap等，在比對(duì)過程中，一般對(duì)比對(duì)參數(shù)進(jìn)行設(shè)置，設(shè)置一條reads最多允許有 s個(gè)堿基錯(cuò)配（mismatch)，例如設(shè)置s < 4,若reads中有超過s個(gè)堿基發(fā)生錯(cuò)配，則視為 該reads無法比對(duì)到（比對(duì)上）參考序列，而所稱的支持某個(gè)SNV的讀段指比對(duì)上該位點(diǎn) 的讀段的相應(yīng)位置的堿基與該位點(diǎn)的一致。所稱的界限是屬于病毒的界限，界限可能是個(gè) 臨界值、數(shù)值范圍或關(guān)系式，但能夠界定出病毒。與（a)-(c)中至少之一對(duì)應(yīng)的界限是利用 感染病毒的樣本的測(cè)序數(shù)據(jù)來確定的，與（a)對(duì)應(yīng)的界限包括平均熵值界限和中位數(shù)熵值 界限至少之一、以及突變位點(diǎn)比例界限，與（b)對(duì)應(yīng)的界限包括平均突變率界限和中位數(shù) 突變率界限至少之一、以及突變位點(diǎn)比例界限，與（c)對(duì)應(yīng)的界限包括突變位點(diǎn)比例界限。 (a)-(c)共包含組裝序列的非重復(fù)的五個(gè)要素，五個(gè)要素對(duì)應(yīng)的界限都可根據(jù)已知感染有 病毒的樣本的核酸測(cè)序數(shù)據(jù)來確定。
[0008] 本發(fā)明的這一方面的病毒鑒定方法，實(shí)現(xiàn)非基于同源序列比對(duì)來鑒定病毒，且適 用于任何新一代高通量測(cè)序產(chǎn)出的數(shù)據(jù)，能有效提高病毒鑒定的檢出率和準(zhǔn)確率。
[0009] 根據(jù)本發(fā)明的一個(gè)實(shí)施例，該方法中的測(cè)序數(shù)據(jù)的第一部分不包括符合以下 (1)-(3)至少之一的讀段：（1)包含接頭序列，和/或平均單堿基錯(cuò)誤率大于0.01，（2)長(zhǎng)度 小于18nt，和/或長(zhǎng)度大于44nt，（3)包含堿基質(zhì)量值小于10的堿基個(gè)數(shù)大于2。過濾掉 被測(cè)序接頭污染的和低質(zhì)量的讀段，讀段的可靠性提高利于準(zhǔn)確鑒定病毒。
[0010] 根據(jù)本發(fā)明的一個(gè)實(shí)施例，該方法中的組裝包括：分別對(duì)所述測(cè)序數(shù)據(jù)的第一部 分中的全部讀段和所述測(cè)序數(shù)據(jù)的第一部分中的至少一部分讀段進(jìn)行第一組裝和第二組 裝，獲得第一組裝結(jié)果和第二組裝結(jié)果；合并所述第一組裝結(jié)果和第二組裝結(jié)果。任選的， 所述第一組裝和/或所述第二組裝為混合組裝。所稱的混合組裝指采用多種Kmer分別組 裝，再合并不同Kmer的組裝結(jié)果。Kmer指長(zhǎng)度為K的一段序列，一條長(zhǎng)度為L(zhǎng)的reads產(chǎn) 生的Kmer數(shù)量為L(zhǎng)-K+1，組裝時(shí)，將每個(gè)讀段分解成其包含的所有長(zhǎng)度為K的固定序列，利 于快速組裝。基于不同數(shù)據(jù)量和采用不同Kmer分別進(jìn)行組裝，再合并各個(gè)組裝結(jié)果，能夠 提尚組裝的精確性。
[0011] 根據(jù)本發(fā)明的一個(gè)實(shí)施例，在獲得組裝序列之后，去除所述組裝序列中的已知序 列。例如，通過與公開數(shù)據(jù)庫(kù)中收錄的物種信息確定的序列進(jìn)行比對(duì)，將比對(duì)上的組裝序列 去除，判定是否比對(duì)上可以依據(jù)所使用的比對(duì)軟件或程序的默認(rèn)設(shè)置或者設(shè)置軟件或程序 的參數(shù)數(shù)值定義。組裝序列中的已知序列由于公開數(shù)據(jù)中已經(jīng)披露了其來源，可直接鑒定 出其來源，無需進(jìn)行后續(xù)步驟來鑒定，去除這些序列減少了數(shù)據(jù)量，利于利用機(jī)器快速運(yùn)行 該病毒鑒定方法。
[0012] 根據(jù)本發(fā)明的一個(gè)實(shí)施例，在確定每條組裝序列的（a)-(c)中的至少之一之前， 去除該組裝序列中深度小于100X的位點(diǎn)。所稱的位點(diǎn)為組成該條組裝序列的各個(gè)位點(diǎn)，包 括SNP位點(diǎn)，也包括非SNP位點(diǎn)，將組裝序列上的深度小于100X的位點(diǎn)排除在任一（a) - (C) 的計(jì)算確定過程，有利于準(zhǔn)確進(jìn)行病毒鑒定。
[0013] 根據(jù)本發(fā)明的一個(gè)實(shí)施例，該方法還包括：設(shè)計(jì)引物，對(duì)判定來自病毒的組裝序列 進(jìn)行RT-PCR延伸，獲得延伸產(chǎn)物；利用延伸產(chǎn)物驗(yàn)證所述組裝序列來自病毒。利用判定來 自病毒的組裝序列進(jìn)行引物設(shè)計(jì)，RT-PCR對(duì)這條或這些組裝序列進(jìn)行延伸和/或連接，以 驗(yàn)證判定的組裝序列的確來自病毒。
[0014] 根據(jù)本發(fā)明的一些實(shí)施例，確定該方法中的界限，包括：獲取至少一個(gè)已知病毒 感染的樣本的RNA測(cè)序結(jié)果，所述測(cè)序結(jié)果包括多個(gè)讀段；對(duì)所述測(cè)序結(jié)果的至少一部分 進(jìn)行組裝，獲得組裝片段；將所述組裝片段與參考序列進(jìn)行第一比對(duì)，進(jìn)行組裝片段物種 注釋，獲得物種注釋結(jié)果；基于所述物種注釋結(jié)果對(duì)所述組裝片段進(jìn)行分類，獲得第一類 組裝片段和第二類組裝片段，所述第一類組裝序列來自病毒，所述第二類組裝片段來自宿 主；分別將所述測(cè)序結(jié)果的至少一部分與所述第一類組裝片段和所述第二類組裝片段進(jìn) 行第二比對(duì)，獲得第二比對(duì)結(jié)果；基于所述第二比對(duì)結(jié)果，確定所述第一類組裝片段和所 述第二類組裝片段上的突變位點(diǎn)，所述突變位點(diǎn)包括SNV，以及確定每條第一類組裝片段 和每條第二類組裝片段的五個(gè)因素，所述五個(gè)因素為平均熵值、中位數(shù)熵值、平均突變率、 中位數(shù)突變率以及突變位點(diǎn)比例，其中，一條組裝片段的平均熵值為其上的SNV的熵值的 平均值，一條組裝片段的中位數(shù)熵值為其上的SNV的熵值的中位數(shù)，一條組裝片段的平均 突變率為其上的SNV的突變率的平均值，一條組裝片段的中位數(shù)突變率為其上的SNV的突 變率的中位數(shù)，一條組裝片段的突變位點(diǎn)比例為其上的SNV的數(shù)目所占的比例，SNV的熵 值=-100* Σ (Pi*logPi)，Pi為該SNV的各種堿基的深度，SNV的突變率=支持該SNV的 讀段數(shù)目/比對(duì)上該SNV的讀段數(shù)目；基于所述第一類組裝片段和第二類組裝片段、以及 每條所述第一類組裝片段和每條所述第二類組裝片段的五個(gè)因素，利用SVM，確定所述五 個(gè)因素中至少之一的界限。支持向量機(jī)（Support Vector Machine，SVM)是比較流行的 用來數(shù)據(jù)分類的方法，它的基本法則可看成，將給定數(shù)據(jù)定為+1或者-1，然后根據(jù)給定數(shù) 據(jù)之外的數(shù)據(jù)和給定數(shù)據(jù)的相似程度，判斷給定數(shù)據(jù)之外的數(shù)據(jù)應(yīng)歸為+1還是-1，例如， 根據(jù)物種注釋結(jié)果，將第一類組裝片段標(biāo)記為病毒（+1)，將第二類組裝片段標(biāo)記為宿主 (-1)，使用 LibSVM 軟件包（Chih-Chung Chang, C.-J.L. (2011). "LIBSVM:a library for support vector machines. 〃ACM Transactions on Intelligent Systems and Technology 2:27:1-27:27.)，使用 SVC (support vector classification)算法來找出規(guī)律即五個(gè)要 素的界限來區(qū)分病毒和宿主序列。在本發(fā)明的一個(gè)實(shí)施例中，要求所述第一類組裝片段的 數(shù)目大于30,利于最后確定出的界限能用于準(zhǔn)確分類鑒定。在本發(fā)明的一個(gè)實(shí)施例中，在確 定界限后，對(duì)確定的界限的分類可靠性進(jìn)行R0C曲線評(píng)估，要求其AUC值不小于0. 7。在本 發(fā)明的一個(gè)實(shí)施例中，要求確定出的界限對(duì)病毒預(yù)測(cè)的準(zhǔn)確率達(dá)90 %。
[0015] 根據(jù)本發(fā)明的一個(gè)實(shí)施例，所述基于第一類組裝片段和第二類組裝片段、以及每 條第一類組裝片段和每條第二類組裝片段的五個(gè)因素，利用SVM，確定所述五個(gè)因素中至少 之一的界限，包括：將所述第一類組裝片段和第二類組裝片段分為多份片段包，每份所述片 段包包含多條第一組裝片段和多條第二組裝片段，利用其中一份片段包中的組裝片段的五 個(gè)因素的至少之一，判斷其它各份片段包中的每條組裝片段為第一組裝片段還是第二組裝 片段，依據(jù)判斷的正確率，確定所述五個(gè)因素的至少之一的界限。
[0016] 依據(jù)本發(fā)明的另一方面，本發(fā)明提供一種鑒定病毒的裝置，包括：數(shù)據(jù)輸入單元， 用于輸入數(shù)據(jù)；數(shù)據(jù)輸出單元，用于輸出數(shù)據(jù)；存儲(chǔ)單元，用于存儲(chǔ)數(shù)據(jù)，其中包括計(jì)算機(jī) 可執(zhí)行程序；處理器，與所述數(shù)據(jù)輸入單元、數(shù)據(jù)輸出單元和存儲(chǔ)單元連接，用于執(zhí)行所述 程序，執(zhí)行所述程序包括完成上述鑒定病毒的方法。
[0017] 依據(jù)本發(fā)明的再一方面，本發(fā)明提供一種計(jì)算機(jī)可讀存儲(chǔ)介質(zhì)，用于存儲(chǔ)供計(jì)算 機(jī)執(zhí)行的程序，本領(lǐng)域普通技術(shù)人員可以理解，在執(zhí)行該程序時(shí)，通過指令相關(guān)硬件可完成 上述鑒定病毒的方法的全部或部分步驟。所稱存儲(chǔ)介質(zhì)可以包括：只讀存儲(chǔ)器、隨機(jī)存儲(chǔ) 器、磁盤或光盤等。
[0018] 本發(fā)明的方法和/或裝置，基于小RNA (small RNAs)，包括病毒介導(dǎo)的小 RNA(virus-derived small interfering RNAs，vsiRNAs)，基于病毒和宿主之間的作用機(jī) 制，利用機(jī)器學(xué)習(xí)和模擬的辦法確定未知序列五個(gè)要素中至少一個(gè)要素的病毒界限，實(shí)現(xiàn) 非基于同源序列比對(duì)來鑒定病毒，且適用于任何新一代高通量測(cè)序產(chǎn)出的數(shù)據(jù)，能有效提 高病毒鑒定的檢出率和準(zhǔn)確率。
【附圖說明】
[0019] 本發(fā)明的上述和/或附加的方面和優(yōu)點(diǎn)從結(jié)合下面附圖對(duì)實(shí)施方式的描述中將 變得明顯和容易理解，其中：
[0020] 圖1是本發(fā)明的一個(gè)實(shí)施例中的病毒鑒定方法包含的步驟的示意圖。
[0021] 圖2是本發(fā)明的一個(gè)實(shí)施例中的病毒鑒定方法的流程圖。
[0022] 圖3是本發(fā)明的一個(gè)實(shí)施例中的測(cè)序飽和度評(píng)估結(jié)果示意圖。
[0023] 圖4是本發(fā)明的一個(gè)實(shí)施例中的不同組裝方式對(duì)最大contig影響的評(píng)估結(jié)果示 意圖。
[0024] 圖5是本發(fā)明的一個(gè)實(shí)施例中的組裝序列的聚類示意圖。
[0025] 圖6是本發(fā)明的一個(gè)實(shí)施例中的訓(xùn)練集評(píng)估的R0C曲線。
[0026] 圖7是本發(fā)明的一個(gè)實(shí)施例中的病毒序列試驗(yàn)驗(yàn)證示意圖。
【具體實(shí)施方式】
[0027] 下面詳細(xì)描述本發(fā)明的實(shí)施例，所述實(shí)施例的示例在附圖中示出，其中，自始至終 相同或類似的標(biāo)號(hào)表示相同或類似的元件或具有相同或類似功能的元件。下面通過參考附 圖描述的實(shí)施例是示例性的，僅用于解釋本發(fā)明，而不能理解為對(duì)本發(fā)明的限制。需要說明 的，本文中所使用的術(shù)語(yǔ)"第一"、"第二"或者"第一部分"等僅為方便描述，不能理解為指 示或暗示相對(duì)重要性，也不能理解為之間有先后順序關(guān)系。在本發(fā)明的描述中，除非另有說 明，"多個(gè)"的含義是兩個(gè)或兩個(gè)以上。在本文中，除非另有明確的規(guī)定和限定，術(shù)語(yǔ)"相連"、 "連接"等術(shù)語(yǔ)應(yīng)做廣義理解，例如，可以是固定連接，也可以是可拆卸連接，或一體地連接； 可以是機(jī)械連接，也可以是電連接；可以是直接相連，也可以通過中間媒介間接相連，可以 是兩個(gè)元件內(nèi)部的連通。
[0028] 如圖1所示，依據(jù)本發(fā)明的一方面，本發(fā)明提供一種鑒定病毒的方法，該方法包 括：S10獲取待測(cè)樣本的RNA測(cè)序數(shù)據(jù)，所述測(cè)序數(shù)據(jù)包括多個(gè)讀段；S20對(duì)所述測(cè)序數(shù)據(jù) 的第一部分進(jìn)行組裝，獲得組裝序列，所述測(cè)序數(shù)據(jù)的第一部分包括不能比對(duì)上ncRNA參 考序列的讀段；S30將所述測(cè)序數(shù)據(jù)的第一部分與所述組裝序列進(jìn)行比對(duì)，獲得比對(duì)結(jié)果； S40基于所述比對(duì)結(jié)果，確定所述組裝序列上的突變位點(diǎn)，所述突變位點(diǎn)包括SNV，以及確 定各條組裝序列的（a)-(c)中的至少之一 ：(a)平均熵值和中位數(shù)熵值至少之一，以及突 變位點(diǎn)比例，（b)平均突變率和中位數(shù)突變率至少之一，以及突變位點(diǎn)比例，（c)突變位點(diǎn) 比例，一條組裝序列的平均熵值為其上的SNV的熵值的平均值，一條組裝序列的中位數(shù)熵 值為其上的SNV的熵值的中位數(shù)，一條組裝序列的平均突變率為其上的SNV的突變率的平 均值，一條組裝序列的中位數(shù)突變率為其上的SNV的突變率的中位數(shù)，一條組裝序列的突 變位點(diǎn)比例為其上的SNV的數(shù)目所占的比例，SNV的熵值=-100* Σ (Pi*logPi)，Pi為該 SNV的各種堿基的深度，SNV的突變率=支持該SNV的讀段數(shù)目/比對(duì)上該SNV的讀段數(shù) 目；S50將所述確定的組裝序列的（a)-(c)至少之一和與其對(duì)應(yīng)的界限比較，判定落入所述 界限的組裝序列來自病毒。其中，所稱的測(cè)序數(shù)據(jù)通過測(cè)序獲得，測(cè)序方法依據(jù)測(cè)序平臺(tái)的 不同可選擇但不限于IIlumina公司的Hisq2000/2500測(cè)序平臺(tái)、Life Technologies公司 的Ion Torrent平臺(tái)和單分子測(cè)序平臺(tái)，測(cè)序方式可以選擇單端測(cè)序，也可以選擇雙末端 測(cè)序，獲得的下機(jī)數(shù)據(jù)是測(cè)讀出來的片段，稱為讀段（reads)。ncRNA(non-coding RNA)為 非編碼RNA，所稱的ncRNA參考序列包括：rRNA參考序列、tRNA參考序列、snRNA參考序列、 snoRNA參考序列和microRNA參考序列中的至少之一，所稱的參考序列指預(yù)先確定的序列， 可以是預(yù)先獲得的待測(cè)樣本所屬或者所包含的生物類別的任意參考模板，例如，若待測(cè)樣 本來源的個(gè)體為人類，參考序列可選擇NCBI數(shù)據(jù)庫(kù)提供的HG19,若待測(cè)樣本來源的個(gè)體宿 主為人類，目標(biāo)核酸是病毒RNA中的ncRNA，參考序列可選擇Rfam數(shù)據(jù)庫(kù)和/或GeneBank 中的人非編碼RNA序列，利于將宿主的已知ncRNA和剩余的ncRNA區(qū)分開，進(jìn)一步地，也可 以預(yù)先配置包含更多參考序列的資源庫(kù)，例如依據(jù)待測(cè)樣本來源個(gè)體的狀態(tài)、地域等因素 選擇或是測(cè)定組裝出更接近的序列作為參考序列。SNV為單核苷酸變異，同SNP，SNV-般 由兩種堿基組成，是一種二等位基因。所稱的深度，即測(cè)序深度，指堿基被測(cè)序或者讀取的 平均次數(shù)。比對(duì)可以利用已知比對(duì)軟件進(jìn)行，例如S0AP、BWA和TeraMap等，在比對(duì)過程中， 一般對(duì)比對(duì)參數(shù)進(jìn)行設(shè)置，設(shè)置一條reads最多允許有s個(gè)堿基錯(cuò)配（mismatch)，例如設(shè) 置s彡4,若reads中有超過s個(gè)堿基發(fā)生錯(cuò)配，則視為該reads無法比對(duì)到（比對(duì)上）參 考序列，而所稱的支持某個(gè)SNV的讀段指比對(duì)上該位點(diǎn)的讀段的相應(yīng)位置的堿基與該位點(diǎn) 的一致。所稱的界限是屬于病毒的界限，界限可能是個(gè)臨界值、數(shù)值范圍或關(guān)系式，但能夠 界定出病毒。與（a)-(c)中至少之一對(duì)應(yīng)的界限是利用感染病毒的樣本的測(cè)序數(shù)據(jù)來確 定的，與（a)對(duì)應(yīng)的界限包括平均熵值界限和中位數(shù)熵值界限至少之一、以及突變位點(diǎn)比 例界限，與（b)對(duì)應(yīng)的界限包括平均突變率界限和中位數(shù)突變率界限至少之一、以及突變 位點(diǎn)比例界限，與（c)對(duì)應(yīng)的界限包括突變位點(diǎn)比例界限。（a)-(c)共包含組裝序列的非 重復(fù)的五個(gè)要素，五個(gè)要素至少之一或者其全部或者部分組合能反映出病毒序列自己的特 性，雖然非病毒序列可能特性并不單一，但至少和病毒序列是有差別的，序列的五個(gè)要素的 至少之一或者其全部或部分的組合能夠用以判斷該序列是否來自病毒，五個(gè)要素對(duì)應(yīng)的界 限都可根據(jù)已知感染有病毒的樣本的核酸測(cè)序數(shù)據(jù)來確定。在本發(fā)明的一個(gè)實(shí)施例中，確 定了組裝序列的（a)-(c)中的至少之二或者全部五個(gè)要素。本發(fā)明的這一方面的病毒鑒定 方法，實(shí)現(xiàn)非基于同源序列比對(duì)來鑒定病毒，且適用于任何新一代高通量測(cè)序產(chǎn)出的數(shù)據(jù)， 能有效提高病毒鑒定的檢出率和準(zhǔn)確率。
[0029] 根據(jù)本發(fā)明的一個(gè)實(shí)施例，該方法中的測(cè)序數(shù)據(jù)的第一部分不包括符合以下 (1)-(3)至少之一的讀段：（1)包含接頭序列，和/或平均單堿基錯(cuò)誤率大于0.01，（2)長(zhǎng)度 小于18nt，和/或長(zhǎng)度大于44nt，（3)包含堿基質(zhì)量值小于10的堿基個(gè)數(shù)大于2。包含接 頭序列指該reads包含測(cè)序接頭，測(cè)讀的不是目標(biāo)區(qū)域，為被測(cè)序接頭污染的讀段。堿基錯(cuò) 誤率和堿基質(zhì)量值為測(cè)序平臺(tái)賦予讀段的值，質(zhì)量值為-l〇*lg(p)，這里，p為測(cè)錯(cuò)的概率， 即單堿基錯(cuò)誤率，當(dāng)一條reads某位置堿基的出錯(cuò)概率為0. 1時(shí)，其質(zhì)量值為10。過濾掉被 測(cè)序接頭污染的和低質(zhì)量的讀段，能使測(cè)序數(shù)據(jù)包含的reads整體質(zhì)量提高，利于減少后 續(xù)分析鑒定對(duì)機(jī)器內(nèi)存的需求，而且讀段的可靠性提高利于準(zhǔn)確鑒定病毒。
[0030] 由于小RNA的高通量測(cè)序產(chǎn)生的數(shù)據(jù)龐大，長(zhǎng)度都小于50nt，并且混雜了不同類 型的小RNA，如：miRNA，siRNA，vsiRNA等，還有屬于宿主的各種RNA片段。對(duì)于這樣的混 合數(shù)據(jù)進(jìn)行組裝，尤其是基于Kmer方式的組裝，單一軟件都有一定的差別。為提高組裝的 精確性，根據(jù)本發(fā)明的一個(gè)實(shí)施例，該方法中的組裝包括：分別對(duì)所述測(cè)序數(shù)據(jù)的第一部分 中的全部讀段和所述測(cè)序數(shù)據(jù)的第一部分中的至少一部分讀段進(jìn)行第一組裝和第二組裝， 獲得第一組裝結(jié)果和第二組裝結(jié)果；合并所述第一組裝結(jié)果和第二組裝結(jié)果。任選的，所 述第一組裝和/或所述第二組裝為混合組裝。組裝可以利用已知序列組裝方法，例如利用 soapdenovo、velvet等。所稱的混合組裝指多個(gè)組裝混合，如采用多種Kmer分別組裝，再合 并不同Kmer的組裝結(jié)果，和/或采用多個(gè)數(shù)據(jù)量分別進(jìn)行組裝。Kmer指長(zhǎng)度為K的一段序 列，一條長(zhǎng)度為L(zhǎng)的reads產(chǎn)生的Kmer數(shù)量為L(zhǎng)-K+1，組裝時(shí)，將每個(gè)讀段分解成其包含的 所有長(zhǎng)度為K的固定序列，利于快速組裝。在本發(fā)明的一個(gè)實(shí)施例中，基于所述測(cè)序數(shù)據(jù)的 第一部分出4份不同等分?jǐn)?shù)據(jù)量分別進(jìn)行組裝，如取25%的測(cè)序數(shù)據(jù)的第一部分?jǐn)?shù)據(jù)量， 50%的測(cè)序數(shù)據(jù)的第一部分的數(shù)據(jù)量，75%的測(cè)序數(shù)據(jù)的第一部分的數(shù)據(jù)量和100%的測(cè) 序數(shù)據(jù)得第一部分?jǐn)?shù)據(jù)量分別進(jìn)行組裝，以分辨不同等分對(duì)組裝出的contigs的影響?；?于不同數(shù)據(jù)量和采用不同Kmer分別進(jìn)行組裝，獲得不同的組裝結(jié)果，能夠提高組裝效果， 提尚組裝的精確性。
[0031] 根據(jù)本發(fā)明的一個(gè)實(shí)施例，在獲得組裝序列之后，去除所述組裝序列中的已知序 列。例如，通過與公開數(shù)據(jù)庫(kù)中收錄的物種信息確定的序列進(jìn)行比對(duì)，將比對(duì)上的組裝序列 去除，判定是否為比對(duì)上可以依據(jù)所使用的比對(duì)軟件或程序的默認(rèn)設(shè)置或者設(shè)置軟件或程 序的參數(shù)數(shù)值定義。組裝序列中的已知序列由于公開數(shù)據(jù)中已經(jīng)披露了其來源，可直接鑒 定出其來源，可以不用進(jìn)行后續(xù)步驟進(jìn)行鑒定，去除這些序列減少了數(shù)據(jù)量，利于利用機(jī)器 快速運(yùn)行該病毒鑒定方法。在本發(fā)明的一個(gè)實(shí)施例中，使用BLASTn和BLASTx將組裝序列 與NCBI的nt庫(kù)和nr庫(kù)分別進(jìn)行比對(duì)，進(jìn)行物種注釋。與Nt庫(kù)的比對(duì)在同源性85%以上 以及序列80%以上能比對(duì)上庫(kù)判定為比對(duì)陽(yáng)性，將比對(duì)陰性的序列進(jìn)行nr庫(kù)比對(duì)，將同源 性50%以上以及序列50%以上比上庫(kù)判定為陽(yáng)性，這樣，比對(duì)結(jié)果可分成三組：比對(duì)上nt 庫(kù)（nt陽(yáng)性）；沒有比對(duì)上nt庫(kù)但比對(duì)上nr庫(kù)（nt陰性但nr陽(yáng)性）；nt和nr都陰性。nt 和nr都陰性的組裝序列為完全未知序列，基于目前根據(jù)同源比對(duì)方法無法確定其來源。
[0032] 根據(jù)本發(fā)明的一個(gè)實(shí)施例，在確定每條組裝序列的（a)-(c)中的至少之一之前， 去除該組裝序列中深度小于100X的位點(diǎn)。所稱的位點(diǎn)為組成該條組裝序列的各個(gè)位點(diǎn)，包 括SNP位點(diǎn)，也包括非SNP位點(diǎn)，將組裝序列上的深度小于100X的位點(diǎn)排除在任一（a) - (c) 的計(jì)算確定過程，有利于準(zhǔn)確進(jìn)行病毒鑒定。
[0033] 根據(jù)本發(fā)明的一個(gè)實(shí)施例，該方法還包括：設(shè)計(jì)引物，對(duì)判定來自病毒的組裝序列 進(jìn)行RT-PCR延伸，獲得延伸產(chǎn)物；利用延伸產(chǎn)物驗(yàn)證所述組裝序列來自病毒。利用判定來 自病毒的組裝序列進(jìn)行引物設(shè)計(jì)，RT-PCR對(duì)這條或這些組裝序列進(jìn)行延伸，任選的對(duì)延伸 產(chǎn)物進(jìn)行順序確定和連接，以驗(yàn)證判定來自病毒的組裝序列來自病毒。
[0034] 據(jù)本發(fā)明的一些實(shí)施例，確定該方法中的界限，包括：獲取至少一個(gè)已知病毒感 染的樣本的RNA測(cè)序結(jié)果，所述測(cè)序結(jié)果包括多個(gè)讀段；對(duì)所述測(cè)序結(jié)果的至少一部分進(jìn) 行組裝，獲得組裝片段；將所述組裝片段與參考序列進(jìn)行第一比對(duì)，進(jìn)行組裝片段物種注 釋，獲得物種注釋結(jié)果；基于所述物種注釋結(jié)果對(duì)所述組裝片段進(jìn)行分類，獲得第一類組 裝片段和第二類組裝片段，所述第一類組裝序列來自病毒，所述第二類組裝片段來自宿主； 分別將所述測(cè)序結(jié)果的至少一部分與所述第一類組裝片段和所述第二類組裝片段進(jìn)行第 二比對(duì)，獲得第二比對(duì)結(jié)果；基于所述第二比對(duì)結(jié)果，確定所述第一類組裝片段和所述第 二類組裝片段上的突變位點(diǎn)，所述突變位點(diǎn)包括SNV，以及確定每條第一類組裝片段和每 條第二類組裝片段的五個(gè)因素，所述五個(gè)因素為平均熵值、中位數(shù)熵值、平均突變率、中位 數(shù)突變率以及突變位點(diǎn)比例，其中，一條組裝片段的平均熵值為其上的SNV的熵值的平均 值，一條組裝片段的中位數(shù)熵值為其上的SNV的熵值的中位數(shù)，一條組裝片段的平均突變 率為其上的SNV的突變率的平均值，一條組裝片段的中位數(shù)突變率為其上的SNV的突變 率的中位數(shù)，一條組裝片段的突變位點(diǎn)比例為其上的SNV的數(shù)目所占的比例，SNV的熵值 =-100* Σ (Pi*logPi)，Pi為該SNV的各種堿基的深度，SNV的突變率=支持該SNV的 讀段數(shù)目/比對(duì)上該SNV的讀段數(shù)目；基于所述第一類組裝片段和第二類組裝片段、以及 每條所述第一類組裝片段和每條所述第二類組裝片段的五個(gè)因素，利用SVM，確定所述五 個(gè)因素中至少之一的界限。支持向量機(jī)（Support Vector Machine，SVM)是比較流行的 用來數(shù)據(jù)分類的方法，它的基本法則可看成，將給定數(shù)據(jù)定為+1或者-1，然后根據(jù)給定數(shù) 據(jù)之外的數(shù)據(jù)和給定數(shù)據(jù)的相似程度，判斷給定數(shù)據(jù)之外的數(shù)據(jù)應(yīng)歸為+1還是-1。一 般的，可以隨機(jī)選出一些數(shù)據(jù)作為模型（訓(xùn)練集），發(fā)現(xiàn)其規(guī)律，利用其它數(shù)據(jù)作為測(cè)試 集來計(jì)算檢測(cè)這個(gè)模型的規(guī)律與訓(xùn)練集的誤差等，從而確定發(fā)現(xiàn)的規(guī)律是否正確。在本 發(fā)明的一個(gè)實(shí)施例，所述基于第一類組裝片段和第二類組裝片段、以及每條第一類組裝片 段和每條第二類組裝片段的五個(gè)因素，利用SVM，確定所述五個(gè)因素中至少之一的界限， 包括：將所述第一類組裝片段和第二類組裝片段分為多份片段包，每份所述片段包包含多 條第一組裝片段和多條第二組裝片段，利用其中一份片段包中的組裝片段的五個(gè)因素的 至少之一，判斷其它各份片段包中的每條組裝片段為第一組裝片段還是第二組裝片段，依 據(jù)判斷的正確率，調(diào)整確定出所述五個(gè)因素的至少之一的界限。具體地，根據(jù)物種注釋結(jié) 果，將每個(gè)樣本的contig標(biāo)記為屬于病毒（+1)和屬于宿主（-1)，使用LibSVM軟件包實(shí) 現(xiàn)基于 SVM 的分類和回歸（Chih-Chung Chang, C.-J.L. (2011). "LIBSVM:a library for support vector machines. 〃ACM Transactions on Intelligent Systems and Technology 2:27:1-27:27.)，使用 SVC (support vector classification)算法來區(qū)分病毒和宿主序 列，為優(yōu)化分類，在建立訓(xùn)練集的時(shí)候，罰分參數(shù)C和RBF(Radical Basis Function)內(nèi)核 參數(shù)γ的設(shè)定將根據(jù)5倍交叉驗(yàn)證來確定，即將用來做訓(xùn)練集的數(shù)據(jù)分5等份，第一次取 第一等份的訓(xùn)練集的數(shù)據(jù)來預(yù)測(cè)剩下4份的數(shù)據(jù)，第二次取第二等份的數(shù)據(jù)預(yù)測(cè)其他4份 數(shù)據(jù)，依次進(jìn)行5次循環(huán)，確定出分類界限。
[0035] 在不同實(shí)施例中，由于構(gòu)建的訓(xùn)練集的不同，例如選取的已知病毒感染的樣本不 同、樣本數(shù)目的不同，使得獲得的組裝片段、物種注釋結(jié)果及組裝片段的五個(gè)要素不同，會(huì) 使確定出的界限不相同，為確保確定出的界限是病毒界限，能夠用于準(zhǔn)確鑒定分類出病毒， 在本發(fā)明的一個(gè)實(shí)施例中，要求所述第一類組裝片段的數(shù)目大于30。在本發(fā)明的一個(gè)實(shí) 施例中，在確定界限后，對(duì)確定的界限的分類效果進(jìn)行R0C曲線評(píng)估，要求其AUC值不小于 0.7。ROC 曲線（receiver operating characteristic curve，接收者操作特征曲線），是 一種二元分類模型，即輸出結(jié)果只有兩種類別的模型。考慮一個(gè)二分問題，即將實(shí)例分成 正類（positive)或負(fù)類（negative)，對(duì)一個(gè)二分問題來說，會(huì)出現(xiàn)四種情況：如果一個(gè)實(shí) 例是正類并且也被預(yù)測(cè)成正類，即為真正類（True positive，TP)，如果實(shí)例是負(fù)類被預(yù)測(cè) 成正類，稱之為假正類（False positive，F(xiàn)P)，相應(yīng)地，如果實(shí)例是負(fù)類被預(yù)測(cè)成負(fù)類，稱 之為真負(fù)類（True negative，TN)，正類被預(yù)測(cè)成負(fù)類則為假負(fù)類（false negative，F(xiàn)N)。 TP:正確肯定的數(shù)目；FN:漏報(bào)，沒有正確找到的匹配的數(shù)目；FP:誤報(bào)，給出的匹配是不正 確的；TN:正確拒絕的非匹配對(duì)數(shù)。在一個(gè)二分類模型中，對(duì)于所得到的連續(xù)結(jié)果，這邊的 連續(xù)結(jié)果指五個(gè)要素的至少之一或者其全部或部分的組合對(duì)組裝序列來源于病毒還是非 病毒的分類結(jié)果，假設(shè)已確定（c)中的突變位點(diǎn)比例的界限，比如為閾值0.5,大于這個(gè)值 的組裝序列劃歸為病毒（正類），小于這個(gè)值則劃到宿主（負(fù)類）。如果減小閥值，減到 0.3,固然能識(shí)別出更多的來源于病毒的序列，也就是提高了識(shí)別出的正類占所有正類的比 例，即TPR(true positive rate，真正類率），但同時(shí)也將更多的負(fù)類當(dāng)作了正類，即提高了 FPR(false positive rate，假正類率）。引入R0C能夠形象化這一變化，R0C曲線可以用于 評(píng)價(jià)一個(gè)分類器，即評(píng)價(jià)這一界限。AUC(Area Under roc Curve)為R0C曲線下方的面積， AUC值介于0. 5到1. 0之間，AUC越大，分類器分類效果越好。
[0036] 依據(jù)本發(fā)明的另一方面，本發(fā)明提供一種鑒定病毒的裝置，包括：數(shù)據(jù)輸入單元， 用于輸入數(shù)據(jù)；數(shù)據(jù)輸出單元，用于輸出數(shù)據(jù)；存儲(chǔ)單元，用于存儲(chǔ)數(shù)據(jù)，其中包括計(jì)算機(jī) 可執(zhí)行程序；處理器，與所述數(shù)據(jù)輸入單元、數(shù)據(jù)輸出單元和存儲(chǔ)單元連接，用于執(zhí)行所述 程序，執(zhí)行所述程序包括完成上述鑒定病毒的方法。
[0037] 依據(jù)本發(fā)明的再一方面，本發(fā)明提供一種計(jì)算機(jī)可讀存儲(chǔ)介質(zhì)，用于存儲(chǔ)供計(jì)算 機(jī)執(zhí)行的程序，本領(lǐng)域普通技術(shù)人員可以理解，在執(zhí)行該程序時(shí)，通過指令相關(guān)硬件可完成 上述鑒定病毒的方法的全部或部分步驟。所稱存儲(chǔ)介質(zhì)可以包括：只讀存儲(chǔ)器、隨機(jī)存儲(chǔ) 器、磁盤或光盤等。
[0038] 本發(fā)明的方法和/或裝置，基于小RNA (smal 1 RNAs)，包括病毒介導(dǎo)的小 RNA(virus-derived small interfering RNAs，vsiRNAs)，基于病毒和宿主之間的作用機(jī) 制，實(shí)現(xiàn)非基于同源序列比對(duì)來鑒定病毒，利用機(jī)器學(xué)習(xí)和模擬的辦法確定未知序列五個(gè) 要素中至少一個(gè)要素的病毒界限，適用于任何新一代高通量測(cè)序產(chǎn)出的數(shù)據(jù)，能有效提高 病毒鑒定的檢出率和準(zhǔn)確率。
[0039] 以下結(jié)合具體實(shí)施例對(duì)本發(fā)明病毒鑒定方法和/或裝置進(jìn)行詳細(xì)的描述。除另有 交待，以下實(shí)施例中涉及的未特別交待的試劑、序列（接頭、標(biāo)簽和引物）、軟件及儀器，都 是常規(guī)市售產(chǎn)品或者開源的，例如購(gòu)買Illumina的轉(zhuǎn)錄組文庫(kù)構(gòu)建試劑盒。
[0040] 實(shí)施例一
[0041] 圖2示意病毒鑒定過程以及界限的確定過程，主要包括：
[0042] 1.樣本選取
[0043] 該實(shí)施例選取了 16個(gè)植物樣本，樣品名稱為：Cooks_footf、GrasslOOf、 PoplarlOOf、TCV_add、TCV、TCV-TYMV_add、TCV-TYMV、TGM-CK、TYMV-2、TYMV、WillowlOOf、 GSM548932、GSM548933、peach_flower、peach_fruit、peach_leaf。其中，TGM-CK 是實(shí)驗(yàn)室 純培養(yǎng)的無菌幼苗，TCV、TCV-TYMV、TYMV是分別在純培養(yǎng)的幼苗基礎(chǔ)上人工主動(dòng)感染病毒， TCV(Turnip crinkle virus)為憲菁皺縮病毒，TYMV(Turnip yellow mosaic virus)為憲 菁黃花葉病毒，樣本TCV感染TCV病毒，樣本TCV-TYMV感染TCV和TYMV兩種病毒，樣本TYMV 感染TYMV病毒，樣本TYMV-2為實(shí)驗(yàn)重復(fù)，樣本TCV_add和TCV-TYMV_add為測(cè)序數(shù)據(jù)重復(fù) (技術(shù)重復(fù)），其他樣本為野外采集樣本，其中樣本GrasslOOf為多種草的混合物。樣本選 擇包括了實(shí)驗(yàn)室的和野外的，有病毒感染的和純凈的，涉及到了實(shí)驗(yàn)重復(fù)和技術(shù)重復(fù)。
[0044] 2.樣本處理
[0045] 將上述樣本組織分別加液氮、然后碾磨，后用Trizol (Invitrogen)提取總RNA。 然后將總RNA用15 %的瓊脂糖凝膠電泳，在maker為50nt左右位置，切約lg瓊脂糖凝 膠，做純化回收。接著在5'和3'加測(cè)序接頭（adaptor)，最后RT-PCR反轉(zhuǎn)成cDNA，HiSeq 2000測(cè)序（Liang, C.，X. Zhang, et al. (2010). "Identification of miRNA from Porphyra yezoensis by high-throughput sequencing and bioinformatics analysis. ^PLoS One 5(5) :el0698.) 〇
[0046] 3.數(shù)據(jù)預(yù)處理
[0047] 測(cè)序產(chǎn)生的序列（讀段，reads)可能會(huì)有接頭污染以及質(zhì)量值低的情況，因此，去 掉有接頭污染的序列，去掉測(cè)序接頭自連的序列以及確保每條序列的平均單堿基錯(cuò)誤率在 0.01以下。最后，丟棄那些長(zhǎng)度在18nt以下以及44nt以上的序列，同時(shí)，丟棄每條序列有 大于兩個(gè)堿基的質(zhì)量值在10以下的序列。過濾完數(shù)據(jù)之后，所有序列通過BLAST 2.2.23 來和非編碼 RNA 數(shù)據(jù)庫(kù)比對(duì)，包括 Rfam(rRNA，tRNA，snRNA，snoRNA，http: //www. sanger. ac. uk/software/Rfam)和 Genbank 的非編碼 RNA 比對(duì)（http://www. ncbi. nlm. nih. gov/)， 比對(duì)條件過濾：允許兩個(gè)堿基錯(cuò)配，然后，所有序列拿來和microRNA數(shù)據(jù)庫(kù)比對(duì)，參考數(shù)據(jù) 版本是miRBase release 18,比對(duì)條件過濾：無錯(cuò)配。
[0048] 4.讀段組裝和比對(duì)
[0049] 將過濾完低質(zhì)量和非編碼RNA的序列數(shù)據(jù)分成4份，分別是25 %的數(shù)據(jù)量，50 %的 數(shù)據(jù)量，75%的數(shù)據(jù)量和100%的數(shù)據(jù)量，依次使用軟件SOAPdenovo-Trans vl. 0進(jìn)行基于 Kmerl5 和 Kmerl7 的混合組裝(http://soap. Renomics. org. cn/SOAPdenovo-Trans. html) (Li, R. , H. Zhu, et al. (2010) · 〃De novo assembly of human genomes with massively parallel short read sequencing. "Genome Res 20(2) :265-272.)。將 4 份的組裝結(jié)果 混合，使用軟件PHRAP進(jìn)行下一步組裝，原理及軟件參數(shù)含義參考（de la Bastide，M. and ff. R. McCombie(2007). ^Assembling genomic DNA sequences with PHRAP. ^Curr Protoc Bioinformatics Chapter 11 :Unitlll4.)，使用的參數(shù)為 50_100overlap match, vector bound 30，max gap 5。將組裝好的長(zhǎng)序列（contigs)進(jìn)行兩步操作，一步是使用BLASTn和 BLASTx比對(duì)到NCBI的nt庫(kù)和nr庫(kù)，進(jìn)行物種注釋。Nt庫(kù)的比對(duì)在同源性85%以上以及 序列80%以上能比上庫(kù)判定為比對(duì)陽(yáng)性，將比對(duì)陰性的序列進(jìn)行nr庫(kù)比對(duì)，將同源性50% 以上以及序列50%以上比上庫(kù)判定為陽(yáng)性，將比對(duì)結(jié)果分成三組：A組，nt陽(yáng)性；B組，nt陰 性但nr陽(yáng)性；C組，nt和nr都陰性。另一步是將過濾完的原始短序列即讀段用軟件B0WTIE v0. 12. 7 (Langmead, B. , C. Trapnell, et al· (2009) ·''Ultrafast and memory-efficient alignment of short DNA sequences to the human genome.''Genome Biol 10 (3) : R25.) 比對(duì)到組裝好的contigs序列上，允許2個(gè)錯(cuò)配。
[0050] 5.突變位點(diǎn)計(jì)算
[0051] 將原始短序列和contig比對(duì)結(jié)果使用軟件Samtools (Li, Η·，B. Handsaker, et al. (2009)·"The Sequence Alignment/Map format and SAMtools. "Bioinformatics 25 (16) : 2078-2079.)進(jìn)行排序，生成一個(gè)pipeup文件，基于這個(gè)文件計(jì)算每條contig的 SNV，并且計(jì)算每個(gè)SNV的熵值S :S = -100*Sum(Pi*logPi)，Pi為每個(gè)SNV位點(diǎn)的各堿基深 度數(shù)值（http://www. fludb. org/brcDocs/documents/IRD_FluPolymorphism. pdf)。每個(gè)突 變位點(diǎn)的覆蓋深度小于100則不參與熵值和突變率的計(jì)算。該位點(diǎn)錯(cuò)配堿基的數(shù)量除以總 的覆蓋深度，取百分?jǐn)?shù)后則為該位點(diǎn)的突變率。突變位點(diǎn)的比例計(jì)算為：該條contig的序 列的突變位點(diǎn)數(shù)除以該條contig所有不小于100覆蓋深度的位點(diǎn)，取百分?jǐn)?shù)。最后，每條 序列都會(huì)生成5個(gè)屬性值，即前述的5個(gè)要素：平均熵值，中位數(shù)熵值，平均突變率，中位數(shù) 突變率，突變位點(diǎn)比例。
[0052] 6.機(jī)器學(xué)習(xí)方法構(gòu)建
[0053] SVM(Support Vector Machine)方法是比較流行的用來數(shù)據(jù)分類的方法，他的基 本法則是，給定數(shù)據(jù)制定為+1或者-1，然后根據(jù)給定數(shù)據(jù)之外的數(shù)據(jù)和給定數(shù)據(jù)的相似 度，判斷歸為+1還是-1。根據(jù)第4步的比對(duì)注釋結(jié)果，將每個(gè)樣本的contigs標(biāo)記為病 毒（+1)和宿主（-1)，使用 LibSVM 軟件包（Chih-Chung Chang，C.-J.L. (2011).〃LIBSVM:a library for support vector machines. 〃ACM Transactions on Intelligent Systems and Technology2:27:1-27:27.)，使用 SVC (support vector classification)算法來區(qū) 分病毒和宿主序列；為優(yōu)化分類，在建立訓(xùn)練集的時(shí)候，罰分參數(shù)C和RBF (Radical Basis Function)內(nèi)核參數(shù)γ的設(shè)定將根據(jù)5倍交叉來確定，即將用來做訓(xùn)練集的數(shù)據(jù)分5等分， 第一次取第一等份的訓(xùn)練集的數(shù)據(jù)來預(yù)測(cè)剩下4份的數(shù)據(jù)，第二次取第二等份的數(shù)據(jù)預(yù)測(cè) 其他4份數(shù)據(jù)，依次進(jìn)行5次循環(huán)，確定出分類界限。
[0054] 7.模型評(píng)價(jià)
[0055] 為了評(píng)估模型的可靠性，進(jìn)行了如下幾個(gè)參數(shù)的評(píng)價(jià)：i)靈敏度（Sensitivity): 正確預(yù)測(cè)為病毒序列的百分比；ii)特異性（Specificity):正確預(yù)測(cè)為宿主序列的數(shù)量 占預(yù)測(cè)為宿主序列的總數(shù)的比例；iii)準(zhǔn)確率（Accuracy):正確預(yù)測(cè)病毒和宿主序列的 比例；iv)二值相關(guān)系數(shù)（MCC，Matthews correlation coefficient)是評(píng)估預(yù)測(cè)質(zhì)量的 值(Matthews, B. ff. (1975). ^Comparison of the predicted and observed secondary structure of T4phage lysozyme. "Biochim Biophys Acta 405(2):442-451. )。MCC 系 數(shù)越靠近1，預(yù)測(cè)結(jié)果越完美，越靠近0,預(yù)測(cè)結(jié)果越隨機(jī)。如果一個(gè)實(shí)例是正類并且也被 預(yù)測(cè)成正類，即為真正類（True positive，TP)，如果實(shí)例是負(fù)類被預(yù)測(cè)成正類，稱之為假 正類（False positive，F(xiàn)P)，相應(yīng)地，如果實(shí)例是負(fù)類被預(yù)測(cè)成負(fù)類，稱之為真負(fù)類（True negative，TN)，正類被預(yù)測(cè)成負(fù)類則為假負(fù)類（false negative，F(xiàn)N)，這4個(gè)評(píng)估參數(shù)可以 通過如下公式算得：
[0056]
[0057]
[0058]
[0059]
[0060] 同時(shí)，為了評(píng)估預(yù)測(cè)模型可靠性，R0C(Receiver Operating Characteristic)曲 線和AUC(Area Under Curve)曲線將用來評(píng)估模型結(jié)果。我們使用軟件包LibSVM里的 plotroc.py腳本進(jìn)行評(píng)估。較佳的，要求預(yù)測(cè)模型（即確定出的界限）的預(yù)測(cè)準(zhǔn)確率達(dá)到 90%，在保證準(zhǔn)確度的前體下，盡量提高靈敏度；和/或要求AUC值不小于0. 7。
[0061] 參照?qǐng)D2,總結(jié)以上，界限的確定和病毒的鑒定包括：1)用已知病毒感染實(shí)驗(yàn)室培 育的無菌幼苗，待出現(xiàn)病毒感染特征后，取有感染特征的葉子組織少許，用液氮冷卻碾磨， 2)提取總RNA。用2%凝膠電泳分離50nt以下的片段，回收凝膠，3)得到small RNA樣本， 4)通過高通量測(cè)序，生成大量的小RNA序列，將小RNA序列通過軟件組裝，5)得到50nt以 上的長(zhǎng)序列（contigs)。6)將長(zhǎng)序列和已知數(shù)據(jù)庫(kù)比對(duì)，得到每條長(zhǎng)序列屬于宿主和病毒 的信息。7)將測(cè)序出來的短序列重新比對(duì)到組裝好的長(zhǎng)序列上，8)在比對(duì)結(jié)果中得到5個(gè) 因素的信息，形成每條長(zhǎng)序列對(duì)應(yīng)的5個(gè)因素的信息。9)用每條長(zhǎng)序列的這5個(gè)因素的信 息和屬于宿主還是病毒的信息，進(jìn)行訓(xùn)練集的構(gòu)建，訓(xùn)練集完成自身評(píng)估后備用。10)同樣 方法處理和測(cè)序出來的未知病毒感染樣本或者野生樣本，得到50nt以上的長(zhǎng)序列，同樣也 生成5因素的信息，11)再用訓(xùn)練集預(yù)測(cè)該樣本，得到每條序列是否是病毒或者宿主，完成 對(duì)未知樣本的預(yù)測(cè)。
[0062] 實(shí)施例二
[0063] 病毒鑒定相關(guān)因素考量與確定，包括測(cè)序數(shù)據(jù)量、組裝方式、五個(gè)要素判斷標(biāo)準(zhǔn)、 預(yù)測(cè)模型的建議和評(píng)估。
[0064] 1.測(cè)序數(shù)據(jù)量確定
[0065] 為盡可能能找到樣本中的病毒，我們期望測(cè)更多的數(shù)據(jù)量，但是測(cè)序量增加成本 也隨之上升，況且，如果測(cè)序量已經(jīng)能覆蓋樣品中的全部序列，測(cè)更多的數(shù)據(jù)也是一種浪 費(fèi)，因此，需要大致確定small RNA的測(cè)序數(shù)據(jù)量。
[0066] 為了評(píng)估這個(gè)量，我們測(cè)了 100M的序列，分別取不同數(shù)據(jù)量中沒有冗余的序列的 數(shù)量比上非冗余序列數(shù)量的比值做為評(píng)估測(cè)序數(shù)據(jù)飽和度的參考。
[0067] 圖3顯示測(cè)序飽和度評(píng)估圖，圖中a線表示非冗余序列條數(shù)比上測(cè)序總條數(shù)的比 值（uniq/total)，如左縱坐標(biāo)軸所示，b線表示沒有冗余的序列的條數(shù)比上非冗余序列的 總數(shù)的比值（Exprl/Uniq)，如右縱坐標(biāo)軸所示。隨著測(cè)序數(shù)據(jù)的增加，uniq/total的比值 在減小，說明越來越多的序列被重復(fù)測(cè)到。但同時(shí)，從1M的序列到16M序列，沒有冗余的序 列比上非冗余序列的比值在增加，說明測(cè)序量越大，不斷有新的序列被測(cè)到，并且新的序列 被測(cè)到的趨勢(shì)要大于已測(cè)到的序列被重復(fù)測(cè)到的概率。從圖3的結(jié)果來看，每個(gè)樣品的測(cè) 序數(shù)據(jù)大于8M較佳。
[0068] 2.組裝方式
[0069] 由于small RNA的高通量測(cè)序產(chǎn)生的數(shù)據(jù)龐大，但是都長(zhǎng)度都小于50nt，并且混 雜了不同類型的small RNA，如：miRNA，siRNA，vsiRNA等，還有屬于宿主的各種RNA片段。 因此，對(duì)于這樣的混合數(shù)據(jù)進(jìn)行組裝，尤其是基于Kmer方式的組裝，單一軟件都有一定的 差別。我們?cè)诖嘶A(chǔ)上，為了提高組裝的精確性，我們嘗試了不同數(shù)據(jù)量和不同數(shù)據(jù)類型的 組裝效果。
[0070] 以下表1顯示基于不同整體分割的序列（讀段）的組裝結(jié)果。其中，A:Exprl，表 示過濾完數(shù)據(jù)后，所有冗余數(shù)為1的序列。B:Filter或者Filtered，表示過濾完數(shù)據(jù)后，去 除已知非編碼RNA的數(shù)據(jù)。C:Full，為過濾完原始數(shù)據(jù)后剩下的質(zhì)量好的全部數(shù)據(jù)。D:No. exprl過濾完數(shù)據(jù)后去除冗余數(shù)為1的短序列。Cut :表示將數(shù)據(jù)累計(jì)等分的數(shù)量，如4 = 1/4+2/4+3/4+4/4的數(shù)據(jù)混合組裝。
[0071] 表 1
[0072]
[0073] 從表1可看出，不同組合的組裝結(jié)果中，取不同的數(shù)據(jù)量對(duì)組裝的序列條數(shù)影響 不大，但是不同類型的數(shù)據(jù)的組裝結(jié)果序列條數(shù)是有差別的。類型D在取不同數(shù)據(jù)量的情 況下，差異最小，雖然在數(shù)據(jù)量少的時(shí)候最長(zhǎng)序列長(zhǎng)度是最長(zhǎng)的，但是數(shù)據(jù)增加基本沒帶來 影響。由于知道的病毒序列基因組在7k多，因此選擇最大長(zhǎng)度在7k左右的。太長(zhǎng)的，可能 有引入組裝錯(cuò)誤，太短的組裝不完整。
[0074] 圖4為不同方式組裝在最大contig長(zhǎng)度水平的評(píng)估結(jié)果，其中，Exprl表示過 濾完數(shù)據(jù)后，所有冗余數(shù)為1的序列。Filtered或Filter表示過濾完數(shù)據(jù)后，去除已知 非編碼RNA的數(shù)據(jù)。Full，表示過濾完原始數(shù)據(jù)后剩下的質(zhì)量好的全部數(shù)據(jù)。No. exprl 表示過濾完數(shù)據(jù)后去除冗余數(shù)為1的短序列。M-reads :將數(shù)據(jù)累計(jì)等分的數(shù)量，如4 = 1/4+2/4+3/4+4/4的數(shù)據(jù)混合組裝。
[0075] 從圖4得知，高通量測(cè)序數(shù)據(jù)的分類組裝法能有效提高針對(duì)small RNA數(shù)據(jù)的組 裝效果。full樣品組裝出的最長(zhǎng)contig增長(zhǎng)趨勢(shì)。no. exprl樣品在缺乏exprl序列的支 持下，其組裝出的最長(zhǎng)contig落后于full樣品在。exprl的數(shù)據(jù)對(duì)于支持contig連接成 更長(zhǎng)的序列起至關(guān)重要的作用。對(duì)于filtered樣品，即過濾非編碼RNA的樣品來說，分成 4份以后最長(zhǎng)contig基本上沒有增加。說明：初級(jí)組裝的序列有可能通過一些非編碼RNA 連接?；谝陨辖Y(jié)果，較佳的，可選取分切4份、全部數(shù)據(jù)組裝作為后續(xù)組裝的標(biāo)準(zhǔn)。
[0076] 3.病毒的5因素判斷標(biāo)準(zhǔn)
[0077] 用原始的短序列（讀段，reads) map (比對(duì)）到組裝好的contig上，每條序列map 最多允許2個(gè)錯(cuò)配。最后每條contig都有5因素屬性。針對(duì)單個(gè)對(duì)每條序列進(jìn)行注釋，并 根據(jù)5因素進(jìn)行距離聚類，對(duì)聚類的方式不作限定，圖5顯示的聚類圖，是將5因素越相近 的contigs歸越近的結(jié)果。
[0078] 圖5顯示出，根據(jù)每條contig的5因素特征，病毒和非病毒序列能聚在一起。黑 色框內(nèi)是病毒序列（V開頭標(biāo)記），其他序列為宿主（非病毒，Η開頭標(biāo)記）序列。
[0079] 由圖5可以看到，V開頭的標(biāo)記序列在一個(gè)分支上，Η開頭的序列在其他分支上，其 他分支根據(jù)具體情況會(huì)有不同分支，這也證明了病毒序列有自己特性，而非病毒序列可能 特性并不單一，但至少和病毒序列是有差別的。
[0080] 4.預(yù)測(cè)模型的建立和評(píng)估
[0081] 從實(shí)施例一中提到的16個(gè)樣本的測(cè)序數(shù)據(jù)中選取了 11個(gè)樣本（Cooks_footf、 GrasslOOf、PoplarlOOf、TCV_add、TCV、TCV-TYMV_add、TCV-TYMV、TGM-CK、TYMV-2、TYMV、 WillowlOOf)構(gòu)建了一個(gè)訓(xùn)練集（命名為trainingll)，為了評(píng)估訓(xùn)練集的準(zhǔn)確性，對(duì)訓(xùn)練 集進(jìn)行了交叉驗(yàn)證。畫了 R0C 曲線（Receiver Operating Characteristic)和 AUC 曲線 (Area Under Curve)，如圖6。對(duì)于訓(xùn)練集來說，AUC值到了 0.8956,說明準(zhǔn)確性比較高。圖 6是訓(xùn)練集評(píng)估的R0C曲線（Receiver Operating Characteristic)圖。圖6中的小表格 是對(duì)訓(xùn)練集進(jìn)行的分割驗(yàn)證。將總的訓(xùn)練集數(shù)據(jù)分成5份，分別獨(dú)立進(jìn)行驗(yàn)證，表格里列出 的是每次的準(zhǔn)確率和準(zhǔn)確的條數(shù)。
[0082] 同時(shí)，為了評(píng)估訓(xùn)練集的可靠性，一些指標(biāo)如Sensitivity (靈敏度，Sn)、 Specificity(特異性，Sp)、Accuracy(準(zhǔn)確率，Acc)、Matthews correlation coefficient (二值相關(guān)系數(shù)，Mcc)和Precision (精確度）用來計(jì)算評(píng)估，其中，上面未提 及的Precision = TP/(TP+FP)，反映了被分類器判定的正例中真正的正例樣本的比重。訓(xùn) 練集重頭分（等分）或者隨機(jī)分成100份，如重頭分割為：1 %、2%、…、100%，然后分別 預(yù)測(cè)自己本身，一百份的每一份都預(yù)測(cè)后計(jì)算那些指標(biāo)值。Sensitivity是預(yù)測(cè)為病毒的 比例，也可以理解為對(duì)于病毒的召回能力，但是這個(gè)值的增加，往往會(huì)帶來Precision的降 低。從訓(xùn)練集的穩(wěn)定性的評(píng)估結(jié)果來看，包括數(shù)據(jù)等分的評(píng)估和隨機(jī)數(shù)據(jù)量的評(píng)估結(jié)果，訓(xùn) 練集trainingll預(yù)測(cè)本身，重頭分割（等分）的所有指標(biāo)的指都穩(wěn)定在80% -100%，而隨 機(jī)分割的穩(wěn)定范圍變大，在30% -100%。此外，precision的值在重頭分割的部分要比在 隨機(jī)分割的部分高，而sensitivity值卻相反。但是，不管這些值如何，在接近100%的地方 這些值都是較高的，因此較佳的，可用100%的數(shù)據(jù)做訓(xùn)練集。
[0083] 實(shí)施例三
[0084] 對(duì)未知序列是否來自病毒的預(yù)測(cè)和驗(yàn)證。
[0085] 1)低質(zhì)量數(shù)據(jù)過濾。將測(cè)序得到的fq文件進(jìn)行去接頭和低質(zhì)量，生成fasta格 式文件，并統(tǒng)計(jì)冗余數(shù)，即，每條序列的ID包含該條序列的重復(fù)數(shù)，如：t00001200。其中 t00001是ID，200是該條序列的重復(fù)數(shù)。
[0086] 2)已知ncRNA過濾。將初步過濾完成的fasta文件和已知的ncRNA數(shù)據(jù)庫(kù)比對(duì)，包 括：http://www. sanger. ac. uk/software/Rfam，GenBank 的非編碼 RNA 數(shù)據(jù)庫(kù)（noncoding RNA database) (http: //www. ncbi. nlm. nih. rov/)，使用 BLAST 2. 2. 23 軟件，要求 e 值 〈0.01。再和miRBase (release 18)數(shù)據(jù)庫(kù)比對(duì)，要求完全一樣比對(duì)上。得到進(jìn)一步過濾的 fasta格式文件。
[0087] 3)數(shù)據(jù)組裝。如果該物種已經(jīng)有基因組序列，那么將生成好的fasta文件使用 BOWTIEvO. 12. 7軟件與基因組進(jìn)行比對(duì)，要求完全匹配。濾除能比上基因組的序列，將剩下 的序列分成4份，分別為總序列條數(shù)的25%，50%，75%和100%，4份數(shù)據(jù)分別進(jìn)行組裝。組 裝軟件使用 SOAPdenovo-Trans vl. 0 (http://soap, genomics, org. cn/SOAPdenovo-Trans. html)，使用Kmer 15和17,將兩個(gè)Kmer組裝結(jié)果混合后去除冗余，定為組裝1。4份數(shù)據(jù)分 別生成組裝1，組裝2,組裝3,組裝4,將4份混合再去冗余，去除50bp以下長(zhǎng)度的序列，最 后使用軟件PHRAP或者minimo進(jìn)行最后的組裝，選擇50 - 100個(gè)堿基長(zhǎng)度作為overlap，生 成最后的組裝結(jié)果，contig的fasta文件。
[0088] 4) Contig和sRNA比對(duì)。將第二步生成的fasta文件和第三步生成的最終Contig 文件使用B0WTIE v0. 12. 7軟件比對(duì)，要求2個(gè)mismatch。比對(duì)完成后，根據(jù)比對(duì)結(jié)果，統(tǒng)計(jì) 每條contig的每個(gè)位點(diǎn)的堿基組成，計(jì)算每個(gè)位點(diǎn)的熵值（參照實(shí)施例一描述），突變率 (單個(gè)位點(diǎn)除了最大頻率堿基的比率外的比率），整條contig最終得到突變位點(diǎn)率（突變 位點(diǎn)的個(gè)數(shù)除以總的位點(diǎn)個(gè)數(shù)的百分比）、平均熵值、中位數(shù)熵值、平均突變率、中位數(shù)突變 率。
[0089] 5)病毒和宿主預(yù)測(cè)。得到contig的5因素信息后，使用R語(yǔ)言的LibSVM軟件包里 的SVC(support vector classification)工具，以我們構(gòu)建好的訓(xùn)練集為模版，對(duì)contig 進(jìn)行分類。結(jié)果中標(biāo)記為1的是病毒，標(biāo)記為-1的為非病毒。具體的，對(duì)樣本中組裝得到 的未知的序列，即，nt和nr比對(duì)都是陰性的contig，進(jìn)行5因素預(yù)測(cè)處理。根據(jù)訓(xùn)練集的 5因素標(biāo)準(zhǔn)，判斷未知contig的5因素綜合結(jié)果是否符合病毒特征，若符合，判斷為病毒，反 之，判斷為非病毒。若未知序列的突變位點(diǎn)率，或者突變位點(diǎn)率、以及平均熵值和中位數(shù)熵 值至少之一，或者突變位點(diǎn)率、以及平均突變率和中位數(shù)突變率的至少之一，或者突變位點(diǎn) 率、平均熵值和中位數(shù)熵值至少之一、以及平均突變率和中位數(shù)突變率的至少之一，或者全 部5個(gè)因素都落在其對(duì)應(yīng)的標(biāo)準(zhǔn)因素的±30%內(nèi)，則判定為符合標(biāo)準(zhǔn)，判定該序列來自病 毒。表2顯示利用建立好的訓(xùn)練集對(duì)未知序列的非同源預(yù)測(cè)結(jié)果。
[0090] 表 2
[0091]
[0092] 為了驗(yàn)證未知序列（contig)是否真的是病毒的，根據(jù)樣本GrasslOOf被預(yù)測(cè)的50 條序列設(shè)計(jì)特異的引物，驗(yàn)證并延伸得到的contig。根據(jù)預(yù)測(cè)的結(jié)果，樣本GrasslOOf中 Contig9、ContiglO、Contigll三條contig被預(yù)測(cè)為病毒，并且都是在nt和nr庫(kù)里比對(duì)陰 性的。分別對(duì)三條序列設(shè)計(jì)引物，RT-PCR驗(yàn)證。圖7顯示驗(yàn)證結(jié)果，如圖7中的A部分，顯 示 congtig9_ll 都有明顯條帶。同時(shí)，對(duì) congtig9-10、congtigl〇-ll、congtig9_ll 分別做 擴(kuò)展延伸驗(yàn)證，最后驗(yàn)證contigl〇-ll和contig9-ll為陽(yáng)性，如圖中的A和B部分所示，驗(yàn) 證三條序列的關(guān)系為congtig9-ll-10,如圖7中的C部分所示，黑色加深的為3730驗(yàn)證區(qū) 域。將三條序列驗(yàn)證、連接后，能鑒定為Bell pepper endornavirus (甜椒內(nèi)源RNA病毒）， 最初組裝得到的同源性很低，不到50%。因此鑒定了病毒在樣本中的存在，驗(yàn)證了預(yù)測(cè)出來 的無法通過同源比對(duì)鑒定的病毒。
[0093] 以上所述僅為本發(fā)明的較佳實(shí)施例，應(yīng)當(dāng)理解，這些實(shí)施例僅用以解釋本發(fā)明，并 不用于限定本發(fā)明。對(duì)于本領(lǐng)域的一般技術(shù)人員，依據(jù)本發(fā)明的思想，可以對(duì)上述具體實(shí)施 方式進(jìn)行變化。
【主權(quán)項(xiàng)】
1. 一種病毒鑒定方法，其特征在于，包括， 獲取待測(cè)樣本的RNA測(cè)序數(shù)據(jù)，所述測(cè)序數(shù)據(jù)包括多個(gè)讀段； 對(duì)所述測(cè)序數(shù)據(jù)的第一部分進(jìn)行組裝，獲得組裝序列，所述測(cè)序數(shù)據(jù)的第一部分包括 不能比對(duì)上ncRNA參考序列的讀段； 將所述測(cè)序數(shù)據(jù)的第一部分與所述組裝序列進(jìn)行比對(duì)，獲得比對(duì)結(jié)果； 基于所述比對(duì)結(jié)果，確定所述組裝序列上的突變位點(diǎn)，所述突變位點(diǎn)包括SNV，以及確 定各條組裝序列的（a)-(c)中的至少之一， (a) 平均熵值和中位數(shù)熵值至少之一，以及突變位點(diǎn)比例， (b) 平均突變率和中位數(shù)突變率至少之一，以及突變位點(diǎn)比例， (c) 突變位點(diǎn)比例， 一條組裝序列的平均熵值為其上的SNV的熵值的平均值， 一條組裝序列的中位數(shù)熵值為其上的SNV的熵值的中位數(shù)， 一條組裝序列的平均突變率為其上的SNV的突變率的平均值， 一條組裝序列的中位數(shù)突變率為其上的SNV的突變率的中位數(shù)， 一條組裝序列的突變位點(diǎn)比例為其上的SNV的數(shù)目所占的比例， SNV的熵值=-100* Σ (Pi*logPi)，Pi為該SNV的各種堿基的深度， SNV的突變率=支持該SNV的讀段數(shù)目/比對(duì)上該SNV的讀段數(shù)目， 將所述確定的組裝序列的（a)-(c)至少之一和與其對(duì)應(yīng)的界限比較，判定落入所述界 限的組裝序列來自病毒。2. 權(quán)利要求1的方法，其特征在于，所述測(cè)序數(shù)據(jù)的第一部分不包括符合（1)-(3)至少 之一的讀段， (1) 包含接頭序列，和/或平均單堿基錯(cuò)誤率大于〇. 01， (2) 長(zhǎng)度小于18nt，和/或長(zhǎng)度大于44nt， (3) 包含堿基質(zhì)量值小于10的堿基個(gè)數(shù)大于2。3. 權(quán)利要求1的方法，其特征在于，所述ncRNA參考序列包括，rRNA參考序列、tRNA參 考序列、snRNA參考序列、snoRNA參考序列和microRNA參考序列。4. 權(quán)利要求1的方法，其特征在于，所述組裝包括， 分別對(duì)所述測(cè)序數(shù)據(jù)的第一部分中的全部讀段和所述測(cè)序數(shù)據(jù)的第一部分中的至少 一部分讀段進(jìn)行第一組裝和第二組裝，獲得第一組裝結(jié)果和第二組裝結(jié)果， 合并所述第一組裝結(jié)果和第二組裝結(jié)果； 任選的， 所述第一組裝和/或所述第二組裝為混合組裝。5. 權(quán)利要求1的方法，其特征在于，在獲得組裝序列之后，去除所述組裝序列中的已知 序列。6. 權(quán)利要求1的方法，其特征在于，在確定每條組裝序列的（a)-(c)中的至少之一之 前，去除該組裝序列中深度小于100X的位點(diǎn)。7. 權(quán)利要求1的方法，其特征在于，還包括， 設(shè)計(jì)引物，對(duì)判定來自病毒的組裝序列進(jìn)行RT-PCR延伸，獲得延伸產(chǎn)物； 利用延伸產(chǎn)物驗(yàn)證所述組裝序列來自病毒。8. 權(quán)利要求1-7任一方法，其特征在于，所述界限的確定，包括， 獲取至少一個(gè)已知病毒感染的樣本的RNA測(cè)序結(jié)果，所述測(cè)序結(jié)果包括多個(gè)讀段； 對(duì)所述測(cè)序結(jié)果的至少一部分進(jìn)行組裝，獲得組裝片段； 將所述組裝片段與參考序列進(jìn)行第一比對(duì)，進(jìn)行組裝片段物種注釋，獲得物種注釋結(jié) 果； 基于所述物種注釋結(jié)果對(duì)所述組裝片段進(jìn)行分類，獲得第一類組裝片段和第二類組裝 片段，所述第一類組裝序列來自病毒，所述第二類組裝片段來自宿主； 分別將所述測(cè)序結(jié)果的至少一部分與所述第一類組裝片段和所述第二類組裝片段進(jìn) 行第二比對(duì)，獲得第二比對(duì)結(jié)果； 基于所述第二比對(duì)結(jié)果，確定所述第一類組裝片段和所述第二類組裝片段上的突變位 點(diǎn)，所述突變位點(diǎn)包括SNV，以及確定每條第一類組裝片段和每條第二類組裝片段的五個(gè)因 素，所述五個(gè)因素為平均熵值、中位數(shù)熵值、平均突變率、中位數(shù)突變率以及突變位點(diǎn)比例， 其中， 一條組裝片段的平均熵值為其上的SNV的熵值的平均值， 一條組裝片段的中位數(shù)熵值為其上的SNV的熵值的中位數(shù)， 一條組裝片段的平均突變率為其上的SNV的突變率的平均值， 一條組裝片段的中位數(shù)突變率為其上的SNV的突變率的中位數(shù)， 一條組裝片段的突變位點(diǎn)比例為其上的SNV的數(shù)目所占的比例， SNV的熵值=-100* Σ (Pi*logPi)，Pi為該SNV的各種堿基的深度， SNV的突變率=支持該SNV的讀段數(shù)目/比對(duì)上該SNV的讀段數(shù)目； 基于所述第一類組裝片段和第二類組裝片段、以及每條所述第一類組裝片段和每條所 述第二類組裝片段的五個(gè)因素，利用SVM，確定所述五個(gè)因素中至少之一的界限； 任選的， 所述第一類組裝片段的數(shù)目大于30 ; 任選的， 所述界限的ROC曲線評(píng)估的AUC值不小于0. 7 ; 任選的， 所述界限的病毒預(yù)測(cè)的準(zhǔn)確率達(dá)90 %。9. 權(quán)利要求8的方法，其特征在于，所述基于第一類組裝片段和第二類組裝片段、以及 每條第一類組裝片段和每條第二類組裝片段的五個(gè)因素，利用SVM，確定所述五個(gè)因素中至 少之一的界限，包括， 將所述第一類組裝片段和第二類組裝片段分為多份片段包，每份所述片段包包含多條 第一類組裝片段和多條第二類組裝片段， 利用其中一份片段包中的組裝片段的五個(gè)因素的至少之一，判斷其它各份片段包中的 每條組裝片段為第一組裝片段還是第二組裝片段，依據(jù)判斷的正確率，調(diào)整確定出所述五 個(gè)因素的至少之一的界限。10. -種鑒定病毒的裝置，其特征在于，包括， 數(shù)據(jù)輸入單元，用于輸入數(shù)據(jù)； 數(shù)據(jù)輸出單元，用于輸出數(shù)據(jù)； 存儲(chǔ)單元，用于存儲(chǔ)數(shù)據(jù)，其中包括計(jì)算機(jī)可執(zhí)行程序； 處理器，與所述數(shù)據(jù)輸入單元、數(shù)據(jù)輸出單元和存儲(chǔ)單元連接，用于執(zhí)行所述程序，執(zhí) 行所述程序包括完成權(quán)利要求1-9任一方法。
【文檔編號(hào)】G06F19/22GK106033502SQ201510125249
【公開日】2016年10月19日
【申請(qǐng)日】2015年3月20日
【發(fā)明人】麻錦敏, 王琿
【申請(qǐng)人】深圳華大基因股份有限公司