用 sourmash 搜索序列时,同一对输入可能得到一个不高的相似度,也可能得到百分之百的包含度。这两个结果并不冲突。前者关心两份序列内容总体有多相像,后者关心某一份查询内容有多少能在另一份数据中找到。
本文安装并运行 sourmash 4.8.14,使用公开参考基因组及其截取片段,实际生成草图、搜索和分解混合输入。通过可以独立复算的集合数量,解释比较方向、采样密度、重复序列和参数兼容性。这样看到一个百分比时,就能先确认它的分母,而不是直接把它当成物种一致率。

草图保留了什么,又省略了什么
序列草图将序列中的短片段转换为哈希值,再按一定规则保留其中一部分,用较小的数据表示原始序列内容。后续比较可以围绕这些保留值的交集、并集和各自数量展开,不必每次重新处理所有碱基。
这种做法适合大量序列的初步比较,也需要接受信息抽样带来的限制。草图中的一个共同值并不是一条完整读段,一组共同值也不是逐碱基比对结果。解释输出前,要先知道使用的短片段长度、采样密度以及是否记录重复出现次数。官方草图参数说明
本次主要使用三十一个碱基的短片段和一千的缩放参数,没有启用丰度跟踪。所有比较都保存实际草图文件,驱动还从文件中读取哈希集合,独立计算交集和并集,核对命令行返回的数值。
输入来自哪里,哪些部分是人为设计的
基础输入为公开参考序列 NC_000913.3,实际长度为四百六十四万一千六百五十二碱基。复现驱动检查长度和完整序列的校验值,避免不同下载文件或不同版本被误当成同一输入。NCBI 参考序列记录
我们从它的开头截取一百万碱基作为部分序列,再截取五千碱基作为短片段。另外固定随机种子生成一条与参考等长的人工序列,用于构成内容已知的混合输入。
混合文件包含完整参考和人工序列两个记录。重复文件则把同一个完整参考写两遍。这里的“混合”只是为了检验算法而构造的输入,没有模拟测序误差、覆盖不均或真实物种丰度,也没有给人工序列赋予分类身份。
这组设计提供了清楚的已知关系。一百万碱基片段完整来自参考,混合文件完整包含参考,重复文件增加了记录数量但没有增加新的序列内容。后续结果可以围绕这些关系逐项检查。
安装中遇到的依赖问题怎样处理
本次在 Python 3.10 的独立环境中安装,实际使用的 sourmash 版本为 4.8.14。首次安装在依赖 screed 的构建阶段报出缺少 pkg_resources,尚未进入草图计算。
我们将构建工具限定到仍支持该依赖的版本范围,单独安装 screed,再安装固定版本的 sourmash。修复后版本命令和下面所有实验均完成。
python3 -m venv sourmash-env
source sourmash-env/bin/activate
pip install "setuptools<81" "setuptools-scm<9" wheel
pip install --no-build-isolation screed==1.1.3
pip install sourmash==4.8.14
sourmash --version
这个处理针对本次实际遇到的构建兼容问题。其他环境若直接安装成功,无须为了复制一个历史错误而调整构建工具。保留准确的异常阶段和最终版本,可以帮助判断问题属于软件安装还是数据分析。
先生成可以逐项核对的草图
下面以完整参考和截取片段为例,分别生成草图。实验驱动还会对其他输入执行同样的操作。
sourmash sketch dna -p k=31,scaled=1000 \
--name-from-first -o full.sig full.fna
sourmash sketch dna -p k=31,scaled=1000 \
--name-from-first -o partial.sig partial.fna
实际完整参考草图包含四千四百七十六个哈希值,片段草图包含九百三十个。二者交集也是九百三十,并集为四千四百七十六。这意味着本例片段保留下来的所有内容都能够在参考草图中找到。
这里不能用原始序列长度直接替代草图中的集合数量。参考中可能存在重复短片段,而且采样只保留部分哈希值。本文的百分比由实际草图集合计算,长度比例只用来帮助理解输入的大致规模。
相似度和包含度使用了不同分母
默认搜索使用的 Jaccard 相似度,可以理解为交集大小除以并集大小。本例九百三十除以四千四百七十六,得到约百分之二十点七八。程序实际返回 0.2077747989,与独立计算一致。
sourmash search partial.sig full.sig \
-t 0 -n 0 -o partial_full_jaccard.csv
如果改为查询包含度,分母变成查询草图自身的大小。同样的交集九百三十,除以查询中的九百三十,结果就是一。程序实际返回百分之百。
sourmash search partial.sig full.sig \
--containment -t 0 -n 0 \
-o partial_full_containment.csv
因此,低整体相似度与完全包含可以同时成立。完整参考还有大量不在片段中的内容,这些内容会扩大并集,降低整体相似度,却不会否定查询片段本身已经全部包含在参考中。
在实际任务中,如果想比较两份基因组整体内容是否相近,可以先考虑整体相似关系。如果想检查一个片段或基因组是否出现在更复杂的数据中,就需要特别关注包含方向及其适用条件。官方搜索命令说明
把查询与参考对调,为什么数字又变小
我们实际交换两个文件的位置,再执行包含度搜索。此时查询是完整参考,被搜索的是一百万碱基片段,结果变成约百分之二十点七八。
两份草图和交集都没有改变,改变的是分母。现在需要回答完整参考有多少内容出现在片段里,因此分母使用完整参考的四千四百七十六个值。
| 查询到参考的方向 | 使用的指标 | 实测结果 |
|---|---|---|
| 片段到完整参考 | Jaccard 相似度 | 20.78% |
| 片段到完整参考 | 查询包含度 | 100.00% |
| 完整参考到片段 | 查询包含度 | 20.78% |
| 完整参考到混合文件 | Jaccard 相似度 | 48.64% |
| 完整参考到混合文件 | 查询包含度 | 100.00% |
| 混合文件到完整参考 | 查询包含度 | 48.64% |
整理结果表时,最好同时保留查询名称、参考名称和指标类型。如果只导出一列叫“相似度”的数字,后续很难发现文件方向发生过交换。比较方向是计算定义的一部分,应与参数一起记录。
混合输入让这种区别更加直观
人工序列草图含四千七百二十七个值,与参考草图在本次条件下没有交集。混合文件因此包含九千二百零三个值,等于两份集合的大小相加。
完整参考与混合文件的 Jaccard 相似度约为百分之四十八点六四,因为参考只占并集的一部分。但以完整参考为查询时,包含度仍为百分之百。反过来以混合文件为查询,则只能由参考解释其中约百分之四十八点六四的内容。
这组三个数字来自同一组输入,分别回答不同问题。不能把百分之百包含直接写成“样本只有这个物种”,也不能把接近一半的整体相似度直接解释为“参考有一半位置发生了突变”。
本次还实际运行混合分解,让程序在完整参考与人工序列两个候选之间解释混合草图。
sourmash gather mixture.sig full.sig unrelated.sig \
--threshold-bp 10000 -o gather.csv
输出中两份输入分别解释约百分之四十八点六四和百分之五十一点三六的查询内容,剩余未解释部分为零。这些比例来自本次未跟踪丰度的草图,不是根据读段覆盖计算的物种相对丰度。两条原始序列虽等长,采样后的集合数量也并非必须完全相等。
重复写两遍参考,会让相似度翻倍吗
我们把同一份完整参考写成两个记录,重新生成草图。实际仍保留四千四百七十六个哈希值,与单份参考一致,比较结果为一。
原因是本次使用的是不跟踪丰度的集合。同一短片段出现更多次,不会因此成为新的集合成员。这有助于比较序列内容,也说明该设置不能表达这条序列在输入里出现了多少遍。
处理真实测序数据时,需要根据任务决定是否跟踪丰度。若关心重复出现次数或覆盖相关信息,就应明确检查相关设置与下游指标。不能一边使用无丰度草图,一边把结果当成读段数量的精确汇总。
缩放参数变大,短片段可能没有留下任何值
我们对一百万碱基片段改用一万的缩放参数,实际保留九十九个哈希值。与一千缩放的完整参考进行兼容比较时,程序在共同采样尺度上计算,包含度仍然为一。
但更稀疏的草图会降低用于估计的证据数量。对五千碱基短片段,本次一千缩放保留四个值,十万缩放则一个都没有留下。这个空草图由采样和输入长度共同造成,不能被解释为该片段没有序列或没有生物学来源。
如果小片段、低覆盖或稀有内容是研究重点,建议先检查实际草图大小,再决定采样设置。输出文件存在并不保证其中保留了足够信息。只对文件大小或任务退出状态做检查,可能漏掉这种情况。
缩放也不是越小越好。更密的草图需要保存和比较更多值,实际选择应结合研究目标和数据规模。本次只是展示一个可复算的取舍,没有用几份小输入推导大规模数据的吞吐量。
不同短片段长度会直接造成不兼容
另一组对照把完整参考的短片段长度改为二十一,仍使用一千的缩放参数,然后用三十一长度的查询去搜索它。实际命令以非零状态退出,并提示没有兼容的草图。
这与“找到一个零相似度结果”不同。程序没有执行能够比较的同类集合匹配,因此不能把缺失结果填成零,再当作序列不相关的证据。
实际合并多个来源的草图库时,应先核对短片段长度、分子类型和采样设置。尤其不要仅凭文件后缀相同就认定它们可以混合搜索。对于不兼容、空草图和低于报告门槛的结果,最好分别记录原因。
还要注意输出中估算长度的单位。本次混合分解报告的查询规模为九百二十万三千,正好来自九千二百零三个哈希值乘以一千的采样缩放量,而两条原始序列的长度之和为九百二十八万三千三百零四。这种差异反映估计方式,不能当成程序实际删除了对应数量的碱基。
如果需要定位具体缺失区域、检查突变位置或确认某段序列的连续匹配,就应进一步使用适合的比对方法。草图结果可以帮助挑选待验证对象,但单独一行集合统计无法提供这些逐位置证据。把筛查与后续验证的输入关联保存,才能从一个候选分数继续追查到具体序列。
使用时可以保留一张怎样的检查表
每次生成草图后,可以保存输入序列版本、完整序列校验值、记录数、总长度、参数和实际哈希数量。这样一份结果异常偏小,能够先回到草图本身检查,而不是立即更换数据库。
每次搜索后,再记录查询方向、指标名称、报告门槛和命中对象。如果需要讨论混合样本,还应检查未解释部分以及多个候选之间的内容重叠。单个高分候选不能自动替代对整体样本组成的评估。
sourmash 的优势是让大量序列内容可以通过较小的草图反复比较,适合候选筛查和内容分解。它的限制也来自同一机制,抽样会省略信息,结果依赖参数与参考覆盖,某些输出仍需要后续比对或其他证据验证。
下载包提供完整输入构造、实际草图数量、九组搜索记录及混合分解结果。驱动还独立复算交并集和包含度,便于将命令输出与具体数字对应起来。下载实跑脚本与本次结果
当项目需要对大量样本建立统一草图库、反复筛查参考并保留可复核结果时,可了解云生信的生物信息学分析计算服务。
