VirStrain：RNA病毒的菌株识别工具

VirStrain: a strain identification tool for RNA viruses

Genome Biology [IF: 13.583]

DOI：https://doi.org/10.1186/s13059-022-02609-x

发表日期：2022-01-31

第一作者: Herui LIAO(廖和睿)¹

通讯作者：Yanni SUN(孙燕妮)(yannisun@cityu.edu.hk)¹

合作作者: Dehan Cai

主要单位：

¹香港城市大学(Department of Electrical Engineering, City University of Hong Kong, Kowloon, China)

摘要

病毒在复制过程中不断变化，这些变化导致了病毒种内的高度多样性。虽然这些变化中有许多变化对于病毒本身是中性或者有害的，但有些变化却可以赋予病毒不同的生物特性，如更好的适应性与更强的传播力。此外，病毒的基因组往往附带相关的元数据，比如其宿主居住地，到访地点等，这些信息有助于推断疫情中的病毒传播路径。因此，病毒的菌株层面分析可以为病毒的定性以及流行病学研究提供重要的信息。

在这项工作中，我们开发了VirStrain，一个以二代测序数据与参考数据库为输入，以 (已知的) 病毒菌株组成为输出的菌株识别工具。我们在多个模拟与真实数据集上测试了VirStrain的性能。结果表示，VirStrain在灵敏度和准确度方面都优于目前较为流行的菌株识别工具或病毒单倍型重构工具。

背景

RNA病毒在复制过程中通常缺乏严格的校正机制，这导致新产生的菌株与原始菌株具有不同的遗传变异。这些变异中有许多对病毒的生存是中性或者有害的。然而，有些变异对病毒的生存能力则是有益的。此外，通过测序得到的病毒菌株基因组往往具有相关的元数据，比如其宿主的感染时间，居住地，性别，种族和药物使用情况等，这些数据对推断病毒的进化以及疫情中病毒的传播非常重要。因此，病毒的菌株组成分析可以为病毒的定性，控制以及疫苗设计提供重要启示。

近年来，随着测序技术的飞速发展，研究人员已经可以以极高的分辨率进行物种组成分析。然而，大部分流行的菌株层面识别工具都是为细菌设计的，且无法直接应用到病毒上。因此对于病毒，研究人员通常选择从头组装出不同病毒菌株/单倍型 (haplotype)的完整基因组，而非进行菌株识别。大量研究以及我们后续的实验已经证实，该策略有两个局限性。首先，大部分病毒单倍型重构工具在重构低丰度菌株或高度相似菌株时 (如SARS-CoV-2的病毒菌株)，效果往往不尽人意。 其次，一些病毒单倍型重构工具，如CliqueSNV，尽管可以重构出高度相似的菌株，但却需要大量的计算资源。因此，利用已有的病毒菌株基因组进行菌株组成分析为菌株层面的病毒研究提供了一种新的可能。

病毒菌株识别的主要挑战来自于菌株基因组间的高度相似性。在真实情况中，具有不同生物学特性的病毒菌株仍可能具有非常高的序列相似性。例如，SARS-CoV-2的Alpha，Beta，Gamma和Delta毒株具有不同的传播性，疾病严重程度和再感染的风险，但是它们的基因组中>99.5%的序列都是相同的。更为详细地说，对于基因组长度为~30000 bp的SARS-CoV-2，上述4种不同毒株间差异的突变数仅有50-70个。由于菌株间的高度相似性，如果直接将测序数据比对到所有的病毒菌株基因组，那么大部分序列可能会比对上许多个菌株基因组，导致很难判断实际存在的菌株。部分现有的菌株识别工具通过特定算法来解析这些比对结果，从而推断可能存在的菌株，然而这需要消耗大量的计算资源。另一部分基于k-mer的方法则十分高效，比如Krakenuniq，通过某些菌株特有的unique k-mer来进行菌株的识别，但是，高度相似性使得许多菌株并不拥有unique k-mer，从而导致Krakenuniq无法识别到这些菌株。此外，对于病毒而言，多菌株感染的情况并不罕见。然而，无论是基于比对还是基于k-mer的工具，都很难区分这些高度相似的菌株。

