关于 16s 序列的注释算法
最近在研究 16s 序列注释的问题，看了几篇文章，稍微整理了下。有些观点我也还需进一步验证，不一定正确，欢迎大家留言讨论。
16S rRNA 扩增子测序已被广泛应用于微生物组研究，其中一个至关重要的步骤就是对相应序列（OTU 或者 ASV）进行分类学注释。虽然已经有大量的注释算法被开发出来，但是我们最普遍用到的仍是朴素贝叶斯分类器（NBC）。追溯其历史，RDP（Ribosomal Database Project）首先使用了NBC 算法[1]来对序列进行分类注释，证明了16S rRNA 序列可以进行属水平分类。很多研究也证明 NBC 在标记基因序列分类上的效果确实十分稳健，比如 QIIME2 内置的 q2-feature-classifier[2] 分类器。
q2-feature-classifier
“Optimizing Taxonomic Classification of Marker-Gene Amplicon Sequences with QIIME 2’s Q2-Feature-Classifier Plugin.” Microbiome 6, no. 1 (May 17, 2018): 90. https://doi.org/10.1186/s40168-018-0470-z.
QIIME2 提供了 3 种不同的方法来进行物种注释，包括两种 alignment-based 的方法（VSEARCH，BLAST+）和一种基于 naive-bayes 机器学习的方法。在模拟数据集中，这些方法的 F-measure 随着分类水平的降低而降低，但在属水平的 F-measure 均大于 0.8（minimum: UCLUST F = 0.81, maximum: naive Bayes bespoke F = 1.00）。朴素贝叶斯的 F-measure、召回率、分类单元检出率、分类单元准确率都高于其他方法，但 Bray-Curtis 距离更小。
在交叉验证中，这些方法的准确性指标相差不大。但对于新种的推断，UCLUST、BLAST+ 和 VSEARCH 的算法显著优于其他算法（P < 0.05）。
q2-feature-classifier 的建模过程进行了两次超参数选择：第一次网格搜索主要用于优化特征提取参数，第二次则是更关注优化 k-mer 长度和置信度参数。
第一次网格搜索的过程，包括 sklearn.feature_extraction.text.HashingVectorizer， sklearn.feature_extraction.text.TfidfTransformer， sklearn.naive_bayes.MultinomialNB。参数见下表。这一步的结论是 TfidfTransformer 步骤没有显着改善分类性能，即 n_features 应设置为 8192，特征向量应使用 L2 归一化进行归一化，并且朴素贝叶斯分类器的 alpha 参数应设置为 0.001。
其实 QIIME2 中的 q2-feature-classifier 插件可以很方便地将 scikit-learn 中的机器学习算法迁移到标记基因序列的分类问题。该插件的默认算法，即是从参考序列中提取 k-mer 计数并训练 scikit-learn 多项式朴素贝叶斯分类器。展开来说，就是用 sklearn.feature_extraction.text.HashingVectorizer 进行特征提取和用 sklearn.naive_bayes.MultinomialNB 进行分类。不同于 RDP 分类器默认的 8-mers，scikit-learn 中的哈希特征提取器则可以使用更长的 k-mers。此外，和大多数 scikit-learn 分类器一样，我们可以在训练多项朴素贝叶斯分类器时设置权重。
所以，这篇文章也提出了物种权重信息可能提高分类准确性的假设，并在一定程度上进行了验证。因为默认的假设是每个分类群出现在每个样本中的可能性相同，这往往和实际情况有所不同。关于这部分结论可参考下文。
“Species Abundance Information Improves Sequence Taxonomy Classification Accuracy.” Nature Communications 10, no. 1 (October 11, 2019): 4643. https://doi.org/10.1038/s41467-019-12669-6.
其他算法
从上文我们可以看到，NBC 的分类性能已经十分优异，同时对算力的要求也不高。那么这就带来了新的思考：
1.更先进的监督学习算法能否超越 NBC 的准确性？2.我们是否已经达到了 16s rRNA 分类性能的上限？
最近的一篇文章就用了随机森林和神经网络算法来进行分类预测，同时也和 NBC 进行了比较。
“Beating Naive Bayes at Taxonomic Classification of 16S rRNA Gene Sequences.” Frontiers in Microbiology 12 (2021): 1664. https://doi.org/10.3389/fmicb.2021.644487.
先说这篇文章的结论，RF 和 CNN 性能接近，但准确性并没有超过加入了权重信息的 NBC。关于 16s rRNA 分类性能的上限，对于150 nt 的 reads 来说，使用加入权重信息的 NBC 几乎已经达到性能极限。对于较长的扩增子（V4 全长） ，性能提升也十分有限。所以，与 RF 和 CNN 分类器相比，对算力要求较低的 NBC 显然对于研究者来说更为友好。
随机森林分类器
和之前的建模流程相似，作者同样使用了两次超参数选择。先在一个相对较小的数据集上进行网格搜索，以选择对模型性能相对重要的超参数。然后在一个更大的数据集上进行第二次网格搜索。初步试验结果表明，max_depth 和 max_features 对分类准确度有显著影响。
这里的分类准确度基于 F-measure 进行计算，用于种水平的分类。使用 Matthews 相关系数 (MCC) 校正不平衡数据集中的 F-measure 的偏差。结果表明，无论使用什么参数，NBC 均优于 RF（Wilcoxon, p < 0.05）。
这幅图也能说明最大树深度 (max_depth) 对分类准确度的影响最大，增加最大深度会导致 F-measure 增加。但要是 max_depth 超过 64 ，模型的性能提升不大。同时，减少纳入的特征数量 (max_features) 也会导致分类准确度略有下降，在较低的置信水平下，这种性能下降较不明显。增加估计量（n_estimator）对分类准确度影响较小。
卷积神经网络
除了使用 RF 建模外，作者还测试了 CNN 的效果。在将 DNA 序列输入网络之前，调用了 word2vec 采用 CBOW（Continuous Bag-Of-Words Model，连续词袋模型）将遗传信息转换为一系列 300 个元素的向量。这里除了能够将 k-mers 转换为数值，而且还携带有关给定序列中任何两个 k-mers 之间相关性/相似性的附加信息。
在大多数测试中，作者使用了一个较为简单的神经网络，包含一个（一维）卷积层，然后是一个全局最大池化层和一个分类层（下图 A，架构 I）。用到的超参数选择范围 （见下表）。
超参数选择的结果表明，增加过滤器的数量或 kernel size 会提高准确度，平均 F-measure 在 0.710 ∓ 0.009 和 0.801 ∓ 0.006 之间。基于这个结果，作者定义了一组参数来测试置信度设置对 CNN 分类性能的影响。结果表明，将置信参数降低至 0.5 可提高性能（平均 F-measure 为 0.849 ∓ 0.005，filters = 512/kernel_size = 7），但要是将置信度进一步降低至 0.25 时，性能并不会进一步提高。进一步的测试也表明，过滤器数量或 kernel size 加倍对分类器精度也几乎没有影响。除此之外，作者还测试了不同网络架构（上图 B 中的架构 II ）和特征提取方法（独热编码）对 CNN 性能的影响。
总而言之，在这篇文章中，CNN 平均 F-measure 的最大值约为 0.85 左右，始终显著低于 NBC（p < 0.05）。
完美分类器
RF 和 CNN 分类器的平庸表现，让人不禁假设 NBC 可能已经接近这个问题的分类准确度上限，因此单独的替代算法不可能超过这个性能。
为了检验这个假设，作者构建了一个完美分类器（Perfect Classifier, PC）来衡量给定分类任务的分类准确率上限。这个分类器只有在两个或多个物种共享完全相同的序列时才会分类失败。当它们共享相同的序列时，算法将随机选择一个匹配的分类结果作为该序列的标签。所以该分类器的性能代表了分类准确度上限。
结果表明，对于 V4 150 nt 的完美分类器，加入权重信息的完美分类器的平均 F-measure 在 0.887 和 0.903 之间，而最佳 NBC 的平均 F-measure 为 0.865 ∓ 0.002 （Wilcoxon, p < 0.05）。对于全长 V4，结果也非常类似，完美分类器的平均 F-measure 得分介于 0.934 和 0.943 之间，而 NBC 为 0.866 ∓ 0.002（Wilcoxon, p < 0.05）。但这里需要注意，NBC 执行的是样本外交叉验证，而完美分类器是训练样本内的交叉验证，所以 F-measure 肯定偏大。因此，这两者的性能其实差距很小，说明 NBC 已接近分类性能极限。
总结
综上所述，16s rRNA 序列的保守型一定程度上限制了分类器的分类性能。目前这些成熟的分类算法已经十分接近分类准确度上限。当然，很多其他变量也会对分类性能有一定的影响，比如 reads 长度、引物选择、标记基因靶标、序列熵等等。参考数据库中的错误注释也导致了这个问题更加复杂。所以分类算法的进一步改进，不太可能仅仅来自新算法的应用，而更可能是某种组合：
1.挖掘 DNA 序列或其他潜在信息中的高位特征；2.利用先前研究中的生态权重信息来提高分类准确性；3.参考序列和分类数据库的进一步优化；4.更长的读长和/或其他标记基因靶标；5.改进的点不仅局限于准确性，还可以包括运行时间、内存使用量等等。
引用链接
[1] NBC 算法: https://doi.org/10.1128/AEM.00062-07
[2] q2-feature-classifier: https://doi.org/10.1186/s40168-018-0470-z