基于多步雙向De Bruijn圖的變長(zhǎng)kmer查詢(xún)的頂點(diǎn)擴(kuò)展方法
【專(zhuān)利摘要】本發(fā)明涉及基因測(cè)序【技術(shù)領(lǐng)域】，提供了一種基于多步雙向De?Bruijn圖的變長(zhǎng)kmer查詢(xún)的頂點(diǎn)擴(kuò)展方法，包括：步驟A：讀取測(cè)序數(shù)據(jù)源文件，構(gòu)造多步雙向De?Bruijn圖；步驟B：在所述多步雙向De?Bruijn圖中對(duì)分叉頂點(diǎn)的變長(zhǎng)kmer進(jìn)行構(gòu)造和統(tǒng)計(jì)；步驟C：在所述多步雙向De?Bruijn圖中基于變長(zhǎng)kmer查詢(xún)的頂點(diǎn)擴(kuò)展。本發(fā)明只選取一些分叉的頂點(diǎn)構(gòu)建非常少的一些變長(zhǎng)kmer，然后對(duì)這些分叉頂點(diǎn)進(jìn)行定向解耦，無(wú)需對(duì)每種kmer長(zhǎng)度都去構(gòu)建一個(gè)De?Bruijn圖，可以方便快速地解決所有長(zhǎng)度小于序列長(zhǎng)度的repeat，最大化contig的長(zhǎng)度和質(zhì)量。
【專(zhuān)利說(shuō)明】基于多步雙向De Bru i jn圖的變長(zhǎng)kmer查詢(xún)的頂點(diǎn)擴(kuò)展方法
【【技術(shù)領(lǐng)域】】
[0001]本發(fā)明涉及基因測(cè)序【技術(shù)領(lǐng)域】，特別是涉及一種基于多步雙向De Bruijn圖的變長(zhǎng)kmer查詢(xún)的頂點(diǎn)擴(kuò)展方法。
【【背景技術(shù)】】
[0002]基因序列分析以算法與數(shù)學(xué)模型為核心，研究?jī)?nèi)容涉及多個(gè)方面，主要包括:基因數(shù)據(jù)的存儲(chǔ)與獲取、序列比對(duì)、測(cè)序與拼接、基因預(yù)測(cè)、生物進(jìn)化與系統(tǒng)發(fā)育分析、蛋白質(zhì)結(jié)構(gòu)預(yù)測(cè)、RNA結(jié)構(gòu)預(yù)測(cè)、分子設(shè)計(jì)與藥物設(shè)計(jì)、代謝網(wǎng)絡(luò)分析、基因芯片、DNA計(jì)算等等。現(xiàn)在生物技術(shù)和計(jì)算機(jī)信息處理技術(shù)的緊密結(jié)合，加快了處理生物信息數(shù)據(jù)的速度，使得在盡量短的時(shí)間內(nèi)對(duì)生物學(xué)意義做出盡量準(zhǔn)確的詮釋?zhuān)涌炝松镄畔W(xué)的發(fā)展。目前，生物信息處理成為當(dāng)前信息【技術(shù)領(lǐng)域】面臨的巨大挑戰(zhàn)之一。
[0003]基因序列分析是對(duì)海量基因序列數(shù)據(jù)進(jìn)行分析，從而提取和挖掘新的生物信息知識(shí)。其中，涉及到計(jì)算機(jī)技術(shù)中的機(jī)器學(xué)習(xí)、模式識(shí)別、書(shū)籍分析與挖掘、組合數(shù)學(xué)、隨機(jī)模型、字符串、圖形算法、分布式計(jì)算、高性能計(jì)算、并行計(jì)算等知識(shí)。其中，全基因組學(xué)的研究是當(dāng)前生物信息學(xué)研究的核心之一。
[0004]基因是人類(lèi)最基本的遺傳密碼，代表著每個(gè)人的生命信息?；蛐蛄猩洗嬖谥z傳位點(diǎn)的細(xì)微差異，這些遺傳密碼的多態(tài)性與人類(lèi)的健康、致病機(jī)理、醫(yī)學(xué)治療有著相當(dāng)密切的關(guān)系。其中，DNA測(cè)序是研究全基因組序列需要完成的基本內(nèi)容之一。
[0005]自1977年Sanger測(cè)序技術(shù)問(wèn)世以來(lái),經(jīng)過(guò)三十多年的發(fā)展,DNA測(cè)序技術(shù)發(fā)展突飛猛進(jìn)，以高通量、短序列為特點(diǎn)的第二代測(cè)序技術(shù)逐漸占領(lǐng)市場(chǎng)，以單分子測(cè)序?yàn)樘攸c(diǎn)的第三代測(cè)序技術(shù)也逐漸出現(xiàn)，它們分別在測(cè)序特點(diǎn)上占有不同的優(yōu)勢(shì)。傳統(tǒng)基因測(cè)序方法的數(shù)據(jù)提取和分析軟件經(jīng)過(guò)近10年來(lái)的研究與開(kāi)發(fā)，目前已經(jīng)較為完善。但是，測(cè)序技術(shù)的發(fā)展，帶來(lái)了測(cè)序數(shù)據(jù)的變化，使得當(dāng)前存在的數(shù)據(jù)處理軟件不能滿足當(dāng)前生物醫(yī)學(xué)研究的需求。
[0006]新一代高通量測(cè)序方法的應(yīng)用，可以在短時(shí)間內(nèi)完成整個(gè)基因組數(shù)據(jù)的測(cè)定。高通量測(cè)序方法的日新月異也同時(shí)對(duì)獲取的基因數(shù)據(jù)的分析處理方法提出了挑戰(zhàn)。在這個(gè)目前炙手可熱的研究領(lǐng)域，迫切需要開(kāi)發(fā)能滿足高通量測(cè)序技術(shù)的海量數(shù)據(jù)處理的生物信息學(xué)平臺(tái)。面對(duì)個(gè)人基因組計(jì)劃及未來(lái)的個(gè)性化醫(yī)療前景，高效低成本的測(cè)序技術(shù)成為必然的趨勢(shì)。同時(shí)，簡(jiǎn)化高效的一站式完備的生物信息學(xué)數(shù)據(jù)分析平臺(tái)等完備的測(cè)序解決方案，也是極為重要、不可或缺的發(fā)展方向。
[0007]然而新一代的高通量測(cè)序方法雖然測(cè)序通量高，但是卻會(huì)引入測(cè)序誤差，同時(shí)測(cè)序樣本本身由于基因突變，測(cè)序不均勻而導(dǎo)致有SNP(Single Nucleotide Polymorphisms,單核苷酸多態(tài)性)的出現(xiàn)，而這些測(cè)序誤差、SNP、測(cè)序不均將會(huì)在基因組組裝時(shí)構(gòu)造的多步雙向De Bruijn圖中引入一些錯(cuò)誤的雙向邊，其中有一部分是自環(huán)雙向邊。而這些錯(cuò)誤的自環(huán)雙向邊在De Bruijn圖中，能夠阻礙圖的收縮，contig無(wú)法擴(kuò)展,最終使得contig的長(zhǎng)度和質(zhì)量都很低。
[0008]綜上所述，新一代的高通量測(cè)序方法產(chǎn)生的短基因片段的組裝導(dǎo)致了大量的測(cè)序錯(cuò)誤，這大大加大了組裝算法的計(jì)算量。大量的測(cè)序錯(cuò)誤，使得組裝錯(cuò)誤率增加，嚴(yán)重影響了組裝結(jié)果。能否有效地解決這個(gè)問(wèn)題，成為評(píng)價(jià)一個(gè)組裝算法優(yōu)劣的關(guān)鍵。
[0009]目前組裝算法的策略主要分為兩類(lèi)，一個(gè)是基于Overlap-Layout-Consensus (OLC)的算法，另外一個(gè)是基于De Bruijn圖的算法。其中基于OLC組裝算法開(kāi)發(fā)的軟件，如SSAKE、VCAKE, SHARCGS等，在基因長(zhǎng)序列組裝中更占有優(yōu)勢(shì)，但并不完全適用于新一代的短序列組裝。與OLC組裝算法不同，De Bruijn算法不再以read為單位組織數(shù)據(jù)，而是以k-mers為單位進(jìn)行數(shù)據(jù)組裝，其優(yōu)點(diǎn)主要有以下幾個(gè)方面:首先，以k-mers為單位進(jìn)行序列組裝，不影響節(jié)點(diǎn)的質(zhì)量，減少了冗余數(shù)據(jù)量。其次，在圖中重復(fù)區(qū)域只出現(xiàn)一次，便于識(shí)別，可以避免錯(cuò)誤的組裝，減小出錯(cuò)率。最后，采取將有重疊區(qū)域映射到同一條弧上的策略，從而簡(jiǎn)化了搜索路徑。目前，很多短序列組裝算法都使用這種框架，如 Velvet、IDBA、SOAPdenovo, ABySS 等。
[0010]其中Velvet有效地利用了 De Brui jn圖，實(shí)現(xiàn)了高效地短序列組裝。Velvet以k-mer為基本單位構(gòu)建De Bruijn圖,利用圖的結(jié)構(gòu)，結(jié)合相應(yīng)的序列特征,簡(jiǎn)化圖的構(gòu)造，最終找到一條最優(yōu)路徑完成組裝過(guò)程。Velvet把焦點(diǎn)集中在錯(cuò)誤的數(shù)據(jù)產(chǎn)生的三種結(jié)構(gòu)上，即tip,bubble,以及erroneous connection。它依照長(zhǎng)度原則和少數(shù)性原則,將長(zhǎng)度小于2k的均去除；利用Tour Bus算法中的深度優(yōu)先搜索策略合并bubble，最后利用覆蓋度閾值法去除了 erroneous connection。該方法也充分利用了 paired-end雙端信息,進(jìn)一步解決repeat問(wèn)題,優(yōu)化了組裝效果。Velvet充分利用圖的結(jié)構(gòu)性質(zhì),簡(jiǎn)化了數(shù)據(jù)冗余,速度較之前的算法有了很大的改進(jìn)。雖然它沒(méi)有在預(yù)處理階段對(duì)序列進(jìn)行糾錯(cuò)，但是其對(duì)錯(cuò)誤的預(yù)防機(jī)制，很大程度上彌補(bǔ)了這方面的缺陷，這使得它更好地應(yīng)用在大型基因組序列的組裝中。
[0011]IDBA也是基于De Bruijn圖，實(shí)現(xiàn)了簡(jiǎn)便且高效的短序列組裝。IDBA以k-mer為基本單位，與以往不同的是，它采用一個(gè)變化的k值域(Kmin-Kmax)，代替使用固定的k值來(lái)得到k-mers的長(zhǎng)度。由于基因組裝以k-mers為單位，通常會(huì)形成很多個(gè)重疊單元，這使得組裝面臨著錯(cuò)誤位置組裝、頂點(diǎn)缺失和覆蓋度低的問(wèn)題。正確的選擇k值的大小成為組裝的一個(gè)關(guān)鍵因素。一些錯(cuò)誤的reads的產(chǎn)生,也導(dǎo)致產(chǎn)生了大量的branching。K值越小，branching問(wèn)題越嚴(yán)重,k值越大,貝U出現(xiàn)的repeat區(qū)域則變少,這直接影響了組裝的質(zhì)量。IDBA采用不固定的k值進(jìn)行組裝，可以很好的解決branching問(wèn)題，從而，提高了組裝的質(zhì)量。另外IDBA通過(guò)刪除低覆蓋率的錯(cuò)誤k-mers而使得IDBA的內(nèi)存使用率明顯降低，同時(shí)也提升了 IDBA的處理速度。
[0012]SOAPdenovo能夠高效高質(zhì)量地完成數(shù)以?xún)|計(jì)的reads的組裝。SOAPdenovo繼承了OLC算法和De Bruijn圖算法的優(yōu)點(diǎn)，使得其組裝質(zhì)量大為提高。SOAP通過(guò)預(yù)置k-mer閾值的方法，采取過(guò)濾、糾錯(cuò)的方式減少了錯(cuò)誤序列的產(chǎn)生。同時(shí)，它借鑒了 Velvet軟件的方法成功處理了 bubble，使得其平均覆蓋度增加。另外，SOAPdenovo利用了雙端信息進(jìn)行重疊區(qū)域匹配，并合并read生成contig片段,生成基于contig的圖結(jié)構(gòu),從而，SOAPdenovo大大簡(jiǎn)化了 contig圖的復(fù)雜性。
[0013]ABySS引進(jìn)并行計(jì)算的思想，搭建了一個(gè)Iinux集群，在集群上建立了一個(gè)分布式的De Bruijn圖結(jié)構(gòu)，將數(shù)據(jù)分布式存儲(chǔ)于每個(gè)節(jié)點(diǎn)上。其采用MPI通信機(jī)制完成節(jié)點(diǎn)之間的相互通信。從構(gòu)建圖、糾錯(cuò)處理到后面的定點(diǎn)融合，最后完成整個(gè)基因組序列的再現(xiàn)，其在運(yùn)行時(shí)間和內(nèi)存消耗方面占有很大的優(yōu)勢(shì)，并且其錯(cuò)誤率極低，在性能方面特別是cluster中單機(jī)內(nèi)存使用上均有很大的提升，正在得到越來(lái)越廣泛的應(yīng)用。
[0014]現(xiàn)有的主要序列組裝軟件，例如SOAPdenovo, Velvet, ABySS, Ray等，都是基于給定長(zhǎng)度的kmer進(jìn)行De Brui jn構(gòu)建,然后進(jìn)行收縮。其優(yōu)化的辦法也只有去選擇最佳的一個(gè)kmer長(zhǎng)度。這種基于固定長(zhǎng)度kmer的組裝策略對(duì)于所有長(zhǎng)度大約kmer長(zhǎng)度的重復(fù)序列無(wú)法解耦。雖然IDBA能夠?qū)Χ喾Nkmer長(zhǎng)度進(jìn)行迭代收縮De Bruijn圖，但是它需要對(duì)每種kmer長(zhǎng)度去將所有的序列進(jìn)行分解存儲(chǔ)和計(jì)算，該策略將會(huì)耗費(fèi)巨大的內(nèi)存和計(jì)算時(shí)間。
[0015]鑒于此，克服該現(xiàn)有技術(shù)所存在的缺陷是本【技術(shù)領(lǐng)域】亟待解決的問(wèn)題。
【
【發(fā)明內(nèi)容】
】
[0016]本發(fā)明要解決的技術(shù)問(wèn)題是提供一種基于多步雙向De Bruijn圖的變長(zhǎng)kmer查詢(xún)的頂點(diǎn)擴(kuò)展方法。
[0017]本發(fā)明采用如下技術(shù)方案:
[0018]—種基于多步雙向De Bruijn圖的變長(zhǎng)kmer查詢(xún)的頂點(diǎn)擴(kuò)展方法,包括:
[0019]步驟A:讀取測(cè)序數(shù)據(jù)源文件，構(gòu)造多步雙向De Bruijn圖；
[0020]步驟B:在所述多步雙向De Brui jn圖中對(duì)分叉頂點(diǎn)的變長(zhǎng)kmer進(jìn)行構(gòu)造和統(tǒng)計(jì)；
[0021]步驟C:在所述多步雙向De Bruijn圖中基于變長(zhǎng)kmer查詢(xún)的頂點(diǎn)擴(kuò)展。
[0022]進(jìn)一步地,所述步驟B中，對(duì)所述多步雙向De Brui jn圖中的頂點(diǎn)所有可能的分叉合并路徑上的k+2長(zhǎng)的變長(zhǎng)kmer構(gòu)造權(quán)重表，選擇權(quán)重最高的一組分叉路徑組合進(jìn)行分叉路徑上的雙向邊合并。
[0023]進(jìn)一步地，所述步驟C中，對(duì)于一個(gè)給定的分叉頂點(diǎn)U，在查詢(xún)頂點(diǎn)u所有的k+2長(zhǎng)的變長(zhǎng)kmer權(quán)重值之后，選擇權(quán)重最高的一組分叉路徑組合進(jìn)行分叉路徑上的雙向邊合并，同時(shí)刪除合并前的被選擇的分叉雙向邊。
[0024]進(jìn)一步地，所述權(quán)重為變長(zhǎng)kmer出現(xiàn)次數(shù)或變長(zhǎng)kmer模糊匹配加權(quán)次數(shù)。
[0025]進(jìn)一步地,所述步驟B進(jìn)一步包括:
[0026]步驟B1:遍歷所述多步雙向De Bruijn圖中的每個(gè)頂點(diǎn)u ；
[0027]步驟B2:統(tǒng)計(jì)頂點(diǎn)u中正向邊的個(gè)數(shù)P和反向邊的個(gè)數(shù)q ;
[0028]步驟B3:若p+q大于等于3且p和q均至少為1，則執(zhí)行步驟B4，否則返回執(zhí)行步驟BI ；
[0029]步驟B4:計(jì)算出頂點(diǎn)u的q個(gè)反向邊的對(duì)偶雙向邊,并將對(duì)偶雙向邊的倒數(shù)第k+1個(gè)字符取出存到入邊字符數(shù)組m ;
[0030]步驟B5:將頂點(diǎn)u的P個(gè)正向邊的第一個(gè)字符存到出邊字符數(shù)組η ；
[0031]步驟Β6:將(m,頂點(diǎn)u的正向字符串，η)所有的組合構(gòu)成的k+2長(zhǎng)的kmer記錄為變長(zhǎng)kmer數(shù)組。
[0032]進(jìn)一步地,所述步驟C進(jìn)一步包括:[0033]步驟Cl:打開(kāi)測(cè)序序列文件，逐個(gè)讀取每條序列；
[0034]步驟C2:將所述變長(zhǎng)kmer數(shù)組逐個(gè)匹配讀入的序列，并對(duì)每個(gè)變長(zhǎng)kmer計(jì)數(shù)；
[0035]步驟C3:遍歷所述多步雙向De Bruijn圖中的每個(gè)頂點(diǎn)u ；
[0036]步驟C4:統(tǒng)計(jì)頂點(diǎn)u中正向邊的個(gè)數(shù)P，反向邊的個(gè)數(shù)q ;
[0037]步驟C5:若p+q大于等于3且p和q均至少為1，則執(zhí)行步驟C6，否則返回執(zhí)行步驟C3 ；
[0038]步驟C6:計(jì)算出頂點(diǎn)u的q個(gè)反向邊的對(duì)偶雙向邊,并將對(duì)偶雙向邊的倒數(shù)第k+1個(gè)字符取出存到入邊字符數(shù)組m ;
[0039]步驟C7:將頂點(diǎn)u的P個(gè)正向邊的第一個(gè)字符存到出邊字符數(shù)組η ；
[0040]步驟C8:查詢(xún)由(m,頂點(diǎn)u的正向字符串，η)的所有組合構(gòu)成的k+2長(zhǎng)的kmer的出現(xiàn)次數(shù)，選擇出現(xiàn)次數(shù)最大的一組正向邊和反向邊進(jìn)行合并擴(kuò)展。
[0041]進(jìn)一步地，所述步驟C2中將所述變長(zhǎng)kmer數(shù)組逐個(gè)匹配讀入的序列包括對(duì)反向序列的匹配。
[0042]進(jìn)一步地,所述步驟A進(jìn)一步包括:
[0043]壓縮存儲(chǔ)步驟，具體為
[0044]All、讀取一個(gè)序列s ;
[0045]A12、將序列s用滑動(dòng)窗口切割為多個(gè)片段t ；
[0046]A13、對(duì)每個(gè)片段t，使用核酸編碼表進(jìn)行編碼，并表示為一個(gè)64位的整數(shù)a ；
[0047]A14、將片段t進(jìn)行反轉(zhuǎn)，使用對(duì)稱(chēng)互補(bǔ)表將反轉(zhuǎn)的片段互補(bǔ)處理，得到互補(bǔ)片段ν,并再次使用步驟A13中的核酸編碼表將互補(bǔ)片段進(jìn)行編碼，并表示為一個(gè)64位的整數(shù)b ；
[0048]A15、取整數(shù)a和整數(shù)b的最大數(shù)，作為片段t和互補(bǔ)片段ν的k分子的標(biāo)志數(shù)；
[0049]A16、重復(fù)步驟A11-A15，直至所有序列完成；
[0050]和De Brui jn圖構(gòu)造步驟,具體為[0051 ] A21、讀取一個(gè)序列s ；
[0052]A22、將序列s用滑動(dòng)窗口切割為多個(gè)片段t，選取一片段t其標(biāo)志數(shù)為cur、并標(biāo)記其前、后的片段的標(biāo)志數(shù)分別為pre、Iat ；
[0053]A23、若t的編碼小于其互補(bǔ)片段編碼，則交換pre，Iat的值；
[0054]A24、在cur的正向位置映射表的相應(yīng)bit位置I來(lái)表示指向pre的邊；
[0055]A25、在cur的反向位置映射表的相應(yīng)bit位置I來(lái)表示指向Iat的邊；
[0056]A26、重復(fù)步驟A22-A25，處理序列s的其他片段t，直至完成序列s的全部片段t，執(zhí)行步驟A27 ；
[0057]A27、讀取一個(gè)新的序列S，重復(fù)步驟A22-A26 ;直至處理完所有的序列，執(zhí)行步驟A28 ；
[0058]A28、完成雙向多步De Bruijn圖的構(gòu)造。
[0059]進(jìn)一步地，所述步驟A12、A22中的滑動(dòng)窗口為長(zhǎng)度為k的滑動(dòng)窗口，其中0〈k〈32且k為奇數(shù)；
[0060]所述步驟A13中的核酸編碼表為{A:00，C:01，G:10，T:11}；
[0061 ] 所述步驟A14中的對(duì)稱(chēng)互補(bǔ)表為{A->T，C->G, G->C, T_>A}；[0062]所述步驟A14具體為，將片段t的字符串進(jìn)行反轉(zhuǎn)，使用對(duì)稱(chēng)互補(bǔ)表將反轉(zhuǎn)的字符串中每個(gè)字符變?yōu)槠浠パa(bǔ)字符，得到互補(bǔ)字符的字符串V，并再次使用步驟A13中的核酸編碼表將字符串ν進(jìn)行編碼，并表示為一個(gè)64位的整數(shù)b。
[0063]進(jìn)一步地，所述步驟A22中，若片段t沒(méi)有之前或之后的片段，則對(duì)pre或者Iat值賦為空或NULL ；
[0064]步驟A24中正向位置映射表為{A:0，C:l，G:2，T:3}，位置查詢(xún)字符為pre的最后一個(gè)字符；
[0065]步驟A25中反向位置映射表為{A:4，C:5, G:6, T:7}，位置查詢(xún)字符為Iat的第一個(gè)字符的互補(bǔ)字符。
[0066]與現(xiàn)有技術(shù)相比，本發(fā)明的有益效果在于:
[0067](I)針對(duì)De Bruijn圖上存在的分叉頂點(diǎn)構(gòu)造長(zhǎng)度為k+2的變長(zhǎng)kmer組合，并在輸入序列中統(tǒng)計(jì)其出現(xiàn)次數(shù)，然后根據(jù)其出現(xiàn)次數(shù)選擇頂點(diǎn)上具有最大權(quán)重的雙向邊合并；而IDBA方法是通過(guò)迭代所有的kmer長(zhǎng)度,需要將每個(gè)kmer長(zhǎng)度的所有可能的kmer都構(gòu)造出來(lái)后，再收縮De Bruijn圖，其方法將導(dǎo)致更大的內(nèi)存消耗和計(jì)算時(shí)間消耗；
[0068](2)選擇最優(yōu)的分叉邊的合并，將頂點(diǎn)上的分叉邊的錯(cuò)誤合并的可能降到最低；
[0069](3)可以顯著提高contig的長(zhǎng)度,也能夠?qū)ontig的質(zhì)量損失降到最??；相比于其他已有方法的提高contig長(zhǎng)度就得犧牲c(diǎn)ontig質(zhì)量，本發(fā)明有了一定程度上的控制和改善。
【【專(zhuān)利附圖】

