武鵬飛2014年清華杯數(shù)學(xué)建模競賽

1、2014 年“杯”數(shù)學(xué)建模競賽承諾書數(shù)學(xué)建模競賽的競賽規(guī)則.仔細閱讀了完全明白，在競賽開始后參賽隊員不能以任何方式（包括、電子郵件、網(wǎng)上等）與隊外的任何人（包括指導(dǎo)教師）研究、與賽題有關(guān)的問題。知道，別人的成果是競賽規(guī)則的, 如果別人的成果或其他公開的資料（包括網(wǎng)上查到的資料），必須按照規(guī)定的參考文獻的表述方式在正文處和參考文獻中明確列出。鄭重承諾，嚴(yán)格遵守競賽規(guī)則，以保證競賽的公正、公平性。競賽規(guī)則的行為，受到嚴(yán)肅處理。參賽的題目是：B 題：組組裝 參賽隊員 (打印并簽名) ：1.武鵬飛2.3.日期： 2014年 4月 28 日一種 DNA摘要：使用一種基于 de序列組裝算法的設(shè)計與測試Br

2、uijn 圖的算法，對中的片段組裝提出新的解決方案，并結(jié)合實例分析該方法的科學(xué)性與合理性。：；序列組裝；實例分析1一、，即分析特定 DNA 分子問題背景組中堿基的排列方式，在生命科學(xué)研DNA究中具有重要意義。目前人們能直接測量的長度較小，所以對較長的 DNA序列只能采用間接測量的方法。一種較為廣泛使用的長 DNA方法稱為“鳥槍法”1，其做法是將組一定份數(shù)，然后隨機打斷為大量短片段，對這些短片段進序，然后利用重合部分信息進行組裝。常用的組裝算法一般基于OLC（Overlap/Layoonsensus）方法、貪婪圖方法、de Bruijn 圖方法等。本文探索了一種基于 de Bruijn 圖的方法

3、，目標(biāo)是降低執(zhí)行算法所需的時間與空間，并盡可能提高組裝序列的連續(xù)性、完整性和準(zhǔn)確性。二、假設(shè)片段測量情況十分復(fù)雜，為簡化分析，一些主要的假設(shè)為：實際的假設(shè)所有認為所有片段的所有堿基是完全隨機并獨立的。均為獨立的A,C,G,T上均勻分布。實際上為組成有效的基因，堿基的分布可能不是均勻或獨立的。假設(shè)片段的也是完全隨機并獨立的。假設(shè)若兩條制的結(jié)果。片段中連續(xù) n 個堿基相同，則這兩條片段是同一位置復(fù)設(shè) n=29，根據(jù)之前的假設(shè)，一條長 120000 的組包含一組 n 長重復(fù)2120000 2 2.510-8，認為可以忽略。片段的概率不大于294實際的組可能存在或短或長的重復(fù)片段，但為了簡化的匹配與片

4、段關(guān)系圖的構(gòu)造，在組成 de Bruijn 圖時忽略重復(fù)性。假設(shè)若兩條等長片段的對應(yīng)位置堿基差異數(shù)與片段長度的比值小于 d，則兩條片段是同一位置的結(jié)果。1.030取 d =，計算與前一假設(shè)類似。因為直接測量的結(jié)果有較多錯誤，需要設(shè)定標(biāo)準(zhǔn)過濾不正確的數(shù)據(jù)。假設(shè)同一位置的正確測量遠多于錯誤測量。否則將無法推斷出正確的序列。三、算法DeBruijn 圖的一般算法如下圖所示2：取所有讀長的長度為 k 的子串（稱為 k-子串），并確定各子串的前驅(qū)/后繼關(guān)系，有向圖。然后簡化 de Bruijn圖，合并無分支的路徑，并根據(jù)一定規(guī)則過濾錯誤的支路。最后遍歷有向圖，求解最優(yōu)路徑。2圖 1de Bruijn 圖

