題:
計算至少與讀取序列的x%匹配的重複kmers序列
hilta007
2018-03-14 16:45:08 UTC
view on stackexchange narkive permalink

在處理fastQ文件時,我希望得到給定長度的所有可能kmers的重複序列,這些重複序列至少佔整個數據集讀取長度的90%。

例如:

的長度為6,帶有kmer“ ATTGGG”和一個數據集,每個數據集包含300個鹼基的讀數。

我想獲取以下kmers:

  ATTGGGATTGGGATTGGGATTGGGATTGGGATTGGG ......(讀取長度的90%,無間隙)ATTGGGATTGGGATTGGGATTGGGATTGGGATTGGGA .....(91%)ATTGGGATTGGGATTGGGATTGGGATTGGGATTGGGAT ....(92%)  

這對於長度為6的kmers的所有可能組合。

我進行了一些研究並測試了許多軟件,例如(KMC2,Jellyfish,Khmer,Dsk, Scturtle),但沒有一個符合我的意願,因為我想獲得重複kmer的出現,並且在一定百分比的讀取覆蓋率之間沒有任何可能的差距。

我已經完成了天真的算法,但t優化顯然不是最好的。有工具可以達到目的嗎?

一 回答:
user172818
2018-03-14 20:03:23 UTC
view on stackexchange narkive permalink

以下javascript可以滿足您的需求。您需要運行node.js。將代碼轉換為Python應該很容易。我沒有仔細測試。謹慎使用。

編輯(對新註釋的響應):程序已更改為計算長度之和。注意,僅計數長度為k * 2或更長的串聯重複序列。例如,按順序 ATTGGGATTGGGATTcGGGATTGGG ,腳本返回15,因為第二個時間不夠長。小於k * 2的長度計數會更慢且更複雜,但通常不是正確的選擇。例如,在 ATTGGGATTGGGcccccccccgaaatcgatagcatcgaGGGATTGcgatc 中,我們不會將第二個 GGGATTG 視為串聯重複。

在性能上,對於長度為$ l的輸入字符串$,此腳本的時間複雜度是$ O(l)$和一個小常數。

如果使用 TRF,則應在序列上運行它,然後測試是否有任何重複的串聯重複序列找到匹配的 GGGAAT

–步驟與我的腳本類似,不同之處在於TRF還會發現不正確的重複。

  //給定重複單元的長度_k_,找到所有2 * k或更長時間的串聯重複函數trf_k(k,str){var streak = 0,a = []; for(var i = k; i < = str.length; ++ i){if(i(i < str.length && str [i] == str [i-k])){++ streak; } else {if(streak > = k)a.push([i-streak-k,i]);條紋= 0;返回一個串聯重複長度之和,其中_kmer_為重複單位函數trf_kmer(kmer,str){var a = trf_k(kmer.length,str),sum = 0;對於(var i = 0; i < a.length; ++ i)如果(str.substring(a [i] [0],a [i] [1])。indexOf(kmer)> = 0)sum + = a [i] [1]-a [i] [0];返回總和;}
每個級別='persistentGggattagattagassagattagattagattagattagattagsuks'; console.log(side_kamer('gagutt',level));  
對不起,我的意思是我想要每個可能重複序列的完美匹配,就您的算法而言,當存在缺口時,仍然算作平方公里。確實,任務是高度定制的。我想我將開始在程序上使用並行性。我會看一下TRF,謝謝。
@hilta007感謝您的解釋(順便說一句,我要刪除我的舊評論)。如果要查找總長度,只需將“ max”更改為“ sum”,但要注意(有關詳細信息,請參見編輯)。我的實施應該相當有效。我不會太擔心性能。一個好的單線程腳本有時可能比一個壞的多線程C / C ++程序要快。
感謝您所做的修改,現在等待TRF的結果和我自己的算法來比較質量/速度。我會牢記您的算法。


該問答將自動從英語翻譯而來。原始內容可在stackexchange上找到,我們感謝它分發的cc by-sa 3.0許可。
Loading...