本發(fā)明屬于生物信息工程領(lǐng)域，涉及生物信息技術(shù)和計(jì)算機(jī)應(yīng)用技術(shù)，具體地說，涉及DNA序列二代測(cè)序的短序列快速比對(duì)分析方法。

背景技術(shù)：

DNA測(cè)序在解讀物種生命的基因序列編碼中，起著最基礎(chǔ)和最廣泛的作用。早在發(fā)現(xiàn)DNA雙螺旋之初，就有人報(bào)道過DAN的測(cè)序技術(shù)，只是流程過于復(fù)雜。隨后不久，在1977年Sanger發(fā)明了末端終止測(cè)序法，具有里程碑的意義。迄今隨著生物信息技術(shù)科學(xué)的發(fā)展，Sanger測(cè)序法已經(jīng)無法滿足研究的需要，于是成本更低，通量更高，速度更快的第二代測(cè)序技術(shù)應(yīng)運(yùn)而生。它的核心思想是邊合成邊測(cè)序。DNA序列比對(duì)是將測(cè)序得到的DNA短序列(reads)與參考基因組進(jìn)行比對(duì)分析，常用以分析物種相似性和同源性，也可以通過計(jì)算機(jī)技術(shù)進(jìn)行基因信息的挖掘和分析。

目前，針對(duì)二代測(cè)序序列比對(duì)的技術(shù)主要有兩個(gè)方向的發(fā)展：一是基于設(shè)計(jì)hash算法映射構(gòu)建hash表庫(kù)思路，這種技術(shù)優(yōu)點(diǎn)比對(duì)速度快，準(zhǔn)確率極高，但占用內(nèi)存和CUP較高。該技術(shù)實(shí)現(xiàn)的軟件代表有SOAP，MAQ等。二是基于BWT轉(zhuǎn)換算法構(gòu)建后綴樹查詢數(shù)據(jù)結(jié)構(gòu)的技術(shù)，該技術(shù)先對(duì)基因組序列循環(huán)移位，然后排序，利用BWT壓縮并建立索引。在比對(duì)時(shí)，利用查找和回溯來定位序列的位置信息。該技術(shù)需要事先構(gòu)建索引文件分步驟操作，比對(duì)速度也不及接近常數(shù)平均時(shí)間的hash表結(jié)構(gòu)。它通常能夠?qū)崿F(xiàn)所有測(cè)序數(shù)據(jù)的常規(guī)比對(duì)，但相對(duì)hash數(shù)據(jù)結(jié)構(gòu)，樹結(jié)構(gòu)要占用龐大的內(nèi)存資源，并且查詢比對(duì)的速度遠(yuǎn)不及平均常數(shù)時(shí)間的hash表庫(kù)。

技術(shù)實(shí)現(xiàn)要素：

有鑒于此，本發(fā)明提供一種二代測(cè)序短序列快速比對(duì)分析方法及裝置，用以解決測(cè)序數(shù)據(jù)的比對(duì)效率低以及內(nèi)存占用高的問題。

一方面，本發(fā)明實(shí)施例提出一種二代測(cè)序短序列快速比對(duì)分析方法，包括：

S1、獲取測(cè)序得到的DNA短序列，并采用第一hash算法和第二hash算法分別映射編碼所述DNA短序列，分別得到第一索引和第二索引；

S2、基于預(yù)設(shè)的index查詢庫(kù)、所述第一索引和第二索引將所述DNA短序列和參考基因組進(jìn)行比對(duì)，其中，所述index查詢庫(kù)由單元結(jié)構(gòu)體數(shù)組構(gòu)成，每個(gè)單元結(jié)構(gòu)體包含有value值和index2值，所述value值為對(duì)應(yīng)的組成所述參考基因組的K-mer片段的位置信息摘要的壓縮，所述index2為所述K-mer片段的第二hash算法的映射結(jié)果，存儲(chǔ)每個(gè)所述單元結(jié)構(gòu)體數(shù)組索引偏移量為對(duì)應(yīng)的K-mer片段的第一hash算法的映射結(jié)果index1，K為片段序列長(zhǎng)度；

