本發(fā)明涉及生物和計(jì)算機(jī)領(lǐng)域，特別涉及一種通過計(jì)算機(jī)尋找樣本的染色體突變位點(diǎn)的分析方法和分析裝置。

背景技術(shù)：

基于高通量測(cè)序數(shù)據(jù)尋找和致病基因、癌癥治療、個(gè)性化用藥相關(guān)的染色體突變位點(diǎn)為臨床應(yīng)用提供了不可估量的前景，例如：癌癥基因組、無(wú)創(chuàng)產(chǎn)前診斷、藥物靶標(biāo)篩選。高通量測(cè)序技術(shù)能否在這些領(lǐng)域取得進(jìn)展，其關(guān)鍵在于分析方法及軟件的創(chuàng)新。

由于測(cè)序技術(shù)的不斷進(jìn)步，測(cè)序價(jià)格的不斷降低，我們得到的數(shù)據(jù)越來(lái)越多，如何快速的挖掘出需要的信息已經(jīng)是亟待解決的問題。目前可以進(jìn)行同樣分析的軟件主要有Samtools、GATK等。而這些軟件是單線程，對(duì)于大量數(shù)據(jù)的分析是很慢的，因此急需尋找一個(gè)新的方式來(lái)縮短分析的時(shí)間。

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

有鑒于此，本發(fā)明基于尋找基因突變的理論和分布式計(jì)算框架提供了一種分析速度更快的數(shù)據(jù)的軟件分析方法和分析裝置。

本發(fā)明的實(shí)施例提供了一種基于高通量測(cè)序數(shù)據(jù)尋找樣本的染色體突變位點(diǎn)的分析方法，所述方法包括：

根據(jù)測(cè)序序列和參考序列的比對(duì)結(jié)果為并行計(jì)算準(zhǔn)備數(shù)據(jù)；

通過并行計(jì)算過濾掉準(zhǔn)備好的數(shù)據(jù)中無(wú)效的堿基；

根據(jù)過濾后保留的堿基確定樣本的染色體突變位點(diǎn)。

所述測(cè)序序列和參考序列的比對(duì)結(jié)果為并行計(jì)算準(zhǔn)備數(shù)據(jù)包括：

對(duì)包含所述比對(duì)結(jié)果的數(shù)據(jù)進(jìn)行切分；

生成包含切分后的數(shù)據(jù)的執(zhí)行實(shí)體。

優(yōu)選地，所述對(duì)包含所述比對(duì)結(jié)果的數(shù)據(jù)進(jìn)行切分包括：

將所述比對(duì)結(jié)果的數(shù)據(jù)的文件切分為文件塊；

每隔預(yù)定時(shí)間間隔將從所述文件塊中讀取到的數(shù)據(jù)切分為組成DStream的RDD。

優(yōu)選地，所述生成包含切分后的數(shù)據(jù)的執(zhí)行實(shí)體包括根據(jù)RDD生成執(zhí)行實(shí)體job。

所述通過并行計(jì)算過濾掉準(zhǔn)備好的數(shù)據(jù)中無(wú)效的堿基包括：

通過執(zhí)行實(shí)體對(duì)準(zhǔn)備好的數(shù)據(jù)進(jìn)行并行計(jì)算：當(dāng)由樣本產(chǎn)生的測(cè)序序列與參考序列對(duì)應(yīng)位點(diǎn)上的堿基不一致時(shí)，對(duì)此測(cè)序序列上的每個(gè)堿基計(jì)算堿基比對(duì)質(zhì)量值；

通過執(zhí)行實(shí)體并行地在所述堿基比對(duì)質(zhì)量值和此堿基的測(cè)序質(zhì)量值中取較小值作為此堿基的最終質(zhì)量值；

通過執(zhí)行實(shí)體并行地進(jìn)行過濾：如果所述最終質(zhì)量值小于此堿基的第一閾值，則過濾掉此堿基。

所述根據(jù)過濾后保留的堿基確定樣本的染色體突變位點(diǎn)包括：

根據(jù)過濾后保留的堿基確定樣本的每個(gè)堿基位點(diǎn)上的突變堿基所占的比率；

如果所述樣本的某一位點(diǎn)上的突變堿基所占的比率大于等于第二閾值，則此位點(diǎn)為樣本的染色體突變位點(diǎn)。

本發(fā)明的實(shí)施例還提供了一種基于高通量測(cè)序數(shù)據(jù)尋找樣本的染色體突變位點(diǎn)的分析裝置，所述裝置包括：

并行準(zhǔn)備模塊，用于根據(jù)測(cè)序序列和參考序列的比對(duì)結(jié)果為并行計(jì)算準(zhǔn)備數(shù)據(jù)；

并行計(jì)算模塊，用于通過并行計(jì)算過濾掉準(zhǔn)備好的數(shù)據(jù)中無(wú)效的堿基；

結(jié)果合并模塊，用于根據(jù)過濾后保留的堿基確定樣本的染色體突變位點(diǎn)。

優(yōu)選地，所述并行準(zhǔn)備模塊包括：

數(shù)據(jù)切分單元，用于對(duì)包含所述比對(duì)結(jié)果的數(shù)據(jù)進(jìn)行切分；