【附圖說(shuō)明】】
`[0070]圖1是本發(fā)明實(shí)施例基于多步雙向De Bruijn圖的變長(zhǎng)kmer查詢(xún)的頂點(diǎn)擴(kuò)展方法流程圖；
[0071]圖2是圖1步驟A中的壓縮存儲(chǔ)步驟流程圖；
[0072]圖3是圖1步驟A中De Bruijn圖構(gòu)造步驟流程圖；
[0073]圖4是步驟B在多步雙向De Bruijn圖中對(duì)分叉頂點(diǎn)的變長(zhǎng)kmer進(jìn)行構(gòu)造和統(tǒng)計(jì)的流程圖；
[0074]圖5是步驟C在多步雙向De Bruijn圖中基于變長(zhǎng)kmer查詢(xún)的頂點(diǎn)擴(kuò)展的流程圖。
【【具體實(shí)施方式】】
[0075]為了使本發(fā)明的目的、技術(shù)方案及優(yōu)點(diǎn)更加清楚明白，以下結(jié)合附圖及實(shí)施例，對(duì)本發(fā)明進(jìn)行進(jìn)一步詳細(xì)說(shuō)明。應(yīng)當(dāng)理解，此處所描述的具體實(shí)施例僅僅用以解釋本發(fā)明，并不用于限定本發(fā)明。
[0076]此外，下面所描述的本發(fā)明各個(gè)實(shí)施方式中所涉及到的技術(shù)特征只要彼此之間未構(gòu)成沖突就可以相互組合。
[0077]本發(fā)明的目的在于設(shè)計(jì)一種基于變長(zhǎng)kmer查詢(xún)的分叉頂點(diǎn)擴(kuò)展方法，它將使DeBruijn圖繼續(xù)收縮，contigs繼續(xù)擴(kuò)展，同時(shí)不會(huì)引入錯(cuò)誤，造成contig質(zhì)量的下降,準(zhǔn)確度降低。
[0078]本發(fā)明提供了一種基于多步雙向De Bruijn圖的變長(zhǎng)kmer查詢(xún)的頂點(diǎn)擴(kuò)展方法，如圖1所示，該方法包括:
[0079]步驟A:讀取測(cè)序數(shù)據(jù)源文件，構(gòu)造多步雙向De Bruijn圖；
[0080]步驟B:在多步雙向De Bruijn圖中對(duì)分叉頂點(diǎn)的變長(zhǎng)kmer進(jìn)行構(gòu)造和統(tǒng)計(jì)；
[0081]步驟C:在多步雙向De Bruijn圖中基于變長(zhǎng)kmer查詢(xún)的頂點(diǎn)擴(kuò)展。
[0082]其中，步驟A可通過(guò)如下方式具體實(shí)現(xiàn):
[0083]壓縮存儲(chǔ)步驟，所需原始數(shù)據(jù)包括第一代，第二代和新一代的測(cè)序儀器產(chǎn)生出來(lái)的FASTA格式文件，將FASTA文件中的序列逐個(gè)切割成k分子并且用二進(jìn)制編碼進(jìn)行壓縮存儲(chǔ)為一個(gè)64位的長(zhǎng)整型k分子的標(biāo)志數(shù)。
[0084]如圖2所示，具體為
[0085]All、讀取一個(gè)序列s ;其中，序列s取自FASTA格式文件；
[0086]A12、將序列s用滑動(dòng)窗口切割為多個(gè)片段t ；
[0087]A13、對(duì)每個(gè)片段t，使用核酸編碼表進(jìn)行編碼，并表示為一個(gè)64位的整數(shù)a ；
[0088]A14、將片段t進(jìn)行反轉(zhuǎn)，使用對(duì)稱(chēng)互補(bǔ)表將反轉(zhuǎn)的片段互補(bǔ)處理，得到互補(bǔ)片段，并再次使用步驟A13中的核酸編碼表將互補(bǔ)片段進(jìn)行編碼，并表示為一個(gè)64位的整數(shù)b ;
[0089]A15、取整數(shù)a和整數(shù)b的最大數(shù)，作為片段t和互補(bǔ)片段ν的k分子的標(biāo)志數(shù)；
[0090]A16、重復(fù)步驟A11-A15，直至所有序列完成。
[0091]通過(guò)上述步驟將兩個(gè)傳統(tǒng)的De Brujin圖中的kmer,轉(zhuǎn)化為一個(gè)64位的k分子的標(biāo)志數(shù)來(lái)存儲(chǔ)。該步驟可以將其他軟件例如velvet、IDBA、S0APdenovo里的兩個(gè)壓縮kmer存儲(chǔ)為一個(gè)壓縮k分子的標(biāo)志數(shù)，并且在得到k分子的標(biāo)志數(shù)后也可以反過(guò)來(lái)求出該k分子的長(zhǎng)度為k的片段t和它的互補(bǔ)片段V。
[0092]步驟A12、A22中的滑動(dòng)窗口為長(zhǎng)度為k的滑動(dòng)窗口，其中0〈k〈32且k為奇數(shù)；步驟A13中的核酸編碼表為{A:00，C:01，G:10，T:11};步驟A14中的對(duì)稱(chēng)互補(bǔ)表為{A->T, C->G, G->C, T->A};步驟A14具體為，將片段t的字符串進(jìn)行反轉(zhuǎn)，使用對(duì)稱(chēng)互補(bǔ)表將反轉(zhuǎn)的字符串中每個(gè)字符變?yōu)槠浠パa(bǔ)字符，得到互補(bǔ)字符的字符串V，并再次使用步驟A13中的核酸編碼表將字符串ν進(jìn)行編碼，并表示為一個(gè)64位的整數(shù)b。
[0093]和De Bruijn圖構(gòu)造步驟，1、使用上述壓縮存儲(chǔ)步驟中計(jì)算k分子的標(biāo)志數(shù),2、將每個(gè)片段以及和它前后相鄰的片段的擴(kuò)展字符作為該k分子和其前后相鄰的片段的對(duì)應(yīng)的k分子的邊并初始化k分子數(shù)據(jù)結(jié)構(gòu)的邊；3、將初始化后的k分子數(shù)據(jù)結(jié)構(gòu)以k分子的標(biāo)志數(shù)為關(guān)鍵值存入hash_map。
[0094]如圖3所示，具體為
[0095]A21、讀取一個(gè)序列s;
[0096]A22、將序列s用滑動(dòng)窗口切割為多個(gè)片段t，選取一片段t其標(biāo)志數(shù)為cur、并標(biāo)記其前、后的片段的標(biāo)志數(shù)分別為pre、Iat ；
[0097]A23、若t的編碼小于其互補(bǔ)片段編碼，則交換pre，Iat的值；
[0098]A24、在cur的正向位置映射表的相應(yīng)bit位置I來(lái)表示指向pre的邊；
[0099]A25、在cur的反向位置映射表的相應(yīng)bit位置I來(lái)表示指向Iat的邊；
[0100]A26、重復(fù)步驟A22-A25，處理序列s的其他片段t，直至完成序列s的全部片段t，執(zhí)行步驟S27 ；
[0101]A27、讀取一個(gè)新的序列S，重復(fù)步驟A22-A26 ;直至處理完所有的序列，執(zhí)行步驟A28;
[0102]A28、完成雙向多步De Bruijn圖的構(gòu)造。
[0103]步驟A22中，若片段t沒(méi)有之前或之后的片段，則對(duì)pre或者Iat值賦為空或NULL ；步驟A24中正向位置映射表為{A:0，C:1，G:2，T:3}，位置查詢(xún)字符為pre的最后一個(gè)字符；步驟A25中反向位置映射表為{A:4，C:5, G:6, T:7}，位置查詢(xún)字符為Iat的第一個(gè)字符的
互補(bǔ)字符。
[0104]步驟B中，對(duì)多步雙向De Bruijn圖中的頂點(diǎn)所有可能的分叉合并路徑上的k+2長(zhǎng)的變長(zhǎng)kmer構(gòu)造權(quán)重表，選擇權(quán)重最高的一組分叉路徑組合進(jìn)行分叉路徑上的雙向邊合并。權(quán)重為變長(zhǎng)kmer出現(xiàn)次數(shù)或變長(zhǎng)kmer模糊匹配加權(quán)次數(shù)。
[0105]本方法中設(shè)定多步雙向De Bruijn圖中的每個(gè)頂點(diǎn)的數(shù)據(jù)結(jié)構(gòu)為:
[0106]
【權(quán)利要求】
1.一種基于多步雙向De Bruijn圖的變長(zhǎng)kmer查詢(xún)的頂點(diǎn)擴(kuò)展方法,其特征在于,包括: 步驟A:讀取測(cè)序數(shù)據(jù)源文件，構(gòu)造多步雙向De Bruijn圖； 步驟B:在所述多步雙向De Bruijn圖中對(duì)分叉頂點(diǎn)的變長(zhǎng)kmer進(jìn)行構(gòu)造和統(tǒng)計(jì)； 步驟C:在所述多步雙向De Bruijn圖中基于變長(zhǎng)kmer查詢(xún)的頂點(diǎn)擴(kuò)展。
2.如權(quán)利要求1所述的方法，其特征在于，所述步驟B中，對(duì)所述多步雙向DeBruijn圖中的頂點(diǎn)所有可能的分叉合并路徑上的k+2長(zhǎng)的變長(zhǎng)kmer構(gòu)造權(quán)重表，選擇權(quán)重最高的一組分叉路徑組合進(jìn)行分叉路徑上的雙向邊合并。
3.如權(quán)利要求1所述的方法，其特征在于，所述步驟C中，對(duì)于一個(gè)給定的分叉頂點(diǎn)U，在查詢(xún)頂點(diǎn)u所有的k+2長(zhǎng)的變長(zhǎng)kmer權(quán)重值之后，選擇權(quán)重最高的一組分叉路徑組合進(jìn)行分叉路徑上的雙向邊合并，同時(shí)刪除合并前的被選擇的分叉雙向邊。
4.如權(quán)利要求2或3所述的方法，其特征在于，所述權(quán)重為變長(zhǎng)kmer出現(xiàn)次數(shù)或變長(zhǎng)kmer模糊匹配加權(quán)次數(shù)。
5.如權(quán)利要求1所述的方法，其特征在于，所述步驟B進(jìn)一步包括: 步驟B1:遍歷所述多步雙向De Bruijn圖中的每個(gè)頂點(diǎn)u ； 步驟B2:統(tǒng)計(jì)頂點(diǎn)u中正向邊的個(gè)數(shù)P和反向邊的個(gè)數(shù)q ； 步驟B3:若p+q大于等于3且P和q均至少為1，則執(zhí)行步驟B4，否則返回執(zhí)行步驟BI ； 步驟B4:計(jì)算出頂點(diǎn)u的q個(gè)反向邊的對(duì)偶雙向邊,并將對(duì)偶雙向邊的倒數(shù)第k+Ι個(gè)字符取出存到入邊字符數(shù)組m ; 步驟B5:將頂點(diǎn)u的P個(gè)正向邊的第一個(gè)字符存到出邊字符數(shù)組η ； 步驟Β6:將(m,頂點(diǎn)u的正向字符串，η)所有的組合構(gòu)成的k+2長(zhǎng)的kmer記錄為變長(zhǎng)kmer數(shù)組。
6.如權(quán)利要求1所述的方法，其特征在于，所述步驟C進(jìn)一步包括: 步驟Cl:打開(kāi)測(cè)序序列文件，逐個(gè)讀取每條序列； 步驟C2:將所述變長(zhǎng)kmer數(shù)組逐個(gè)匹配讀入的序列，并對(duì)每個(gè)變長(zhǎng)kmer計(jì)數(shù)； 步驟C3:遍歷所述多步雙向De Bruijn圖中的每個(gè)頂點(diǎn)u ； 步驟C4:統(tǒng)計(jì)頂點(diǎn)u中正向邊的個(gè)數(shù)P，反向邊的個(gè)數(shù)q ； 步驟C5:若p+q大于等于3且P和q均至少為1，則執(zhí)行步驟C6，否則返回執(zhí)行步驟C3 ； 步驟C6:計(jì)算出頂點(diǎn)u的q個(gè)反向邊的對(duì)偶雙向邊,并將對(duì)偶雙向邊的倒數(shù)第k+Ι個(gè)字符取出存到入邊字符數(shù)組m ; 步驟C7:將頂點(diǎn)u的P個(gè)正向邊的第一個(gè)字符存到出邊字符數(shù)組η ； 步驟C8:查詢(xún)由(m,頂點(diǎn)u的正向字符串,η)的所有組合構(gòu)成的k+2長(zhǎng)的kmer的出現(xiàn)次數(shù)，選擇出現(xiàn)次數(shù)最大的一組正向邊和反向邊進(jìn)行合并擴(kuò)展。
7.如權(quán)利要求6所述的方法，其特征在于，所述步驟C2中將所述變長(zhǎng)kmer數(shù)組逐個(gè)匹配讀入的序列包括對(duì)反向序列的匹配。
8.如權(quán)利要求1所述的方法，其特征在于，所述步驟A進(jìn)一步包括: 壓縮存儲(chǔ)步驟，具體為Al 1、讀取一個(gè)序列s ； A12、將序列s用滑動(dòng)窗口切割為多個(gè)片段t ; A13、對(duì)每個(gè)片段t，使用核酸編碼表進(jìn)行編碼，并表示為一個(gè)64位的整數(shù)a ； A14、將片段t進(jìn)行反轉(zhuǎn)，使用對(duì)稱(chēng)互補(bǔ)表將反轉(zhuǎn)的片段互補(bǔ)處理，得到互補(bǔ)片段V，并再次使用步驟A13中的核酸編碼表將互補(bǔ)片段進(jìn)行編碼，并表示為一個(gè)64位的整數(shù)b ;A15、取整數(shù)a和整數(shù)b的最大數(shù)，作為片段t和互補(bǔ)片段ν的k分子的標(biāo)志數(shù)； A16、重復(fù)步驟A11-A15,直至所有序列完成； 和De Brui jn圖構(gòu)造步驟,具體為 A21、讀取一個(gè)序列s ； A22、將序列s用滑動(dòng)窗口切割為多個(gè)片段t，選取一片段t其標(biāo)志數(shù)為cur、并標(biāo)記其前、后的片段的標(biāo)志數(shù)分別為pre、Iat ； A23、若t的編碼小于其互補(bǔ)片段編碼，則交換pre，Iat的值； A24、在cur的正向位置映射表的相應(yīng)bit位置I來(lái)表示指向pre的邊； A25、在cur的反向位置映射表的相應(yīng)bit位置I來(lái)表示指向Iat的邊； A26、重復(fù)步驟A22-A25，處理序列s的其他片段t，直至完成序列s的全部片段t，執(zhí)行步驟A27 ； A27、讀取一個(gè)新的序列S，重復(fù)步驟A22-A26 ;直至處理完所有的序列，執(zhí)行步驟A28 ； A28、完成雙向多步De Bruijn圖的構(gòu)造。`
9.如權(quán)利要求8所述的方法，其特征在于，所述步驟A12、A22中的滑動(dòng)窗口為長(zhǎng)度為k的滑動(dòng)窗口，其中0〈k〈32且k為奇數(shù)； 所述步驟A13中的核酸編碼表為{A:00，C:01，G:10，T:11}； 所述步驟A14中的對(duì)稱(chēng)互補(bǔ)表為{A->T，C->G，G->C，T->A}； 所述步驟A14具體為，將片段t的字符串進(jìn)行反轉(zhuǎn)，使用對(duì)稱(chēng)互補(bǔ)表將反轉(zhuǎn)的字符串中每個(gè)字符變?yōu)槠浠パa(bǔ)字符，得到互補(bǔ)字符的字符串V，并再次使用步驟A13中的核酸編碼表將字符串ν進(jìn)行編碼，并表示為一個(gè)64位的整數(shù)b。
10.如權(quán)利要求8所述的方法，其特征在于，所述步驟A22中，若片段t沒(méi)有之前或之后的片段，則對(duì)pre或者Iat值賦為空或NULL ； 步驟A24中正向位置映射表為{A:0，C:l, G:2, Τ: 3}，位置查詢(xún)字符為pre的最后一個(gè)字符； 步驟A25中反向位置映射表為{A:4，C:5, G:6, T:7}，位置查詢(xún)字符為Iat的第一個(gè)字符的互補(bǔ)字符。
【文檔編號(hào)】G06F19/22GK103699819SQ201310670752
【公開(kāi)日】2014年4月2日 申請(qǐng)日期:2013年12月10日 優(yōu)先權(quán)日:2013年12月10日
【發(fā)明者】孟金濤, 張慧琳, 彭豐斌, 魏彥杰, 馮圣中 申請(qǐng)人:深圳先進(jìn)技術(shù)研究院