S3、根據(jù)比對(duì)的結(jié)果，若比對(duì)上，則獲取與對(duì)應(yīng)的DNA短序列比對(duì)上的K-mer片段的value值，并根據(jù)所述與對(duì)應(yīng)的DNA短序列比對(duì)上的K-mer片段的value值確定出所述對(duì)應(yīng)的DNA短序列所在染色體號(hào)和所在染色體上的位點(diǎn)。

優(yōu)選地，所述第一hash算法為XDDHash算法，所述第二hash算法為MWFHash算法。

優(yōu)選地，所述S2，包括：

對(duì)于每一個(gè)DNA短序列，查找首地址為該DNA短序列的第一索引的內(nèi)存塊，若查找到首地址為該DNA短序列的第一索引的內(nèi)存塊，則判斷該內(nèi)存塊的數(shù)據(jù)存儲(chǔ)位置是否為1個(gè)；

若為1個(gè)，則確定出該DNA短序列與該內(nèi)存塊的數(shù)據(jù)存儲(chǔ)位置處存儲(chǔ)的結(jié)構(gòu)體對(duì)應(yīng)的K-mer片段比對(duì)上。

優(yōu)選地，所述方法還包括：

若不為1個(gè)，則從所述首地址開始，采用移位線性探測(cè)的方法獲取當(dāng)前地址處存儲(chǔ)的結(jié)構(gòu)體的index2值，判斷該index2值是否與該DNA短序列的第二索引相同；

若該index2值與該DNA短序列的第二索引相同，則確定出該DNA短序列與該當(dāng)前地址處存儲(chǔ)的結(jié)構(gòu)體對(duì)應(yīng)的K-mer片段比對(duì)上。

優(yōu)選地，所述index查詢庫(kù)的裝填因子為0.607。

優(yōu)選地，在所述S2之前，所述方法還包括：

S40、依據(jù)先驗(yàn)數(shù)據(jù)模型，對(duì)所述參考基因組的24條染色體按照歷史比對(duì)上的DNA序列的數(shù)量從多到少進(jìn)行排序；

S41、根據(jù)排序的結(jié)果，按照從前到后的順序，對(duì)當(dāng)前染色體進(jìn)行滑窗式截取，得到連續(xù)的K-mer片段，并標(biāo)記各片段所在的位置信息；

S42、去除包含簡(jiǎn)并堿基的K-mer片段和重復(fù)的K-mer片段；

S43、對(duì)于得到的每一個(gè)K-mer片段，對(duì)該K-mer片段分別進(jìn)行第一hash算法映射和第二hash算法映射，分別得到index1值和index2值，并計(jì)算該K-mer片段的位置信息摘要的壓縮value值，生成包含所述value值和index2值的結(jié)構(gòu)體；

S44、將index1值作為索引構(gòu)建所述index查詢庫(kù)。

優(yōu)選地，在所述DNA短序列和參考基因組進(jìn)行比對(duì)時(shí)，采用多線程技術(shù)實(shí)現(xiàn)正鏈和反鏈短序列的并行比對(duì)。

另一方面，本發(fā)明實(shí)施例提出一種二代測(cè)序短序列快速比對(duì)分析裝置，包括：

映射單元，用于獲取測(cè)序得到的DNA短序列，并采用第一hash算法和第二hash算法分別映射編碼所述DNA短序列，分別得到第一索引和第二索引；