執(zhí)行實(shí)體生成單元，用于生成包含切分后的數(shù)據(jù)的執(zhí)行實(shí)體。

所述數(shù)據(jù)切分單元，用于：

將包含所述比對(duì)結(jié)果的數(shù)據(jù)的文件切分為文件塊；

每隔預(yù)定時(shí)間間隔將從所述文件塊中讀取到的數(shù)據(jù)切分為組成DStream的RDD。

所述執(zhí)行實(shí)體生成單元用于根據(jù)RDD生成執(zhí)行實(shí)體job。

優(yōu)選地，所述并行計(jì)算模塊包括：

質(zhì)量計(jì)算單元，用于：

通過執(zhí)行實(shí)體對(duì)準(zhǔn)備好的數(shù)據(jù)進(jìn)行并行計(jì)算：當(dāng)由樣本產(chǎn)生的測(cè)序序列與參考序列對(duì)應(yīng)位點(diǎn)上的堿基不一致時(shí)，對(duì)此測(cè)序序列上的每個(gè)堿基計(jì)算堿基比對(duì)質(zhì)量值，

通過執(zhí)行實(shí)體并行地在所述堿基比對(duì)質(zhì)量值和此堿基的測(cè)序質(zhì)量值中取較小值作為最終質(zhì)量值；

過濾單元：用于通過執(zhí)行實(shí)體并行地進(jìn)行過濾：當(dāng)所述最終質(zhì)量值小于此堿基的第一閾值時(shí)，過濾掉此堿基。

優(yōu)選地，所述結(jié)果合并模塊包括：

比率計(jì)算單元，用于根據(jù)過濾后保留的堿基確定樣本的每個(gè)堿基位點(diǎn)上的突變堿基所占的比率；

突變位點(diǎn)確定單元，用于當(dāng)所述樣本的某一位點(diǎn)上的突變堿基所占的比率大于等于第二閾值時(shí)，確定此位點(diǎn)為樣本的染色體突變位點(diǎn)。

采用了本發(fā)明的技術(shù)方案以后，可以采用多機(jī)器彈性并行計(jì)算來(lái)加快分析速度。Spark采用內(nèi)存迭代計(jì)算，從緩存中存入和讀取的時(shí)間很短，幾乎可以忽略，而另一種可選方案需要向HDFS讀寫數(shù)據(jù)則需要耗費(fèi)大量時(shí)間。本發(fā)明的技術(shù)方案與其他軟件的解決方案相比，分析速度大幅提升。

附圖說明

圖1為本發(fā)明第一實(shí)施例提供的基于高通量測(cè)序數(shù)據(jù)尋找樣本的染色體突變位點(diǎn)的分析方法的示意圖；

圖2為本發(fā)明第二實(shí)施例提供的基于高通量測(cè)序數(shù)據(jù)尋找樣本的染色體突變位點(diǎn)的分析方法的示意圖；

圖3為本發(fā)明第三實(shí)施例提供的尋找樣本的染色體突變位點(diǎn)的分析方法的示意圖；

圖4為隱馬爾可夫模型(HMM)的示意圖；

圖5為本發(fā)明實(shí)施例提供的基于高通量測(cè)序數(shù)據(jù)尋找樣本的染色體突變位點(diǎn)的分析裝置的結(jié)構(gòu)框圖。

具體實(shí)施方式

本發(fā)明提供了一種基于高通量測(cè)序數(shù)據(jù)尋找樣本的染色體突變位點(diǎn)的分析方法和分析裝置。下面三個(gè)實(shí)施例對(duì)本發(fā)明進(jìn)行說明。

請(qǐng)參閱圖1，在本發(fā)明第一實(shí)施例中，首先根據(jù)堿基比對(duì)結(jié)果為并行計(jì)算進(jìn)行準(zhǔn)備，然后通過并行計(jì)算在準(zhǔn)備好的數(shù)據(jù)中過濾掉無(wú)效的堿基，最后根據(jù)過濾后保留的堿基確定樣本的染色體突變位點(diǎn)。

S101、根據(jù)堿基比對(duì)結(jié)果為并行計(jì)算準(zhǔn)備數(shù)據(jù)。

將包含測(cè)序序列和參考序列堿基比對(duì)結(jié)果的文件切分成多個(gè)文件塊，保存在文件系統(tǒng)中。

S102、通過并行計(jì)算過濾掉準(zhǔn)備好的數(shù)據(jù)中無(wú)效的堿基。

多臺(tái)計(jì)算機(jī)同時(shí)從文件系統(tǒng)中的多個(gè)文件塊中讀取數(shù)據(jù)，通過并行計(jì)算過濾掉其中無(wú)效的堿基。

S103、根據(jù)過濾后保留的堿基確定樣本的染色體突變位點(diǎn)。

請(qǐng)參閱圖2，在本發(fā)明的第二實(shí)施例中，通過使用hadoop框架進(jìn)行并行計(jì)算尋找樣本的染色體突變位點(diǎn)：首先將包含堿基比對(duì)結(jié)果的大文件切分為多個(gè)小數(shù)據(jù)塊保存到HDFS中，再利用多臺(tái)計(jì)算機(jī)并行地從所述多個(gè)小數(shù)據(jù)塊中讀取數(shù)據(jù)，通過并行計(jì)算過濾掉無(wú)效的堿基。將多臺(tái)計(jì)算機(jī)上并行計(jì)算的結(jié)果合并，確定樣本的染色體突變位點(diǎn)。