5、算法當(dāng)數(shù)據(jù)量較大時，枚舉并比較所有的 k-子串耗時較長。注意到上圖中兩個 k-子串相匹配時，將這兩個子串后移一位，得到兩個新的子串，它們繼續(xù)匹配的概率也是相當(dāng)大的。為此設(shè)計有向圖生成算法如下：首先創(chuàng)建一個散列表 Table（鍵不重復(fù)）。（散列表是一種類似于字典的結(jié)構(gòu)，可以向其中快速地或查找數(shù)據(jù)，但占用空間較大。如果限定鍵不相同，的序列 s2，從左到右遍歷其 k-則第二次相同的數(shù)據(jù)會失敗。）對所有子串，并嘗試將其Table（附上子串的位置信息）。如果失敗，則表明該 k-子串與之前的子串重復(fù)，從 Table這個子串所屬的序列 s1 和偏移量。這時停止錄入，在 s1 和 s2 中將該 k-子串向兩側(cè)

6、擴展，尋找最大的匹配子串 overlap。從 s2 中刪除 overlap 段。此時需考慮以下幾種情況：如果 s2 在 overlap 的右側(cè)仍有字符，則創(chuàng)建一個新的序列 s3，其中包3含 overlap 的右側(cè)(k-1)個字符，以及 s2 右側(cè)的剩余部分。將 s3序列庫中，稍候與的序列一同處理。如果 s2 以 overlap 結(jié)尾，則刪除 overlap 后不做任何操作。如果 s2 在 overlap 的左側(cè)有字符（一般可認為是測錯的字符），則應(yīng)“回退”s2 至與 s1前的狀態(tài)，即不斷從 Table 中刪除 s2 的最右側(cè) k-子串并刪除 s2 的最右側(cè)字符，直到 s2 只含 overlap

7、 的前(k-1)個字符。如果 s2 以 overlap 開頭，則在右側(cè)處理（可能生成 s3）之后，直接刪除全部 s2，并將所有對 s2 的應(yīng)位置。（前驅(qū)/后繼，見下）轉(zhuǎn)移至 s1 的對同時，應(yīng)序列中寫入前驅(qū)/后繼關(guān)系：s2s1s3（意為：s2/s1 之后某位置可用 s1/s3 替換），并標(biāo)明位置。根據(jù) overlap 的長度，適當(dāng)增加 s1 的權(quán)重 (見下)。經(jīng)過如上操作，可以保證至當(dāng)前為止的所有 k-子串都有唯一的歸屬。例如：s1 = AACA ATTs2 =G則處理后生成：（設(shè) k=3）s1 s2 s3=AACCGTCCA GCCTCATT這樣，無論是 AAC、GCC、TCA 等邊緣位置還

8、是 overlap 中的 CGT 等，都有唯一確定的從屬序列和位置偏移量。另外，這種方式最大限度地保證了 s1 和 s2 的左側(cè)不被改變。實際處理時 s1是已經(jīng)完整錄入的序列，s2 是從左側(cè)處理至當(dāng)前狀態(tài)的序列，這樣的保證可以減少對已完成處理的子串進行回退的不必要操作。當(dāng)所有的讀入序列和上述 s3 序列均被處理之后，一個各元素長度不均等的圖便建成了。理想情況下（如果隨機性的假設(shè)成立），整個圖是連通的，而且每個片段都只有一個前驅(qū)與一個后繼（除 2 端點）。在實際情況中，序列間有重復(fù)，而且錯誤的序列也會被編入 de Bruijn 圖中。這使得圖中有不止一條可能的完整路徑，而且各個位置都可能有完全不

9、正確、甚至錯誤地與其他片段相連的分支。因此，對圖中每個元素循環(huán)，觀察它的后繼。如果存在 2 個后繼使連接好的序列近似相同（根據(jù)假設(shè)，堿基差異數(shù)與片段長度1.030的比值小于 d =），則認為這兩個后繼是等價的。此時根據(jù)假設(shè)，認為權(quán)重越大正確率越高，須比較各后繼的權(quán)重，以權(quán)重較高的后繼為標(biāo)準(zhǔn)，修正權(quán)重較低的后繼。此時相對準(zhǔn)確率較高的de Bruijn 圖已構(gòu)造完成。下一步應(yīng)是尋找最優(yōu)的路徑，使其通過盡可能多的支路。但以問題二中數(shù)據(jù)測試時發(fā)現(xiàn)，需要枚舉的數(shù)據(jù)量過大，不太可能短時間內(nèi)在使用的測試系統(tǒng)上完成。因此，本文中的程序使用了簡化的處理方式：根據(jù)假設(shè)，認為大段 重復(fù)對結(jié)果的影響是可以忽略的。這樣