比對(duì)單元，用于基于預(yù)設(shè)的index查詢庫(kù)、所述第一索引和第二索引將所述DNA短序列和參考基因組進(jìn)行比對(duì)，其中，所述index查詢庫(kù)由單元結(jié)構(gòu)體數(shù)組構(gòu)成，每個(gè)單元結(jié)構(gòu)體包含有value值和index2值，所述value值為對(duì)應(yīng)的組成所述參考基因組的K-mer片段的位置信息摘要的壓縮，所述index2為所述K-mer片段的第二hash算法的映射結(jié)果，存儲(chǔ)每個(gè)所述單元結(jié)構(gòu)體的數(shù)組索引偏移量為對(duì)應(yīng)的K-mer片段的第一hash算法的映射結(jié)果index1，K為片段序列長(zhǎng)度；

計(jì)算單元，用于根據(jù)比對(duì)的結(jié)果，若比對(duì)上，則獲取與對(duì)應(yīng)的DNA短序列比對(duì)上的K-mer片段的value值，并根據(jù)所述與對(duì)應(yīng)的DNA短序列比對(duì)上的K-mer片段的value值確定出所述對(duì)應(yīng)的DNA短序列所在染色體號(hào)和所在染色體上的位點(diǎn)。

本發(fā)明具有如下有益效果：

1、采用了hash算法，比對(duì)速度快，接近常數(shù)時(shí)間O(1)，提高項(xiàng)目產(chǎn)品時(shí)間效率；

2、采用二次hash映射技術(shù)，無需直接保存reads序列，節(jié)省大量?jī)?nèi)存資源，提高項(xiàng)目產(chǎn)品資源效益；

3、采用先驗(yàn)數(shù)據(jù)模型重排hash索引，預(yù)先比對(duì)出現(xiàn)頻率高的reads序列，加速比對(duì)；

4、采用多線技術(shù)比對(duì)并行執(zhí)行比對(duì)任務(wù)，縮短了比對(duì)時(shí)間。

附圖說明

圖1為本發(fā)明二代測(cè)序短序列快速比對(duì)分析方法一實(shí)施例的流程示意圖；

圖2為本發(fā)明用來做一次hash函數(shù)映射的XDDHash算法的時(shí)間性能圖；

圖3為本發(fā)明用來做一次hash函數(shù)映射的XDDHash算法的沖突量性能圖；

圖4為本發(fā)明用來做二次hash函數(shù)映射的MWFHash算法的時(shí)間性能圖；

圖5為本發(fā)明用來做二次hash函數(shù)映射的MWFHash算法的沖突量性能圖；

圖6為本發(fā)明構(gòu)建的hash表的裝填因子和平均探測(cè)次數(shù)的關(guān)系圖；

圖7為本發(fā)明方法在設(shè)計(jì)成軟件程序產(chǎn)品后，比對(duì)樣本的數(shù)據(jù)量和花費(fèi)時(shí)間的關(guān)系圖；

圖8為本發(fā)明二代測(cè)序短序列快速比對(duì)分析裝置一實(shí)施例的結(jié)構(gòu)示意圖。

具體實(shí)施方式

為使本發(fā)明實(shí)施例的目的、技術(shù)方案和優(yōu)點(diǎn)更加清楚，下面將結(jié)合本發(fā)明實(shí)施例中的附圖，對(duì)本發(fā)明實(shí)施例中的技術(shù)方案進(jìn)行清楚地描述，顯然，所描述的實(shí)施例是本發(fā)明一部分實(shí)施例，而不是全部的實(shí)施例?；诒景l(fā)明中的實(shí)施例，本領(lǐng)域普通技術(shù)人員在沒有做出創(chuàng)造性勞動(dòng)前提下所獲得的所有其他實(shí)施例，都屬于本發(fā)明保護(hù)的范圍。

參看圖1，本實(shí)施例公開一種二代測(cè)序短序列快速比對(duì)分析方法，包括：

S1、獲取測(cè)序得到的DNA短序列，并采用第一hash算法和第二hash算法分別映射編碼所述DNA短序列，分別得到第一索引和第二索引；