具體地，在本發(fā)明的第二實(shí)施例中：

第一閾值是指當(dāng)測(cè)序序列某一堿基的最終質(zhì)量值小于此值時(shí)，認(rèn)為此堿基無(wú)效，并將其過濾掉，不作為尋找樣本的染色體突變位點(diǎn)的依據(jù)；

第二閾值是指對(duì)于樣本的某一確定位點(diǎn)，如果其突變堿基所占的比率大于等于此值，則確定此位點(diǎn)為樣本的染色體突變位點(diǎn)。

S201、將包含堿基比對(duì)結(jié)果的大數(shù)據(jù)文件切分為多個(gè)小數(shù)據(jù)塊。

從數(shù)據(jù)源輸入的文件為bam或者sam文件，其中包含了測(cè)序序列與參考序列堿基的比對(duì)結(jié)果的大量數(shù)據(jù)。

多個(gè)客戶端通過網(wǎng)絡(luò)將原始的bam或sam大文件上傳到HDFS(HadoopDistributed File System)，HDFS將此大文件根據(jù)一定的規(guī)則切分成小數(shù)據(jù)塊并保存。切分規(guī)則如下：比如輸入的大文件為3G，在spark streaming中將一個(gè)block設(shè)置為128M，那么輸入的大文件總共將被切分為3*1024/128＝24塊。

S202、Hadoop讀取HDFS中的小數(shù)據(jù)塊，啟動(dòng)job，job為每一個(gè)小文件塊生成實(shí)際的執(zhí)行實(shí)體map task和reduce task。

參與并行計(jì)算的多臺(tái)計(jì)算機(jī)上的Hadoop并行地讀取HDFS中的小數(shù)據(jù)塊，并啟動(dòng)一個(gè)job，job為每一個(gè)小數(shù)據(jù)塊生成一個(gè)map task和一個(gè)reduce task。Map task和reduce task是并行計(jì)算的執(zhí)行實(shí)體。

S203、在map task中，通過并行計(jì)算堿基比對(duì)質(zhì)量值(BAQ)和最終質(zhì)量值來(lái)過濾掉測(cè)序序列位點(diǎn)上無(wú)效的堿基，并將過濾后保留下來(lái)的堿基保存到HDFS中。

Hadoop在參與并行計(jì)算的多臺(tái)計(jì)算機(jī)上并行地執(zhí)行步驟S202中生成的map task和reduce task，先運(yùn)行map task，然后運(yùn)行reduce task。

為了尋找樣本的染色體突變位點(diǎn)，需要通過計(jì)算堿基比對(duì)質(zhì)量值(BAQ)和最終質(zhì)量值過濾掉測(cè)序序列上的無(wú)效堿基。這一部分并行計(jì)算是在map task中進(jìn)行的。有關(guān)計(jì)算堿基比對(duì)質(zhì)量值BAQ和最終質(zhì)量值以及過濾掉測(cè)序序列位點(diǎn)上無(wú)效堿基的算法描述請(qǐng)參見本說明書步驟S303至S306的相關(guān)描述。

S204、在reduce task中，根據(jù)保留下來(lái)的堿基通過并行計(jì)算確定樣本的每一個(gè)位點(diǎn)上的突變堿基所占的比率，并根據(jù)此比率確定樣本的染色體突變位點(diǎn)。

在reduce task中，首先從HDFS中讀取通過步驟S203過濾后保留下來(lái)的堿基。然后，對(duì)于樣本上的每一個(gè)位點(diǎn)，累加由測(cè)序數(shù)據(jù)產(chǎn)生的ATCG各堿基的次數(shù)(通過多臺(tái)計(jì)算機(jī)對(duì)其上的小數(shù)據(jù)塊進(jìn)行并行計(jì)算得出)，從而得到此位點(diǎn)上的突變堿基所占的比率。

如果樣本某一位點(diǎn)上的突變堿基所占的比率大于等于第二閾值，則確定此位點(diǎn)是樣本的染色體突變位點(diǎn)；否則確定此位點(diǎn)不是樣本的染色體突變位點(diǎn)。通過這種方法確定樣本的所有染色體突變位點(diǎn)。

請(qǐng)參閱圖3，在本發(fā)明的第三實(shí)施例中，通過使用并行計(jì)算框架spark的分析方法基于高通量測(cè)序數(shù)據(jù)尋找樣本的染色體突變位點(diǎn)：首先將包含堿基比對(duì)結(jié)果的大文件切分為多個(gè)小文件塊，然后spark s treaming每隔預(yù)定時(shí)間間隔將從所述文件塊中讀取到的數(shù)據(jù)切分為RDD，并將其封裝生成job。在多臺(tái)計(jì)算機(jī)并行地執(zhí)行多個(gè)job的過程中，通過計(jì)算堿基比對(duì)質(zhì)量值(BAQ)過濾掉無(wú)效的堿基。最后將多臺(tái)計(jì)算機(jī)針對(duì)樣本上同一位點(diǎn)的并行計(jì)算結(jié)果進(jìn)行合并，確定樣本的染色體突變位點(diǎn)。

