#sketch #min-hash #bioinformatics

已删除 finch_lib

基因组数据的最小感知独立排列局部敏感哈希('MinHashing')的实现,以及用于操作的命令行工具

0.4.0 2021年1月22日

#11 in #min-hash

MIT 许可证

190KB
4.5K SLoC

Finch

CI DOI

Finch 是基因组数据的最小感知独立排列局部敏感哈希("MinHashing")的实现。此仓库提供了一个库和命令行界面,用于重新实现 One Codex 的 现有内部聚类/序列搜索工具(并添加了新功能/扩展!)在 Rust 中。

入门指南

安装

您可以从源代码构建 Finch,需要 Rust >= 1.43。Rust 的 Cargo 包管理器(有关 Cargo 安装说明,请参阅 rustup)可以自动构建和安装 Finch,使用 cargo install finch。如果您需要 Python 绑定,您必须采取额外步骤(请参阅 Python 支持)。或者,下载预构建的二进制文件

示例用法

要开始使用,我们首先为几个 FASTA 或 FASTQ 文件计算草图。这些草图是底层基因组数据的紧凑、采样表示,允许 finch 快速估计数据集之间的距离。使用 finch sketch 命令进行草图文件

finch sketch example.fastq example2.fastq

生成的草图文件(example.fastq.skexample2.fastq.sk)然后可以与其他 finch 命令(以及其他 MinHash 实现1)一起使用。请注意,Finch 的所有命令都可以接受草图或原始序列文件。如果传递后者,finch 将在飞行中绘制草图。然而,在飞行中生成的草图不会保存,因此如果您打算多次使用草图,应调用 finch sketch

一旦绘制了草图,就可以比较多个测序运行以确定它们有多相似

finch dist example.fastq.sk example2.fastq.sk

这将打印包含一些关键距离统计数据的结果(JSON 格式),包括 containmentjaccard 相似度分数以及 mashDistance 距离估计

[
  {
    "commonHashes": 30,
    "containment": 0.03,
    "jaccard": 0.015228426395939087,
    "mashDistance": 0.1669789474914277,
    "query": "example2.fastq",
    "reference": "example1.fastq",
    "totalHashes": 1000
  }
]

在这种情况下,这些文件之间的估计距离为 ~0.17,包含度为 0.03(即两个FASTQ共享其min-mers的3%)。请注意,使用较大的 --n-hashes 参数重新计算sketch可以提供对高度相似数据集的额外分辨率。

接下来,我们可能想找到与我们的示例FASTQ 在整个RefSeq中最接近的基因组。为此,我们只需将示例查询文件作为第一个参数,包含RefSeq中所有基因组的sketch文件作为第二个参数(有关与finch一起工作的预计算RefSeq数据库,请参阅示例数据部分)。

finch dist ./example.fastq.sk ./refseq_sketches_21_1000.sk --max-dist 0.2

在这里,我们还将最大距离设置为0.2,以便过滤掉不太相关的基因组(距离为0表示完全相同的基因组)。设置最大值可以确保只返回相关结果 -- 省略此参数将返回到RefSeq中所有基因组的距离。

注意:以下每个命令都进一步详细说明,并且可以通过将--help标志传递给每个命令来获取更多信息,例如:finch dist --help

设计目标

我们对Finch有以下三个主要设计目标

  1. 支持跟踪k-mer/minmer 计数

  2. 改进了用于原始序列数据(即FASTQs)的错误过滤;

  3. 改进性能。

支持计数

Finch是第一个支持跟踪每个唯一k-mer/minmer丰度的MinHash实现之一。这允许对数据进行有用的快速解释(例如,使用finch hist可视化测序深度)以及为更高级的过滤(在此实现)和计数感知距离度量(未来研究的领域)提供支持。由于Finch还支持新的Mash JSON兼容格式,我们可以将此计数信息保存下来以供下游处理和比较。

⚠️  注意,尽管Finch支持常见的JSON交换格式,但由于哈希种子值不兼容,仍可能存在不兼容性(Finch使用种子0,而Mash使用42;有关更多讨论,请参阅此问题)。

过滤

Finch的一个关键设计目标是能够准确且稳定地绘制细菌分离株的可变深度测序(即不需要高质量的组装)。虽然Mash包括指定固定绝对截止值(或使用布隆过滤器过滤掉样本中的唯一minmers)的能力,但实际上,这种过滤机制在变化的测序深度、基因组大小和测序错误配置文件上的性能表现不佳。这里实现的Finch过滤方法扩展了原始Mash论文中的几个想法,以及关于聚类临床C. diff分离株的先前工作