10、 de Bruijn 圖一定是如圖 1 所示的網(wǎng)狀結(jié)構(gòu)。同時，由獨立性的假設(shè)，在每個分支點決策的優(yōu)劣，對后續(xù)各分支點決策是沒有影響的。所以，從頭枚舉所有的支路，但在每個分支片段緩存決策的結(jié)果，這樣可以將指數(shù)級的時間復(fù)雜度降為線性。4CCGTC CCGTC其中：“從頭”的起始點并不能被明確的找出（根據(jù)測試，圖中有大量沒有前驅(qū)的片段，這可能是錯誤造成的）。這時只有枚舉所有可能的起始點。為便于分析，將整個子串空間劃分為若干個連通空間，使每個子空間中各子串都是以有限次前驅(qū)/后繼方式相連通的，而子空間的任意兩子串均不連通。對每個子空間，枚舉從所有點開始的最優(yōu)路徑（緩存各分支點決策）。為了簡化計算，將判斷

11、決策與路徑優(yōu)劣的標(biāo)準(zhǔn)定為其包含子串的個數(shù)。因為這里的子串都是原的變化也不大。段直接或切分得到，子串的長度并不長，連接偏移量一旦獲得了最優(yōu)的路徑，便可以根據(jù) de Bruijn 圖生成時各片段的連接位置信息，將路徑上所有片段連接成一條長鏈。對所有子空間的最優(yōu)路徑進行連接，就得到了此算法的（多條）輸出。四、測試、敏感度分析本程序的 Test 函數(shù)生成偽隨機數(shù)序列用于測試（隨機數(shù)是可重現(xiàn)的）。測試生成長度為 n 的序列，并截取 c 次長為 88 的子序列作為讀長。同時設(shè)置一個隨機分布 dist3，當(dāng)這個值隨機產(chǎn)生 1 時，某一序列的某一堿基將被隨機的 ACGT之一替換（模擬錯誤）。測試發(fā)現(xiàn)，當(dāng)序列完

12、全準(zhǔn)確時（dist3 不產(chǎn)生 1），即使原始序列多達 200000個，截取 100000 次，拼接成的序列也幾乎與原序列完全相同。當(dāng) dist3 產(chǎn)生 1概率為千分之一時，結(jié)果也比較準(zhǔn)確（除開頭結(jié)尾部分，這些位置很可能測點極少）。但是，dist3 錯誤率上升至 1%時，結(jié)果開始變得不穩(wěn)定，拼接總長度也開始變短。在這種情況下，每個子序列堿基錯誤個數(shù)期望近似為 0.88 個，即幾乎每個子序列都有錯誤，導(dǎo)致算法的準(zhǔn)確度下降。對問題二的計算，本算法的效果也并不十分理想。計算產(chǎn)生了數(shù)百條無法繼續(xù)連接的序列，其中最長的也僅有約 40000 個堿基，距離總長 120000 還是有很大差距的。這可能是由于真實

13、情況與假設(shè)的偏離（重復(fù)區(qū)段被視為同一區(qū)域，某些位置錯誤率過高）導(dǎo)致的。所以，本算法受重復(fù)與錯誤的影響還是很大的。為最大限度地節(jié)省時間，將計算時間控制在測試系統(tǒng)可接受的范圍內(nèi)，后期對 de Bruijn 圖的處理并不是全局最優(yōu)的，也沒有對重復(fù)的情形做特殊的處理。本算法能夠基于前面相同的 k-子串發(fā)現(xiàn)很多錯誤，但更隱蔽、短距離內(nèi)多個或重復(fù)多次發(fā)生的錯誤，還是會對算法起到不利的影響，降低處理效率與準(zhǔn)確度。五、結(jié)語本文嘗試構(gòu)造了一種較為直接快速的基于 de Bruijn 圖的拼接算法，并實現(xiàn)其程序，對隨機構(gòu)造的序列以及問題二中給出的序列進行了測試與分析。這一算法執(zhí)行較快，但準(zhǔn)確率等方面仍有改進的空間。