为了解决该难题，VirStrain利用了一种新的特征: unique k-mer组合来进行菌株识别。一些高度相似的菌株虽然没有unique k-mer，但是它们可能拥有unique k-mer组合，利用该特征，就有可能区分开高度相似的病毒菌株。图1展示了如何利用unique k-mer组合来区分5条高度相似的序列。

图1 使用unique k-mer组合来区分五条高度相似的序列。每个序列都有一个unique k-mer组合

由图1可知，虽然右边每个k-mer都不是unique的，但是，这些k-mer的组合对于每条序列来说，却是unique的，该信息完全可以被用于菌株识别。所以，为了从菌株基因组中找到这些unique k-mer组合，我们开发了一种全新的贪婪覆盖算法，可从菌株基因组的多比对结果中提取到每个菌株独特的SNV位点组合，进而获得这些unique k-mer组合。

在这项工作中，我们开发了VirStrain，用于从二代测序数据中快速且准确地识别可能存在的病毒菌株及预测其相对丰度。VirStrain通过利用新的特征: unique k-mer组合，以及利用该特征进行菌株识别的迭代菌株搜索算法，提升了病毒菌株识别的准确度与分辨率。此外，VirStrain是一款不依赖于序列组装的工具，所以，当因数据质量过差或测序深度太低而无法组装出完整的病毒菌株基因组时，用户依然可以使用VirStrain来识别样本中可能存在的菌株。VirStrain的输出包含两部分，一部分是文本格式的报告，包含了识别到的病毒菌株，及其对应的相对丰度，元数据等。另一部分则是可交互的Html网页，主要显示每个识别到的菌株的基因组位点覆盖度。

结果

VirStrain在模拟数据集上的表现

Performance on the simulated datasets

在该实验中，我们主要比较了VirStrain与其它流行的菌株识别工具 (Kraken2、Krakenuniq、Centrifuge、Pathoscope2、Sigma) 及病毒单倍型重构工具 (ShoRAH、PredictHaplo、aBayesQR、CliqueSNV、HaROLD、TenSQR) 在三种病毒 (SARS-CoV-2、H1N1、HIV) 的模拟数据集上的表现。这三种病毒的菌株基因组均从公共数据库下载，并输入到所有的菌株识别工具用以建立参考数据库。这里用到的模拟数据分成两种类型，一种是只包含一株菌株，测序深度为100X的单菌株数据，另一种是包含两株菌株的多菌株数据，而多菌株数据中的两株菌是按照测序深度100X:10X模拟生成的。对于这三种病毒，每种类型的模拟数据都包含100组，所以总共有10032=600组模拟数据。为了确保公平性，每个模拟数据中的菌株，均存在于所有工具的参考数据库中。对于菌株识别工具而言，我们主要测试其是否能够从模拟数据中识别到正确的菌株。而对病毒单倍型重构工具而言，我们主要测试其重构的菌株是否和数据库中的正确菌株最为相似，如果是，则视为正确的识别，反之则是错误的识别。在这里，我们使用准确率 (Accuracy)
作为主要的评估指标。准确率表示有多少个菌株被正确的识别。

如图2所示，VirStrain在所有的测试数据集上都有非常高的准确率。尤其是在SARS-CoV-2的多菌株模拟数据集上，VirStrain的准确率比第二名的Sigma高出了40%，足以凸显VirStrain在区分高度相似菌株上的优势。对于这三种病毒，我们不难发现，大部分工具对于HIV的模拟数据都有较好的表现，而在SARS-CoV-2的模拟数据上表现很差。这是因为三种病毒具有不同的菌株序列相似性。 HIV菌株间的相似度相对较低，故更容易识别。SARS-CoV-2的菌株之间则普遍有相当高的序列相似度。而对于H1N1，则是部分菌株极为相似但部分菌株差异较大

图2 十二种工具在模拟数据集上的菌株识别准确率

12种工具的精度比较。单菌株数据集有 100 组模拟读长，多菌株数据集有 100 组模拟读长。对于每组多菌株模拟读长，有两个菌株分别具有 100X 和 10X 的覆盖率。

