mudskipper

快速入门

# convert a bulk RNA-Seq genomic BAM to a transcriptomic BAM for quantification with Salmon
mudkipper bulk --gtf annotation.gtf --alignment genomic.bam --out transcriptomic.bam
# convert a single-cell RNA-Seq genomic SAM to a transcriptomic RAD for quantification with alevin-fry
mudkipper sc --gtf annotation.gtf --alignment genomic.sam --out transcriptomic_dir

# build and store a GTF index; useful if you want to convert multiple BAM/SAM files
mudskipper index --gtf annotation.gtf --dir-index gtf_index
# run mudskipper using the pre-built GTF index
mudkipper bulk --index gtf_index --alignment genomic.bam --out transcriptomic.bam
mudkipper sc --index gtf_index --alignment genomic.sam --out transcriptomic_dir

目录

• 快速入门
• 介绍
• 从源代码构建
• 界面
  • 大量RNA-Seq读对齐投影
  • 单细胞RNA-Seq读对齐投影
  • 构建和存储GTF区间树
• 用例
  • 从基因组对齐进行大量RNA-Seq样本的转录量计算
  • 从基因组对齐进行单细胞RNA-Seq样本的转录量计算
• 局限性

介绍

mudskipper 是一个将基因组BAM/SAM文件转换为转录组BAM/RAD文件的工具。更具体地说，它根据GTF格式的给定转录注释将每个对齐条目的基因组坐标投影到转录组坐标，如果基因组对齐与多个转录重叠，则生成多个对齐条目。

虽然 mudskipper 被设计为通用的对齐转换/投影工具，但当前 mudskipper 的重点是允许从基因组对齐进行转录量计算，而无需将短读段重新映射到转录组（例如，使用 salmon 对大量RNA-Seq样本和 alevin-fry 对单细胞RNA-Seq样本）。开发 mudskipper 的动机是因为存在许多需要将短RNA-Seq读段与参考基因组对齐的工具。另一方面，转录量计算工具通常期望将短RNA-Seq读段与转录组对齐。 mudskipper 允许用户使用基因组对齐而不是从头开始整个过程（通过避免直接将短RNA-Seq读段映射到转录组）来进行转录量计算。

算法细节： mudskipper 解析输入的 GTF 文件，并使用 coitrees crate 构建一个区间树，存储 GTF 文件中所有转录本的外显子坐标。然后，它处理给定输入 BAM/SAM 文件的每个基因组比对。请注意，RNA-Seq 读取的基因组比对可能是剪接的（即，CIGAR 字符串中可能有 N）。因此，mudskipper 使用生成的区间树将每个基因组比对投影到满足以下条件的所有转录本

1. 转录本在基因组上的起始和终止坐标标记的区域完全包含比对
2. 转录本的外显子坐标与基因组比对的剪接相匹配

对于这样的每个转录本，创建一个 BAM 记录来存储到该转录本的正确比对。注意 上述第一个条件意味着不会报告到基因间或内含子区域的投影。这也意味着 mudskipper 目前不报告支持读取，这些读取以重叠方式对转录本进行比对。将来可以放宽此要求以支持重叠比对，并可能允许任意投影。

从源代码构建

mudskipper 是用 Rust 编写的，可以按以下方式构建：[要求： Rust 工具链；最低 rustc 版本 1.51.0]

cargo build --release

这将创建一个位于 target/release/mudskipper 的二进制可执行文件。您可以为此文件添加到环境 PATH 变量以方便使用。这可以通过以下方式暂时完成：

export PATH=`pwd`/target/release/:$PATH

界面

大量RNA-Seq读对齐投影

mudskipper bulk [OPTIONS] --alignment <FILE> (--gtf <FILE>|--index <DIR>) --out <FILE>

OPTIONS:
    -a, --alignment <FILE>      Input SAM/BAM file
    -g, --gtf <FILE>            Input GTF/GFF file
    -i, --index <DIR>           Index directory containing parsed GTF files
    -o, --out <FILE>            Output file name
    -s, --max-softclip <INT>    Max allowed softclip length [default: 50]
    -t, --threads <INT>         Number of threads for processing bam files [default: 1]
    -r, --rad                   Output in RAD format instead of BAM
    -h, --help                  Prints help information
    -V, --version               Prints version information