具體地，在本發(fā)明的第三實(shí)施例中：

第一閾值是指當(dāng)測(cè)序序列堿基的最終質(zhì)量值小于此值時(shí)，認(rèn)為此堿基無(wú)效，并將其過濾掉，不作為尋找樣本的染色體突變位點(diǎn)的依據(jù)；

第二閾值是指對(duì)于樣本的某一確定位點(diǎn)，如果突變堿基所占的比率大于等于此值，則確定此位點(diǎn)為樣本的染色體突變位點(diǎn)。

Spark：是UC Berkeley AMPLab開發(fā)的一種計(jì)算框架。

Spark streaming：是建立在Spark上的實(shí)時(shí)計(jì)算框架，通過它提供豐富的API，基于內(nèi)存的高速執(zhí)行引擎，用戶可以結(jié)合流式、批處理和交互式查詢應(yīng)用。

DStream是指離散流(Discretized Stream)，它是Spark Streaming提供的基本的抽象，代表一個(gè)實(shí)時(shí)的數(shù)據(jù)流。

RDD是指彈性分布式數(shù)據(jù)集(Resilient Distributed Datasets)，它是可容錯(cuò)的并行數(shù)據(jù)結(jié)構(gòu)，使用戶能夠顯式地在內(nèi)存中保存中間的運(yùn)算結(jié)果，通過控制RDD的分區(qū)來(lái)優(yōu)化數(shù)據(jù)的布局，并使用豐富的轉(zhuǎn)換算子進(jìn)行操作。

Job是指在spark streaming中，封裝了RDD用于進(jìn)行并行計(jì)算的執(zhí)行實(shí)體。

S301、將包含堿基比對(duì)結(jié)果的大數(shù)據(jù)文件切分為多個(gè)小文件塊。

從數(shù)據(jù)源輸入的文件為bam或者sam文件，其中包含了測(cè)序序列與參考序列比對(duì)結(jié)果的大量數(shù)據(jù)。

多個(gè)客戶端通過網(wǎng)絡(luò)上傳原始的bam或sam大文件，spark streaming將此大文件根據(jù)一定的規(guī)則切分成小文件塊保存到kafka中。切分規(guī)則如下：比如輸入的大文件為3G，在spark streaming中將一個(gè)block設(shè)置為128M，那么輸入的大文件總共將被切分為3*1024/128＝24塊。

S302、參與并行計(jì)算的多臺(tái)計(jì)算機(jī)并行地將從多個(gè)小文件塊中讀取到的數(shù)據(jù)切分為RDD，并根據(jù)RDD生成并行計(jì)算的執(zhí)行實(shí)體job。

參與并行計(jì)算的多臺(tái)計(jì)算機(jī)上的spark streaming并行地以預(yù)定時(shí)間間隔(比如3秒)從kafka中保存的多個(gè)小文件塊中讀取數(shù)據(jù)，并在此過程中將讀取到的數(shù)據(jù)切分為RDD，供spark streaming后續(xù)的并行計(jì)算使用。參與并行計(jì)算的所有計(jì)算機(jī)按時(shí)間片(即預(yù)定時(shí)間間隔)讀取到的所有數(shù)據(jù)組成了代表持續(xù)輸入的堿基比對(duì)結(jié)果的DStream。Spark streaming將對(duì)DStream的操作轉(zhuǎn)換為DStream的有向無(wú)環(huán)圖(DAG)。對(duì)于每個(gè)時(shí)間片(即預(yù)定時(shí)間間隔)，spark streaming根據(jù)DStream DAG生成RDD DAG,并根據(jù)RDD DAG生成相應(yīng)的執(zhí)行實(shí)體job。后續(xù)的步驟S303至S306都是在這里的多個(gè)job中進(jìn)行并行計(jì)算的。

S303、確定需要計(jì)算堿基比對(duì)質(zhì)量值(BAQ)的堿基。

Spark streaming在參與并行計(jì)算的多臺(tái)計(jì)算機(jī)上并行地執(zhí)行步驟S302中生成的多個(gè)job。

為了尋找樣本的染色體突變位點(diǎn)，需要過濾掉測(cè)序序列一部分位點(diǎn)上的堿基：測(cè)序數(shù)據(jù)是基于高通量測(cè)序獲取的，而高通量測(cè)序的錯(cuò)誤率在0.01％左右，所以需要去掉低質(zhì)量的數(shù)據(jù)(即無(wú)效的堿基)，不將其作為尋找樣本的染色體突變位點(diǎn)的依據(jù)。

過濾過程首先需要判斷是否需要對(duì)堿基計(jì)算BAQ值。確定方法如下：如果測(cè)序序列與參考序列所有對(duì)應(yīng)位點(diǎn)上的堿基完全一致，則將此測(cè)序序列所有位點(diǎn)的堿基作為尋找樣本的染色體突變位點(diǎn)的依據(jù)，無(wú)需計(jì)算它們的BAQ值；如果測(cè)序序列某一位點(diǎn)與參考序列對(duì)應(yīng)位點(diǎn)上的堿基不一致，則需要對(duì)此測(cè)序序列所有位點(diǎn)的堿基計(jì)算BAQ值。

