基于C程序的DNA序列的kmerindex數(shù)據(jù)查找

1、論文題目：基于c程序的dna序列的k-mer index數(shù)據(jù)查找姓名 學(xué)院 年級(jí) 專業(yè) 學(xué)號(hào) 聯(lián)系電話 數(shù)學(xué)分析 高等代數(shù) 高等數(shù)學(xué) 線性代數(shù) 概率統(tǒng)計(jì) 數(shù)學(xué)實(shí)驗(yàn) 數(shù)學(xué)模型 cet4 cet6 電氣工程學(xué)院2013 級(jí)電氣工程及其自動(dòng)化918996560電氣工程學(xué)院2013 級(jí)電氣工程及其自動(dòng)化918885578540電氣工程學(xué)院2013 級(jí)電氣工程及其自動(dòng)化858989554基于c程序的dna序列的k-mer index數(shù)據(jù)查找 摘要dna 是生命體的基本遺傳物質(zhì)，其組成和序列變化創(chuàng)造了形形色色的生命世界。快速、準(zhǔn)確地獲取生物體的遺傳信息對(duì)于生命科學(xué)的研究具有重要意義1?，F(xiàn)需要給定一種數(shù)據(jù)索

2、引方法 ，利用一種查詢算法查詢百萬(wàn)條序列中是否存在相應(yīng)的片段，如果存在，則輸出相應(yīng)片段所在的位置。針對(duì)問題一，運(yùn)用karp-rabin算法，在c程序環(huán)境下編寫字符串匹配算法。具體做法是將堿基序列映射成四進(jìn)制的數(shù)串，對(duì)給定的k，構(gòu)造合適的哈希函數(shù)，將四進(jìn)制數(shù)串內(nèi)每個(gè)長(zhǎng)度為k的子數(shù)串譯為唯一的十進(jìn)制數(shù)，按順序放進(jìn)索引數(shù)組（哈希表）。查找相同的字符串等價(jià)于判斷相應(yīng)的hash值是否相同。此法可以大大提高建立索引和查詢的時(shí)間。針對(duì)問題二，對(duì)不同k值經(jīng)過大量多次的嘗試，一般來說建立索引約10秒，查詢約0.3秒（visual c+ 6.0運(yùn)行環(huán)境）。針對(duì)問題三和問題四，本算法使用for循環(huán)函數(shù)，先計(jì)算出建立

3、索引與使用索引的算法中每一個(gè)語(yǔ)句的執(zhí)行次數(shù)，然后再相加，最后依據(jù)去低階項(xiàng)，去掉常數(shù)項(xiàng)，去掉高階項(xiàng)的常參的原則得到時(shí)間復(fù)雜度。然后根據(jù)算法臨時(shí)存儲(chǔ)的空間大小來計(jì)算空間復(fù)雜度，觀察有無(wú)臨時(shí)存儲(chǔ)大小以及臨時(shí)存儲(chǔ)大小與輸入變量的關(guān)系。經(jīng)分析，本算法在建立索引時(shí)的時(shí)間復(fù)雜度為o（n*m），空間復(fù)雜度為o（n）；在使用索引時(shí)的時(shí)間復(fù)雜度為o（n*m），空間復(fù)雜度為o（1）。針對(duì)問題五，首先分析了整形數(shù)據(jù)最大值的限制，k最大只能取14。為了擴(kuò)展k的最大值，我們提出將字符串分成幾個(gè)不相交的完備的字串，其長(zhǎng)度不超過14，分別比較每一個(gè)字串，然后對(duì)結(jié)果取交集。并給出了這種方法下所需要的內(nèi)存與k的關(guān)系，得出理論上可

