写在前面
前述,咱们已经完成了:
- 重测序读段回帖到基因组上,得到了 SAM 文件
- 对 SAM 文件进行排序,得到 Position-Sorted 的 BAM 文件
于是接下来,就是 SNP Calling (也就是变异检测)之前的步骤....即 “Mark Duplicates”。
功能界面
原本写成两个插件,但想想我还是把这两块合并,主要是他们之间的代码基本是一样的。有区别的,那么只是调用的命令有不同。
看一看 Mark Duplicates 的界面
于是,问题解决。下面来个示例
输出文件也简单
写在最后
Emmm,时间差不多了。剩下的功能就不多了:
- Call SNP (不含 Indels)
- Call SNP (含 Indels)
至于,为什么此处两块要分开?因为 Indels 的检测相对麻烦。一般可能会使用 GATK 等流程。但现在这个流程,说实话,用起来很不容易。参考已公布的评测结果,会选用一个「性价比」最好的策略....
说的太多,整了三个插件了。也还没看看到底自己比对的是啥。用 IGV 看看。
Emmm... 看起来 1M paried-reads 还是太少了。不过就比对来说,在我的 PC 上,就需要 30min了。
按照 黄瓜 基因组为 230Mb 来计算,测序方式是 PE150,于是可以得到,覆盖 30X 至少需要 23M Reads,配对起来,那么就大概是 12M PE,于是需要 6 个小时才能搞定一个样品的比对。但事实上,很可能会测到比如 50X 甚至 100X。所以在比对的时间上,或许会变成 12h,甚至 24 小时一个BSA测序样品。Anyway,我计算的只是 4.0GHz 6 个线程的情况下。假设你有个不错的PC,可以有 18 个线程 比如,那么你又可能把时间缩短为 8 个小时.... 这个咋说呢,「过夜培养」就好了。