必需参数

-a, --比对<FILE>

使用此参数可以指定输入比对文件。目前支持 BAM/SAM 格式。此 BAM/SAM 文件应包含短 RNA-Seq 读取与参考基因组比对。

✏️ 此文件中存储的比对可能是剪接的。

-g, --gtf<FILE>

可以使用此参数传递 GTF 格式的基因注释。在这种情况下，区间树是在运行时从 GTF 文件构建的。或者，可以使用 --index 参数。这意味着 --gtf 和 --index 是互斥的。

✏️ 确保该 GTF 文件与短读取对齐的参考基因组版本相对应。如果 GTF 文件中缺少某些目标序列，则将对这些目标序列的比对自动删除。

-i, --索引<DIR>

此参数指定了先前从 GTF 文件创建的预构建区间树的路径（有关如何创建此类索引的更多详细信息，请参阅 构建和存储 GTF 区间树）。如果您希望对多个 BAM/SAM 文件运行 mudskipper，则这很有帮助。请注意，--gtf 和 --index 是互斥的。

-o, --输出<FILE>

输出比对文件的路径。默认情况下，输出比对文件为 BAM 格式。如果传递了 --rad，则输出比对文件将为 RAD 格式。

可选参数

-r, --rad

传递此参数以将比对输出为 RAD 格式而不是 BAM。默认情况下不设置此参数。

-t, --线程<INT>

用于读取和写入 BAM 文件的线程数。默认情况下，此值设置为 1。

-s, --最大-软剪辑<INT>

删除软剪辑基数量超过 INT 的比对。默认情况下，此值设置为 50。

单细胞RNA-Seq读对齐投影

mudskipper sc [OPTIONS] --alignment <FILE> (--gtf <FILE>|--index <DIR>) --out <FILE/DIR>

OPTIONS:
    -a, --alignment <FILE>       Input SAM/BAM file
    -g, --gtf <FILE>             Input GTF/GFF file
    -i, --index <DIR>            Index directory containing parsed GTF files
    -o, --out <FILE/DIR>         Output file name (or directory name when --rad is passed)
    -s, --max-softclip <INT>     Max allowed softclip length [default: 50]
    -t, --threads <INT>          Number of threads for processing bam files [default: 1]
    -r, --rad                    Output in RAD format instead of BAM
    -m, --rad-mapped <FILE>      The name of output rad file; Only used with --rad [default: map.rad]
    -u, --rad-unmapped <FILE>    The name of file containing the number of unmapped reads for each barcode; Only used with --rad [default: unmapped_bc_count.bin]
    -c, --corrected-tags         Output error-corrected cell barcode and UMI
    -h, --help                   Prints help information
    -V, --version                Prints version information

必需参数

-a, --比对<FILE>

使用此参数可以指定输入比对文件。目前支持 BAM/SAM 格式。此 BAM/SAM 文件应包含短 RNA-Seq 读取与参考基因组比对。

✏️ 此文件中存储的比对可能是剪接的。

-g, --gtf<FILE>

此参数可用于传递 GTF 格式的基因注释。在这种情况下，将动态地从 GTF 文件构建区间树。或者，可以使用 --index 参数。这意味着 --gtf 和 --index 是互斥的。

✏️ 确保该 GTF 文件与短读取对齐的参考基因组版本相对应。如果 GTF 文件中缺少某些目标序列，则将对这些目标序列的比对自动删除。

-i, --索引<DIR>

此参数指定了先前从 GTF 文件创建的预构建区间树的路径（有关如何创建此类索引的更多详细信息，请参阅 构建和存储 GTF 区间树）。如果您希望对多个 BAM/SAM 文件运行 mudskipper，则这很有帮助。请注意，--gtf 和 --index 是互斥的。

-o, --输出<FILE/DIR>

输出 BAM 格式比对文件的路径。如果传递了 --rad，则此参数指定包含 RAD 格式以及其他由 alevin-fry 执行转录定量所需的文件的输出目录。

可选参数

-r, --rad