4、以在8g的內(nèi)存下支持所有k值的查找。關(guān)鍵詞：dna 數(shù)據(jù)索引 查詢算法 算法復(fù)雜度 karp-rabin 哈希函數(shù)c程序 字符串匹配一、 問題重述這個(gè)問題來自 dna序列的k-mer index問題。給定一個(gè)dna序列，這個(gè)系列只含有4個(gè)字母atcg，如 s =“ctgtactgtat”。給定一個(gè)整數(shù)值k，從s的第一個(gè)位置開始，取一連續(xù)k個(gè)字母的短串，稱之為k-mer（如k= 5，則此短串為ctgta），然后從s的第二個(gè)位置，取另一k-mer（如k= 5，則此短串為tgtac），這樣直至s的末端，就得一個(gè)集合，包含全部k-mer 。 如對(duì)序列s來說，所有5-mer為ctgta，tgtac，gt

5、act，tactg，actgt，tgtat通常這些k-mer需一種數(shù)據(jù)索引方法，可被后面的操作快速訪問。例如，對(duì)5-mer來說，當(dāng)查詢ctgta，通過這種數(shù)據(jù)索引方法，可返回其在dna序列s中的位置為1，6。問題現(xiàn)在以文件形式給定 100萬(wàn)個(gè) dna序列，序列編號(hào)為1-1000000，每個(gè)基因序列長(zhǎng)度為100 。（1）要求對(duì)給定k， 給出并實(shí)現(xiàn)一種數(shù)據(jù)索引方法，可返回任意一個(gè)k-mer所在的dna序列編號(hào)和相應(yīng)序列中出現(xiàn)的位置。每次建立索引，只需支持一個(gè)k值即可，不需要支持全部k值。（2）要求索引一旦建立，查詢速度盡量快，所用內(nèi)存盡量小。（3）給出建立索引所用的計(jì)算復(fù)雜度，和空間復(fù)雜度分析。（

6、4）給出使用索引查詢的計(jì)算復(fù)雜度，和空間復(fù)雜度分析。（5）假設(shè)內(nèi)存限制為8g，分析所設(shè)計(jì)索引方法所能支持的最大k值和相應(yīng)數(shù)據(jù)查詢效率。（6）按重要性由高到低排列，將依據(jù)以下幾點(diǎn)，來評(píng)價(jià)索引方法性能 索引查詢速度索引內(nèi)存使用8g內(nèi)存下，所能支持的k值范圍建立索引時(shí)間二、問題分析針對(duì)問題一，要求對(duì)給定k，給出并實(shí)現(xiàn)一種數(shù)據(jù)索引方法，可返回任意一個(gè)k-mer所在的dna序列編號(hào)和相應(yīng)序列中出現(xiàn)的位置。以及每次建立索引，只需支持一個(gè)k值即可，不需要支持全部k值。首先利用c程序的中的fp文件讀取函數(shù)將目標(biāo)文件中的一億個(gè)單字符數(shù)據(jù)讀出，用于建立數(shù)據(jù)庫(kù)，由于給定的是100萬(wàn)個(gè)dna序列，每個(gè)基因序列長(zhǎng)度為1

7、00。而c程序中靜態(tài)數(shù)組的范圍限制遠(yuǎn)遠(yuǎn)小于100萬(wàn)，于是定義全局變量建立動(dòng)態(tài)數(shù)組，動(dòng)態(tài)分配儲(chǔ)存空間，即建立行為100萬(wàn)，列為100的1000000100的動(dòng)態(tài)數(shù)組，也就是索引的建立過程。采用karp-rabin算法2，可令a=0,c=1,g=2,t=3，將堿基序列映射成一個(gè)四進(jìn)制數(shù)串。對(duì)于給定的k，可以構(gòu)造這樣一個(gè)集合，即這個(gè)集合囊括了四進(jìn)制數(shù)串中所有長(zhǎng)度為k的子數(shù)串。構(gòu)造相應(yīng)的哈希函數(shù)，將每個(gè)四進(jìn)制子數(shù)串轉(zhuǎn)換成唯一的十進(jìn)制數(shù)字。待查kmer也做相同的處理，然后遍歷整個(gè)集合里的十進(jìn)制數(shù)，記錄下與待查kmer對(duì)應(yīng)十進(jìn)制數(shù)相同的位置。由于轉(zhuǎn)換的是唯一的十進(jìn)制數(shù)字，所以只要相等，便找到了相應(yīng)的k-m