S2、基于預(yù)設(shè)的index查詢庫(kù)、所述第一索引和第二索引將所述DNA短序列和參考基因組進(jìn)行比對(duì)，其中，所述index查詢庫(kù)由單元結(jié)構(gòu)體數(shù)組構(gòu)成，每個(gè)單元結(jié)構(gòu)體包含有value值和index2值，所述value值為對(duì)應(yīng)的組成所述參考基因組的K-mer片段的位置信息摘要的壓縮，所述index2為所述K-mer片段的第二hash算法的映射結(jié)果，存儲(chǔ)每個(gè)所述單元結(jié)構(gòu)體的數(shù)組索引偏移量為對(duì)應(yīng)的K-mer片段的第一hash算法的映射結(jié)果index1，K為片段序列長(zhǎng)度；

S3、根據(jù)比對(duì)的結(jié)果，若比對(duì)上，則獲取與對(duì)應(yīng)的DNA短序列比對(duì)上的K-mer片段的value值，并根據(jù)所述與對(duì)應(yīng)的DNA短序列比對(duì)上的K-mer片段的value值確定出所述對(duì)應(yīng)的DNA短序列所在染色體號(hào)和所在染色體上的位點(diǎn)。

查詢時(shí)序列位置信息的獲取如下：

內(nèi)存索引index獲取：通過第一hash算法和第二hash算法，將樣本中序列信息轉(zhuǎn)換為第一索引和第二索引，取得內(nèi)存中value值。

序列位置獲?。?/p>

DNA短序列所在染色體號(hào)K-mer_number＝value％26(求模運(yùn)算，即求余)；

DNA短序列所在染色體上位點(diǎn)K-mer_location＝value/26(求整除)。

本發(fā)明實(shí)施例中，采用二次hash算法，一次hash的映射結(jié)果用來查詢定位，二次hash的映射結(jié)果用來沖突篩選定位，能夠加快比對(duì)速度，且無須直接保存龐大的reads序列信息即可統(tǒng)計(jì)出比對(duì)結(jié)果，降低對(duì)內(nèi)存資源的浪費(fèi)。

可選地，在本發(fā)明二代測(cè)序短序列快速比對(duì)分析方法的另一實(shí)施例中，所述第一hash算法為XDDHash算法，所述第二hash算法為MWFHash算法。

一次映射的XDDhash壓縮算法，用來索引關(guān)鍵字信息摘要。其算法偽函數(shù)如下：

hash[i+1]＝(hash[i]*seed+str[i])^0x00000000FFFFFFFF(i＝0,1,2,…,N-2)

其中，hash[0]設(shè)定為63，乘法因子種子seed取33，str[i]為reads序列的第i+1個(gè)堿基字符的ASCⅡ值，N為待查詢比對(duì)reads序列的長(zhǎng)度。P^0x00000000FFFFFFFF表示對(duì)P取其低32位二進(jìn)制值。做循環(huán)遞歸運(yùn)算直到求得hash[N-1]即為第一索引。

經(jīng)理論和程序檢測(cè)的實(shí)踐驗(yàn)證，該算法具有較好的性能：沖突量少，壓縮映射速度快。其性能統(tǒng)計(jì)見說明書附圖圖2和3，其中圖2為用來做一次hash函數(shù)映射的XDDHash算法的時(shí)間性能圖，反映XDDHash算法隨著reads數(shù)量變化耗時(shí)關(guān)系圖，圖3為用來做一次hash函數(shù)映射的XDDHash算法的沖突量性能圖，反映XDDHash算法隨著reads數(shù)量變化產(chǎn)生的沖突量的關(guān)系圖。

二次映射的MWFhash壓縮算法，用來分離索引沖突值，對(duì)奇偶位采用不同的移位因子算法如下，

奇數(shù)位索引算法(i為奇數(shù)值時(shí))：

hash[i+1]＝(hash[i]*seedodd+i*str[i])^0x00000000FFFFFFFF

偶數(shù)位索引算法(i為偶數(shù)值時(shí))：

hash[i+1]＝(hash[i]*seedeven+i*str[i])^0x00000000FFFFFFFF