传递此参数以将比对输出为 RAD 格式而不是 BAM。默认情况下不设置此参数。

-m, --rad-比对<FILE>

指定输出 rad 文件的名称。当传递 --rad 时，使用此选项。在这种情况下，该文件将存储在由 --out 指定的目录中。默认情况下，此设置为 map.rad。

-u, --rad-未比对<FILE>

指定包含每个条码未比对读数的文件名称。当传递 --rad 时，使用此选项。在这种情况下，该文件将存储在由 --out 指定的目录中。默认情况下，此设置为 unmapped_bc_count.bin。

-t, --线程<INT>

用于读取和写入 BAM 文件的线程数。默认情况下，此值设置为 1。

-s, --最大-软剪辑<INT>

删除软剪辑基数量超过 INT 的比对。默认情况下，此值设置为 50。

-c, --校正-标签

默认情况下，mudskipper sc 期望在输入 BAM/SAM 文件中看到 CR 和 UR 标签，分别存储细胞条码和 UMI。如果传递此参数，则必须存在 CB 和 UB 标签，分别存储校正的细胞条码和校正的 UMI。

构建和存储GTF区间树

mudskipper index --gtf <FILE> --dir-index <DIR>

OPTIONS:
    -g, --gtf <FILE>         Input GTF/GFF file
    -d, --dir-index <DIR>    Output index directory name
    -h, --help               Prints help information
    -V, --version            Prints version information

当使用 mudskipper bulk 或 mudskipper sc 时，如果传递基因注释文件，mudskipper 将动态构建区间树。此区间树用于查询与给定基因组区域重叠的转录组坐标。您可以通过构建区间树并将其存储以供以后使用来避免动态构建区间树。如果您希望在不同注释下对表示对同一参考基因组对齐的多个 BAM/SAM 文件运行 mudskipper，这将非常有用。

必需参数

-g, --gtf<FILE>

指定 GTF 格式的基因/转录注释文件。

-d, --目录-索引<DIR>

将存储区间树文件的目录路径。

✏️ 如果该目录不存在，则将其创建。

用例

从基因组对齐进行大量RNA-Seq样本的转录量计算

如果我们有一个包含 RNA-Seq 读比对参考基因组的大样本 BAM/SAM 文件，我们无法使用它进行转录定量，例如使用像 salmon 这样的许多最先进的工具。尽管 salmon 具有非常快的嵌入式映射模块，但我们仍然可以通过使用 mudskipper bulk 命令将基因组比对投影到转录组坐标来避免重新映射。这使得我们可以在基于对齐的模式下使用 salmon。更具体地说，mudskipper 将输出 BAM 文件中的转录组对齐，可以通过 -a 或 --alignments 输入选项传递给 salmon quant。

从基因组对齐进行单细胞RNA-Seq样本的转录量计算

mudskipper 还可以促进单细胞RNA测序样本的转录定量，这些样本是通过基因组对齐获得的。 alevin-fry 是一个处理单细胞RNA测序样本的工具，它以RAD格式输出 salmon alevin 的结果作为输入。使用 mudskipper，如果我们有一个包含单细胞RNA测序样本与参考基因组对齐的BAM/SAM文件，我们可以避免运行 salmon alevin，而是可以直接从BAM文件中提取所需信息——以正确的格式。为此，我们使用带有 --rad 选项的 mudskipper sc 命令生成所需的RAD文件。要处理此RAD文件，可以使用正常的 alevin-fry generate-permit-list、alevin-fry collate 和 alevin-fry quant 命令。

局限性

mudskipper 目前仍处于开发初期，有很多改进的空间。到目前为止，mudskipper 只被测试用于转录定量。目前，它有以下已知限制：

• 不报告嵌合对齐。这意味着目前只有当所有段都落在同一目标上时，才会报告对齐。
• 它只报告完全包含在转录中的投影对齐。换句话说，它目前不报告任何悬空对齐。[#10]
• 输出BAM的一些字段和可选标记可能没有正确更新。[#13]
• 对于单细胞样本，它丢弃了在其条形码中包含 N 的读的配对。[#15]

欢迎提交错误报告和建议。