S304、測(cè)序序列與參考序列堿基不一致時(shí)計(jì)算測(cè)序序列所有位點(diǎn)的堿基比對(duì)質(zhì)量值(BAQ)。

請(qǐng)參閱圖4，如果測(cè)序序列(read序列)和參考序列(reference序列)不完全一致，則對(duì)read序列和reference序列建立一個(gè)隱馬爾可夫模型(HMM)，并對(duì)read序列中狀態(tài)為M的堿基計(jì)算其堿基比對(duì)質(zhì)量值(BAQ)，該值代表此堿基比對(duì)結(jié)果錯(cuò)誤的可能性。

在HMM中，M代表堿基匹配(包括完全一致和錯(cuò)配兩種情況)，I代表插入，D代表缺失，S代表開始，E代表結(jié)束。S指向所有的I和M，所有的M和I指向E。

BAQ算法分為前向遞推和后向遞推。算法中用到的變量如下：

前向遞推：

(a)值的初始化，其中k＝1,…,L；i＝0,1。

f_0,S＝1

(b)遞推，其中i＝2,…,l；k＝1,…,L。

對(duì)于i＝l+1，有：

其中，a_ij由以下狀態(tài)轉(zhuǎn)換矩陣確定：

后向遞推：

(a)值的初始化，其中k＝1,…,L。

b_l+1,E＝1

(b)遞推，其中i＝l-1,…,1；k＝L,…,1。

當(dāng)i＝1時(shí)，δ_i1＝1，其他則為0。

對(duì)于i＝0，有：

BAQ通過下式得到：

其中，b_ij由以下狀態(tài)轉(zhuǎn)換矩陣確定：

S305、確定已計(jì)算過BAQ的堿基的最終質(zhì)量值。

對(duì)于已計(jì)算過BAQ的堿基，在BAQ和此堿基的測(cè)序質(zhì)量值(測(cè)序質(zhì)量值可自bam或sam文件中獲得)中取較小值作為最終質(zhì)量值。

S306、過濾掉無(wú)效的堿基。

在測(cè)序序列所有的堿基中，對(duì)于已計(jì)算過BAQ的堿基，如果其最終質(zhì)量值小于第一閾值，則認(rèn)為此堿基無(wú)效，并將其過濾掉，不作為尋找樣本的染色體突變位點(diǎn)的依據(jù)；否則認(rèn)為此堿基有效，將其保留并作為尋找樣本的染色體突變位點(diǎn)的依據(jù)。

S307、根據(jù)保留下來(lái)的堿基計(jì)算樣本的每個(gè)位點(diǎn)上的突變堿基所占的比率。

對(duì)于樣本上的每一個(gè)位點(diǎn)，將所有測(cè)序序列在該位點(diǎn)產(chǎn)生的ATCG各堿基的次數(shù)(通過對(duì)多臺(tái)計(jì)算機(jī)上的多個(gè)RDD進(jìn)行并行計(jì)算得出)累加，得到該位點(diǎn)上的突變堿基所占的比率。

S308、確定樣本的染色體突變位點(diǎn)。

如果樣本某一位點(diǎn)上的突變堿基所占的比率大于等于第二閾值，則確定此位點(diǎn)是樣本的染色體突變位點(diǎn)；否則確定此位點(diǎn)不是樣本的染色體突變位點(diǎn)。通過這種方法確定樣本所有的染色體突變位點(diǎn)。

請(qǐng)參閱圖5，本發(fā)明還提供了一種尋找樣本的染色體突變位點(diǎn)的分析裝置，所述裝置包括：

并行準(zhǔn)備模塊510，用于根據(jù)測(cè)序序列和參考序列的比對(duì)結(jié)果為并行計(jì)算準(zhǔn)備數(shù)據(jù)；

并行計(jì)算模塊520，用于通過并行計(jì)算過濾掉準(zhǔn)備好的數(shù)據(jù)中無(wú)效的堿基；

結(jié)果合并模塊530，用于根據(jù)過濾后保留的堿基確定樣本的染色體突變位點(diǎn)。

所述并行準(zhǔn)備模塊510包括：

數(shù)據(jù)切分單元5101，用于對(duì)包含所述比對(duì)結(jié)果的數(shù)據(jù)進(jìn)行切分；

執(zhí)行實(shí)體生成單元5102，用于生成包含切分后的數(shù)據(jù)的執(zhí)行實(shí)體。

優(yōu)選地，所述數(shù)據(jù)切分單元5101用于：

將包含所述比對(duì)結(jié)果的數(shù)據(jù)的文件切分為文件塊；

每隔預(yù)定時(shí)間間隔將從所述文件塊中讀取到的數(shù)據(jù)切分為組成DStream的RDD。

優(yōu)選地，所述執(zhí)行實(shí)體生成單元5102用于根據(jù)RDD生成執(zhí)行實(shí)體job。

所述并行計(jì)算模塊520包括：

質(zhì)量計(jì)算單元5201，用于：

