本發(fā)明屬于拷貝數(shù)變異技術(shù)領(lǐng)域，尤其涉及一種基于狀態(tài)轉(zhuǎn)移模型的新一代測序拷貝數(shù)變異仿真方法。
背景技術(shù)：
拷貝數(shù)變異是由基因組發(fā)生重排而導(dǎo)致的，一般指長度為1kb以上的基因組大片段的拷貝數(shù)增加或者減少，主要表現(xiàn)為亞顯微水平的缺失(deletion)和重復(fù)(insertion)。拷貝數(shù)變異是基因組結(jié)構(gòu)變異(Structural variation，SV)的重要組成部分，它的位點突變率遠高于SNP(Single nucleotide polymorphism)，是人類疾病的重要致病因素之一。按照發(fā)生場所的不同，可以將拷貝數(shù)變異分為生殖細胞拷貝數(shù)變異(Copy number variation，CNV)和體細胞拷貝數(shù)變異(Copy number alternation，CNA)，顧名思義就是他們的發(fā)生場所分別為生殖細胞和體細胞。CNV具有遺傳效應(yīng)，CNA沒有遺傳效應(yīng)，這是由它們所處細胞的機制所決定的?？截悢?shù)變異的仿真就是設(shè)計仿真算法，用程序?qū)截悢?shù)變異的過程進行模擬。近年來，基因組測序領(lǐng)域發(fā)展迅速，這非常有助于對許多生物系統(tǒng)的理解。在過去的五年中，計算機生物學(xué)家和生物信息學(xué)專家針對發(fā)現(xiàn)、分析和解釋不同的基因組變異的高通量測序數(shù)據(jù)，提出了新的、更好的和更有效的檢測拷貝數(shù)變異的工具。在使用檢測工具時，可靠的模擬數(shù)據(jù)集是必不可少的，模擬數(shù)據(jù)的獲得是測試新開發(fā)檢測工具的第一步。雖然目前已經(jīng)有很多可用的拷貝數(shù)變異仿真工具，但是這些工具的功能都不是很全面，要么就是只可以模擬CNV和CNA中的一個功能，要么就是沒有一個可信的狀態(tài)轉(zhuǎn)移模型，要么就是只有序列生成部分。因此，開發(fā)一個有效的關(guān)于CNV和CNA的模擬器和序列生成器是必要的，它要能夠模擬拷貝數(shù)變異且考慮到真實生物樣品的錯誤率。不同的下一代測序儀所生成的reads的length和error profile也不同，目前最流行的測序數(shù)據(jù)是從Illumina測序平臺所產(chǎn)生的，它采用了化學(xué)方法來進行序列合成并生成reads，要開發(fā)的高效模擬器正是基于Illumina測序平臺的。正是由于Illumina平臺產(chǎn)生的數(shù)據(jù)受歡迎且應(yīng)用廣泛這個特點，任何其它的測序平臺通過提供一個特定的錯誤配置文件就可以對其進行使用。目前可用的仿真軟件可以生成基于特定平臺的相關(guān)錯誤配置文件的reads，也可以跨平臺生成reads。已經(jīng)存在的一些仿真軟件都有各自的優(yōu)點，但是同時它們也存在著一些缺陷?，F(xiàn)有的仿真軟件最大的缺陷是不能能同時仿真CNV和CNA，下面針對一些仿真軟件的性質(zhì)和功能分別加以說明。SInc是用C語言開發(fā)的一款仿真軟件，它是開源的，擁有CLI接口，也有自己的error model，但是它存在的問題沒有仿真CNA變異；MetaSim是用JAVA語言開發(fā)的一款仿真軟件，它擁有CLI和GUI接口，不是開源的，可以仿真pair-end數(shù)據(jù)，它的缺陷是只有序列生成部分而沒有變異仿真部分，即沒有將quality value賦值給reads；FlowSim是用Haskell語言編寫的，它擁有CLI接口，是一款開源仿真軟件，它的缺陷是沒有變異仿真部分且不能仿真Illumina平臺的數(shù)據(jù)，即不能仿真paie-end數(shù)據(jù)；GenFrag仿真軟件是開源的，有CLI接口，它的缺點是沒有變異仿真部分且它的erroe model過分簡單；DwgSim仿真軟件是由變異仿真和序列生成兩個部分構(gòu)成的，它有CLI接口且是開源的，它的缺陷是不能模擬真實數(shù)據(jù)。當然，這些仿真軟件共有的一個缺陷是不能仿真CNA變異。
現(xiàn)有拷貝數(shù)變異仿真軟件的實現(xiàn)方法存在存在以下問題：沒有將生殖細胞和體細胞的兩種拷貝數(shù)變異集成在一起，即只可以仿真單個的CNV，不能同時仿真CNA,這就導(dǎo)致了仿真軟件的功能不全面，仿真出來的數(shù)據(jù)比較片面，不是完整的拷貝數(shù)變異后的數(shù)據(jù)，限制了用戶的使用；沒有合適的model來確定拷貝數(shù)變異的各個狀態(tài)，因為拷貝數(shù)變異中各個狀態(tài)之間的轉(zhuǎn)換是符合某種轉(zhuǎn)換機制的，要是沒有相應(yīng)的model的話，仿真出來的數(shù)據(jù)和真實數(shù)據(jù)的差異較大，仿真結(jié)果的可信度將受到影響；沒有將拷貝數(shù)變異的兩種形式變異仿真和序列生成集成在一塊，一般存在這種問題的仿真軟件大都只有序列生成的部分，沒有變異仿真的部分，即它能生成最終的fq文件，但是將生物變異的部分省略，這樣顯然是不符合實際情況的，因為拿到的真實樣本不一定是完全沒有發(fā)生拷貝數(shù)變異的，相反發(fā)生變異的比例還很大，所以加上變異仿真這一步是相當有必要的。
技術(shù)實現(xiàn)要素：
本發(fā)明的目的在于提供一種基于狀態(tài)轉(zhuǎn)移模型的新一代測序拷貝數(shù)變異仿真方法，旨在解決為拷貝數(shù)變異檢測提供合適的模擬數(shù)據(jù)的問題。
本發(fā)明是這樣實現(xiàn)的，一種基于狀態(tài)轉(zhuǎn)移模型的新一代測序拷貝數(shù)變異仿真方法，所述基于狀態(tài)轉(zhuǎn)移模型的新一代測序拷貝數(shù)變異仿真方法采用拷貝數(shù)變異仿真算法；在仿真算法中增加狀態(tài)轉(zhuǎn)移模型和序列生成部分；
所述拷貝數(shù)變異包括CNV和CNA；
基于Illumina測序平臺的Profile文件的生成，核心步驟是將fq文件的reads說明部分的ASCii碼轉(zhuǎn)換成堿基的quality value，相應(yīng)方法是對應(yīng)字符的ASCii碼減去33；
將仿變異真后的fa文件和生成的profile文件作為輸入，設(shè)置read length，利用多線程和序列生成算法，生成并輸出最終的fq文件。
進一步，所述CNV仿真算法和包括：
(a)確定發(fā)生CNV變異的位置、尺寸、類型；
(b)根據(jù)a中確定的CNV變異的參數(shù)執(zhí)行CNV變異，并打印變異參數(shù)的記錄文件和變異后的fa文件。
進一步，所述CNV狀態(tài)轉(zhuǎn)移模型為：
Normal:
Paa＝Pa Pnn＝Pn Pdd＝Pd
Pa＝Paa*Pnn*Pdd/(2-Paa*Pnn*Pdd)
Pd＝(1-Pa)*Pnn
Pn＝1-Pa-Pd
Insertion:
Paa＝Pa Pnn＝Pn Pdd＝Pd
Pd＝Paa*Pnn*Pdd/(2-Paa*Pnn*Pdd)
Pd＝(1-Pd)*Paa
Pa＝1-Pn-Pd
Deletion:
Paa＝Pa Pnn＝Pn Pdd＝Pd
Pn＝Paa*Pnn*Pdd/(2-Paa*Pnn*Pdd)
Pd＝(1-Pn)*Pdd
Pd＝1-Pa-Pn。
進一步，所述CNA仿真算法包括：
(a)確定發(fā)生CNA變異的位置、尺寸、類型；
(b)根據(jù)a中確定的CNA變異的參數(shù)，執(zhí)行CNA變異，并打印變異參數(shù)的記錄文件和變異后的fa文件。
進一步，所述CNA狀態(tài)轉(zhuǎn)移模型為：
Normal:
Paa＝Pa Pnn＝Pn Pdd＝Pd
Pa＝Paa*Pnn*Pdd/(2-Paa*Pnn*Pdd)
Pd＝(1-Pa)*Pnn
Pn＝1-Pa-Pd
Insertion:
Paa＝Pa Pnn＝Pn Pdd＝Pd
Pd＝Paa*Pnn*Pdd/(2-Paa*Pnn*Pdd)
Pd＝(1-Pd)*Paa
Pa＝1-Pn-Pd
Deletion:
Paa＝Pa Pnn＝Pn Pdd＝Pd
Pn＝Paa*Pnn*Pdd/(2-Paa*Pnn*Pdd)
Pd＝(1-Pn)*Pdd
Pd＝1-Pa-Pn。
本發(fā)明的另一目的在于提供一種應(yīng)用所述基于狀態(tài)轉(zhuǎn)移模型的新一代測序拷貝數(shù)變異仿真方法的CNV和CNA模擬器。
本發(fā)明的另一目的在于提供一種應(yīng)用所述基于狀態(tài)轉(zhuǎn)移模型的新一代測序拷貝數(shù)變異仿真方法的CNV和CNA序列生成器。
本發(fā)明提供的基于狀態(tài)轉(zhuǎn)移模型的新一代測序拷貝數(shù)變異仿真方法，能夠解決現(xiàn)有拷貝數(shù)變異的仿真程序只考慮了CNV或CNA的情況，將CNA和CNV集成在同一個仿真軟件中，并擁有自己獨特的拷貝數(shù)變異(包括CNV和CNA)仿真算法；在CNV和CNA仿真算法的基礎(chǔ)上還增加了狀態(tài)轉(zhuǎn)移模型，有了這個model之后，整個仿真的過程就更加具有可信度。沒有model的時候，拷貝數(shù)變異類型是由程序設(shè)定按一定比例產(chǎn)生的，一般是發(fā)生缺失的比例與發(fā)生插入的比例為4:1：有狀態(tài)轉(zhuǎn)移模型的時候，下一狀態(tài)發(fā)生變化的類型和上一狀態(tài)有關(guān)，至于有什么關(guān)系，主要取決于model的設(shè)定。下面分別是沒有狀態(tài)轉(zhuǎn)移模型和有狀態(tài)轉(zhuǎn)移模型時候拷貝數(shù)變異記錄文件的對比圖，由對比圖可以發(fā)現(xiàn)，沒有狀態(tài)轉(zhuǎn)移模型的時候，變異類型狀態(tài)之間的轉(zhuǎn)換服從穩(wěn)定的4:1的比例，但是加上狀態(tài)轉(zhuǎn)移模型之后，變異狀態(tài)之間的比例不一定服從4:1，這更符合實際數(shù)據(jù)，因為真實數(shù)據(jù)發(fā)生變異的情況受環(huán)境等多種因素影響，不可能服從一個穩(wěn)定的變化比率，這需要反復(fù)訓(xùn)練真實數(shù)據(jù)來獲得。在CNA和CNV及其model的基礎(chǔ)之上，增加了序列生成部分，使得這個仿真軟件不僅擁有變異仿真功能，還擁有序列生成功能。
本發(fā)明采用狀態(tài)轉(zhuǎn)移模型，將狀態(tài)轉(zhuǎn)移模型加在拷貝數(shù)變異的仿真算法中，在此基礎(chǔ)上，還應(yīng)加上序列生成的功能；在實現(xiàn)仿真拷貝數(shù)變異的過程中，不僅仿真了生殖細胞的拷貝數(shù)變異，還仿真了體細胞的變異。
本發(fā)明在仿真CNV和CNA的基礎(chǔ)上，還加入了序列生成的部分，模擬了序列生成部分中可能發(fā)生的錯誤機制，增加了數(shù)據(jù)的真實性。綜合了CNV和CNA于一個仿真軟件中，方便了用戶的使用和研究，并且仿真算法中加入了狀態(tài)轉(zhuǎn)移模型，使得生物變異的仿真更加具有可信度。除此之外，本發(fā)明還加入了序列生成部分，考慮了真實測序中可能發(fā)生的錯誤，生成錯誤配置文件，進而生成最終的reads?？偠灾?，和其它仿真軟件相比，本拷貝數(shù)變異仿真軟件功能更加完整，數(shù)據(jù)更加貼近真實數(shù)據(jù)。下面三幅圖是本發(fā)明最終的成果圖，由圖可以看到本仿真軟件有仿真CNV和CNA的功能，有生成錯誤配置文件的功能，有生成最終的reads的功能及一些其它的附加功能，功能完善齊全，還有完整的GUI(圖形用戶界面)，使用起來方便快捷。
附圖說明
圖1是本發(fā)明實施例提供的基于狀態(tài)轉(zhuǎn)移模型的新一代測序拷貝數(shù)變異仿真方法流程圖。
具體實施方式
為了使本發(fā)明的目的、技術(shù)方案及優(yōu)點更加清楚明白，以下結(jié)合實施例，對本發(fā)明進行進一步詳細說明。應(yīng)當理解，此處所描述的具體實施例僅僅用以解釋本發(fā)明，并不用于限定本發(fā)明。
下面結(jié)合附圖對本發(fā)明的應(yīng)用原理作詳細的描述。
如圖1所示，本發(fā)明實施例提供的基于狀態(tài)轉(zhuǎn)移模型的新一代測序拷貝數(shù)變異仿真方法包括以下步驟：
S101:CNV仿真算法和狀態(tài)轉(zhuǎn)移模型的設(shè)計；
S102:CNA仿真算法和狀態(tài)轉(zhuǎn)移模型的設(shè)計；
S103：基于Illumina測序平臺的Profile文件的生成，核心步驟是將fq文件的reads說明部分的ASCii碼轉(zhuǎn)換成堿基的quality value，相應(yīng)方法是對應(yīng)字符的ASCii碼減去33；
S104:將仿變異真后的fa文件和生成的profile文件作為輸入，設(shè)置合適的read length，利用多線程和序列生成算法，生成并輸出最終的fq文件。
下面結(jié)合具體實施例對本發(fā)明的應(yīng)用原理作進一步的描述。
本發(fā)明的實施例在拷貝數(shù)變異生物特性及仿真算法的基礎(chǔ)之上，建立狀態(tài)轉(zhuǎn)移模型，經(jīng)過反復(fù)對仿真數(shù)據(jù)進行訓(xùn)練，設(shè)置合適的狀態(tài)轉(zhuǎn)移模型，對仿真算法進行改進。
本發(fā)明實施例的技術(shù)方案
(1)CNV仿真算法和狀態(tài)轉(zhuǎn)移模型的設(shè)計
CNV仿真算法:
(a)確定發(fā)生CNV變異的位置、尺寸、類型；
(b)根據(jù)a中確定的CNV變異的參數(shù)執(zhí)行CNV變異，并打印變異參數(shù)的記錄文件和變異后的fa文件。
CNV狀態(tài)轉(zhuǎn)移模型：
Normal:
Paa＝Pa Pnn＝Pn Pdd＝Pd
Pa＝Paa*Pnn*Pdd/(2-Paa*Pnn*Pdd)
Pd＝(1-Pa)*Pnn
Pn＝1-Pa-Pd
Insertion:
Paa＝Pa Pnn＝Pn Pdd＝Pd
Pd＝Paa*Pnn*Pdd/(2-Paa*Pnn*Pdd)
Pd＝(1-Pd)*Paa
Pa＝1-Pn-Pd
Deletion:
Paa＝Pa Pnn＝Pn Pdd＝Pd
Pn＝Paa*Pnn*Pdd/(2-Paa*Pnn*Pdd)
Pd＝(1-Pn)*Pdd
Pd＝1-Pa-Pn
(2)CNA仿真算法和狀態(tài)轉(zhuǎn)移模型的設(shè)計
CNA仿真算法:
(a)確定發(fā)生CNA變異的位置、尺寸、類型；
(b)根據(jù)a中確定的CNA變異的參數(shù)，執(zhí)行CNA變異，并打印變異參數(shù)的記錄文件和變異后的fa文件。
CNA狀態(tài)轉(zhuǎn)移模型：
Normal:
Paa＝Pa Pnn＝Pn Pdd＝Pd
Pa＝Paa*Pnn*Pdd/(2-Paa*Pnn*Pdd)
Pd＝(1-Pa)*Pnn
Pn＝1-Pa-Pd
Insertion:
Paa＝Pa Pnn＝Pn Pdd＝Pd
Pd＝Paa*Pnn*Pdd/(2-Paa*Pnn*Pdd)
Pd＝(1-Pd)*Paa
Pa＝1-Pn-Pd
Deletion:
Paa＝Pa Pnn＝Pn Pdd＝Pd
Pn＝Paa*Pnn*Pdd/(2-Paa*Pnn*Pdd)
Pd＝(1-Pn)*Pdd
Pd＝1-Pa-Pn
(3)Profile文件的生成
本發(fā)明的測序數(shù)據(jù)是基于Illumina測序平臺，所以生成profile文件時所用的fq文件也是由Illumina測序平臺產(chǎn)生。Profile文件其實是統(tǒng)計序列中某個堿基出現(xiàn)的次數(shù)，是根據(jù)它的quality value來判斷它出現(xiàn)的次數(shù)的，所以在生成profile文件之前應(yīng)該將fq文件中reads的序列說明部分的ASCii碼轉(zhuǎn)換成堿基的quality value，相應(yīng)的方法是對應(yīng)字符的ASCii碼減去33。
(4)Reads的生成
將仿變異真后的fa文件和生成的profile文件作為輸入，設(shè)置合適的read length，利用多線程和序列生成算法，生成并輸出最終的文件。
(5)算法的性能評估
對設(shè)計的仿真算法和狀態(tài)轉(zhuǎn)移模型進行評估和改進，對序列生成過程進行優(yōu)化，形成最終的仿真軟件。
以上所述僅為本發(fā)明的較佳實施例而已，并不用以限制本發(fā)明，凡在本發(fā)明的精神和原則之內(nèi)所作的任何修改、等同替換和改進等，均應(yīng)包含在本發(fā)明的保護范圍之內(nèi)。