Lib.rs

Randomly subsample sequencing reads or alignments.

Hall, M. B. (2022). Rasusa: Randomly subsample sequencing reads to a specified coverage. Journal of Open Source Software, 7(69), 3941, https://doi.org/10.21105/joss.03941

目录

• 动机
• 安装
  • cargo
  • conda
  • 容器
  • homebrew
  • 发布二进制文件
  • 本地构建
• 使用方法
  • 基本使用 - 读段
  • 基本使用 - 比对
  • 必需参数
  • 可选参数
  • 完整使用方法
• 基准测试
  • 单个长读段输入
  • 配对端输入
• 贡献
• 引用
  • Bibtex

动机

我找不到一个满足我要求的子采样 fastq 读段的工具。我发现的所有策略都存在缺陷，因为它们要么只想对读段进行子采样以达到一定数量或百分比，要么如果它们对覆盖率进行子采样，则假设所有读段的大小都相同（例如 Illumina）。由于我主要处理长读段数据，因此当我想要将文件子采样到特定覆盖率时，由于读段的长度从未被考虑，这就会成为一个问题。 rasusa 解决了这一缺点。

我一直使用的一个权宜之计是使用 filtlong。它足够简单，我只是计算出我需要多少碱基来达到（理论上的）样本覆盖率。例如，我有一个来自 E. coli 样本的 fastq 文件，包含 500 万个读段，我想将其子采样到 50x 覆盖率。我只需要将样本基因组的预期大小（460 万个碱基对）乘以我想要的覆盖率，我就有了我的目标碱基 - 2.3 亿个碱基对。在 filtlong 中，我可以这样做

target=230000000
filtlong --target_bases "$target" reads.fq > reads.50x.fq

然而，从技术上讲，这并不是filtlong的预期功能；它是一个质量过滤工具。最终你得到的是在（理论上的）50x覆盖率下“得分最高”的读取子集。“最高得分”的读取。根据你的情况，这可能是你想要的。然而，你倾向于数据集中最好的/最长的读取——这不是你整个数据集的公平代表。还有可能偏袒产生更长/更高质量读取的基因组区域。De Maio 等人甚至发现，通过随机子采样纳米孔读取，你可以实现比过滤后更好的基因组组装。

因此，根据你的情况，你可能需要一个无偏见的读取子集。如果是这种情况，rasusa可以满足你的需求。

安装

总结：预编译的二进制文件

curl -sSL rasusa.mbh.sh | sh
# or with wget
wget -nv -O - rasusa.mbh.sh | sh

你可以像这样传递选项到脚本中

$ curl -sSL rasusa.mbh.sh | sh -s -- --help
install.sh [option]

Fetch and install the latest version of rasusa, if rasusa is already
installed it will be updated to the latest version.