通過執(zhí)行實(shí)體對(duì)準(zhǔn)備好的數(shù)據(jù)進(jìn)行并行計(jì)算：當(dāng)由樣本產(chǎn)生的測(cè)序序列與參考序列對(duì)應(yīng)位點(diǎn)上的堿基不一致時(shí)，對(duì)此測(cè)序序列上的每個(gè)堿基計(jì)算堿基比對(duì)質(zhì)量值，

通過執(zhí)行實(shí)體并行地在所述堿基比對(duì)質(zhì)量值和此堿基的測(cè)序質(zhì)量值中取較小值作為此堿基的最終質(zhì)量值；

過濾單元5202，用于通過執(zhí)行實(shí)體并行地進(jìn)行過濾：當(dāng)所述最終質(zhì)量值小于此堿基的第一閾值時(shí)，過濾掉此堿基。

所述結(jié)果合并模塊530包括：

比率計(jì)算單元5301，用于根據(jù)過濾后保留的堿基確定樣本的每個(gè)堿基位點(diǎn)上的突變堿基所占的比率；

突變位點(diǎn)確定單元5302，用于當(dāng)所述樣本的某一位點(diǎn)上的突變堿基所占的比率大于等于第二閾值時(shí)，確定此位點(diǎn)為樣本的染色體突變位點(diǎn)。

在本發(fā)明的第一實(shí)施例中：

并行準(zhǔn)備模塊510根據(jù)堿基的比對(duì)結(jié)果為并行計(jì)算準(zhǔn)備數(shù)據(jù)。并行準(zhǔn)備模塊510將包含測(cè)序序列和參考序列堿基的比對(duì)結(jié)果的文件切分成多個(gè)文件塊，保存在文件系統(tǒng)中。并行準(zhǔn)備模塊510將從文件系統(tǒng)中的多個(gè)文件塊中讀取數(shù)據(jù)為并行計(jì)算進(jìn)行準(zhǔn)備。并行計(jì)算模塊520通過并行計(jì)算過濾掉準(zhǔn)備好的數(shù)據(jù)中無(wú)效的堿基。結(jié)果合并模塊530根據(jù)過濾后保留的堿基確定樣本的染色體突變位點(diǎn)。

在本發(fā)明的第二實(shí)施例中：

并行準(zhǔn)備模塊510中的數(shù)據(jù)切分單元5101將包含堿基的比對(duì)結(jié)果的大數(shù)據(jù)文件切分為多個(gè)小數(shù)據(jù)塊。

從數(shù)據(jù)源輸入的文件為bam或者sam文件，其中包含了測(cè)序序列與參考序列堿基比對(duì)結(jié)果的大量數(shù)據(jù)。

多個(gè)客戶端通過網(wǎng)絡(luò)將原始的bam或sam大文件上傳到HDFS(Hadoop Distributed File System)，并行準(zhǔn)備模塊510中的數(shù)據(jù)切分單元5101將此大文件根據(jù)一定的規(guī)則切分成小數(shù)據(jù)塊并保存。切分規(guī)則如下：比如輸入的大文件為3G，在spark streaming中將一個(gè)block設(shè)置為128M，那么輸入的大文件總共將被切分為3*1024/128＝24塊。

并行準(zhǔn)備模塊510中的執(zhí)行實(shí)體生成單元5102讀取HDFS中的小數(shù)據(jù)塊，啟動(dòng)job，job為每一個(gè)小文件塊生成實(shí)際的執(zhí)行實(shí)體map task和reduce task。

并行準(zhǔn)備模塊510中的執(zhí)行實(shí)體生成單元5102讀取HDFS中的小數(shù)據(jù)塊，并啟動(dòng)一個(gè)job，job為每一個(gè)小數(shù)據(jù)塊生成一個(gè)map task和reduce task。Map task和reduce task是并行計(jì)算的執(zhí)行實(shí)體。

在map task中，并行計(jì)算模塊520中的質(zhì)量計(jì)算單元5201通過并行計(jì)算堿基比對(duì)質(zhì)量值(BAQ)和最終質(zhì)量值來(lái)過濾掉測(cè)序序列位點(diǎn)上無(wú)效的堿基，并將過濾后保留下來(lái)的堿基保存到HDFS中。

并行計(jì)算模塊520并行地執(zhí)行map task和reduce task，先運(yùn)行map task，然后運(yùn)行reduce task。

為了尋找樣本的染色體突變位點(diǎn)，并行計(jì)算模塊520中的質(zhì)量計(jì)算單元5201需要并行地計(jì)算堿基比對(duì)質(zhì)量值(BAQ)和最終質(zhì)量值，并行計(jì)算模塊520中的過濾單元5202根據(jù)最終質(zhì)量值過濾掉測(cè)序序列位點(diǎn)上無(wú)效的堿基。這一部分并行計(jì)算是在map task中進(jìn)行的。有關(guān)計(jì)算BAQ和最終質(zhì)量值以及過濾掉測(cè)序序列位點(diǎn)上無(wú)效的堿基的算法描述請(qǐng)參見本說明書步驟S303至S306的相關(guān)描述。

在reduce task中，結(jié)果合并模塊530中的比率計(jì)算單元5301通過并行計(jì)算根據(jù)保留下來(lái)的堿基計(jì)算樣本的每一個(gè)位點(diǎn)上的突變堿基所占的比率，結(jié)果合并模塊530中的突變位點(diǎn)確定單元5302根據(jù)此比率確定樣本的突變位點(diǎn)。