对于多菌株感染的情况，不同菌株的相对丰度也是一项非常重要的指标。因此，我们进一步比较了不同工具在100组SARS-CoV-2多菌株模拟数据上预测相对丰度的性能。由图2可知，大部分工具在SARS-CoV-2多菌株模拟数据集上的表现都非常差，因此，我们只选取了那些准确率大于50%的工具进行比较，即VirStrain、Sigma、Centrifuge、HaROLD、TenSQR。

图3中的结果显示，相较于其它工具，VirStrain所预测的相对丰度更接近于正确答案。 除VirStrain外，剩余工具在很多测试数据集中都未能检测到低丰度菌株 (10X)，且给出了偏离真实答案较远的相对丰度预测。

图3 五种工具在100组SARS-CoV-2多菌株模拟数据中所预测的菌株丰度

100 个模拟 SARS-CoV-2 双菌株数据集的预测丰度分布气泡图。每个圆圈的中心代表一个工具输出的两种菌株的相对丰度。当一个工具在多个数据集上产生相同的丰度分布时，我们使用一个圆圈表示相同的输出，其大小代表这些数据集的数量。“Truth”是指每个数据集中两个菌株的真实相对丰度，通过标准化测序深度（100X 和 10X）计算得出。它的圆圈包含 100 个数据集(样本)。 许多圆的中心 x 坐标为 0，这意味着这些工具仅输出一个菌株。

VirStrain可以检测新型菌株的最近亲属

VirStrain detects the closest relative for novel strains

病毒的突变是非常迅速的，因此，很容易诞生之前从未有过的新型菌株。对于新型菌株，由于其参考基因组不在已有的参考数据库中，所以当使用菌株识别工具进行识别时，我们期望其能找到数据库中与新型菌株最为相似的菌株，这种菌株我们称之为最近亲属 (closest relative)。为了测试不同工具在检测最近亲属上的性能，我们挑选了53株所有工具都能正确识别的SARS-CoV-2菌株，并基于这些菌株的原始基因组模拟了具有5、7、9、11、13个随机点突变的突变基因组，我们将这些新获得的突变基因组标记为M5、M7、M9、M11、M13。最终，基于这些不存在于参考数据库中的突变基因组，我们生成了53 5=265组单菌株模拟数据集，再加上包含原始菌株 (标记为Raw) 的53组单菌株模拟数据集，该实验总共包含265+53=318组模拟数据。对于每个突变基因组，我们使用了MegaBLAST来搜索参考数据库中与之最为相似的菌株，并把该菌株视作其最近亲属。不同菌株识别工具在该数据集上检测最近亲属的准确率如*图4所示。

图4 五种工具在318组模拟数据集中检测最近亲属的准确率

5 种工具在 318 个模拟数据集中检测最近亲的准确度比较。。本实验中使用的 53 株菌株可以被图 3 中的所有工具正确识别。“Raw”表示来自参考基因组的数据集，M5、M7、M9、M11和M13表示来自突变株的模拟数据集。每组共有 53 个数据集，每个数据集包含一个菌株。

图4中的结果表明，VirStrain和Sigma在所有测试数据集上都能找到正确的最近亲属，而其它工具在某些数据集中则未能找到正确的最近亲属。 尽管在检测最近亲属上都有非常高的准确率，VirStrain所需的运行时间却远远小于基于序列比对的Sigma，我们会在下一章节展开讨论。

运行时间比较

Running time comparison

在该实验中，我们比较了VirStrain与其它菌株识别工具在不同数据集上的运行时间。在运行时，所有测试的工具都使用了8个线程。为了评估所有工具在大数据集上的计算效率，我们特意使用了一个包含已知SARS-CoV-2菌株的真实宏基因组测序数据集 (SRR10971381)，如果某工具能正确识别到该数据中存在的菌株，则表中显示为“Y”，否则为“N”。如图5所示，VirStrain相较于那些计算高效的工具如Krakenuniq和Centrifuge拥有更高的准确率，而相比于那些基于序列比对的工具如Sigma，VirStrain则具有更短的运行时间。

图5 五种工具在不同数据集上的运行时间

Sim_single_strain 和 Sim_multiple_strain 分别代表模拟的单菌株和多菌株数据集。对于真实数据，识别结果用 Y 和 N 表示，其中 Y 表示正确识别，N 表示错误识别。由于Sigma运行时间较长，没有识别结果。

