寫在前面
前述，咱們已經(jīng)完成了：
- 重測(cè)序讀段回帖到基因組上，得到了 SAM 文件
- 對(duì) SAM 文件進(jìn)行排序，得到 Position-Sorted 的 BAM 文件
于是接下來，就是 SNP Calling （也就是變異檢測(cè)）之前的步驟....即 “Mark Duplicates”。
功能界面
原本寫成兩個(gè)插件，但想想我還是把這兩塊合并，主要是他們之間的代碼基本是一樣的。有區(qū)別的，那么只是調(diào)用的命令有不同。
看一看 Mark Duplicates 的界面
于是，問題解決。下面來個(gè)示例
輸出文件也簡(jiǎn)單
寫在最后
Emmm，時(shí)間差不多了。剩下的功能就不多了：
- Call SNP (不含 Indels)
- Call SNP (含 Indels)
至于，為什么此處兩塊要分開？因?yàn)?Indels 的檢測(cè)相對(duì)麻煩。一般可能會(huì)使用 GATK 等流程。但現(xiàn)在這個(gè)流程，說實(shí)話，用起來很不容易。參考已公布的評(píng)測(cè)結(jié)果，會(huì)選用一個(gè)「性價(jià)比」最好的策略....
說的太多，整了三個(gè)插件了。也還沒看看到底自己比對(duì)的是啥。用 IGV 看看。
Emmm... 看起來 1M paried-reads 還是太少了。不過就比對(duì)來說，在我的 PC 上，就需要 30min了。
按照 黃瓜 基因組為 230Mb 來計(jì)算，測(cè)序方式是 PE150，于是可以得到，覆蓋 30X 至少需要 23M Reads，配對(duì)起來，那么就大概是 12M PE，于是需要 6 個(gè)小時(shí)才能搞定一個(gè)樣品的比對(duì)。但事實(shí)上，很可能會(huì)測(cè)到比如 50X 甚至 100X。所以在比對(duì)的時(shí)間上，或許會(huì)變成 12h，甚至 24 小時(shí)一個(gè)BSA測(cè)序樣品。Anyway，我計(jì)算的只是 4.0GHz 6 個(gè)線程的情況下。假設(shè)你有個(gè)不錯(cuò)的PC，可以有 18 個(gè)線程 比如，那么你又可能把時(shí)間縮短為 8 個(gè)小時(shí).... 這個(gè)咋說呢，「過夜培養(yǎng)」就好了。