在reduce task中，首先結(jié)果合并模塊530中的比率計(jì)算單元5301從HDFS中讀取過濾后保留下來(lái)的堿基。然后，對(duì)于樣本上的每一個(gè)位點(diǎn)，比率計(jì)算單元5301累加由測(cè)序數(shù)據(jù)產(chǎn)生的ATCG各堿基的次數(shù)(通過多臺(tái)計(jì)算機(jī)對(duì)其上的小數(shù)據(jù)塊進(jìn)行并行計(jì)算得出)，從而得到此位點(diǎn)上的突變堿基所占的比率。

結(jié)果合并模塊530中的突變位點(diǎn)確定單元5302按照以下方法確定樣本的染色體突變位點(diǎn)：如果樣本某一位點(diǎn)上的突變堿基所占的比率大于等于第二閾值，則確定此位點(diǎn)是樣本的染色體突變位點(diǎn)；否則確定此位點(diǎn)不是樣本的染色體突變位點(diǎn)。通過這種方法確定樣本的所有染色體突變位點(diǎn)。

在本發(fā)明的第三實(shí)施例中：

并行準(zhǔn)備模塊510中的數(shù)據(jù)切分單元5101將包含堿基比對(duì)結(jié)果的大數(shù)據(jù)文件切分為多個(gè)小文件塊。

從數(shù)據(jù)源輸入的文件為bam或者sam文件，其中包含了測(cè)序序列與參考序列的堿基比對(duì)結(jié)果的大量數(shù)據(jù)。

多個(gè)客戶端通過網(wǎng)絡(luò)上傳原始的bam或sam大文件，并行準(zhǔn)備模塊510中的數(shù)據(jù)切分單元5101將此大文件根據(jù)一定的規(guī)則切分成小文件塊保存到kafka中。切分規(guī)則如下：比如輸入的大文件為3G，一個(gè)block為128M，那么輸入的大文件總共將被切分為3*1024/128＝24塊。

并行準(zhǔn)備模塊510中的數(shù)據(jù)切分單元5101將從多個(gè)小文件塊中讀取到的數(shù)據(jù)切分為RDD，并行準(zhǔn)備模塊510中的執(zhí)行實(shí)體生成單元5201根據(jù)RDD生成并行計(jì)算的執(zhí)行實(shí)體job。

并行準(zhǔn)備模塊510中的數(shù)據(jù)切分單元5101以預(yù)定時(shí)間間隔(比如3秒)從kafka中保存的多個(gè)小文件塊讀取數(shù)據(jù)，并在此過程中將讀取到的數(shù)據(jù)切分為RDD。數(shù)據(jù)切分單元5101按時(shí)間片(即預(yù)定時(shí)間間隔)讀取到的所有數(shù)據(jù)組成了代表持續(xù)輸入的堿基比對(duì)結(jié)果的DStream。數(shù)據(jù)切分單元5101將對(duì)DStream的操作轉(zhuǎn)換為DStream的有向無(wú)環(huán)圖(DAG)。對(duì)于每個(gè)時(shí)間片(即預(yù)定時(shí)間間隔)，數(shù)據(jù)切分單元5101根據(jù)DStream DAG生成RDD DAG。并行準(zhǔn)備模塊510中的執(zhí)行實(shí)體生成單元510根據(jù)RDD DAG生成相應(yīng)的執(zhí)行實(shí)體job。

并行計(jì)算模塊520中的質(zhì)量計(jì)算單元5201確定需要計(jì)算堿基比對(duì)質(zhì)量值(BAQ)的堿基。

并行計(jì)算模塊520在參與并行計(jì)算的多臺(tái)計(jì)算機(jī)上并行地執(zhí)行前面步驟生成的多個(gè)job。

為了尋找樣本的染色體突變位點(diǎn)，需要過濾掉測(cè)序序列一部分位點(diǎn)上的堿基：測(cè)序數(shù)據(jù)是基于高通量測(cè)序獲取的，而高通量測(cè)序的錯(cuò)誤率在0.01％左右，所以需要去掉低質(zhì)量的數(shù)據(jù)(即無(wú)效的堿基)，不將其作為尋找樣本的染色體突變位點(diǎn)的依據(jù)。

過濾過程首先需要由并行計(jì)算模塊520中的質(zhì)量計(jì)算單元5201對(duì)測(cè)序序列的一部分位點(diǎn)的堿基計(jì)算堿基比對(duì)質(zhì)量值(BAQ)。質(zhì)量計(jì)算單元5201確定需要計(jì)算BAQ的堿基的方法如下：如果測(cè)序序列與參考序列所有對(duì)應(yīng)位點(diǎn)上的堿基完全一致，則將此測(cè)序序列所有位點(diǎn)的堿基作為尋找樣本的染色體突變位點(diǎn)的依據(jù)，無(wú)需計(jì)算它們的BAQ值；如果測(cè)序序列某一位點(diǎn)與參考序列對(duì)應(yīng)位點(diǎn)上的堿基不一致，則需要對(duì)此測(cè)序序列所有位點(diǎn)的堿基計(jì)算BAQ值。