8、er所在的dna序列編號(hào)和相應(yīng)序列中出現(xiàn)的位置，然后以txt文件形式輸出所有找到的dna序列編號(hào)以及相應(yīng)的位置。針對(duì)問題二，要求索引一旦建立，查詢速度盡量快，所用內(nèi)存盡量小。首先要求查詢速度非?？?，所以對(duì)算法的要求就比較高，要求算法的時(shí)間復(fù)雜度與空間復(fù)雜度都應(yīng)該比較低，也就是在c語(yǔ)言的程序中盡量少一些嵌入式循環(huán)，本次的索引算法中最多只有2個(gè)for循環(huán)嵌入使用，所以能夠保證時(shí)間復(fù)雜度較低，故查詢速度較快。本算法采用動(dòng)態(tài)數(shù)組，并且最后的索引使用也與k值的大小有關(guān)，k值越大，使用的內(nèi)存越少，所以能夠在一定程度上保證所用的內(nèi)存較小。針對(duì)問題三，此處給出建立索引的時(shí)間復(fù)雜度的分析。一個(gè)算法的空間復(fù)雜度只

9、考慮在運(yùn)行過程中為局部變量分配的存儲(chǔ)空間的大小，它包括為參數(shù)表中形參變量分配的存儲(chǔ)空間和為在函數(shù)體中定義的局部變量分配的存儲(chǔ)空間兩個(gè)部分。所以本算法的空間復(fù)雜度第一要考慮所使用的堆?？臻g的大小，它等于一次調(diào)用所分配的臨時(shí)存儲(chǔ)空間的大小乘以被調(diào)用的次數(shù)。本算法使用for循環(huán)函數(shù)，先計(jì)算出建立索引的算法中每一個(gè)語(yǔ)句的執(zhí)行次數(shù)，然后再相加，最后依據(jù)去低階項(xiàng)，去掉常數(shù)項(xiàng)，去掉高階項(xiàng)的常參的原則得到時(shí)間復(fù)雜度。然后根據(jù)算法臨時(shí)存儲(chǔ)的空間大小來計(jì)算空間復(fù)雜度，觀察有無(wú)臨時(shí)存儲(chǔ)大小以及臨時(shí)存儲(chǔ)大小與輸入變量的關(guān)系。針對(duì)問題四，與建立索引的復(fù)雜度的分析相似，先找到 for循環(huán)函數(shù)，計(jì)算出使用索引的算法中每一個(gè)

10、語(yǔ)句的執(zhí)行次數(shù)，然后再相加，最后依據(jù)去低階項(xiàng)，去掉常數(shù)項(xiàng)，去掉高階項(xiàng)的常參的原則得到時(shí)間復(fù)雜度。然后根據(jù)算法臨時(shí)存儲(chǔ)的空間大小來計(jì)算空間復(fù)雜度，觀察有無(wú)臨時(shí)存儲(chǔ)大小以及臨時(shí)存儲(chǔ)大小與輸入變量的關(guān)系。針對(duì)問題五，假設(shè)內(nèi)存限制為8g，分析所設(shè)計(jì)索引方法所能支持的最大k值和相應(yīng)數(shù)據(jù)查詢效率。內(nèi)存消耗最大的應(yīng)該是程序當(dāng)中的兩個(gè)超大數(shù)組，一個(gè)是存放文件的數(shù)組，另一個(gè)是存放索引表（hash表）的數(shù)組。通過對(duì)兩個(gè)數(shù)組所占內(nèi)存的分析，我們得出了算法所需要的內(nèi)存與k的函數(shù)關(guān)系。顯然，k若越大，則對(duì)應(yīng)的十進(jìn)制數(shù)也越大，在整型數(shù)據(jù)有最大值的條件下，我們通過理論分析給出了k最大能取14。為了擴(kuò)展k的值，我們考慮將長(zhǎng)度