Finch包括两类自动过滤(默认情况下用于FASTQs)来应对这一挑战

  1. 自适应,基于计数的过滤:我们的默认实现“过度sketch”一个输入文件(即包含比用户指定的多n倍的minmers)并动态确定一个适当的基于计数的过滤截止值,以去除低丰度的k-mers(这些主要由随机错误组成)。作为后备,此实现永远不会排除可能被认为是错误的minmers的上限比例(k倍指定的--err-filter错误率,默认为21%)。

  2. 链过滤:我们还应用了一个“链过滤”来移除在单一方向中发现的k-mer;在实践中,这些通常是一些测序适配器或其他不具代表性的样本(可能错误地降低两个草图之间的相似度)。

实际上,这些优化显著提高了可变深度隔离群和完成组装之间的距离估计的稳定性和准确性。以下图表显示了高质量细菌组装与代表各种测序深度的FASTQ文件草图的距离:无过滤(cutoff=0)、“天真”过滤(cutoff=1)和Finch的适应性、可变截止过滤。我们使用包含度量,结果越接近100%越好

Accuracy versus sequencing depth for different filtering schemes

值得注意的是,这种策略在测序深度方面表现良好,而固定截止值0和1从未达到接近100%的预期包含(前者在较低的深度上甚至失败得非常严重)。它还对重复错误(即计数大于等于2的k-mer/minmers)具有鲁棒性,这些错误可能在测序深度低至10-20X时开始引起问题。

速度/性能

除了Rust提供的某些内存安全性保证外,我们还在现有实现上看到了相当大的速度提升。理想情况下,散列应该是MinHashing中的速率限制步骤。分析表明finch大约花费了其时间的三分之一在murmurhash3_x64_128上,因此我们应该在这个理论极限的十倍之内。

Mash Mash(过滤后的) Sourmash Finch Finch(过滤后的)
时间 238秒 276秒 518秒 99秒 104秒
最大内存(MB) 1.2 501.1 13.9 21.8 60.0

注意:这些基准在2015年早期Macbook Pro上运行。基准是对4.8Gb FASTQ文件进行草图绘制(SRR5132341.fastq),草图大小为n=10,000。Finch过滤后的结果使用默认选项,而Mash使用-b 500Mb分配一个500Mb的布隆过滤器。基准使用了以下Mash(23776db)和sourmash(5da5ee7)的提交。

用法

Finch支持四种主要操作,其中许多操作具有类似参数。

共享参数

Finch可以接受许多参数来控制草图绘制的方式(这些选项也适用于disthistinfo命令)。

  • -n <N> / --n-hashes <N>控制草图的整体大小(更高的值在比较中提供更好的分辨率)。默认值为1000
  • -k <K> / --kmer-length <K>设置要散列的kmer的大小(更高的值使比较在分类学上更加特定)。默认值为21
  • --seed <S>设置散列的种子。只有在直接导出草图与其他版本Mash算法的比较时才应更改此设置,后者使用非零默认种子。默认值为0

草图绘制后,将执行过滤,并通过几个选项进行控制

  • -f / --filter / --no-filter 用于确定是否应用过滤(如果没有指定,默认对FASTQ文件进行过滤,对FASTA文件不进行过滤)。
  • --min-abun-filter <MIN> / --max-abun-filter <MAX> 设置所有kmer必须存在的绝对最小和最大丰度(包含)。最小过滤将覆盖任何自适应错误过滤猜测(以下)。
  • --err-filter <ERR_VALUE> 是默认的自适应过滤方案。从概念上讲,ERR_VALUE 应该大约是测序仪的错误率或更高。请注意,高错误率测序数据(例如长读测序)通常与MinHashing表现不佳。
  • --strand-filter <V> 设置链过滤截止值。

请注意,如果从oversketch中留下的kmer不足以满足sketch大小,sketch将失败。有两种选项可能会有所帮助

  • --oversketch <N> 可以用于增加oversketch的大小(通常是200倍),并增加过滤版本足够大的可能性。
  • --no-strict 将允许Finch使用小于指定大小的sketch继续进行。请谨慎使用。

finch sketch

finch sketch 将读取FASTA或FASTQ文件并生成一个“sketch”,可以用于进一步。

如果正在读取的文件是sketch且已设置过滤标志(-f/--filter),则将对sketch重新应用丰度过滤器,以允许post-sketch过滤。但是,链过滤只能应用于原始FAST(A/Q)文件,因为在保存sketch时丢失了链级数据。

默认情况下,sketch将在它sketch的FAST(A/Q)文件旁边生成以.sk结尾的文件。这可以通过传递-O(大写O)将其写入标准输出或-o <file>(小写o)将其写入文件来覆盖。要从标准输入读取,请使用文件名-;这允许将文件流式传输到finch,例如cat testfile.fq | finch sketch -o testfile.sk -