當(dāng)測(cè)序序列與參考序列堿基不一致時(shí)，并行計(jì)算模塊520中的質(zhì)量計(jì)算單元5201計(jì)算測(cè)序序列所有位點(diǎn)的堿基比對(duì)質(zhì)量值(BAQ)。

請(qǐng)參閱圖4，如果測(cè)序序列(read序列)和參考序列(reference序列)不完全一致，則質(zhì)量計(jì)算單元5201對(duì)read序列和reference序列建立一個(gè)隱馬爾可夫模型(HMM)，并對(duì)read序列中狀態(tài)為M的堿基計(jì)算其堿基堿基比對(duì)質(zhì)量值(BAQ)，該值代表此次堿基比對(duì)結(jié)果錯(cuò)誤的可能性。

在HMM中，M代表堿基匹配(包括完全一致和錯(cuò)配兩種情況)，I代表插入，D代表缺失，S代表開始，E代表結(jié)束。S指向所有的I和M，所有的M和I指向E。

BAQ算法分為前向遞推和后向遞推。算法中用到的變量如下：

前向遞推：

(a)值的初始化，其中k＝1,…,L；i＝0,1。

f_0,S＝1

(b)遞推，其中i＝2,…,l；k＝1,…,L。

對(duì)于i＝l+1，有：

其中，a_ij由以下狀態(tài)轉(zhuǎn)換矩陣確定：

后向遞推：

(a)值的初始化，其中k＝1,…,L

b_l+1,E＝1

(b)遞推，其中i＝l-1,…,1；k＝L,…,1

當(dāng)i＝1時(shí)，δ_i1＝1，其他則為0。

對(duì)于i＝0，有：

BAQ通過下式計(jì)算得到：

其中，b_ij由以下狀態(tài)轉(zhuǎn)換矩陣確定：

并行計(jì)算模塊520中的質(zhì)量計(jì)算單元5201確定已計(jì)算過BAQ的測(cè)序序列所有堿基的最終質(zhì)量值。

對(duì)于已計(jì)算過BAQ的測(cè)序序列的堿基，質(zhì)量計(jì)算單元5201在BAQ和此堿基的測(cè)序質(zhì)量值(測(cè)序質(zhì)量值可自bam或sam文件中獲得)中取較小值作為最終質(zhì)量值。

并行計(jì)算模塊520中的過濾單元5202過濾掉無(wú)效的堿基。

并行計(jì)算模塊520中的過濾單元5202按照以下方法過濾掉測(cè)序序列上無(wú)效的堿基：在測(cè)序序列所有的堿基中，對(duì)于已計(jì)算過BAQ的堿基，如果其最終質(zhì)量值小于第一閾值，則認(rèn)為此堿基無(wú)效，并將其過濾掉，不作為尋找樣本的染色體突變位點(diǎn)的依據(jù)；否則認(rèn)為此堿基有效，將其保留并作為尋找樣本的染色體突變位點(diǎn)的依據(jù)。

結(jié)果合并模塊530中的比率計(jì)算單元5301根據(jù)保留下來(lái)的堿基計(jì)算樣本的每個(gè)位點(diǎn)上的突變堿基所占的比率。

對(duì)于樣本上的每一個(gè)位點(diǎn)，結(jié)果合并模塊530中的比率計(jì)算單元5301累加所有測(cè)序序列在該位點(diǎn)產(chǎn)生的ATCG各堿基的次數(shù)(通過對(duì)多臺(tái)計(jì)算機(jī)上的多個(gè)RDD進(jìn)行并行計(jì)算得出)，從而得到此位點(diǎn)上的突變堿基所占的比率。

結(jié)果合并模塊530中的突變位點(diǎn)確定單元5301按照以下方法確定每一個(gè)位點(diǎn)是否為突變位點(diǎn)：如果樣本某一位點(diǎn)上的突變堿基所占的比率大于等于第二閾值，則確定此位點(diǎn)是樣本的染色體突變位點(diǎn)；否則確定此位點(diǎn)不是樣本的染色體突變位點(diǎn)。突變位點(diǎn)確定單元5301通過這種方法確定樣本所有的染色體突變位點(diǎn)。

本領(lǐng)域普通技術(shù)人員可以理解，實(shí)現(xiàn)上述本發(fā)明實(shí)施例中的基于高通量測(cè)序數(shù)據(jù)尋找樣本的染色體突變位點(diǎn)的分析方法和分析裝置可以通過程序指令相關(guān)的硬件來(lái)完成，所述的程序可以存儲(chǔ)于可讀取存儲(chǔ)介質(zhì)中，該程序在執(zhí)行時(shí)執(zhí)行上述方法中的對(duì)應(yīng)步驟。

以上所述僅是本發(fā)明的優(yōu)選實(shí)施方式，應(yīng)當(dāng)指出，對(duì)于本技術(shù)領(lǐng)域的普通技術(shù)人員來(lái)說，在不脫離本發(fā)明原來(lái)的前提下，還可以做出若干改進(jìn)和潤(rùn)飾，這些改進(jìn)和潤(rùn)飾也應(yīng)視為本發(fā)明的保護(hù)范圍。