11、超過14的分成多個(gè)部分，每個(gè)部分限制最長(zhǎng)為14。這樣，判斷兩個(gè)字符串相等的條件成為要每個(gè)部分都相等，對(duì)應(yīng)在算法里就是對(duì)每個(gè)部分的比較結(jié)果取交集。顯然，這種方法會(huì)導(dǎo)致使用更多的大數(shù)組，相應(yīng)會(huì)占更多的內(nèi)存。三、符號(hào)說明符號(hào)說明fp文件d:bdata放置dna序列文件的位置solexa_100_170_1.fa前半部分dna序列文件solexa_100_170_2.fa后半部分dna序列文件shuju.txt存放查找出的kmerndna總序列數(shù)m每個(gè)序列的長(zhǎng)度n宏定義1000000total 將4進(jìn)制轉(zhuǎn)換為10進(jìn)制后的存儲(chǔ)數(shù)組p讀取文件時(shí)的二維存儲(chǔ)數(shù)組k需查找的kmer序列的長(zhǎng)度內(nèi)存函數(shù)擴(kuò)展內(nèi)存函數(shù)

12、四、模型建立與求解4.1 問題一的模型建立與求解4.1.1文件處理對(duì)一百萬(wàn)條dna序列的讀取采用c語(yǔ)言的文件操作語(yǔ)句循環(huán)讀入。在每條序列循環(huán)讀取的時(shí)候，首先判斷是否為“eof”，若為空，則表示文件讀取完畢，結(jié)束，若不為“eof”，則繼續(xù)讀取dna序列，若讀取的是“atcg”中的任意一個(gè)，則以具體的堿基序列存入二維數(shù)組中，若不是“atcg”任意一個(gè)，繼續(xù)從文件中讀取下一個(gè)字符。算法流程圖如下：將該字符放進(jìn)一個(gè)二維字符數(shù)組結(jié)束是該字符是否為“atcg”中的一個(gè)否否從文件中讀取一個(gè)字符，判斷是否為eof開始打開文件是4.1.2 建立數(shù)據(jù)索引在建立索引時(shí)，首先從存放字符的數(shù)組中讀取一個(gè)字符，判斷是否為

13、該行（100個(gè)）的最后一個(gè)字母，如果是，則跳到下一行，如果不是，則分別令可令a=0,c=1,g=2,t=3，即轉(zhuǎn)化為四進(jìn)制數(shù)串。然后將長(zhǎng)度為k的字串按遞推的順序放進(jìn)索引數(shù)組。遞推的方式為將整個(gè)長(zhǎng)度為k的字符串舍去第一個(gè)字符，再加上下一個(gè)字符，成為一個(gè)新的長(zhǎng)度為k的字符串。這種遞推直到這一行的最后一個(gè)字符為止。如 s =“ctgtactgtat”。則字串為ctgta，tgtac，gtact，tactg，actgt，tgtat。對(duì)每一個(gè)字串，可以構(gòu)造這樣的哈希函數(shù)：其中sj指的是字串里從右往左數(shù)第j個(gè)字母對(duì)應(yīng)的數(shù)值（對(duì)應(yīng)規(guī)則為a=0,c=1,g=2,t=3）。每一個(gè)dna鏈只用計(jì)算一次哈希函數(shù)，后

14、面的k-mer對(duì)應(yīng)的hash值（十進(jìn)制數(shù)值）可以由以下的遞推公式給出：其中sj+k是下一個(gè)字符對(duì)應(yīng)的數(shù)值，sj是當(dāng)前第一個(gè)字符對(duì)應(yīng)的數(shù)值，這兩個(gè)字符的距離剛好是k。 如此一來，就已經(jīng)將每個(gè)k-mer與一個(gè)唯一的十進(jìn)制數(shù)對(duì)應(yīng)起來了。然后將這個(gè)十進(jìn)制數(shù)存放到一個(gè)二維數(shù)組total1000000100里面去。每個(gè)k-mer在數(shù)組中的位置由如下規(guī)則確定：行號(hào)是所在dna鏈的編號(hào)，列號(hào)是這個(gè)k-mer的第一個(gè)字符在該列dna鏈中的位置（從左往右數(shù)）。這樣就建立好了索引數(shù)組，程序中用total表示。算法流程圖如下：讀取下一個(gè)字符 結(jié)束是否所有字符轉(zhuǎn)換為了相應(yīng)的hash值是否是是否是第a個(gè)字符 （a=k）跳