如果编辑它们的src/mash/hash.h并将哈希值设置为0或通过设置--seed 42手动覆盖Finch的种子值,则sketch应与原始Mash实现兼容。

finch dist

finch dist 将计算不同sketch之间的Jaccard距离。如果列出的文件是FASTA/Q而不是sketch,Finch将自动使用与sketch中相同的命令行参数将它们sketch到内存中,或者对于第一个文件之后的文件,使用用于sketch第一个文件的参数。

距离和包含性将仅从查询到一组参考的一组查询中进行计算;默认情况下,列表中的第一个sketch将用作查询,所有其他sketch将用作参考。可以通过传递--queries <sketch_1>,<sketch_2手动覆盖此行为,并可以使用其他sketch作为参考。此外,通过传递--pairwise选项将计算所有sketch之间的距离。

由于不同的计数算法和停止标准,距离可能略不同于原始Mash程序和较旧版本finch的计算结果。通过设置-old-dist标志将回退到较旧版本的Finch计算;自版本0.3起,已取消对Mash精确距离计算的支持。

finch hist

finch hist将为每个提供的sketch输出JSON格式的直方图。直方图是每个深度的minmers数量列表,例如,对于一个具有两个minmers的sketch,一个深度为1(第一个位置),一个深度为3(第三个位置),其格式如下:{"sketch_name": [1, 0, 1]}

⚠️  您可以使用以下命令与Matplotlib一起快速获取直方图:finch hist test.fastq.sk | python -c 'import json; import matplotlib.pyplot as plt; import sys; v = json.loads(sys.stdin.read()).values()[0]; plt.plot(range(1, len(v)+1), v); plt.show()'。请注意,finch hist的JSON输出在未来版本中可能会改变,例如,变为更紧凑的{count: value}格式或类似格式。

finch info

finch info将为每个提供的sketch输出格式化的GC%,覆盖度等信息列表。

⚠️  请注意,从这些值返回的结果是近似的,并且用于计算的算法仍然比较粗糙,可能会发生变化。

示例数据

我们已使用此脚本于2017年3月27日绘制了NCBI RefSeq集合的草图(作为草图),并为每个细菌和病毒基因组提供了单独的sketch的tar包。链接:[k=21 和 n=1,000](https://static.onecodex.com/public/finch-rs/refseq_sketches_21_1000.sk.gz),[k=31 和 n=1,000](https://static.onecodex.com/public/finch-rs/refseq_sketches_31_1000.sk.gz),[k=21 和 n=10,000](https://static.onecodex.com/public/finch-rs/refseq_sketches_21_10000.sk.gz),和[k=31 和 n=10,000](https://static.onecodex.com/public/finch-rs/refseq_sketches_31_10000.sk.gz)。

参考文献 & 注释

除了这个版本外,Mash算法还有其他几个实现,它们应该与这个版本兼容/可比,尤其是:

  • Mash - 首个实现和理论论文
  • SourMash - Python中的较新实现;提供了一些实验性功能

注释

Python支持##e

您可以使用以下命令将Finch编译成Python库:

pip install maturin
cd lib
maturin install --cargo-extra-args="--features=python" --release --strip
# or maturin develop, etc
# to cross-compile linux wheels:
cp ../README.md .
# and edit the readme key of config.toml
docker run --rm -v $(pwd):/io konstin2/maturin:master build --cargo-extra-args="--features=python" --release --strip

然后,例如,要计算两个大肠杆菌之间的相似性:

from cli import sketch_file

sketch_one = sketch_file('WIS_EcoB_v2.fas')
sketch_two = sketch_file('WIS_Eco10798_DRAFTv1.fas')
cont, jacc = sketch_one.compare(sketch_two)

Cap'n Proto

src/serialization文件中有一个finch.capnp,以及MinHash模式(https://github.com/marbl/Mash/blob/54e6d66b7720035a2605a02892cad027ef3231ef/src/mash/capnp/MinHash.capnp)的输出

两者都是在安装capnpcargo install capnpc之后生成的,使用以下命令

capnp compile -orust:finch-lib/src/serialization/ --src-prefix=finch-lib/src/serialization/ finch-lib/src/serialization/finch.capnp
capnp compile -orust:finch-lib/src/serialization/ --src-prefix=finch-lib/src/serialization/ finch-lib/src/serialization/mash.capnp

然后搜索和替换以修复crate::路径crate::serialization::

贡献

可以通过GitHub问题报告问题或改进建议。我们很高兴接受和/或指导任何新增或修复(最好以拉取请求的形式提交)。有关我们的行为准则,请参阅CODE_OF_CONDUCT.md

依赖

~5–7MB
~140K SLoC