VirStrain用于检测真实数据中的SARS-CoV-2菌株

VirStrain detects SARS-CoV-2 strains from real sequencing data

为了评估VirStrain在检测SARS-CoV-2菌株中的表现，我们在32个真实数据上进行了实验。这32个样本都是从被SARS-CoV-2感染的人或动物中采样，其中有8个样本的病毒菌株是没有参考基因组的，其余样本的菌株都存在参考基因组。对于有菌株基因组的样本，我们使用MegaBLAST来搜索其在参考数据库中的最近亲属，若没有菌株基因组，则标记为“Unknown”。但由于某些菌株的参考基因组不完整，所以MegaBLAST的结果会出现错误比对或比对上多条序列的情况。此外，我们也记录每个样本附带的地区来源信息，并将该信息与VirStrain识别到的菌株的宿主地区信息进行比较。由图6可知，对于大部分样本，VirStrian识别到的菌株和MegaBLAST搜索到的最近亲属相同，且对应的宿主地区信息都和其附带的地区信息保持一致。这表明VirStrain的准确度已经接近以基因组为输入的高准确度比对工具。而在部分参考基因组质量较差或没有参考基因组的样本中，VirStrain所识别的菌株其宿主地区信息仍和样本附带的信息保持一致。

图6 VirStrain在32个真实SARS-CoV-2测序数据中的菌株识别结果

X 轴：输入的重叠群的长度； Y轴：准确率。

图6中存在两个比较有趣的样本，即加粗的SAMN14560168和SAMN14643484。SAMN14560168是柬埔寨的首位SARS-CoV-2患者，他在入院前曾到访过中国。VirStrain的鉴定结果显示，它的最近亲属是MT470175.1，其宿主来自中国。以上结果表明，这名柬埔寨患者有可能是在中国感染后返回柬埔寨，这与该患者的旅行史相一致。另一个有趣的病例，SAMN14643484，来自以色列。VirStrain的识别结果显示其感染的菌株宿主来自钻石公主号邮轮。而查询公共数据库记录信息后我们发现，该以色列病人确实是钻石公主号的乘客，并在那里感染了SARS-CoV-2。这些实验结果表明VirStrain可以为流行病学研究提供重要信息。

此外，VirStrain也可以被应用于检测SARS-CoV-2的多菌株感染。最近的一项研究报告了一个案例，一名患者同时被两种不同的SARS-CoV-2毒株感染。据该研究描述，两个样本 (Sample1 和 Sample2) 来自同一个病人，采样时间相隔八天，且两个样本中均包含相同但丰度不同的两种SARS-CoV-2菌株。因此，我们将VirStrain、 HaROLD、 CluqueSNV、 TenSQR和aBayesQR应用于这两个样本，来检验是否这些工具能够检测到SARS-CoV-2的多菌株感染。

图7中，我们比较了原始研究 (标记为Report) 与五个工具得到的结果。结果表明，VirStrain是唯一一个能够在两个样本中都识别到两种菌株的工具，并且所预测的丰度和原始研究所报道的保持一致。该实验表明，VirStrain可用于检测多菌株感染，即使是面对高度相似且丰度较低的菌株，VirStrain仍能有较好的表现。

图7 五种工具在两个真实SARS-CoV-2多菌株感染样本中的菌株识别结果

原始研究中报告的丰度(简称“报告”)与5个工具预测的丰度之间的比较。Strain1 属于进化枝 20C，strain2 属于由 nextstrain 定义的进化枝 20B。根据最初的研究，strain1和strain2之间只有8个突变。由于内存使用，TenSQR无法分析“Sample2”，所以丰度为“0%”。类似地，aBayesQR在使用8个线程分析Sample2超过7天后终止，所以丰度也是“0%”

VirStrain从已知组成的测序数据中检测出5株HIV

VirStrain identifies 5 strains from HIV mock data

该实验中，我们将VirStrain应用于一个已知组成的真实测序数据中。 该数据包含五种不同的HIV菌株 (JRCSF、89.6、NL43、YU2和HXB2) 并且每种菌株的相对丰度都各不相同。识别结果如图8所示，VirStrain识别到了所有的HIV菌株且预测的相对丰度非常接近真实答案。该实验结果表明VirStrain在样本中存在超过两株菌株的情况下，依然有较好表现。