15、到下一行從存放字符的數(shù)組中讀取一字符開始將這個(gè)四進(jìn)制數(shù)轉(zhuǎn)化為十進(jìn)制；存入索引數(shù)組是否是這一行最后一個(gè)字符是否a=0t=3g=2c=1用遞推公式計(jì)算下一個(gè)hash值否是否是第k個(gè)字符 是否4.1.3 使用數(shù)據(jù)索引計(jì)算輸入字符串的kmer值，與索引數(shù)組的total中的字符串依次比較，如果匹配，則說明找到了相應(yīng)字符串在100萬(wàn)個(gè)序列中的位置，然后將位置輸出到txt文件中，利用循環(huán)遍歷整個(gè)索引數(shù)組，找出所有匹配的字符串的位置。算法流程圖如下：否j=0，i+結(jié)束是i是否大于1000000j是否大于（100-）否是kmer值與索引數(shù)組totalij是否匹配開始i=0,j=0 計(jì)算輸入字符串的kmer值是否

16、輸出（i,j）到一個(gè)txt文件中j+是4.2問題三的模型建立與求解4.2.1建立索引所用的計(jì)算復(fù)雜度，和空間復(fù)雜度分析4.2.1.1建立索引時(shí)的時(shí)間復(fù)雜度分析：建立索引部分程序如下：該部分個(gè)語(yǔ)句需要運(yùn)行的次數(shù)如下圖所示：其中n1=1000000，m1=100，代表dna的1000000個(gè)序列，m代表每個(gè)序列的長(zhǎng)度，h為輸入的kmer值，該部分算法的時(shí)間復(fù)雜度等于各個(gè)語(yǔ)句的執(zhí)行次數(shù)相加，即o（n1+n1+n1*h+n1*h+n1*h+n1*h+n1*（m1-h）+n1*（m1-h）+n1+n1）=o（2n1*m1+2n1*h+2n1*h+2n1）o表示去低階項(xiàng)，去掉常數(shù)項(xiàng)，去掉高階項(xiàng)的常參得到時(shí)

17、間復(fù)雜度o,即o（n1*m1）.4.2.1.2建立索引時(shí)的空間復(fù)雜度分析：臨時(shí)存儲(chǔ)空間大小的部分為：空間復(fù)雜度是對(duì)一個(gè)算法在運(yùn)行過程中臨時(shí)占用存儲(chǔ)空間大小的量度，建立索引過程中數(shù)組hash是作為臨時(shí)占用儲(chǔ)存空間的數(shù)組，由于在for循環(huán)中有a=1）于是（k=1）（單位：mb）即k每增大1，所占用的內(nèi)存就會(huì)減少3.8147mb。n（k）max=476.8374mb4.3.2 k的最大取值分析限制的因素：限制k的主要原因是整型數(shù)據(jù)的最大值只有231-1；而在建立哈希表的過程中要用到指數(shù)函數(shù)pow（4，k-1） ；k的最大值可以由以下的表達(dá)式確定：即，可以求得。4.3.3對(duì) k15時(shí)索引情況的討論事實(shí)

18、上，對(duì)于pow函數(shù)（指數(shù)函數(shù)）來講，如果計(jì)算結(jié)果溢出，其輸出將自動(dòng)歸為0，在這里，最多能計(jì)算pow（4,15） 。因此，計(jì)算機(jī)在運(yùn)行時(shí)，當(dāng)建立k15的索引時(shí)（哈希表），計(jì)算機(jī)并不會(huì)報(bào)錯(cuò)，甚至可以實(shí)現(xiàn)查找,不過這種查找可能不會(huì)是精確的。這種不精確性體現(xiàn)在下面兩方面：1、由于pow函數(shù)溢出后輸出為0引起的。2、由于在建立哈希表時(shí)使用了遞推公式，可能會(huì)使數(shù)據(jù)溢出。4.3.4假設(shè)內(nèi)存限制為8g，分析所設(shè)計(jì)索引方法所能支持的最大k值和相應(yīng)數(shù)據(jù)查詢效率。4.3.4.1改進(jìn)后的最大k值分析我們可以擴(kuò)展k的取值。我們定義從一個(gè)字符串中從第i個(gè)位置到第j個(gè)位置的字串為“字串（i，j）”顯然，比較兩個(gè)字符串是否相