其中，hash[0]設(shè)定為0，奇數(shù)移位因子種子seedodd取值5，偶數(shù)移位因子種子seedeven取值7。做循環(huán)遞歸運(yùn)算得到hash[N-1]值,即為第二索引。

經(jīng)理論和程序檢測(cè)的實(shí)踐驗(yàn)證，該算法具有較好的性能：沖突量少，壓縮映射速度快。其性能統(tǒng)計(jì)見說明書附圖圖4和5，圖4為用來做二次hash函數(shù)映射的MWFHash算法的時(shí)間性能圖，反映MWFHash算法隨著reads數(shù)量變化所耗時(shí)的關(guān)系圖，圖5為用來做二次hash函數(shù)映射的MWFHash算法的沖突量性能圖，反映MWFHash算法隨著reads數(shù)量變化產(chǎn)生的沖突量的關(guān)系圖。

本發(fā)明實(shí)施例中，在hash算法編碼映射reads序列時(shí)，直接編碼取其低32位二進(jìn)制值簡(jiǎn)化了算法時(shí)間復(fù)雜度。

可選地，在本發(fā)明二代測(cè)序短序列快速比對(duì)分析方法的另一實(shí)施例中，所述S2，包括：

對(duì)于每一個(gè)DNA短序列，查找首地址為該DNA短序列的第一索引的內(nèi)存塊，若查找到首地址為該DNA短序列的第一索引的內(nèi)存塊，則判斷該內(nèi)存塊的數(shù)據(jù)存儲(chǔ)位置是否為1個(gè)；

若為1個(gè)，則確定出該DNA短序列與該內(nèi)存塊的數(shù)據(jù)存儲(chǔ)位置處存儲(chǔ)的結(jié)構(gòu)體對(duì)應(yīng)的K-mer片段比對(duì)上。

可選地，在本發(fā)明二代測(cè)序短序列快速比對(duì)分析方法的另一實(shí)施例中，還包括：

若不為1個(gè)，則從所述首地址開始，采用移位線性探測(cè)的方法獲取當(dāng)前地址處存儲(chǔ)的結(jié)構(gòu)體的index2值，判斷該index2值是否與該DNA短序列的第二索引相同；

若該index2值與該DNA短序列的第二索引相同，則確定出該DNA短序列與該當(dāng)前地址處存儲(chǔ)的結(jié)構(gòu)體對(duì)應(yīng)的K-mer片段比對(duì)上。

本發(fā)明實(shí)施例中，依據(jù)樣本短序列hash映射值的特點(diǎn)，利用改進(jìn)過的線性探測(cè)的方法解決hash值的沖突問題。具體是說，在索引向下移位探測(cè)時(shí)，只需要索引自加運(yùn)算即可，降低了比對(duì)搜索時(shí)的時(shí)間復(fù)雜度。

可選地，在本發(fā)明二代測(cè)序短序列快速比對(duì)分析方法的另一實(shí)施例中，所述index查詢庫(kù)的裝填因子為0.607。

所述index查詢庫(kù)實(shí)質(zhì)是hash表，Hash表裝填因子K的選?。?/p>

裝填因子：hash表中實(shí)際存放的記錄數(shù)與hash表大小的比值，是衡量hash表性能的一個(gè)重要指標(biāo)。其值越大，hash值沖突數(shù)越多，導(dǎo)致hash表的查詢速度降低。其值過小，浪費(fèi)內(nèi)存資源。我們基于XDDHash算法和MWFHash算法設(shè)計(jì)程序測(cè)試統(tǒng)計(jì)裝填因子和平均探測(cè)次數(shù)(次數(shù)越多，耗時(shí)越多)之間的關(guān)系圖見圖6所示，它反映所構(gòu)建表的查詢比對(duì)的效率性能，本發(fā)明設(shè)計(jì)的裝填因子K為0.607，可以觀察到其平均查找次數(shù)僅在3次左右，其接近常數(shù)時(shí)間查找。

