任意分子标记的组装

  缘起：课题组需要组装一批物种的matK、rbcL与ITS分子标记，用于后续的生态学研究。

  对于matK、rbcL与ITS这类常用分子标记，根据课题组的经验，都是使用GetOrganelle组装的。前两者是成功组装并注释叶绿体后获取的。ITS则是GetOrganelle附带可组装的对象。于是我以为可以很轻松的做完。但结果是无论如何调整GetOrganelle软件的参数，这包括尝试不同的组装轮数（-R的设置），使用同属物种的参考基因组作为种子与基因标签（seed与label），仍然止步于59.38%的叶绿体组装率与37.5%的ITS组装率。值得一提的是，除了一次性成功的组装，其后的努力对组装质量的提高没有任何效果（通过Bandage查看），这很不尽如人意并使人疑惑。
  对此，我首先考虑了测序问题，我们获得的单样品数据量都是10G。这批样品的采集笔者全程参与，因此可以排除原始样品DNA过度降解的可能（至少不会大比重降解）。至于建库测序，是基于DNBSEQ的BGI方法，我看测序结果报告Q20达99以上，Q30达97以上，故可排除测序误差（这套流程是全自动机械化过程）。其次，考虑到样品涉及大量被子植物的不同科属（均为木本），那么是否可能有些物种的基因组太大导致测序深度不够，使得叶绿体与ITS的组装困难呢？虽然这是有可能的，但我检查了那些失败的组装，物种的基因组并不大（小于1G，几百M，至少在我检查的几个如此）。而且这很可能不是组装失败的主要原因，因为包括一个我非常熟悉且知其拥有组装高成功率的物种也失败了（Quercus jenseniana）。

测序公司的DNA跑胶图

  于是我只能认为是DNA提取存在重大嫌疑。我看了公司的DNA提取报告，最直截了当的琼脂糖凝胶电泳结果大抵还行，但是颜色很寡淡。使用Nanodrop与Qubit对DNA浓度与含量进行检测，结果显示：前者所有物种的C (ng/μl)在15.0–192.9之间，260/280在1.67–2.00之间，260/230在0.75–1.74之间；后者C (ng/μl)在2.16–46.6之间，总量(μg)在0.13–2.80之间。如果单纯看这份报告。首先，毫不意外的是Nanodrop与Qubit检测结果差很多，通常认为前者所得的定性结果误差很大，后者的精度更高。那么有19个物种（59.38%）的总量＜1 μg，这非常低。对于DNA总量，在华大基因官网（及其他一些公司的官网）中的送样建议是≥1 μg（对应浓度≥16 ng/μL）。再看Nanodrop检测所能说明的问题，260/280多数在纯DNA的范畴内（1.8–2.0），仅1个不在；但260/230全部低于标准（1.8–2.2），表明DNA提取过程中有很多杂质污染，包括蛋白质、胍盐、苯酚、Trizol、EDTA、乙醇等。综合以上分析，DNA提取的质量应该较差，不仅浓度低还有一定程度的污染。在所有未能组装成功叶绿体的样品中，DNA总量不足1 μg的有7个，占比53.85%，且占所有DNA总量不足1 μg样品的36.84%。以上表明，DNA总量低不一定导致组装失败，并且组装失败还有其它原因。尽管如此，较低的DNA提取总量与浓度，会对测序数据总体质量产生重要影响，如果是做对数据质量要求很高的分析时更为如此。最近，我了解到不同的DNA提取策略、文库制备与样品的存储都会对测序结果产生影响，并且DNA提取对测序的影响最为深远。那么，除却提取浓度与污染，是否是DNA试剂盒影响了我们的结果呢？目前不得而知。
  为了获取到matK、rbcL、ITS，我只好另辟蹊径。我采用今年刚发布的新软件Geneminer实现。使用Geneminer，给以参考序列，可以组装任意分子标记。

# 单参考。-rtfa 指定参考序列
geneminer.py -1 skimming_data1.fq.gz -2 skimming_data2.fq.gz -rtfa shallow_ref/ITS.fasta -o ITS_out
# 多参考。
geneminer.py -1 skimming_data1.fq.gz -2 skimming_data2.fq.gz -rtfa shallow_ref/ -o out1 # -rtfa指定到装有多个参考序列的目录。
# 其他常用参数：
-t：线程
-b：基于参考序列组装完成后向序列侧翼延伸多少碱基，默认75 bp。例如-b 200则双端各延伸200 bp。

Geneminer组装结果

注意：
  （1）Geneminer的算法组装特别依赖参考序列（参考基因组），如果亲缘关系较远，则不能获得结果，以至于我在使用同属的物种时也曾发生1次组装失败的案例，但换用更近缘的物种时就能成功。
  （2）如果-rtfa指定文件夹，那么Geneminer会对所有参考进行评估，并找出适用的参考得出结果，因此耗费时间较长。如果有两个参考都能适配样品序列做组装，将输出两条结果。
  （3）Geneminer不依赖其他软件，使用其自开发的算法，其前身是Easy353。就开发者的测试来看，在计算速度、硬盘占用较aTRAM、HybPiper优秀。从对被子植物353个通用低拷贝核基因的组装表现来看，Geneminer的精度、准确度、覆盖度最高。
  由于我没试过HybPiper，无法评价。但在使用Geneminer做组装之前，我试图使用GetOrganelle做组装。

get_organelle_from_reads.py -s seed/Rhua_matK.fasta --genes seed/Rhua_matK.fasta \
-1 resequence1.fq.gz -2 resequence2.fq.gz -o target_dir -R 20 -t 8 -k 21,45,65,85,105 \
-F anonym  
# -s：参考序列，--genes：与-s指定完全相同的文件（因为软件设计限制，必须与-s协同使用，当-F anonym时）。

  GetOrganelle可以得到组装结果，但序列非常长，远大于参考的长度，显然目标序列蕴含在内了，这给后续在NCBI里上传序列带来麻烦。类似地，GetOrganelle组装的ITS也很长，一般的被子植物包含ITS1/2序列的全长只有600–700 bp。我试验了将这么长的序列拿去与长度适当的ITS Mafft，似乎存在对齐偏差。
  由此可见，GetOrganelle组装叶绿体是不二之选，但对分子标记还得靠更专业的软件。

---

参考资料：

[1] 超详细Nanodrop结果判读！（下）——A260/A280与A260/A230比值偏高怎么办？
[2] 干货分享|一文全面解读迷之又迷的260/280、260/230
[3] 核酸定量哪家强？Nanodrop vs. Qubit
[4] 如何使用酶标仪对DNA和RNA 进行快速准确的检测
[5] Nature：不同提取方法对实验样本DNA的影响
[6] Costea P I, Zeller G, Sunagawa S, et al. Towards standards for human fecal sample processing in metagenomic studies[J]. Nature biotechnology, 2017, 35(11): 1069-1076.
[7] Xie P, Guo Y, Teng Y, et al. GeneMiner: A tool for extracting phylogenetic markers from next‐generation sequencing data[J]. Molecular Ecology Resources, 2024, 24(3): e13924.
[8] Zhang Z, Xie P, Guo Y, et al. Easy353: A tool to get Angiosperms353 genes for phylogenomic research[J]. Molecular Biology and Evolution, 2022, 39(12): msac261.
[9] Genemimer-tutorial online
[10] GetOrganelle-How to assemble a target organelle genome using my own reference?