19、等，等價(jià)于他們每個(gè)字串都分別向相等，特別地，也等價(jià)于”字串（1，i）“和”字串（i+1，k）“這兩個(gè)字串分別相等。 于是，我們采用這樣的方法，將字符串分成”字串（1,14）“、”字串（15,28）“.”字串（85，98）“、”字串（99,100）“這樣的8個(gè)字串，即建立一共8個(gè)索引數(shù)組，然后根據(jù)不同的k值決定要使用幾個(gè)索引數(shù)組。對(duì)每個(gè)待查的k-mer也采用相同的辦法計(jì)算hash值，將其與索引數(shù)組進(jìn)行對(duì)比。然后對(duì)最后的結(jié)果取交集，可以正確找出匹配項(xiàng)的位置。顯然，對(duì)給定的k，一共需要（ ）個(gè)索引數(shù)組，因此，這樣看來，k越大，對(duì)內(nèi)存的需求也越大。最后，所需要的內(nèi)存與k的關(guān)系可以由以下的公式給出：將n

20、（k）代入上式可以得3.8147*（101-k）+95.3674對(duì)上式求導(dǎo)數(shù)，得:故這是一個(gè)單調(diào)遞增的函數(shù)，最大值在k=100處取得，于是 n*（k）max=n*（100）=793.4568mb4.3.4.2效率分析：為了擴(kuò)展k的最大值，我們無(wú)疑是犧牲了效率，因?yàn)樗饕龜?shù)組被建立了多次，但如果不是因?yàn)檎螖?shù)據(jù)溢出，這些操作是完全沒有必要的。如若原來時(shí)的效率為a，則當(dāng)時(shí)的效率為。五、模型一般性測(cè)試1. 查找k=5，k-mer=atcga的dna序列程序運(yùn)行圖在原始fa文件中查找對(duì)比來驗(yàn)證正確性，選取其中2個(gè)位置查找驗(yàn)證：在文件1（solexa_100_170_1）中，第7條dna序列中read_1

21、70_7_1 random_genome_10000000 1211051 100 gaggactcagcccctcctagcgataggattgtatgtagccagcttggaacagaatatcgaggctaacgcgatgtcccgcactgagcccatggcactgtat在文件2（solexa_100_170_2）中，第499987條dna序列中（所有序列中的第999987條）read_170_499987_2 random_genome_10000000 7315838 100 aatccgattccttatctacgatatgctatcgatacgacaaattgtcaatca

22、aagggcgtggacgacctagtgcgtgaagtgcccatgaacgaaggttact1. 查找k=10，k-mer=atcgtgcata的dna序列程序運(yùn)行圖：在文件1（solexa_100_170_1）中，第12727條dna序列中read_170_12727_1 random_genome_10000000 8805313 100 gaaagggcccatcatcgtgcatatccgcaaaagaacccaacatagtacaaaattcgatgaggatagtccaaatgttatgaaccatactagaaaaagttat在文件2（solexa_100_170_2）中，

23、第486940條dna序列中（所有序列中的第986840條）read_170_486940_2 random_genome_10000000 9507192 100 tgaacgcatctttggaggaaatatgctccttaagagtggccataccatgacgtcccagcgtggaatacgatcgtgcataatgatcgtcctgatgtagtga3.查找k=15，k-mer=的dna序列在文件1（solexa_100_170_1）中，第21823條dna序列中read_170_21823_1 random_genome_10000000 3021063 100 ctgtggg

24、tagacggatacaactacaagatcggaatgcccttcgtgaggcacaggctgtttcagtgtgcggtcaacgattgaggattgccggcttcggaat在文件2（solexa_100_170_2）中，第407297條dna序列中（所有序列中的第907297條）read_170_407297_2 random_genome_10000000 3021120 100 ctggcgtcgcatagttaggctccgcggcaatactccttgcatgagtgtatgaagtacctgtgggtagacggatacaactacaagatcggaatgcccttcg

