组会论文 | Microbiome 利用二代测序数据进行高分辨率菌株水平微生物组成分析

题目：High-resolution strain-level microbiome composition analysis from short reads
DOI：https://doi.org/10.1186/s40168-023-01615-w
期刊：Microbiome

本文开发的软件为：StrainScan

1. 尽管已有大量的基于二代测序的微生物组成分析工具，但是它们并没有优化以应对菌株水平分析中的挑战：高度相似的菌株基因组和一个样本中存在一个物种下的多个菌株。
2. 本文作者提出了一个菌株水平分析工具StrainScan，通过应用基于树的kmers索引结构在菌株水平识别准确率和计算复杂性取得平衡。
3. 在大量的模拟和真实测序数据进行测试，StrainScan在细菌菌株识别方面的准确率和分辨率均优于一些流行的菌株分析工具。
   
   

现有的工具对微生物物种水平的识别已经达到非常高的准确率，但想要进一步进行菌株水平的识别仍然存在比较大的一些挑战。

1. 样本中存在高度相似的菌株
   尽管有菌株分析工具，但它们可能需要来自同一群体的多个样本，只输出优势菌株，或者对菌株之间的相似性构成限制。
2. 菌株水平识别的分辨率
   参考数据库越大，分辨率越高；尽管有些菌株有很高的相似性，但没有已知的相似性临界值，低于该临界值可以忽略遗传差异。即使菌株含有高ANI（>99.9%），也可以被不同的噬菌体感染，表现出不同的防御或吸附机制。
3. 低丰度菌株的识别
   许多菌株水平分析工具需要10×覆盖度才能返回准确的识别。
4. 菌株识别时间
   当参考数据库很大时候，基于比对的方法需要的计算成本很高。

Microbiome | 香港城市大学孙燕妮组‍发表高分辨率细菌菌株识别工具

廖和睿，公众号：宏基因组Microbiome | 香港城市大学孙燕妮组发表高分辨率细菌菌株识别工具

本文作者已经在简洁凝练地分享过文章的结果部分，下面我在此分享下该工具的方法。

图1. StrainScan的整体工作流程

StrainScan接受用户提供的二代测序数据 (reads) 和用户感兴趣的某种细菌的菌株基因组，然后输出该样本中包含的一株或者多株细菌及其相对丰度。因为样本里的菌株往往和参考菌株不同，所以StrainScan会输出和真实菌株最相似的菌株。

算法简介：StrainScan利用了一种全新的数据结构 (层次聚类树) 以及更多的基因组特征来实现高效，准确且高分辨率的菌株识别。如图1所示，对于输入的细菌菌株基因组，StrainScan首先会根据他们的相似度进行层次聚类 (图1a)，然后根据聚类结果，基于基因组k-mer构建一棵层次聚类树 (图1b)。该树的每一个叶子节点即代表一个菌株簇。而对于每一个簇，当簇内菌株数目多于一株时，StrainScan会进一步提取该簇内部所有菌株的基因组特征用于簇内菌株识别 (图1c)。值得一提的是，这里提取的特征不仅包含普遍使用的strain-specific k-mer，也包含了部分菌株共享的group-specific k-mer以及所有菌株拥有的joint k-mer。这些特征的共同使用则可以提升菌株识别的分辨率。接着，将这些提取的特征保存在一个矩阵中供后续识别使用 (图1d)。于是，当一个样本输入时，StrainScan会首先利用层次聚类树快速锁定样本中可能存在的簇 (图1e)。然后，StrainScan会在识别到的簇中利用先前建立的特征矩阵进一步识别可能存在的菌株。最终，StrainScan会利用弹性网回归模型来预测菌株的测序深度及相对丰度 (图1f)。总结来说，在StrainScan的整个工作流程中，我们新开发的层次聚类树通过快速簇识别提升了菌株识别的计算效率与分辨率，而后续的簇内部识别方法则进一步提升了菌株识别的分辨率。

参考数据库构建

给定参考基因组：

一、cluster search tree(CST)