可選地，在本發(fā)明二代測(cè)序短序列快速比對(duì)分析方法的另一實(shí)施例中，在所述S2之前，還包括：

S40、依據(jù)先驗(yàn)數(shù)據(jù)模型，對(duì)所述參考基因組的24條染色體按照歷史比對(duì)上的DNA序列的數(shù)量從多到少進(jìn)行排序；

S41、根據(jù)排序的結(jié)果，按照從前到后的順序，對(duì)當(dāng)前染色體進(jìn)行滑窗式截取，得到連續(xù)的K-mer片段，并標(biāo)記各片段所在的位置信息；

S42、去除包含簡(jiǎn)并堿基的K-mer片段和重復(fù)的K-mer片段；

在具體實(shí)施例中，包含簡(jiǎn)并堿基的K-mer片段在參考基因組hg19代表未知序列，構(gòu)建index庫(kù)時(shí)不可用，去掉后也可以節(jié)省內(nèi)存占用。

S43、對(duì)于得到的每一個(gè)K-mer片段，對(duì)該K-mer片段分別進(jìn)行第一hash算法映射和第二hash算法映射，分別得到index1值和index2值，并計(jì)算該K-mer片段的位置信息摘要的壓縮value值，生成包含所述value值和index2值的結(jié)構(gòu)體；

S44、將index1值作為索引構(gòu)建所述index查詢庫(kù)。

value值是K-mer片段序列的位置信息摘要的壓縮，壓縮函數(shù)如下：

Value＝K-mer_number*26+K-mer_location，

其中，K-mer_number為K-mer片段序列所在染色體號(hào)，K-mer_location為K-mer片段序列所在染色體上的位點(diǎn)。

Index2值和value值存儲(chǔ)于hash結(jié)構(gòu)體中：

Typedef struct Node{

unsigned int value；

unsigned int index2；

}HashType

這樣構(gòu)建的hash表的偽表達(dá)式可以簡(jiǎn)單寫為：

Hash[index1]＝HashType element；

就實(shí)現(xiàn)在內(nèi)存index1處存儲(chǔ)HashType類型的值。

本發(fā)明實(shí)施例中，依據(jù)先驗(yàn)數(shù)據(jù)模型，優(yōu)先排序出現(xiàn)頻率高的短序列，以提高查詢速度，這個(gè)具體方法是說，在處理有重復(fù)沖突的hash索引時(shí)，優(yōu)先查找實(shí)際測(cè)序樣本中出現(xiàn)率高的reads序列，減少索引移位搜索次數(shù)，避開不必要的時(shí)間開銷。

可選地，在本發(fā)明二代測(cè)序短序列快速比對(duì)分析方法的另一實(shí)施例中，在所述DNA短序列和參考基因組進(jìn)行比對(duì)時(shí)，采用多線程技術(shù)實(shí)現(xiàn)正鏈和反鏈短序列的并行比對(duì)。

在具體應(yīng)用中，可以實(shí)時(shí)監(jiān)視存儲(chǔ)測(cè)序得到的DNA序列的下機(jī)目錄，只要有新下機(jī)樣本數(shù)據(jù)，就進(jìn)行比對(duì)分析出結(jié)果數(shù)據(jù)。

實(shí)施例1基于二代測(cè)序短序列快速比對(duì)分析方法應(yīng)用于無創(chuàng)產(chǎn)前檢測(cè)的數(shù)據(jù)分析

本發(fā)明實(shí)現(xiàn)無創(chuàng)產(chǎn)前檢測(cè)的模塊流程首先是利用perl語言和shell命令預(yù)處理人類參考基因組hg19.fa，處理成index查詢庫(kù)所需的格式，包括三列：

第一列：36-mer DNA序列片段；

第二列：片段所在染色體號(hào)；

第三列：片段所在染色體上的位點(diǎn)。

然后其核心方法是基于C語言算法和數(shù)據(jù)結(jié)構(gòu)來實(shí)現(xiàn)。