25、六 模型的優(yōu)缺點(diǎn)及改進(jìn)方向6.1模型的優(yōu)缺點(diǎn)6.1.1模型的優(yōu)點(diǎn)：優(yōu)點(diǎn)：1.在程序算法中使用了小的數(shù)據(jù)類型，并且通過減少運(yùn)算的強(qiáng)度使得索引建立時(shí)間短，檢索速度快。2.算法使用索引所占內(nèi)存空間小，空間復(fù)雜度低，因?yàn)橹皇褂昧薿(1)。6.1.2缺點(diǎn)：1.該算法是基于4進(jìn)制與10進(jìn)制間的轉(zhuǎn)換，當(dāng)k14時(shí)，數(shù)據(jù)溢出，會(huì)覆蓋掉其他內(nèi)存地址的數(shù)據(jù)，導(dǎo)致檢索缺失。參考文獻(xiàn)1 郝甜甜 李強(qiáng)飛 李國(guó)治 陳禹翰 鄧衛(wèi)東 畜牧與飼料科學(xué) 2014年 第3期 云南農(nóng)業(yè)大學(xué)動(dòng)物科學(xué)技術(shù)學(xué)院 云南昆明 6502012 何建強(qiáng) 對(duì)karp-rabin串匹配隨機(jī)算法的改進(jìn) 廣西科學(xué)院學(xué)報(bào) 第18卷第4期 2002年11月附錄

26、實(shí)際使用的軟件名稱：visual c+ 6.0附錄dna序列的k-mer index數(shù)據(jù)查找程序#include #include #include #include #define n 1000000double totaln100;char pn100;int main()printf(準(zhǔn)備運(yùn)行程序.n);int o,n; int i,j; char l100;o=1000100;n=105;file *fp1;file *fp2;if(fp1=fopen(d:bdatasolexa_100_170_1.fa , r)=null)/文件路徑一定要正確，請(qǐng)將數(shù)據(jù)fa文件放置在d盤bdata文

27、件夾中。printf(打開文件失敗，請(qǐng)確認(rèn)文件所在路徑正確);getchar();return(1);if(fp2=fopen(d:bdatasolexa_100_170_2.fa , r)=null) /文件路徑一定要正確，請(qǐng)將數(shù)據(jù)fa文件放置在d盤bdata文件夾中。printf(打開文件失敗，請(qǐng)確認(rèn)文件所在路徑正確);getchar();return(1);j=0;i=0;int k;char ch1,ch2;printf(- 讀取文件中. -n);while( ch1!=eof )ch1=fgetc(fp1);switch(ch1)caset:k=0;break;casea:k=0;b

28、reak;caseg:k=0;break;casec:k=0;break;default:k=1;if(k=1)continue;if(i=100)i=0;j+;if(j=500000)break;pji=ch1 ;i+;fclose(fp1);while( ch2!=eof )ch2=fgetc(fp2);switch(ch2)caset:k=0;break;casea:k=0;break;caseg:k=0;break;casec:k=0;break;default:k=1;if(k=1)continue;if(i=100)i=0;j+;pji=ch2 ;i+;if(j=1000000)b

29、reak;fclose(fp2);int h,m;printf(請(qǐng)輸入您想要查找的k值:);scanf(%d,&h);printf(- 正在建立索引. -n);clock_t begin, end;double cost;begin = clock();int f=1000100,g=105; int b,q; int a,r=0;double hash100,s100;for(q=0;q=1000000;q+)totalq0=0;for(b=0;bh;b+)switch(pqb)casea:hashb=0;break;caseg:hashb=2;break;caset:hashb=3;break;casec:hashb=1;break; sb=hashb;totalq0=totalq0+hashb*pow(4,(h-b-1);for(a=0;a100-h;a+,b+)switch(pqb)casea:hasha=0;break;caseg:hasha=2;break;caset:hasha=3;break;casec:hasha=1;break;sa+h=hasha;totalqa+1=hasha+4*(totalqa-sa*pow(4,(h-1);end = cloc