CST定义：完整二叉树的拓扑结构；每个节点都有一组kmers集合。

1. CST拓扑结构
   （1） 使用Dashing计算菌株基因组之间的kmer jaccard similarity
   （2） 基于距离矩阵结果进行层次聚类(single linkage)，得到树状图的拓扑结构
   （3） 选择cutoff，默认0.95（约为99.89% ANI），将树状图切割为许多cluster，这些cluster由一个或多个菌株组成。
   （4） 保留树的拓扑结构，不保留距离信息，每个cluster作为树的叶子节点，形成一棵完整的二叉树。
2. 每个节点分配kmers

对于节点v，以该节点为根的子树Tv，分配给节点v的kmers遵循两个标准：Tv叶子节点内大部分菌株共享；kmers是唯一分配给Tv中的菌株。

（1）叶子节点最初的kmers集合定义Lv 

（2）递归兄弟节点之间的共享kmers移向父节点 

（3）在节点中移除出现超过一次的kmers

最后，每个节点v都得到一组唯一的kmers集合，表示为Kv，具体定义如下，可以看上图B更直观地理解。

但是，有一个很大的问题是，每个节点得到的一组kmers有时候会很少，这就非常容易匹配，这样子会导致假阳性。这些节点分别被定义为弱节点和强节点，cutoff为1000，通过比较reads分布与kmer分布得到的。
进一步对这些弱节点进行增强，在基于广度优先搜索的策略下，只需要确保相同深度节点不存在重复的kmer即可，因此对于弱节点，可以容忍其包含PREv部分的节点的kmers。

第二步：基因组特征kmers矩阵

如下图，使用具有区分能力的kmers：

（1）strain-specific kmers 来自菌株特异性的区域 

（2）group-specific kmers  可能来自某些菌株常见的SV 

（3）joint kmers 包含所有基因组中存在的核心基因组区域的SNV和indel

使用Sibeliaz软件提取局部共线性块并对结果进行处理得到kmers矩阵X=M×N。这个kmers矩阵行（M）是kmers，列（N）是菌株基因组。只有0和1两个值，反映相应的kmer是否存在于菌株基因组中。

菌株识别

给定二代测序数据：

1. 在CST搜索cluster
   （1）首先提取CST中所有的kmer，使用Jellyfish软件快速计算二代测序数据的kmers并与之匹配。
   （2）kmer match count映射回CST，每个节点分配一维向量Cv
   （3）广度优先，逐层检查
   计算匹配的kmers分数fracv和平均kmers匹配计数abundv，然后基于二项分布检验来决定是否遍历节点v的后代。二项检验目标是区分kmer match是由于测序错误引起还是真正菌株中真正的匹配。
   （4）如果只有一个叶子节点被识别，可以直接回溯节点v到根节点；如果找到多个叶子节点，首先会找有具有识别节点的最大子树，然后再回溯。

最后计算匹配的kmers加权平均分数Fv和加权平均kmers匹配计数Av，这里计算了路径的所有节点。Fv>0.4被认为这个cluster存在。至此完成了第一步识别，也就是找到cluster。

2. 通过kmers矩阵确定共存菌株 这一步中指的kmers与上一步是不同的。
   （1）使用Jellyfish计算矩阵中的kmers在测序数据中的计数→向量y
   （2）对于每个菌株基因组，计算其得分score fj=X[:,j]*j
   （3）根据所有菌株的得分进行排序，取第一个菌株，对于这个菌株中出现的kmers重置为0，迭代，直到一个阈值（kmers总数）。

最后得到一个新的矩阵X’=M×N’，行（M）是kmers，列（N’）是已识别的菌株。

预测测序深度和相对丰度

使用弹性网络回归模型预测测序深度，矩阵X′，y是kmers在测序数据中计数，因此β就是相当于真实的测序深度，残差平方和和两个正则化项组成损失函数，通过最小化损失函数得到最优的β’，也就是预测的测序深度。对β’标准化后可以得到相对丰度。

感谢阅读！

由于作者水平有限，

论文解读难免存在不准确，欢迎批评指正。