14、例如，使用 2 個二進制位（而不是字符）代表一個堿基以節(jié)省內(nèi)存使用與 I/O 時間3，優(yōu)化圖的處理以適應(yīng)較復(fù)雜的數(shù)據(jù)/較好的執(zhí)行條件，圖生成后做進一步與處理以降低錯誤率與遍歷次數(shù)，等等。5附錄一：程序源代碼使用 C+11 編寫，在 Windows 8.1 上用 Visual Studio 2013 編譯通過。使用了boost 開源庫4。#include #include #include #include #include #include #include #include #include #include #include #include #include#include #inclu

15、destd:ofstream fout(out.txt, std:ios:out| std:ios:trunc);struct xSequencing voidAdd(conststd:string&src)Contigs.push_back(*newxContig(src); std:list Pros() unsigned contig_sz = Contigs.front().Str.size(); unsigned sz = contig_sz / 3;xContig& last_contig = Contigs.back(); for (auto p = Contigs.begin(

16、);) auto autonext = std:next(p);is_last = (&*p = &last_contig);ProsContig(*p, sz); /* may delete p */if (is_last) break;elsep = next;6std:list result;while (!Contigs.empty() boost:rusive:list group;MoveContigToGroup(group, Contigs.front(); if (group.size() 1) result.push_back(xPathResolver(group).Re

17、solve();group.clear_and_dise(std:default_delete();return std:move(result);private:struct xContig;voidMoveContigToGroup(boost: xContig& contig) if (!contig.Flag) return;rusive:list&out,Contigs.erase(Contigs.iterator_to(contig); out.push_back(contig);contig.Flag = false;for (auto& p : contig.Prevs) Mo

18、veContigToGroup(out, *p); for (auto& p : contig.Nexts) MoveContigToGroup(out, *p.);void ProsContig(xContig& contig, unsigned sz) xStr& str = contig.Str;for (auto p = str.begin(); p second.Offset,-second.Contig != value.Contig) -second.Contig-Overlap(this,ins_res.value.Contig, value.Offset, sz);break

19、;7private:typedef std:string xStr;template sic xRoe(x return (val (sizeof(x)*CHAR_BIT - 1);struct xStrHash std:size_t operator()(const xStr& s) conststd:size_t for (auto cv = Rov = 0;: s) e(Roe(v);for (unsigned i = 0; i = 3; +i) v = i;break;if(c=ACGTi)return std:hash()(v);struct xContig;struct xCont

20、igLocation xContig* Contig; unsigned Offset;struct xContigNextInfo unsigned Offset1, Offset2;public:xContigNextInfo(unsigned Offset1, unsigned: Offset1(Offset1), Offset2(Offset2) ;Offset2)struct xContig : boost:rusive:list_base_hook xContig(xStr s) : Str(std:move(s) void Overlap(xSequencing* parent,

21、 unsigned offset1, xContig*other,8unsigned offset2, unsigned sz) xStr& str1 = this-Str, str2 = other-Str;xContig* new_contig = nullptr;unsigned_diff_offset2;/ compare rightunsigned compare_sz = std:min(str1.size() - offset1, str2.size() - offset2);_diff_offset2 = offset2 + compare_sz;for (unsigned i

22、 = offset1 + sz, j = offset2 + sz; i str1.size() & j str2.size(); +i, +j) if (str1i != str2j) _diff_offset2 = j; break;if (_diff_offset2 Contigs.push_back(*new xContig(str2.substr(_diff_offset2 + 1 - sz);new_contig = &parent-Contigs.back();Nexts.insert( new_contig, offset2,sz - 1 );new_contig-Prevs.

23、insert(this);_diff_offset2+offset1-unsigned_overlap_offset2;/ compare leftunsigned min_offset = std:min(offset1, offset2);_overlap_offset2 = offset2 - min_offset;unsigned i = min_offset;while (i- 0) / i goes to 0 (from min_offset - 1) if (str1i + offset1 - min_offset != str2i +min_offset) offset2-_o

24、verlap_offset2=i+offset2-9min_offset + 1;break;for (unsigned i =_overlap_offset2; i Table.erase(str2.substr(i, sz);if (_overlap_offset2 0) other-Str.erase( unsigned_overlap_offset2_overlap_offset2 + sz - 1);_overlap_offset1=+ offset1 - offset2; other-Nexts.insert( this, _overlap_offset2 +sz- 1,_over

25、lap_offset1 + sz - 1 ); Prevs.insert(other); else other-Eaten(parent, this, offset1 - offset2);Weight +=_diff_offset2 -_overlap_offset2;if (new_contig) parent-ProsContig(*new_contig, sz);private:void Eaten(xSequencing* parent, xContig* eater, unsigned offset)for (auto* prev : Prevs) auto p = prev-Ne

26、xts.find(this);if (p != prev-Nexd() auto value = p-second; prev-Nexts.insert(value.Offset2 + offset );eater,value.Offset1,eater-Prevs.insert(prev); prev-Nexts.erase(p);eater-Weight += Weight;auto iter = parent-Contigs.iterator_to(*this);parent-Contigs.erase_and_dis std:default_delete();e(iter,10priv

27、ate:friend xSequencing; xStr Str;std:set Prevs; std:map Nexts; bool Flag = true;unsigned Weight = 0;std:pair Deci= nullptr, 0 ;struct xPathResolver xPathResolver(boost:Group(Group) xStr Resolve() rusive:list&Group):for (auto p = Group.begin(); p != Group.end(); std:liststd:pair strs; for (auto q0 =

28、p-Nexts.begin(); q0 != p-Nex+p) d(); +q0)strs.push_back( q0-, Merge(&*p,*q0) );if (strs.size() = 2) strs.sort(const std:pairxContig*,const std:pair& r)xStr& l,return l.-Weight r.-Weight;);for std:next(strs.begin();(autoq1=strs.begin(),q2=q2 != strs.end(); +q1, +q2) if (IsSimilar(q1-second, q2-second

29、) if (q2-second.size() second.size() q2-second.resize(q1-second.size();std:copy(q1-second.begin(), q1-second.end(), q2-second.begin();auto info = p-Nexts.at(q2-);-Str.begin()q2-Str.replace(q2-+ info.Offset2,q2-Str.end(),q2-second.substr(info.Offset1);unsigned max_depth = 0;11xContig*= nullptr;for (a

30、uto& c : Group) if (!c.Flag) DecidePath(&c);if (c.Deci.second max_depth) max_depth = c.Deci= &c;.second;return MergePath();void RemoveFlag(xContig* p) for (;) p-Flag = false; std:pair&pa= p-Deci;if (!(p = pa.) break;xStr MergePath(xContig* p) auto s = unsigned for (;)p-Str; offset = 0;std:pair&pa= p

31、-Deci;if (!pa. auto& c2 s.repla) break;*p-Nexts.find(pa.=);+.begin()+offsetc2.second.Offset1,s.end(),c2. offset += p = pa.-Str.substr(c2.second.Offset2); c2.second.Offset1 - c2.second.Offset2;return s;void DecidePath(xContig* contig) if (contig-Flag+) contig-Deci return;= nullptr, 100 ;12if (contig-

32、Deci.second) else if (contig-Nexts.empty() contig-Deci else unsigned max_v xContig* deci= nullptr, 100 ;= 0;= nullptr;for (auto& next_pa : contig-Nexts)if (!next_pa.-Flag) DecidePath(next_pa. auto v = next_pa. if (v max_v) );-Deci.second+ 1;max_v deci=v;= next_pa.;contig-Deci if (max_v = deci, max_v

33、 ;0) contig-Deci= nullptr, 100 ;contig-Flag = false;xStrMerge(xContig*c1,conststd:pair& c2) auto s = s.replac2. return s;c1-Str;.begin() +c2.second.Offset1, s.end(),-Str.substr(c2.second.Offset2);bool IsSimilar(const xStr& l, const xStr& r) auto sz = std:min(l.size(), unsigned c = 0;for (unsigned i

34、= 0; i sz; if (li != ri)+c;return double(c) / sz = 1.0r.size();+i)/ 30;private:boost:;rusive:list&Group;13public:xSequencing() Contigs.clear_and_dise(std:default_delete();private:boost:rusive:list Contigs;std:unordered_map Table;void Test() std:minstd_rand rd;std:uniform_distributiondist(0, 3);unsigned n = 200000; unsigned c = 100000;std:uniform_ std:uniform_distribution_distributiondist