图8 VirStrain在HIV真实数据中预测的菌株丰度

Ground truth 和 VirStrain 之间 HIV 模拟数据的丰度比较。按降序排列的真实平均丰度分别为 31.35%、20.04%、18.11%、17.65% 和 12.86%。

扩展：VirStrain用于识别DNA病毒菌株

Extension: VirStrain identifies strains of DNA viruses from mix-strain sequencing data

由于RNA病毒有大量高度相似的菌株基因组，所以我们之前的实验主要是在RNA病毒上评估VirStrain。但事实上VirStrain也可以被应用到DNA病毒的菌株识别，我们通过三个实验来证实了VirStrain在DNA病毒的菌株识别上仍具有较高的准确率。

首先，我们将VirStrain与其它两个计算高效的工具Centrifuge和Krakenuniq应用到了两个包含HBV菌株的样本。这两个样本都包含相同的两株HBV菌株，且具有不同的相对丰度。从图9中可知，三个工具中仅有VirStrain正确识别了所有的菌株。

图9 三种工具在两个HBV真实数据中的菌株识别结果

表中的“-”表示该菌株未被识别，因此丰度未知。

紧接着，我们又使用了10组已知组成的HCMV真实数据测试了VirStrain的性能。这十组数据是由不同HCMV菌株按照不同比例混合然后测序获得的。比如图10中的“TA-1-1”表示该数据是由TB40/E (T) 和AD169 (A) 两种菌株按照1：1的比例混合测序所得。由于HaROLD在其它研究中也被用于重构HCMV菌株基因组，所以该实验中我们将VirStrain与HaROLD进行了比较。

图10 VirStrain和HaROLD在10个HCMV真实数据中的菌株识别结果

表中的“-”表示没有对这些数据集进行HaROLD测试，因为它们只包含一个菌株。这三个菌株分别为TB40/E (T)、AD169 (A)和Merlin (M)。每个样品名称以两个复合菌株的首字母缩写词开头，后跟菌株比。

如图10，在这10个样本中，VirStrain实现了100%的识别准确率。对于两个存在低丰度菌株的样本“TA-1-50”和“TM-1-50”，VirStrain仍能很好地检测出低丰度菌株。对于“TM”组的三个样本，HaROLD只能重构出Merlin菌株的基因组，这与我们之前实验中的结论保持一致，即HaROLD倾向于低估样本中存在的菌株数量。

最后，我们将VirStrain应用到了一个研究HCMV母婴传播的队列，以展示VirStrain在检测DNA病毒多菌株感染方面的潜力。该队列中包含的样本是从五位来自肯尼亚，感染HCMV的母亲及其婴儿身体不同部位，在不同时间点采取的。样本的来源主要有三个：母亲的母乳 (BM)，宫颈 (CV) 以及婴儿的血斑 (BS)。该队列中共包含五个家庭，家庭 (Family) 代表一个母亲和她的婴儿。 为了研究HCMV在母婴间的传播，我们首先使用VirStrain去识别每个样本中存在的HCMV菌株，接着根据菌株间的进化距离进行聚类，把在同一类的菌株划分成同一个菌株“基因型”(genotype) 并用相同的颜色标记。最后，绘制出每个样本HCMV菌株基因型的相对丰度进行分析。

如图11所示，五位母亲的母乳中HCMV的菌株多样性明显高于宫颈和婴儿的血斑。此外，我们可以发现12号和123号家庭的婴儿最初感染的HCMV基因型和母亲羊水中的HCMV基因型相同，而后来感染的HCMV基因型里则包含了只在母亲母乳里才有的HCMV基因型，这说明婴儿后续的HCMV感染有可能是通过母亲的母乳传播给婴儿的。更为有趣的是，VirStrain在母乳中识别到的一些HCMV菌株(MK422176.1, MK290742.1等)，是从来自赞比亚，也感染了HCMV的母亲的母乳中分离出来的。同样，在宫颈中发现的一些菌株(KJ361966.1, KR534203.1等)也是来自羊水中。这些结果表明，VirStrain可以为多菌株感染提供更为全面的菌株识别信息。