本發(fā)明設(shè)計(jì)過程和運(yùn)行需要硬件和軟件環(huán)境：Linux系統(tǒng)；3個(gè)核心以上；35內(nèi)存以上；Linux平臺(tái)下C庫(kù)；Gcc編譯器；Gdb調(diào)試軟件。

輸入樣本文件格式：

Illumina測(cè)序平臺(tái)NIPT原始測(cè)序數(shù)據(jù)，為fastq格式，每4行為一條read信息：

第一行：以@開頭，標(biāo)示一條序列，包括測(cè)序儀名字，標(biāo)記位置，單/雙端測(cè)序以及過濾情況和引物接頭等相關(guān)信息描述；

第二行：reads序列；

第三行：以‘+’開頭，后面是序列標(biāo)示符、描述信息，或者什么也不加；

第四行：是reads質(zhì)量信息，和第二行的序列相對(duì)應(yīng)。

運(yùn)行以及輸出結(jié)果

(1)利用linux下命令nohup和&轉(zhuǎn)后臺(tái)，可以一直掃描下機(jī)目錄，只要有新批次樣本數(shù)據(jù)出現(xiàn)就可以出結(jié)果。運(yùn)行命令：

nohup./hashtb_pth.c&；

(2)輸入樣本文件數(shù)據(jù)：

我們將7.4M大小樣本文件Human_A1600169-R169_good_1.fq數(shù)據(jù)和Human_A1600171-R171_good_1.fq存入下機(jī)目錄。輸入文件：

Human_A1600169-R169_good_11.fq、Human_A1600169-R169_good_12.fq、Human_A1600169-R169_good_13.fq、Human_A1600169-R169_good_14.fq、Human_A1600169-R169_good_15.fq、Human_A1600169-R169_good_1.fq和Human_A1600171-R171_good_1.fq，其中，前6個(gè)是同一樣本數(shù)據(jù)的copy，用以檢測(cè)程序的穩(wěn)定性，最后一個(gè)用來檢測(cè)程序處理差異樣本的能力。

(3)比對(duì)結(jié)果

結(jié)果目錄中產(chǎn)生的文件：

Human_A1600169-R169_good_1、Human_A1600169-R169_good_11、Human_A1600169-R169_good_12、Human_A1600169-R169_good_13、Human_A1600169-R169_good_14、Human_A1600169-R169_good_15和Human_A1600171-R171_good_1，

樣本Human_A1600169-R169_good_1.fq的6個(gè)copy的結(jié)果數(shù)據(jù)完全一致，反映基于本發(fā)明方法實(shí)現(xiàn)的軟件在即使連續(xù)處理多個(gè)樣本數(shù)據(jù)時(shí)是穩(wěn)定可靠的。

由結(jié)果數(shù)據(jù)可以簡(jiǎn)單計(jì)算出：

樣本唯一比對(duì)率為：72.772％，

樣本比對(duì)率為：74.443％。

對(duì)差異樣本樣本Human_A1600171-R171_good_1.fq結(jié)果數(shù)據(jù)后4行的結(jié)果做簡(jiǎn)單計(jì)算：

樣本唯一比對(duì)率：71.531％，

樣本比對(duì)率：73.487％。

結(jié)果和Human_A1600169-R169_good_1.fq不同，說明軟件程序能夠即時(shí)連續(xù)處理同一批次目錄下的不同樣本數(shù)據(jù)。

(4)比對(duì)用時(shí)

處理7個(gè)樣本用時(shí)：77.850秒。

統(tǒng)計(jì)樣本reads量和比對(duì)時(shí)間之間的關(guān)系圖，見圖7，圖7為本發(fā)明方法在設(shè)計(jì)成軟件程序產(chǎn)品后，比對(duì)樣本的數(shù)據(jù)量和花費(fèi)時(shí)間的關(guān)系圖，這是由真實(shí)樣本實(shí)例測(cè)試得出，它綜合反映了本發(fā)明方法由C編程語言實(shí)現(xiàn)后的效果。它是采用無創(chuàng)產(chǎn)前基因檢測(cè)樣本的原始數(shù)據(jù)做的測(cè)試。

