写在前面
快一年了，写了重测序三兄弟插件，几乎所有用户都可以用 TBtools 轻松完成 BSAseq 数据分析。当然，我自己基本也没怎么折腾过，毕竟当时只是为了上个课。前几天将家人送回老家后，开始进入浑浑噩噩的状态。大体是每天也不知道干啥或者不干啥，尽管确实也在干活。索性，就还是简单整理一下一些插件功能，也再测试测试看看。于是决定，干脆先重现一下一些论文分析结果。比如这篇 Nature Plant 论文，定位了番茄分支数目的决定位点论文（分支越多，果实越多嘛，很好理解）。
于是下载了数据
SRR8307487 suppressed_A
SRR8307489 branched_A
于是我下载了基因组和测序数据
看起来不错，数据挺大的。可以开干，整体步骤如下：
- 读段回帖
- 比对结果排序
- 标记PCR重复
1.读段回帖
首先，需要测序得到的reads比对到基因组上，使用 TBtools 的 「BWA-mem2 插件」即可。具体界面如下：
过了大概 24 小时，第一个样品比对结束了，并产生了一个 10Gb 的 BAM 文件。这里似乎速度比较慢。有点怀疑数据是否是有接头。回头我可以跑一个 FastQC 看看。当然，逻辑上没啥问题，毕竟一共是 32Gb+ 的reads。
整体时间差不多。也要考量本地电脑 6 个线程确实慢。手上我的CPU也是 2016年的。
2. 比对结果排序
这一步很快，大体20分钟不到。
3. 标记PCR重复
重测序数据分析，最大的问题就是PCR重复会影响SNP的检测，所以需要标记出来
这一步也很快
4. 检测基因组变异 Call SNPs
开始检测变异
点击开始，然后报错
调整后，重新 Start，于是开始跑
我们用的是 bcftools，速度还是很快的
5. 过滤低质量 SNP
几乎所有变异检测的软件或流程出来的结果都是需要过滤的。因为每个人认为可靠的 SNP 的标准是不同的。比如有些人觉得测序深度要10X，有些人觉得5X就足够了。在TBtools的这个Pipeline中，我们直接用默认的。逻辑上对BSAseq影响不大。
这一步逻辑上非常快。事实上，整体流程瓶颈还是在比对。等待几分钟即可。
6. 进行QTL位点检测
使用相对简单，不过还是有两三个参数要手动给一下，具体在设置输入的 VCF 文件时，会有预览窗口，通过预览文本窗口的信息，就可以直接复制黏贴设置，整体如下
不会很久，不过一般也要几分钟
复现结果如下
与原文结果比较
注意：图片美化上自然还有很多办法....；所以结果挺好
具体输出的表格信息也可以看一下定位的区间
尴尬了，区间没对应上。当然，第一反应是他没有问题，同时咱们的流程也没问题。
随后查看了下 2019年05月，过去三年多了，基因组版本更新了，现在我们用的是新版本。看看原来的注释在我们用的版本的注释几何？
结果没问题，偏离一点点。当然了。精准的结果我们完全可以通过指定某个染色体，重新跑一次，结合 deltaSNP 和 G’value 来看。此外，我们其实只用了论文中的一部分数据，每个混池都只选了10+个材料，完全可以再下载两组数据，合并好然后跑一下。另外还注意到，原文作者对群体做了分层，
OK，复现结果。我相信，只要合并数据，完全就可以解决这个问题，感兴趣的朋友可以试试
SRR8307487 suppressed_A
SRR8307488 suppressed_B
# 下面两组
SRR8307489 branched_A
SRR8307490 branched_B
如果你复现了，有结果了，欢迎投稿，有稿酬
写在最后
至于为什么要来复现这个图稿，自然我仍然希望 TBtools的「插件模式」可以给更多朋友打开数据分析大门，自己动手，丰衣足食。逻辑上现在只要你有群体，其实可能2000块钱不到，就可以搞定两个F2或者F1群体的混池测序，为何不试试呢？也没多少钱。
最后，大伙用系列插件时，请引用Wrapper对应的原软件，上述包括但不仅限于：
bwa-mem2, samtools, bcftools, qtlseqr
至于 TBtools 最近我们观测的最高引用频率是 每天在Web of Science 被引 10 次，所以是否引用，大伙可以看心情。
网友评论