图11 每个样本中母乳、宫颈 和血斑的菌株丰度，以及每个家族的菌株簇的最大似然系统发育树

左侧为每个家庭母亲母乳 (BM)，宫颈 (CV) 与婴儿血斑 (BS) 不同时间点HCMV菌株基因型丰度变化，右侧为每个家庭中所识别到的HCMV菌株的进化树。对于同一家庭中识别到的菌株，颜色相同则代表为同一基因型。树中显示的值表示自举值。使用 iTOL 生成系统发育树的可视化。

讨论

在这项工作中，我们提出了VirStrain，通过使用新发明的贪婪覆盖算法，VirStrain可从菌株基因组的多比对结果中提取一种新的特征：unique k-mer组合，再结合利用该特征的迭代菌株搜索算法，VirStrain能够从二代测序数据中精准且高效地识别高度相似的病毒菌株并预测其相对丰度。在大量模拟及真实数据实验中，VirStrain都比现有的菌株识别工具或病毒单倍型重构工具有更高的准确率。尤其是在对高度相似且丰度较低的多菌株识别上，VirStrain相比第二好的工具提升了40%的准确率。对于不在参考数据库中的新型菌株，VirStrain仍能够返回参考数据库中与之最为相似的最近亲属。此外，VirStrain在计算上还十分高效，分析数据量为19.5G的宏基因组数据，仅需不到4分钟。综上所述，该工具可为流行病学研究，病毒多菌株感染检测，精准用药等方面提供十分有价值的信息。我们的下一步工作，是实现网页版VirStrain，使得没有编程背景的用户也能十分轻松的使用VirStrain进行菌株识别。此外，我们也开发了针对细菌的菌株组成分析工具StrainScan，欢迎访问https://doi.org/10.1101/2022.01.26.477962了解详情。

作者简介

第一作者：廖和睿，香港城市大学博士研究生，主要研究方向是宏基因组菌株层面识别工具开发，宏基因组特征在疾病诊断中的应用，宏基因组在线分析网站开发等。目前发表文章有1篇Microbiome、 1篇 Genome Biology等。欢迎关注作者个人Github主页: https://github.com/liaoherui

通讯作者：孙燕妮，香港城市大学电机工程系(Electrical Engineering)副教授，博士生导师。在美国圣路易斯的华盛顿大学(Washington University in Saint Louis)取得计算机系博士学位后在密西根州立大学计算机系担任助理教授和副教授(with tenure)。2018年加入香港城市大学。主要研究方向是生物信息学，序列分析，宏基因组学，和病毒基因组学。具体的研究课题，发表的论文，以及实验室的位置请参加作者个人主页：https://yannisun.github.io/

Reference

Herui Liao, Dehan Cai & Yanni Sun. (2022). VirStrain: a strain identification tool for RNA viruses. Genome Biology 23, 38, doi: https://doi.org/10.1186/s13059-022-02609-x

猜你喜欢

• 10000+: 菌群分析
  宝宝与猫狗 提DNA发Nature 实验分析谁对结果影响大  Cell微生物专刊 肠道指挥大脑
• 系列教程：微生物组入门 Biostar 微生物组  宏基因组
• 专业技能：生信宝典 学术图表 高分文章 不可或缺的人
• 一文读懂：宏基因组 寄生虫益处 进化树
• 必备技能：提问 搜索  Endnote
• 文献阅读 热心肠 SemanticScholar Geenmedical
• 扩增子分析：图表解读 分析流程 统计绘图
• 16S功能预测   PICRUSt  FAPROTAX  Bugbase Tax4Fun
• 在线工具：16S预测培养基 生信绘图
• 科研经验：云笔记  云协作 公众号
• 编程模板: Shell  R Perl
• 生物科普: 肠道细菌 人体上的生命 生命大跃进  细胞暗战 人体奥秘  

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外5000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。

学习扩增子、宏基因组科研思路和分析实战，关注“宏基因组”

点击阅读原文，跳转最新文章目录阅读
https://mp.weixin.qq.com/s/5jQspEvH5_4Xmart22gjMA