(5)結(jié)果驗(yàn)證和比較

對(duì)照結(jié)果我們采用主流可靠的BWA軟件來實(shí)現(xiàn)：

我們利用BWA對(duì)樣本Human_A1600169-R169_good_1.fq比對(duì)。流程如下：

1.利用參考基因組構(gòu)建index庫(kù)，運(yùn)行

bwa index-a bwtsw hg19.fa，

其中，hg19.fa為參考基因組，由于參考基因組hg19.fa大于2GB，-a參數(shù)采用bwtsw。

2.樣本序列和參考序列index庫(kù)比對(duì)，運(yùn)行

bwa aln-f Human_A1600169-R169_good_1.sai hg19.faHuman_A1600169-R169_good_1.fq

生成sai文件，包含有SA coordinates信息摘要。

3.將sai文件轉(zhuǎn)換為sam輸出，運(yùn)行

bwa samse-f Human_A1600169-R169_good_1.sam hg19.fa

Human_A1600169-R169_good_1.sai

Human_A1600169-R169_good_1.fastq，

其中，sam文件包含比對(duì)結(jié)果的各種信息。

4.利用perl語言提取比對(duì)結(jié)果數(shù)據(jù)并統(tǒng)計(jì)

樣本唯一比對(duì)率：72.772％，

樣本比對(duì)率：74.443％。

本軟件BWA和結(jié)果一致，但bwa整個(gè)流程需要半小時(shí)并占用大量資源，同時(shí)bwa不能即時(shí)連續(xù)處理多個(gè)樣本，單樣本處理也不能直接給出分析NIPT所需的統(tǒng)計(jì)結(jié)果，要另外寫腳本統(tǒng)計(jì)分析。

針對(duì)無創(chuàng)產(chǎn)前檢測(cè)項(xiàng)目二者比較(單個(gè)樣本)：

參看圖8，本實(shí)施例公開一種二代測(cè)序短序列快速比對(duì)分析裝置，包括：

映射單元1，用于獲取測(cè)序得到的DNA短序列，并采用第一hash算法和第二hash算法分別映射編碼所述DNA短序列，分別得到第一索引和第二索引；

比對(duì)單元，用于基于預(yù)設(shè)的index查詢庫(kù)、所述第一索引和第二索引將所述DNA短序列和參考基因組進(jìn)行比對(duì)，其中，所述index查詢庫(kù)由單元結(jié)構(gòu)體數(shù)組構(gòu)成，每個(gè)單元結(jié)構(gòu)體包含有value值和index2值，所述value值為對(duì)應(yīng)的組成所述參考基因組的K-mer片段的位置信息摘要的壓縮，所述index2為所述K-mer片段的第二hash算法的映射結(jié)果，存儲(chǔ)每個(gè)所述單元結(jié)構(gòu)體的數(shù)組索引偏移量為對(duì)應(yīng)的K-mer片段的第一hash算法的映射結(jié)果index1，K為片段序列長(zhǎng)度；

計(jì)算單元3，用于根據(jù)比對(duì)的結(jié)果，若比對(duì)上，則獲取與對(duì)應(yīng)的DNA短序列比對(duì)上的K-mer片段的value值，并根據(jù)所述與對(duì)應(yīng)的DNA短序列比對(duì)上的K-mer片段的value值確定出所述對(duì)應(yīng)的DNA短序列所在染色體號(hào)和所在染色體上的位點(diǎn)。

雖然結(jié)合附圖描述了本發(fā)明的實(shí)施方式，但是本領(lǐng)域技術(shù)人員可以在不脫離本發(fā)明的精神和范圍的情況下做出各種修改和變型，這樣的修改和變型均落入由所附權(quán)利要求所限定的范圍之內(nèi)。