Options
        -V, --verbose
                Enable verbose output for the installer

        -f, -y, --force, --yes
                Skip the confirmation prompt during installation

        -p, --platform
                Override the platform identified by the installer [default: apple-darwin]

        -b, --bin-dir
                Override the bin installation directory [default: /usr/local/bin]

        -a, --arch
                Override the architecture identified by the installer [default: x86_64]

        -B, --base-url
                Override the base URL used for downloading releases [default: https://github.com/mbhall88/ssubmit/releases]

        -h, --help
                Display this help message

cargo

先决条件：[href="https://www.rust-lang.net.cn/tools/install" rel="noopener ugc nofollow"]rust工具链（至少v1.74.1）

cargo install rasusa

conda

先决条件：[href="https://docs.conda.org.cn/projects/conda/en/latest/user-guide/install/" rel="noopener ugc nofollow"]conda（以及bioconda通道正确设置）

conda install rasusa

感谢Devon Ryan（[href="https://github.com/dpryan79" rel="noopener ugc nofollow">@dpryan79帮助调试bioconda配方。

容器

Docker镜像托管在[.href="https://quay.io/repository/mbhall88/rasusa" rel="noopener ugc nofollow"]quay.io。对于0.3.0和更早的版本，镜像托管在[.href="https://hub.docker.com/r/mbhall88/rasusa" rel="noopener ugc nofollow"]Dockerhub。

singularity

先决条件：[href="https://sylabs.io/guides/3.4/user-guide/quick_start.html#quick-installation-steps" rel="noopener ugc nofollow"]singularity

URI="docker://quay.io/mbhall88/rasusa"
singularity exec "$URI" rasusa --help

上面的命令将使用最新版本。如果你想指定一个版本，请使用标签（或提交）如下。

VERSION="0.8.0"
URI="docker://quay.io/mbhall88/rasusa:${VERSION}"

docker

先决条件：[href="https://docs.docker.net.cn/v17.12/install/" rel="noopener ugc nofollow"]docker

docker pull quay.io/mbhall88/rasusa
docker run quay.io/mbhall88/rasusa --help

你可以在quay.io仓库上找到所有可用的标签。注意：0.4.0之前的版本托管在Docker Hub。

本地构建

先决条件：[href="https://www.rust-lang.net.cn/tools/install" rel="noopener ugc nofollow"]rust工具链

git clone https://github.com/mbhall88/rasusa.git
cd rasusa
cargo build --release
target/release/rasusa --help
# if you want to check everything is working ok
cargo test --all

使用方法

基本使用 - 读段

子采样fastq读取

rasusa reads --coverage 30 --genome-size 4.6mb in.fq

上述命令将输出子采样文件到stdout。

或者，如果你有配对的Illumina

rasusa reads --coverage 30 --genome-size 4g -o out.r1.fq -o out.r2.fq r1.fq r2.fq

有关上述选项和更多选项的更多详细信息，请参阅以下内容。

基本使用 - 比对

子采样比对

rasusa aln --coverage 30 in.bam | samtools sort -o out.bam

这将将对齐中的每个位置子采样到30x覆盖率。

必需参数

运行rasusa reads需要三个必需选项。

输入

此位置参数指定包含你想要子采样的读取或比对文件的文件（s）。对于reads命令，文件必须是有效的fasta或fastq格式，并且可以是压缩的（使用例如gzip的工具）。对于aln命令，文件必须是有效的索引SAM/BAM文件。
如果向reads传递两个文件，则rasusa将假设它们是配对端读取。

Bash巫师提示🧙：让glob为你工作 r*.fq

覆盖率

-c，--coverage

如果为reads提供了--bases，则不需要。

此选项用于确定最小覆盖度，以便对读段进行子采样。对于reads命令，它可以指定为一个整数（100）、十进制/浮点数（100.7）或者以上两种的前缀加上一个'x'（100x）。对于aln命令，它仅限于整数。

由于在reads命令中确定实现所需覆盖度所需基数的算法，实际覆盖度最终可能会略高于请求的值。例如，如果最后包含的读段非常长。日志消息应告知您最终的覆盖度。

对于aln命令，覆盖度是每个对齐位置应存在的最小读段数。如果一个位置读段数少于请求的数量，则该位置的读段将全部包含。此外，还将有（小的）读段数超过请求数量的区域 - 通常局限于读段对齐结束的地方。这是因为当所选读段的对齐结束时，下一个读段将基于它跨越上一个对齐的末尾进行选择。在选择此下一个对齐时，我们优先选择起点最接近上一个对齐末尾的对齐，以确保与上一个对齐的最小重叠。请参阅下面的IGV截图以获得视觉示例。

基因组大小

-g，--genome-size

对于aln不适用

如果为reads提供了--bases，则不需要。

输入的基因组大小也是必需的。它用于确定实现所需覆盖度所需多少个碱基。这可以像您喜欢的那么精确或粗糙。
基因组大小可以通过多种方式传入。作为一个普通的整数（1600），或者使用度量单位后缀（1.6kb）。所有度量单位后缀都可以有一个可选的'b'后缀，可以是小写、大写或混合大小写。因此，“Kb”、“kb”和“k”都会被推断为“kilo”。有效的度量单位后缀包括

• 碱基（b） - 乘以1
• 千（k） - 乘以1,000
• 兆（m） - 乘以1,000,000
• 吉（g） - 乘以1,000,000,000
• 太（t） - 乘以1,000,000,000,000

或者，可以提供一个FASTA/Q索引文件，基因组大小将设置为该文件中所有参考序列的总和。

可选参数

输出

-o，--output

reads

注意：如果将配对Illumina数据传递给reads，则需要此参数。

默认情况下，rasusa将子采样文件输出到stdout（如果只给出一个文件）。如果您希望指定输出文件路径，则使用此选项。

对于Illumina配对文件，必须使用--output两次指定输出 - -o out.r1.fq -o out.r2.fq

输出文件的顺序假定为与输入相同。
注意：输出始终与输入格式相同。您不能以fastq为输入，并要求以fasta为输出。

rasusa reads还将尝试自动推断是否需要压缩输出文件（s）。它通过检测支持的扩展名来完成此操作

• .gz：将使用gzip压缩输出
• .bz 或 .bz2：将输出使用 bzip2 进行压缩
• .lzma：将输出使用 xz LZMA 算法进行压缩

aln

对于 aln 命令，如果写入 stdout，输出文件格式将与输入相同，否则将根据文件扩展名推断。

注意：输出对齐最可能不会排序。您可以使用 samtools sort 对输出进行排序。例如：

rasusa aln -c 5 in.bam | samtools sort -o out.bam

输出压缩/格式

-O，--output-type

reads

使用此选项手动设置用于输出文件（夹）的压缩算法。它将覆盖从输出路径自动检测到的任何格式。

有效选项有

• g：gzip
• b：bzip2
• l 或 x：xz LZMA 算法
• z：zstd
• u：无压缩

aln

使用此选项手动设置输出文件格式。默认情况下，将使用与输入相同的格式，或者如果提供了 --output 路径扩展名，则将猜测格式。有效选项有

• b 或 bam：BAM
• c 或 cram：CRAM
• s 或 sam：SAM

注意：此选项的所有值均不区分大小写。

压缩级别

-l，--compress-level

如果压缩输出，则使用此压缩级别。默认情况下，此值设置为输出压缩类型的默认值。

目标碱基数

-b，--bases

reads 仅

显式设置所需子样本中的碱基数。此选项接受的数字格式与 基因组大小 相同。

注意：如果提供了此选项，则基因组大小和覆盖范围不是必需的，或者如果提供了它们，则将忽略。

读取数

-n，--num

reads 仅

显式设置子样本中的读取数。此选项接受的数字格式与 基因组大小 相同。

当提供配对读取作为输入时，此选项将采样这么多总读取对。例如，当传递 -n 20 r1.fq r2.fq 时，两个输出文件将各有 20 个读取，并且读取 ID 在两者中将是相同的。

注意：如果提供了此选项，则基因组大小和覆盖范围不是必需的。

读取分数

-f，--frac

reads 仅

显式设置子样本中总读取的分数。此选项给出的值可以是浮点数或百分比 - 即，-f 0.5 和 -f 50 都将取一半的读取。

注意：如果提供了此选项，则基因组大小和覆盖范围不是必需的。

随机种子

-s，--seed

此选项允许您指定随机子采样器使用的随机种子。通过显式设置此参数，您可以使输入的子采样可重复。只有当您将来可能希望再次对相同的输入文件进行子采样并希望获得相同的读集时，才需要传递此参数。然而，如果您忘记使用此选项，系统生成的种子将打印到日志输出中，允许您在未来使用它。

详细程度

-v

添加此可选标志将使日志更加详细。默认情况下，日志将生成“info”级别或以上的消息（更多信息请参阅此处）。如果启用了详细程度，您还将获得“debug”级别的日志消息。

完整使用方法

$ rasusa --help
Randomly subsample reads or alignments

Usage: rasusa [OPTIONS] <COMMAND>

Commands:
  reads  Randomly subsample reads
  aln    Randomly subsample alignments to a specified depth of coverage
  cite   Get a bibtex formatted citation for this package
  help   Print this message or the help of the given subcommand(s)

Options:
  -v             Switch on verbosity
  -h, --help     Print help
  -V, --version  Print version

reads 命令

$ rasusa reads --help
Randomly subsample reads

Usage: rasusa reads [OPTIONS] <FILE(S)>...

Arguments:
  <FILE(S)>...
          The fast{a,q} file(s) to subsample.

          For paired Illumina, the order matters. i.e., R1 then R2.

Options:
  -o, --output <OUTPUT>
          Output filepath(s); stdout if not present.

          For paired Illumina pass this flag twice `-o o1.fq -o o2.fq`

          NOTE: The order of the pairs is assumed to be the same as the input - e.g., R1 then R2. This option is required for paired input.
          
  -g, --genome-size <size|faidx>
          Genome size to calculate coverage with respect to. e.g., 4.3kb, 7Tb, 9000, 4.1MB

          Alternatively, a FASTA/Q index file can be provided and the genome size will be set to the sum of all reference sequences.

          If --bases is not provided, this option and --coverage are required

  -c, --coverage <FLOAT>
          The desired depth of coverage to subsample the reads to

          If --bases is not provided, this option and --genome-size are required

  -b, --bases <bases>
          Explicitly set the number of bases required e.g., 4.3kb, 7Tb, 9000, 4.1MB

          If this option is given, --coverage and --genome-size are ignored

  -n, --num <INT>
          Subsample to a specific number of reads

          If paired-end reads are passed, this is the number of (matched) reads from EACH file. This option accepts the same format as genome size - e.g., 1k will take 1000 reads

  -f, --frac <FLOAT>
          Subsample to a fraction of the reads - e.g., 0.5 samples half the reads

          Values >1 and <=100 will be automatically converted - e.g., 25 => 0.25

  -s, --seed <INT>
          Random seed to use

  -v
          Switch on verbosity

  -O, --output-type <u|b|g|l|x|z>
          u: uncompressed; b: Bzip2; g: Gzip; l: Lzma; x: Xz (Lzma); z: Zstd

          Rasusa will attempt to infer the output compression format automatically from the filename extension. This option is used to override that. If writing to stdout, the default is uncompressed

  -l, --compress-level <1-21>
          Compression level to use if compressing output. Uses the default level for the format if not specified

  -h, --help
          Print help (see a summary with '-h')

  -V, --version
          Print version

aln 命令

$ rasusa aln --help
Randomly subsample alignments to a specified depth of coverage

Usage: rasusa aln [OPTIONS] --coverage <INT> <FILE>

Arguments:
  <FILE>
          Path to the indexed alignment file (SAM/BAM/CRAM) to subsample

Options:
  -o, --output <FILE>
          Path to the output subsampled alignment file. Defaults to stdout (same format as input)

          The output is not guaranteed to be sorted. We recommend piping the output to `samtools sort`

  -O, --output-type <FMT>
          Output format. Rasusa will attempt to infer the format from the output file extension if not provided

  -c, --coverage <INT>
          The desired depth of coverage to subsample the alignment to

  -s, --seed <INT>
          Random seed to use

      --step-size <INT>
          When a region has less than the desired coverage, the step size to move along the chromosome to find more reads.

          The lowest of the step and the minimum end coordinate of the reads in the region will be used. This parameter can have a significant impact on the runtime of the subsampling process.

          [default: 100]

  -h, --help
          Print help (see a summary with '-h')

  -V, --version
          Print version

基准测试

“时间飞逝如箭；苍蝇飞舞如香蕉。”
—— 安东尼·G·奥廷格

真正的问题是：rasusa 是否会无端地浪费您在地球上的宝贵时间？

为了进行此基准测试，我将使用hyperfine。

我使用的数据来自

Bainomugisa, Arnold, 等人。 "完全高质量 MinION 纳米孔结核分枝杆菌北京亚型菌株的组装，鉴定了重复 PE/PPE 基因区域的新的变异。" 微生物基因组学 4.7 (2018)。

注意，这些基准测试仅针对 reads，因为没有其他工具能够复制 aln 的功能。

单个长读段输入

下载并重命名 fastq

URL="ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR649/008/SRR6490088/SRR6490088_1.fastq.gz"
wget "$URL" -O - | gzip -d -c > tb.fq

文件大小为 2.9G，共有 379,547 条读数。
我们使用与动机中概述的相同策略，与 filtlong 进行基准测试。

TB_GENOME_SIZE=4411532
COVG=50
TARGET_BASES=$(( TB_GENOME_SIZE * COVG ))
FILTLONG_CMD="filtlong --target_bases $TARGET_BASES tb.fq"
RASUSA_CMD="rasusa reads tb.fq -c $COVG -g $TB_GENOME_SIZE -s 1"
hyperfine --warmup 3 --runs 10 --export-markdown results-single.md \
     "$FILTLONG_CMD" "$RASUSA_CMD" 

结果

总结： rasusa 运行速度比 filtlong 快 21.77 ± 0.29 倍。

配对端输入

下载后，使用pyfastaq解交织 fastq。

URL="ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR648/008/SRR6488968/SRR6488968.fastq.gz"
wget "$URL" -O - | gzip -d -c - | fastaq deinterleave - r1.fq r2.fq

每个文件的大小为 179M，共有 283,590 条读数。
对于此基准测试，我们将使用 seqtk。我们还将测试 seqtk 的 2-pass 模式，因为这类似于 rasusa reads。

NUM_READS=140000
SEQTK_CMD_1="seqtk sample -s 1 r1.fq $NUM_READS > /tmp/r1.fq; seqtk sample -s 1 r2.fq $NUM_READS > /tmp/r2.fq;"
SEQTK_CMD_2="seqtk sample -2 -s 1 r1.fq $NUM_READS > /tmp/r1.fq; seqtk sample -2 -s 1 r2.fq $NUM_READS > /tmp/r2.fq;"
RASUSA_CMD="rasusa reads r1.fq r2.fq -n $NUM_READS -s 1 -o /tmp/r1.fq -o /tmp/r2.fq"
hyperfine --warmup 10 --runs 100 --export-markdown results-paired.md \
     "$SEQTK_CMD_1" "$SEQTK_CMD_2" "$RASUSA_CMD"

结果

总结： rasusa reads 比 seqtk (1-pass) 快 1.84 倍，比 seqtk (2-pass) 快 1.77 倍。

因此，rasusa reads 比 seqtk 快，但不需要固定数量的读数 - 允许您避免进行数学计算以确定需要下采样到特定覆盖率的读数数量。🤓

贡献

如果您想帮助改进 rasusa，您非常受欢迎！

要接受更改，它们必须通过 CI 和覆盖率检查。这包括

• 使用 rustfmt 格式化代码。这可以通过在项目目录中运行 cargo fmt 来完成。
• 没有编译错误/警告。您可以通过运行以下命令来检查：cargo clippy --all-features --all-targets -- -D warnings
• 代码覆盖率没有降低。如果您在推送更改之前想检查覆盖率，我使用的是 kcov。

引用

如果您在研究中使用了 rasusa，非常感谢您能引用它。

Hall, M. B. (2022). Rasusa: Randomly subsample sequencing reads to a specified coverage. Journal of Open Source Software, 7(69), 3941, https://doi.org/10.21105/joss.03941

Bibtex

您可以通过运行 rasusa cite 获取以下引用。

@article{Hall2022,
  doi = {10.21105/joss.03941},
  url = {https://doi.org/10.21105/joss.03941},
  year = {2022},
  publisher = {The Open Journal},
  volume = {7},
  number = {69},
  pages = {3941},
  author = {Michael B. Hall},
  title = {Rasusa: Randomly subsample sequencing reads to a specified coverage},
  journal = {Journal of Open Source Software}
}