QTLseq数据分析，鼠标点点，复现 Nature Plant 番茄论文结果

写在前面

快一年了，写了重测序三兄弟插件，几乎所有用户都可以用 TBtools 轻松完成 BSAseq 数据分析。当然，我自己基本也没怎么折腾过，毕竟当时只是为了上个课。前几天将家人送回老家后，开始进入浑浑噩噩的状态。大体是每天也不知道干啥或者不干啥，尽管确实也在干活。索性，就还是简单整理一下一些插件功能，也再测试测试看看。于是决定，干脆先重现一下一些论文分析结果。比如这篇 Nature Plant 论文，定位了番茄分支数目的决定位点论文（分支越多，果实越多嘛，很好理解）。

于是下载了数据

SRR8307487  suppressed_ASRR8307489  branched_A

于是我下载了基因组和测序数据

看起来不错，数据挺大的。可以开干，整体步骤如下：

1. 读段回帖

2. 比对结果排序

3. 标记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_ASRR8307488  suppressed_B# 下面两组SRR8307489  branched_ASRR8307490  branched_B