利用dna測序片段數(shù)據(jù)確定個(gè)體單體型問題

1snps和單體型為了徹底消除人類遺產(chǎn)信息，人類基金會(huì)于20世紀(jì)90年代開始了人類疾病研究計(jì)劃。2000年6月，人類疾病樣本的圖紙被繪制，2003年4月人類疾病的地圖基本完成。在此之前，人類疾病的共同特征被揭示，但不同的人具有不同的外表和身體、不同的抗逆能力和不同的藥物敏感性。遺傳上，不同的個(gè)體（除同生卵胎外）的基因不完全相同。兩組之間的dna差異約為組的0.1%。單個(gè)氨基酸含量的synps是人類染色體特定點(diǎn)的堿基變化。將synps廣泛分布在人類疾病基因中，預(yù)計(jì)將有數(shù)百個(gè)synps。單核苷酸多態(tài)性是一個(gè)物種中不同個(gè)體表型的主要遺傳來源.識(shí)別SNPs對(duì)基因的精確定位、了解基因功能很有幫助,對(duì)遺傳病等疾病的診斷和藥物研究有重要作用,SNPs可用于個(gè)體識(shí)別、親子鑒定,亦可用于人類各群體的遺傳關(guān)系分析.Stephens等利用單體型研究人類313個(gè)基因中的3899個(gè)SNPs,然后進(jìn)行連鎖不平衡分析,其結(jié)果支持了人類群體在近代擴(kuò)張的說法.Horikawa等根據(jù)SNPs進(jìn)行關(guān)聯(lián)分析在墨西哥裔美國人中把2型糖尿病基因定位在2號(hào)染色體長臂,并發(fā)現(xiàn)CAPN10基因的3個(gè)SNPs和2型糖尿病相關(guān).一個(gè)SNP位點(diǎn)指的是在一個(gè)物種的基因組DNA序列中不同個(gè)體可能出現(xiàn)不同堿基的位置.在一條染色體SNP位點(diǎn)上的堿基序列叫做單體型(haplotype).人類等二倍體生物的染色體是成對(duì)存在的,都有一對(duì)單體型.圖1中的一對(duì)單體型分別是“ATACG”和“GCATG”.單體型在SNPs的上述應(yīng)用中扮演著重要的角色,不幸的是在當(dāng)前的實(shí)驗(yàn)技術(shù)下,把個(gè)體的一對(duì)染色體分開,獨(dú)立測定每一條染色體上的單體型既費(fèi)錢又費(fèi)時(shí)間,因此利用計(jì)算機(jī)技術(shù)來確定個(gè)體的單體型具有極其重要的現(xiàn)實(shí)意義.本文研究如何利用個(gè)體的DNA測序片斷數(shù)據(jù),根據(jù)不同的優(yōu)化準(zhǔn)則確定該個(gè)體的單體型.本文第2節(jié)介紹個(gè)體單體型問題的抽象模型和當(dāng)前的研究現(xiàn)狀;第3節(jié)闡述基于無空隙片斷數(shù)據(jù)的參數(shù)化算法,分析其復(fù)雜度;第4節(jié)給出實(shí)驗(yàn)結(jié)果;最后進(jìn)行總結(jié)和展望.2dna片段和snp矩陣人類等雙倍體生物的DNA序列是按染色體成對(duì)出現(xiàn)的.對(duì)于任意一個(gè)SNP位點(diǎn)來說,2條染色體上的堿基可以是相同的,這種現(xiàn)象叫純合(homozygous);也可以是不同的,這種現(xiàn)象叫做雜合(heterozygous).這樣一條染色體在SNP位點(diǎn)的投影序列即單體型(haplotype)可以用字符集{A,B}上的字符序列來表示,不必用真正的堿基字符,其中‘A’通常表示人群在該位點(diǎn)上常見的SNP值,這樣圖1中的一對(duì)單體型則可表示為一對(duì)字符串“ABAAB”和“BAABB”.基于直接測定單體型的技術(shù)困難性和單體型在遺傳分析上的重要性,Lancia等提出了直接利用個(gè)體的DNA測序片斷數(shù)據(jù)來確定該個(gè)體的一對(duì)單體型,這個(gè)計(jì)算問題就是個(gè)體單體型問題(individualhaplotypingproblem).對(duì)基因組進(jìn)行測序時(shí),因?yàn)榧夹g(shù)上的限制,只能直接對(duì)較短的DNA片斷進(jìn)行測序,這些片斷來自一對(duì)染色體的不同單體,測序過程中不可避免地會(huì)發(fā)生一些錯(cuò)誤.給定某個(gè)個(gè)體一組已測序的DNA片斷數(shù)據(jù),個(gè)體單體型問題就是要去掉最小數(shù)量的數(shù)據(jù),以便能發(fā)現(xiàn)一對(duì)單體型與剩下的數(shù)據(jù)兼容.去掉最小數(shù)量的數(shù)據(jù)可以是去掉最少的片斷或最少的SNP位點(diǎn).n個(gè)SNP位點(diǎn)按在染色體上的次序從左到右記作S:{1,2,…,n},m個(gè)片斷記作F:{1,2,…,m}.任意SNP位點(diǎn)應(yīng)該被某些DNA片斷覆蓋,任意片斷在它所覆蓋的SNP位點(diǎn)的取值為{A,B,-},其中‘-’為空值,表示片斷在該位點(diǎn)的取值未知.DNA片斷的數(shù)據(jù)集可以表示為在{A,B,-}上的一個(gè)m×n的矩陣,叫做SNP矩陣M.圖2是一個(gè)8×8的SNP矩陣.SNP矩陣的列表示SNP位點(diǎn),行表示片段在對(duì)應(yīng)的SNP位點(diǎn)上的取值,Mij表示第i個(gè)片斷在第j個(gè)SNP位點(diǎn)上的取值.為了行文簡潔,下面引入與SNP矩陣M相關(guān)的幾個(gè)定義.定義1.如果(?k(Mik≠‘-’)∧(k≤j))∧(?r(Mir≠‘-’)∧(j≤r))),則稱行i覆蓋列j.行i覆蓋列j就是行i在列j上取值非空,或在列j前至少有一列,在列j后也至少有一列,使得行i在這兩列上的取值均非空.在圖2中行2覆蓋列2～7.如果某一行覆蓋某一列,但是該行在該列上的取值為空,則稱該行在該列上有個(gè)洞(hole).圖2中行2在列5上有個(gè)洞.如果M的所有行都沒有洞,則稱M為無空隙的SNP矩陣.定義2.如果兩行在某一列上的值都不是空值,且這兩行在該列上的值不相等,那么這兩行在該列上沖突.如果兩行在所有的列上均不沖突,則這兩行兼容.如果測序過程沒有任何錯(cuò)誤,代表來自于同一染色體的片段的行必定兩兩兼容.定義3.對(duì)于M的兩列j1和j2,如果M的所有行都可以劃分到兩個(gè)子集中,使得劃分在同一個(gè)子集中的任意兩行在列j1和j2上均不沖突,則列j1和j2兼容,滿足上述條件的劃分叫做在列j1和j2上兼容的劃分;如果不存在著在列j1和j2上兼容的劃分,則列j1和j2沖突.在圖2中,如果要在列3上不沖突,行2,3,6必須劃分到同一子集中,行4,5則必須劃分到另一子集中;但是這樣劃分的話,行6與行2在列4上沖突,由此可見,列3和4沖突.對(duì)于列1和列2而言,把行2和3分到一個(gè)子集中,行4和5分到另一個(gè)子集中的任意劃分都是在列1和2上兼容的劃分,所以列1和列2兼容.定理1.對(duì)于一個(gè)SNP矩陣M,只有‘A’或‘B’值的任意列不會(huì)和其它列沖突.對(duì)于既有‘A’又有‘B’值的列j1,j2,列j1,j2沖突當(dāng)且僅當(dāng)存在著兩行i1、i2:i1≠i2,有Mi1j1,Mi1j2,Mi2j1,Mi2j2均不為空值,且這4個(gè)值之中有3個(gè)相同、1個(gè)不同.證明.不失一般性,令列j1只有‘A’值,對(duì)其它任意列j2,M中的所有行均可如下劃分:在列j2上取A值的行劃分到第一個(gè)集合中,在列j2上取B值的行劃分到第2個(gè)集合中,其它的行隨機(jī)劃分到這兩個(gè)集合中,容易驗(yàn)證,這樣劃分到同一個(gè)集合中的任意兩行在列j1和j2上不會(huì)沖突.因此如果一列只有‘A’或‘B’值,該列不會(huì)和其它任何列沖突.對(duì)于既有‘A’又有‘B’值的列j,如果要使劃分在同一個(gè)子集中的行在列j上不沖突,則在列j上取值為‘A’的所有行必須劃分在同一子集中,在列j上取值為‘B’的行必須劃分在另一個(gè)子集中.對(duì)于既有‘A’又有‘B’值的兩列j1,j2,如果存在著兩行i1,i2,有Mi1j1,Mi1j2,Mi2j1,Mi2j2均不為空值,且這4個(gè)值之中有3個(gè)相同、1個(gè)不同,不失一般性,令Mi1j1=Mi2j1=Mi1j2=‘A’、Mi2j2=‘B’,那么按列j1上的取值,行i1和i2應(yīng)該劃分到同一個(gè)子集中;而按列j2上的取值,行i1和i2應(yīng)該劃分到不同的子集中,這相互矛盾,說明不存在列j1和j2上兼容的劃分.因此列j1,j2沖突.反過來,對(duì)于既有‘A’又有‘B’值的兩列j1,j2,如果列j1,j2不沖突,那么存在著一個(gè)劃分使M的所有行劃分到兩個(gè)子集中,使得相同子集中的任意兩行在列j1和j2上不沖突.這樣對(duì)任意兩行i1,i2,如果有Mi1j1,Mi1j2,Mi2j1,Mi2j2均不為空值,即不是‘A’值就是‘B’值,那么如果這兩行被劃分到同一個(gè)子集中,則這兩行在列j1上的值相等,在列j2上的值也應(yīng)相等;如果這兩行被劃分到不同的子集中,則這兩行在列j1上的值不相等,在列j2上的值也應(yīng)不相等.這樣這4個(gè)值之中就不可能有3個(gè)相同、1個(gè)不同的情況.定理得證.證畢.定義4.如果SNP矩陣M的所有行可以分成2個(gè)不相交的子集,每個(gè)子集中的所有行都相互兼容,則M是可行的.如果M是可行的,則很容易找到一個(gè)劃分,使劃分在同一個(gè)子集中的所有行都兼容,進(jìn)而通過同一個(gè)子集中的行很容易重建與這些行兼容的單體型.如果測序過程沒有任何錯(cuò)誤,測序片段數(shù)據(jù)對(duì)應(yīng)的SNP矩陣必定是可行的,那么通過DNA測序片段數(shù)據(jù)很容易得出個(gè)體的一對(duì)單體型.可是由于測序過程的錯(cuò)誤是不可避免的,因此實(shí)驗(yàn)室測得的片段數(shù)據(jù)對(duì)應(yīng)的SNP矩陣通常是不可行的.Lancia等最先討論對(duì)SNP矩陣進(jìn)行處理使其可行的計(jì)算模型,引入了下面的個(gè)體單體型優(yōu)化問題:問題1.MFR(MinimumFragmentRemoval,最少片斷刪除).給定一個(gè)SNP矩陣M,刪除最少的行使M可行.問題2.MSR(MinimumSNPRemoval,最少SNP刪除).給定一個(gè)SNP矩陣M,刪除最少的列使M可行.不進(jìn)行任何處理,圖2所示的SNP矩陣是不可行的.當(dāng)去掉第6行后,行1,4,5和8相互兼容,構(gòu)成一個(gè)子集,進(jìn)而得出一個(gè)單體型為“BABBABAB”;行2,3和7相互兼容,構(gòu)成另一個(gè)子集,進(jìn)而得出另一個(gè)單體型為“ABAABABA”.或者不去掉行,則去掉列4后,行1,4,5和8相互兼容,構(gòu)成一個(gè)子集,進(jìn)而得出一個(gè)單體型為“BAB-ABAB”;行2,3,6和7相互兼容,構(gòu)成另一個(gè)子集,進(jìn)而得出另一個(gè)單體型為“ABA-BABA”(去掉的列所代表SNP位點(diǎn)的值為空值).Lancia等證明了在DNA片斷數(shù)據(jù)有洞(hole)的情況下,MFR是NP-hard;如果片斷上洞的個(gè)數(shù)超過一個(gè),則MSR是NP-hard;并且在這些情況下,MSR和MFR問題都是APX-hard.在DNA片斷數(shù)據(jù)中無洞的情況下,Lancia等證明了MFR和MSR是多項(xiàng)式時(shí)間可解的.Bafna等對(duì)這些問題進(jìn)行了深入的研究,對(duì)于m個(gè)DNA片斷,n個(gè)SNP位點(diǎn)的無空隙的SNP矩陣的MSR和MFR問題,提出了時(shí)間復(fù)雜度為O(mn2)和O(m2n+m3)、空間復(fù)雜度為O(mn+n2)和O(mn+m3)的多項(xiàng)式算法.人有24條不同的染色體,染色體的平均長度約為150000000個(gè)堿基1,SNP的分布密度約為0.1%,這樣一條染色體上的SNP位點(diǎn)數(shù)約為150000.當(dāng)從shotgun全基因組測序所得到的序列數(shù)據(jù)推出個(gè)體的單體型時(shí),能直接測序的片段長度約為1000個(gè)堿基,一個(gè)堿基大約被10個(gè)不同的片段覆蓋,因此得到的片斷數(shù)約為1500000(150000,000×10/1000).在這種情況下,即使DNA片斷數(shù)據(jù)無空隙,解決MSR和MFR的這些算法的時(shí)間復(fù)雜度仍然太大,推斷染色體一個(gè)區(qū)域的單體型時(shí)情況也是如此.因而根據(jù)目前基因組測序的技術(shù)現(xiàn)狀,針對(duì)DNA測序數(shù)據(jù)和個(gè)體單體型問題的特點(diǎn)提出更高效的算法具有重要的現(xiàn)實(shí)意義.3參數(shù)化條件:k1,k1,k目前最流行的SNP探測方法是DNA直接測序,這種方法48h可分析近百萬個(gè)堿基對(duì),雜合的SNPs探測率達(dá)95%以上,SNP聯(lián)盟(SNPConsortium)用這種方法測得了上百萬SNPs.當(dāng)前DNA測序的主導(dǎo)方法是Sanger雙脫氧鏈終止法.采用Sanger雙脫氧鏈終止法測序,一次能測定的DNA序列的長度僅為800～1200個(gè)堿基.各大測序中心使用的第三代測序儀如MageBACE等可測片斷長度約為1000堿基2,而SNPs的平均分布密度約為1/1000,雖然SNPs在整個(gè)染色體上的分布很不均勻,從已有的數(shù)據(jù)來看,一個(gè)長度1000bp的片斷上的SNP位點(diǎn)是極其有限的,通常在10個(gè)以內(nèi).另外出于測序的時(shí)間和代價(jià)的考慮,基因組測序的DNA片斷對(duì)SNP位點(diǎn)的覆蓋度(覆蓋一個(gè)SNP位點(diǎn)的片段數(shù))是有限的.在目前的shotgun實(shí)驗(yàn)測序中,覆蓋度約為10.當(dāng)需要測定某個(gè)體在一條染色體上的單體型時(shí),覆蓋一個(gè)SNP位點(diǎn)的DNA片斷數(shù)遠(yuǎn)小于shotgun法測得的片斷總數(shù).因此,在DNA測序?qū)嶒?yàn)中,一個(gè)片斷覆蓋的SNP位點(diǎn)數(shù)(小于10)和覆蓋一個(gè)SNP位點(diǎn)的片斷數(shù)(約為10)與通常要探測的SNP位點(diǎn)數(shù)(n)及DNA片斷總數(shù)(m)相比,均是很小的數(shù).基于以上事實(shí),我們提出以下參數(shù)化條件.定義5.(k1,k2)參數(shù)化條件:k1,k2是正整數(shù),(k1,k2)參數(shù)化條件定義為片段覆蓋的SNP位點(diǎn)數(shù)不超過k1,覆蓋任意SNP位點(diǎn)的片段數(shù)不超過k2.對(duì)一個(gè)無空隙的SNP矩陣M而言,(k1,k2)參數(shù)化條件等價(jià)為矩陣M的每一行最多只有一塊連續(xù)非‘-’字符序列,該字符序列長度不超過k1個(gè);矩陣M每一列最多有k2個(gè)行的取值為非‘-’字符.下面本文深入研究無空隙的SNP矩陣M的MSR和MFR算法.對(duì)于任意無空隙的SNP矩陣M,很容易得出它的(k1,k2)參數(shù)值,因此為了行文簡潔,除了特別聲明外,以下的SNP矩陣都是指滿足(k1,k2)參數(shù)化條件的無空隙SNP矩陣.首先對(duì)MSR問題和MFR問題參數(shù)化.問題3.P_MSR(ParameterizedMinimumSNPRemoval).給定一個(gè)SNP矩陣M,滿足(k1,k2)參數(shù)化條件,P_MSR問題是要求刪除最少的列,也就是保留最多的列使M可行.問題4.P_MFR(ParameterizedMinimumFragmentRemoval).給定一個(gè)SNP矩陣M,滿足(k1,k2)參數(shù)化條件,P_MFR問題是要求刪除最少的行,也就是保留最多的行使M可行.本文中,一個(gè)SNP矩陣M的P_MSR解指的是使M可行能保留下來的最多的列數(shù),用P_MSR(M)來表示;一個(gè)SNP矩陣M的P_MFR解指的是使M可行能保留下的最多的行數(shù),用P_MFR(M)來表示.函數(shù)left和right分別用來表示M中的行覆蓋的最左邊和最右邊的列號(hào),即行i覆蓋的最左邊的列是left(i),覆蓋的最右邊的列是right(i).M的前i行構(gòu)成的SNP矩陣記作M(i,:).在求解問題3和問題4時(shí),先對(duì)SNP矩陣M進(jìn)行如下預(yù)處理:1.排序.調(diào)整M中各行的次序,使各行按其left值進(jìn)行非降序排列.同時(shí)對(duì)于任意列j,計(jì)算出覆蓋該列的行的序號(hào)的有序集,記作rowset(j).對(duì)于一個(gè)滿足(k1,k2)參數(shù)化條件的SNP矩陣M,排序后如圖3所示.2.去掉冗余列.對(duì)M的每一列j,如果rowset(j)中的所有行在該列的取值中沒有出現(xiàn)‘A’或沒有出現(xiàn)‘B’,那么該列就是冗余列,標(biāo)記該列.最后去掉所有標(biāo)記的列,調(diào)整剩下列的序號(hào),修改對(duì)應(yīng)行的left和right函數(shù)值.3.去掉冗余行:對(duì)M的每一行i,如果該行在剩下的列上的取值均為空,則該列就是冗余行,標(biāo)記該行.最后去掉所有標(biāo)記的行,修改受影響的列的rowset集.令M采用如下數(shù)據(jù)結(jié)構(gòu);對(duì)每一行i,記下該行覆蓋的最左邊和最右邊的列號(hào),即left(i)和right(i),記錄該行從列l(wèi)eft(i)到列right(i)的各列上的值.預(yù)處理時(shí)間復(fù)雜度分析:排序所需的時(shí)間為O(mlogm),對(duì)排好序的行掃描一遍可得到各列的rowset有序集,這樣步1所需時(shí)間為O(mlogm+mk1);去掉冗余列所需時(shí)間為O(nk2),修改left和right值所需時(shí)間為O(mk1),這樣步2所需的時(shí)間為O(nk2+mk1);去掉冗余行所需時(shí)間為O(mk1),修改rowset有序集所需時(shí)間為O(nk2),這樣步3所需的時(shí)間為O(mk1+nk2).所以整個(gè)預(yù)處理所需的時(shí)間為O(mlogm+nk2+mk1).定理2.一個(gè)滿足(k1,k2)參數(shù)化條件的SNP矩陣M經(jīng)過以上預(yù)處理后仍然滿足(k1,k2)參數(shù)化條件.證明.排序、去掉冗余列和冗余行的預(yù)處理顯然不會(huì)增加某一行覆蓋的列數(shù)和覆蓋某一列的行數(shù).證畢.定理3.一個(gè)SNP矩陣M經(jīng)過以上預(yù)處理后得到的SNP矩陣記作M′,去掉的冗余列集和冗余行集分別記作X和Y,則M的P_MSR解等于M′的P_MSR解與X的列數(shù)之和,M的P_MFR解等于M′的對(duì)應(yīng)解與Y的行數(shù)之和.證明.首先證明M的P_MSR解等于M′的P_MSR解與X的列數(shù)之和.因?yàn)槿サ舻牧性谒械男猩系娜≈抵缓小瓵’或‘B’,所以M的任意一行均不會(huì)在去掉的列上和其它行沖突.令S是M′的一些或全部列構(gòu)成的集合,容易看出,如果保留S中的所有列可使M′可行,那么保留S∪X中的所有列也可使M可行.由此可證M的P_MSR解等于M′的P_MSR解與X的列數(shù)之和.同理Y中的行不會(huì)與其它任何行沖突,因此M′的P_MFR解中所保留下的行加上Y中的行所形成的SNP矩陣肯定是可行的,因此可以證明M的P_MFR解等于M′的P_MFR解與Y的行數(shù)之和.證畢.從定理3可知,求解一個(gè)SNP矩陣的P_MSR和P_MFR解可以歸結(jié)為求其預(yù)處理后的SNP矩陣的對(duì)應(yīng)解.對(duì)于行列數(shù)為0的SNP矩陣,易知其P_MSR和P_MFR解均為0.為了敘述的簡潔,下面討論的SNP矩陣M均指預(yù)處理后的SNP矩陣,該矩陣具有如下特點(diǎn):行列數(shù)m,n>0;任何一個(gè)列總有一些行的值是‘A’,還有一些行的值為‘B’;對(duì)于行i≤j,有l(wèi)eft(i)≤left(j).3.1基于pmsr的無孔隙snp矩陣m的最優(yōu)模型由文獻(xiàn)的定理3可知一個(gè)無空隙的SNP矩陣M是可行的當(dāng)且僅當(dāng)其任意兩列均不沖突,因此P_MSR問題等價(jià)于保留最多的列,使保留的列中的任意兩列均不沖突.令C是M一個(gè)列的子集,如果C中的列彼此兼容,且C中的列的最大列號(hào)為j,則C稱為列j的兼容列集.對(duì)于列j的任意兼容列集C,容易看出保留C中的所有列可使M可行.令K(j)表示列j的最大兼容列集.顯然:K(1)={1}(1)由文獻(xiàn)的引理2可知對(duì)于一個(gè)無空隙的SNP矩陣M的任意3列j1<j2<j3,如果j1與j2兼容,j2與j3兼容,則一定有j1與j3兼容.因此對(duì)于任意兩列r、j:r<j,如果列r和j兼容,則K(r)∪{j}一定是列j的兼容列集,由此容易證明Κ(j)={j}∪maxr∶1≤r<j,列r和j兼容Κ(r),其中對(duì)集合進(jìn)行的max運(yùn)算是取元素最多的集合,即maxr∶1≤r<j,列r和j兼容Κ(r)為滿足1≤r<j且列r和j兼容的所有K(r)中的最大集合,如果沒有滿足這些條件的K(r),則為空集?,下文也是如此.對(duì)于列j和它前面的列r,下面判斷列j是否和r相容:Case1.r≤j-k1:因?yàn)镸滿足(k1,k2)參數(shù)化條件,所以M的任一行覆蓋的列數(shù)不會(huì)超過k1,這樣就不可能存在著兩行i1、i2:i1≠i2,Mi1j1,Mi1j2,Mi2j1,Mi2j2均不為空值.根據(jù)定理1,列r與j兼容.Case2.r>j-k1:因?yàn)镸經(jīng)過了預(yù)處理,所以任意一列上既有‘A’又有‘B’.根據(jù)定理1,只需對(duì)在列r與j上取值均非空的行進(jìn)行觀察,這些行的個(gè)數(shù)不大于k2.如果這些行在這兩列上的取值只有“AA”、“BB”兩種類型,或只有“AB”、“BA”兩種類型,那么列r與j兼容;否則就不兼容.令A(yù)(j)表示K(1)到K(j)的最大值,如果j<0,規(guī)定A(j)=?;令OK(j)={r|j>r>j-k1且列r與j兼容}.由上面的討論,易知下面的公式成立:Κ(j)={j}∪max(A(j-k1),maxr∈ΟΚ(j)(Κ(r)))(2)由P_MSR問題的定義可知:P_MSR(M)=max(|A(n-k1)|,|K(n-k1+1)|,…,|K(n)|)(3)在式(1)～(3)的基礎(chǔ)上,可得到求解無空隙SNP矩陣M的P_MSR問題的動(dòng)態(tài)規(guī)劃算法,如圖4所示.定理4.對(duì)于一個(gè)m×n滿足(k1,k2)參數(shù)化條件的無空隙SNP矩陣M,P_MSR算法是正確的,加上預(yù)處理其時(shí)間復(fù)雜度為O(nk1k2+mlogm+mk1),空間復(fù)雜度為O(mk1+nk1).證明.P_MSR算法的正確性由式(1)～(3)的正確性來保證,而式(1)～(3)的正確性在引入的時(shí)候已經(jīng)得到了證明.下面分析其復(fù)雜度.時(shí)間復(fù)雜度分析.步1中掃描M得到k1和k2所需的時(shí)間為O(mk1);步2.2.1中判斷列r與j是否兼容所需時(shí)間為O(k2),步2.2.1最多執(zhí)行nk1次,由此可知步2的時(shí)間復(fù)雜度為O(nk1k2);步3的時(shí)間復(fù)雜度為O(1);步4的時(shí)間復(fù)雜度為O(k1).這樣加上預(yù)處理的時(shí)間整個(gè)算法的時(shí)間復(fù)雜度為O(nk1k2+mlogm+mk1).空間復(fù)雜度分析.SNP矩陣所需的空間為O(mk1),計(jì)算K(j)只需要A(j-k1),K(j-k1+1),…,K(j-1)的值,因此A和K均可采用大小為k1的循環(huán)隊(duì)列,需要的空間為O(nk1),所以整個(gè)算法的空間復(fù)雜度為O(mk1+nk1).定理得證.證畢.3.2多兼容行集的要求考慮m×n的SNP矩陣M的前i行構(gòu)成的矩陣M(i,:),令R為M(i,:)的行的子集,如果保留R中的所有行能使M(i,:)可行,那么R就稱為M(i,:)的一個(gè)兼容行集.如果R是M(i,:)的一個(gè)兼容行集,則必定可以把R劃分成兩個(gè)子集R1和R2,使得劃分在同一子集的行彼此兼容.對(duì)于Rj:j=1,2,找出具有以下特性的行rj作其代表:(1)rj∈Rj;(2)對(duì)于任意行r∈Rj,有right(r)<right(rj),或者right(r)=right(rj)且r≤rj.在“right(r1)<right(r2)”或“right(r1)=right(r2)且r1<r2”的情況下,R叫做行r1,r2代表的M(i,:)的兼容行集;否則R叫做行r2,r1代表的M(i,:)的兼容行集.劃分得來的子集在極端情況下可能是空集,為了上述概念的完整性,規(guī)定空集的代表行為0,而且left(0)=-1,right(0)=-2,行0與所有其它行均兼容.顯然空集也應(yīng)是M(i,:)的一個(gè)兼容行集,所以進(jìn)一步規(guī)定空集是行0,0代表的M(i,:)的兼容行集.令R(d,k,i)表示以行d,k為代表的M(i,:)的最大兼容行集.易知R(0,0,1)=?;R(0,1,1)={1}(4)令R(*,k,i)表示所有滿足right(d)<left(i+1)的最大R(d,k,i),即R(*,k,i)=maxd:right(d)<left(i+1)R(d,k,i).令R(*,*,i)表示所有滿足right(k)<left(i+1)的最大R(*,k,i),即R(*,*,i)=maxk:right(k)<left(i+1)R(*,k,i).這樣就有R(*,0,1)=?;R(*,1,1)={1};R(*,*,1)={{1},right(1)<left(2)?,否則(5)為了使R(*,k,i)、R(*,*,i)在i=m時(shí)有意義,引入行m+1,規(guī)定left(m+1)=n+1.這樣,根據(jù)P_MFR的定義,顯然有P_MFR(M)=|R(*,*,m)|(6)即R(*,*,m)中的行數(shù).定理5.對(duì)于一個(gè)預(yù)處理后的無空隙SNP矩陣M,行i+1與R(d,k,i)中以行d(或k)為代表的子集中所有行兼容當(dāng)且僅當(dāng)行i+1與行d(或k)兼容.證明.如果行i+1與R(d,k,i)中以行d(或k)為代表的子集中所有行兼容,那么顯然行i+1與行d(或k)兼容.下面證明如果行i+1與行d兼容,則行i+1與R(d,k,i)中以行d為代表的子集中所有行兼容:對(duì)于R(d,k,i)中以行d為代表的子集中任意一行r,必有right(r)≤right(d),且行r與d在列max(left(r),left(d))到列right(r)中的任一列上取值必定非空(無空隙),而且相等.由于M是經(jīng)過排序的,且r,d≤i+1,所以max(left(r),left(d))≤left(i+1).由于M是無空隙的SNP矩陣,且行i+1與行d兼容,所以行i+1與行d在列l(wèi)eft(i+1)到列min(right(i+1),right(d))中的任一列上取值必定相等,而且非空.這樣行i+1與行r在列l(wèi)eft(i+1)到列min(right(i+1),right(r))中的任一列上取值必定相等,而在其它的列,總有一行在該列的值為空,因此行i+1與行r在任意列上均不沖突,行i+1與行r兼容.同樣可以證明如果行i+1與行k兼容,則行i+1與R(d,k,i)中以行k為代表的子集中所有行兼容.定理得證.證畢.為了行文簡潔,滿足下列2個(gè)條件的行k的集合記作Sk(i):(1)k≤i;(2)right(k)≥left(i).滿足下列4個(gè)條件的(d,k)的集合記作Sdk(i):(1)d,k<i;(2)right(d)≥left(i);(3)right(k)≥left(i);(4)“right(d)<right(k)”或“right(d)=right(k)且d<k”.令C(i)={r|r∈Sk(i),且行r與i兼容}.下面討論對(duì)于行i:2≤i≤m,在R(*,*,i-1)、所有的R(*,k,i-1):k∈Sk(i)和所有的R(d,k,i-1):(d,k)∈Sdk(i)都已知的條件下,如何求出R(*,*,i)、所有可能的R(*,k,i)(k∈k(i+1))和所有的R(d,k,i)((d,k)∈Sdk(i+1)).首先對(duì)于所有的(d,k)∈Sdk(i),計(jì)算R(d,k,i):R(d,k,i-1)顯然是一個(gè)以行d,k為代表的M(i,:)的兼容行集,而M(i,:)比M(i-1,:)多的行只有行i,易知R(d,k,i)比R(d,k,i-1)最多只會(huì)增加一個(gè)元素,增加的元素只可能是行i.在滿足“i與d兼容,right(i)<right(d)”或“i與k兼容,right(i)<right(k)”的條件下,根據(jù)定理5,{i}∪R(d,k,i-1)顯然是一個(gè)以行d,k為代表的M(i,:)的兼容行集,所以有R(d,k,i)={i}∪R(d,k,i-1)(7)否則,必定有R(d,k,i)=R(d,k,i-1)(8)這是因?yàn)槿缦率聦?shí):條件“i與d兼容,right(i)≤right(d)”得不到滿足,則行i如果劃分到行d代表的子集中,那么或者d所在的子集中的行不再相互兼容,或者d所在的子集的代表不再是d而應(yīng)該是i.而條件“i與k兼容,right(i)≤right(k)”得不到滿足,那么或者k所在的子集中的行不再相互兼容,或者k所在的子集的代表不再是k而應(yīng)該是i.這兩個(gè)條件都得不到滿足,則行i必定不能在以d,k為代表的M(i,:)的兼容行集之中.對(duì)于所有的r∈Sk(i):如果right(r)≤right(i),有下面的公式成立(證明見定理6):R(r,i,i)={i}∪max(R(*,r,i-1),maxd:d∈C(i),(d,r)∈Sdk(i)R(d,r,i-1),maxk:k∈C(i),(r,k)∈Sdk(i),right(k)≤right(i)R(r,k,i-1))(9)否則,即right(r)>right(i),有下面的公式成立(證明見定理6):R(i,r,i)={i}∪max(R(*,r,i-1),maxd:d∈C(i),(d,r)∈Sdk(i),right(d)≤right(i)R(d,r,i-1))(10)對(duì)于所有的k∈Sk(i),計(jì)算R(*,k,i):令I(lǐng)表示maxd:right(d)<left(i)R(d,k,i),Ιdk(i)表示Sdk(i)∪{(i,k)|k∈Sk(i),right(k)>right(i)}.根據(jù)R(*,k,i)的定義:R(*,k,i)=max(Ι,maxd:left(i)≤right(d)<left(i+1)R(d,k,i))=max(Ι,maxd:(d,k)∈Ιdk(i),right(d)<left(i+1)R(d,k,i)).對(duì)于所有的d≤left(i),如果k與i兼容且right(k)>right(i),有R(d,k,i)=R(d,k,i-1)∪{i}(理由與式(7)相同),Ι=maxd:right(d)<left(i)R(d,k,i-1)∪{i}={i}∪R(*,k,i-1)；否則有R(d,k,i)=R(d,k,i-1)(理由與式(8)相同),Ι=maxd:right(d)<left(i)R(d,k,i-1)=R(*,k,i-1).因此有下面公式成立:如果k與i兼容且right(k)>right(i):R(*,k,i)=max(R(*,k,i-1)∪{i},maxd:(d,k)∈Ιdk(i),right(d)<left(i+1)R(d,k,i))(11)否則,必定有R(*,k,i)=max(R(*,k,i-1),maxd:(d,k)∈Ιdk(i),right(d)<left(i+1)R(d,k,i))(12)最后計(jì)算R(*,i,i)和R(*,*,i)(證明見定理6):R(*,i,i)=max({i}∪R(*,*,i-1),{i}∪maxk:k∈C(i),right(k)≤right(i)R(*,k,i-1),maxd:d∈Sk(i),right(d)≤right(i),right(d)<left(i+1)R(d,i,i))(13)R(*?*?i)=max(R(*?*?i-1),maxk:k∈Sk(i)∪{i},right(k)<left(i+1)R(*,k,i))(14)由于Sk(i+1)與Sk(i)的差集或者是空或者是{i},Sdk(i+1)與Sdk(i)的差集最多是{(d,i)|d∈Sk(i),right(d)≤right(i)}∪{(i,k)|k∈Sk(i),right(i)≤right(k)},所以通過式(7)～(14)可以求出R(*,*,i)、所有可能的R(*,k,i)(k∈Sk(i+1))和所有的R(d,k,i)((d,k)∈Sdk(i+1)).根據(jù)式(4)～(14),可得出下面的P_MFR動(dòng)態(tài)規(guī)劃算法,如圖5所示.定理6.對(duì)于一個(gè)預(yù)處理后滿足(k1,k2)參數(shù)化條件的無空隙m×nSNP矩陣M,P_MFR算法是正確的,加上預(yù)處理,其時(shí)間復(fù)雜度為O(mk22+mk1k2+mlogm+nk2),空間復(fù)雜度為O(mk1+mk22).證明.P_MFR算法的正確性取決于式(4)～(14)的正確性.式(4)～(8)、(11)和(12)的證明在公式的引入時(shí)已經(jīng)給出,下面證明其它公式.在式(9)中,顯然只有當(dāng)right(r)≤right(i)時(shí),R(r,i,i)才會(huì)有意義.令I(lǐng)1表示R(*,r,i-1),I2表示maxd:d∈C(i),(d,r)∈Sdk(i)R(d,r,i-1),I3表示maxk:k∈C(i),(r,k)∈Sdk(i),right(k)≤right(i)R(r,k,i-1).首先證明{i}∪I1,{i}∪I2,{i}∪I2都是以r,i為代表的M(:,i)的兼容行集.根據(jù)定義,I1即R(*,r,i-1)存在著一個(gè)劃分,使得劃分在同一個(gè)子集中的行彼此兼容,且其中一個(gè)子集的代表是行r,另一個(gè)子集的代表為行d且right(d)≤left(i).顯然行i與行d沒有共同覆蓋的列,所以行i與d兼容,進(jìn)而根據(jù)定理5,行i與d所代表的子集中的所有行均兼容.這樣行i加入到d所在的子集后,行i將取代d作為該子集的代表,這樣{i}∪I1就是一個(gè)以r,i為代表的M(:,i)的兼容行集.對(duì)于任意(d,r)∈Sdk(i),同樣存在著R(d,r,i-1)的一個(gè)劃分,使得劃分在同一個(gè)子集中的行彼此兼容,且其中一個(gè)子集的代表行是r,另一個(gè)子集的代表為行d.如果d∈C(i),即行i與d兼容,那么根據(jù)定理5,行i與d所代表的子集中的所有行均兼容.行i加入到d所在的子集后,由于right(d)≤right(r)≤right(i),行i將取代d作為該子集的代表,顯然{i}∪R(d,r,i-1)是一個(gè)以r,i為代表的M(:,i)的兼容行集,因此{(lán)i}∪I2也是.同理,對(duì)于任意(r,k)∈Sdk(i),R(r,k,i-1)存在著一個(gè)劃分,使得劃分到同一個(gè)子集中的行彼此兼容,且其中一個(gè)子集的代表是行r,另一個(gè)為k.如果k∈C(i),即行i與k兼容,根據(jù)定理5,行i與k代表的子集中的所有行均兼容.行i加入到k所在的子集后,如果right(k)≤right(i),行i將取代k作為該子集的代表行,這樣{i}∪R(r,k,i-1)就是一個(gè)以r,i為代表的M(:,i)的兼容行集,因此{(lán)i}∪I3也是.現(xiàn)在證明{i}∪max(I1,I2,I2)是一個(gè)以r,i為代表的M(:,i)的最大兼容行集.采用反證法,假設(shè){i}∪max(I1,I2,I3)不是M(:,i)的一個(gè)以r,i為代表的最大兼容行集,這就是說R(r,i,i)中的行數(shù)比{i}∪I1、{i}∪I2和{i}∪I3中的行數(shù)都要多.根據(jù)定義,R(r,i,i)中的行可以劃分成兩個(gè)子集,同一個(gè)子集的行互相兼容,且這兩個(gè)子集中的行的代表分別是r,i.令i所在的子集去掉行i后的代表行為l,顯然l與i兼容.如果right(l)≤left(i),則R(r,i,i)去掉行i后得到的行集R(r,i,i)-{i}是以l,r為代表的M(:,i-1)的一個(gè)兼容行集,如果假設(shè)成立,則R(r,i,i)-{i}比I1大,這與R(*,r,i-1)的定義矛盾.如果right(l)≥left(i),那么或者(l,r)∈Sdk(i)或者(r,l)∈Sdk(i)(因?yàn)閞∈Sk(i)).如果(l,r)∈Sdk(i),則R(r,i,i)-{i}是以l,r為代表的M(:,i-1)的一個(gè)兼容行集.如果假設(shè)成立,則R(r,i,i)-{i}比I2大,這與R(l,r,i-1)的定義矛盾;如果(r,l)∈Sdk(i),則R(r,i,i)-{i}是以r,l為代表的M(:,i-1)的一個(gè)兼容行集.如果假設(shè)成立,則R(r,i,i)-{i}比I3大,這與R(r,l,i-1)的定義矛盾.由此式(9)對(duì)于所有滿足right(r)≤right(i)的(r,i)都成立,用同樣的方法可以證明式(10)成立.式(13)的證明:令I(lǐng)1表示maxd:right(d)<left(i)R(d,i,i)?I2表示maxd:d∈Sk(i),right(d)≤right(i),right(d)<left(i+1)R(d,i,i).根據(jù)定義,R(*,,i,i)=max(I1,I2),這樣證明式(13)只需證明Ι1=max({i}∪R(*?*?i-1),{i}∪maxk:k∈C(i),right(k)≤right(i)R(*,k,i-1))={i}∪max(R(*?*?i-1),maxk:k∈C(i),right(k)≤right(i)R(*,k,i-1)).對(duì)于滿足條件“right(d)≤left(i)”的任意d,用證明式(9)的方法可以證明下式成立:R(d,i,i)={i}∪max(R(*,d,i-1),maxk:k∈C(i),right(k)≤right(i)R(d,k,i-1))?因此Ι1={i}∪max(maxd:right(d)<left(i)R(*,d,i-1),maxk:k∈C(i),right(k)≤right(i)(maxd:right(d)<left(i)R(d,k,i-1)))={i}∪max(R(*?*?i-1),maxk:k∈C(i),right(k)≤right(i)R(*,k,i-1)).式(13)成立.式(14)的證明:令I(lǐng)1代表maxk:right(k)<left(i)R(*,k,i)?I2代表maxk:k∈Sk(i)∪{i},right(k)<left(i+1)R(*,k,i).根據(jù)定義,R(*,*,i)=max(I1,I2).對(duì)于滿足right(k)≤left(i)的任意R(d,k,i)而言,有right(d)≤right(k)≤left(i),因此有R(d,k,i)=R(d,k,i-1)(理由同式(8)),所以Ι1=maxk:right(k)<left(i)(maxd:right(d)<left(i)R(d,k,i))=maxk:right(k)<left(i)(maxd:right(d)<left(i)R(d,k,i-1))=maxk:right(k)<left(i)R(*,k,i-1)=R(*?*?i-1).式(14)成立.算法的時(shí)間復(fù)雜度:算法的主體部分是步3.由于M經(jīng)過了預(yù)處理,那么k≤i意味著left(k)≤left(i),這時(shí)如果還有right(k)≥left(i),則可以肯定行k覆蓋列l(wèi)eft(i),所以Sk(i)中的元素個(gè)數(shù)|Sk(i)|不會(huì)超過覆蓋列l(wèi)eft(i)的行數(shù).由于M滿足(k1,k2)參數(shù)化條件,所以|Sk(i)|≤k2.這樣步3.1循環(huán)不超過(m-1)k2次,在步3.1.3中,判斷兩行是否兼容只需檢查兩行在它們共同覆蓋的列上是否沖突,步3.1.3執(zhí)行一次所需時(shí)間為O(k1),所以步3.1所需時(shí)間為O(mk1k2);可以看出Sdk(i)中的元組個(gè)數(shù)不會(huì)超過k22,所以步3.3和步3.5循環(huán)不超過(m-1)k22次,所需時(shí)間為O(mk22);剩下的步3.4、步3.6和步3.7所需時(shí)間為O(mk1).因此加上預(yù)處理,整個(gè)算法的時(shí)間復(fù)雜度為O(mk22+mk1k2+mlogm+nk2).算法的空間復(fù)雜度:SNP矩陣所需的空間為O(mk1),記錄所有的R(d,k,i)需要的空間O(mk22)(只需記錄相鄰的兩組R值,即i-1和i),所以算法的空間復(fù)雜度為O(mk1+mk22).定理得證.證畢.4模擬數(shù)據(jù)生成器測試與分析我們用C++語言實(shí)現(xiàn)了MFR(從文獻(xiàn)作者的源程序Fasthare中移植過來)、MSR、P_MFR和P_MSR算法(源程序可通過E-mail向本文作者索取),在一臺(tái)Linux服務(wù)器(4個(gè)IntelXeon3.6GHzCPU,4GBRAM)上對(duì)MSR和P_MSR、MFR和P_MFR算法的運(yùn)行時(shí)間(Runningtime)和單體型重建率(Reconstructionrate)進(jìn)行了比較.單體型重建率指的是算法重建出的單體型中正確的SNP位點(diǎn)數(shù)與總的SNP位點(diǎn)數(shù)的比值.實(shí)驗(yàn)中的單體型采用2種方式得到,第1種與文獻(xiàn)相同,采用來自公開數(shù)據(jù)庫的真實(shí)的單體型,本文實(shí)驗(yàn)采用的真實(shí)單體型數(shù)據(jù)來自于國際人類基因組單體型圖計(jì)劃2006年7月發(fā)布的數(shù)據(jù)文件genotypes_chr1_CEU_r21_nr_fwd_phased.gz3,該文件中包含了CEPH樣本(祖籍是北歐或西歐的美國猶他州人)中60個(gè)個(gè)體的單體型,每個(gè)單體型有SNP位點(diǎn)193333個(gè),本文實(shí)驗(yàn)隨機(jī)選擇一個(gè)個(gè)體指定長度的一對(duì)單體型.第2種與文獻(xiàn)一樣用計(jì)算機(jī)模擬生成,即首先隨機(jī)生成指定長度的單體型,根據(jù)指定的兩個(gè)單體型的差異率來隨機(jī)生成另一個(gè)單體型,本文采用差異率與文獻(xiàn)一樣,為20%.由于原始的DNA片段測序數(shù)據(jù)很難得到,在得到一對(duì)單體型的基礎(chǔ)上,上述文獻(xiàn)均根據(jù)指定的參數(shù)利用計(jì)算機(jī)來隨機(jī)生成片段數(shù)據(jù)集.實(shí)驗(yàn)室中,Sanger雙脫氧鏈終止法的DNA測序誤差約為1%,片段的覆蓋度約為10.為了使模擬生成的片段數(shù)據(jù)能很好地反映真實(shí)情況,與文獻(xiàn)一樣,本文采用著名的shotgun測序模擬數(shù)據(jù)生成器Celsim.下面的實(shí)驗(yàn)采用的測序誤差為1%,單體型長度,即SNP位點(diǎn)數(shù),在一個(gè)區(qū)間內(nèi)變化,生成片段的最小長度為3,片段的最大長度和覆蓋度在沒有特別說明的情況下分別為7和10,生成的片段數(shù)則按(單體型長度×片段覆蓋度/片段平均長度)設(shè)置.模擬數(shù)據(jù)生成器的詳細(xì)情況請(qǐng)參照文獻(xiàn).圖6～圖11的每一個(gè)點(diǎn)均為100次重復(fù)測試的平均值.下面各圖中,圖(a)均為在真實(shí)單體型數(shù)據(jù)上的實(shí)驗(yàn)結(jié)果,圖(b)均為在模擬的單體型數(shù)據(jù)上的實(shí)驗(yàn)結(jié)果.因?yàn)樵谀M的單體型數(shù)據(jù)上的實(shí)驗(yàn)結(jié)果與在真實(shí)單體型數(shù)據(jù)上的實(shí)驗(yàn)結(jié)果基本相同,所以下面主要討論在真實(shí)單體型數(shù)據(jù)上的實(shí)驗(yàn)結(jié)果.每個(gè)子圖中,左邊的Y軸表示算法的重構(gòu)精度,右邊的Y軸表示算法的運(yùn)行時(shí)間.圖6和圖7顯示算法性能隨SNP位點(diǎn)數(shù)(對(duì)應(yīng)SNP矩陣的列數(shù))變化的情況.圖6(a)中,當(dāng)位點(diǎn)數(shù)為500時(shí),P_MSR和MSR算法的單體型重構(gòu)率分別是84.62%和84.54%,運(yùn)行時(shí)間為0.0001s和0.191s;當(dāng)位點(diǎn)數(shù)增加到3500時(shí),P_MSR和MSR算法的單體型重構(gòu)率分別是82.76%和83.18%,運(yùn)行時(shí)間分別是0.004s和63.82s.圖7(a)中,當(dāng)位點(diǎn)數(shù)為80時(shí),P_MFR和MFR算法的單體型重構(gòu)率分別是91.40%和91.55%,運(yùn)行時(shí)間分別是0.009s和0.67s;當(dāng)位點(diǎn)數(shù)增加到260時(shí),P_MFR和MFR算法的單體型重構(gòu)率分別是88.68%和88.10%,運(yùn)行時(shí)間分別是0.055s和146.05s.從圖6和圖7可以看出,隨著SNP位點(diǎn)數(shù)的增加,算法的單體型重構(gòu)精度有下降趨勢;MFR和MSR運(yùn)行時(shí)間顯著增長,這是因?yàn)楫?dāng)SNP位點(diǎn)數(shù)增加,覆蓋度不變時(shí),片斷數(shù)也隨著增加,所以MFR和MSR運(yùn)行時(shí)間的增長速度是位點(diǎn)數(shù)增長速度的3次方,而P_MFR和P_MSR的時(shí)間則基本成線性增長.在片段數(shù)和其它參數(shù)保持不變的條件下,圖8和圖9通過改變片段的最大長度和SNP位點(diǎn)數(shù)來比較各算法的性能.圖8的片段數(shù)保持4000不變.圖8(a)中,在片段最大長度為6、SNP位點(diǎn)數(shù)為1800時(shí),P_MSR和MSR算法的單體型重構(gòu)率分別是84.6%和84.7%,運(yùn)行時(shí)間分別是0.001s和9.5s;當(dāng)片段最大長度為9、SNP位點(diǎn)數(shù)為2400時(shí),P_MSR