Lib.rs

一个高度并行化的工具，用于分析按碱基级别的指标。

如果缺少指标，或性能不足。请在 issues 中提交错误/功能工单。

为什么？

为什么在有许多其他工具的情况下选择 perbase？perbase 利用 Rust 的并发系统来自动并行处理输入区域。这导致运行时间快几个数量级，并且可以扩展到您可用的计算资源。此外，perbase 力求比其他工具更准确。例如：perbase 将 DEL 计入深度，bam-readcount 则不，perbase 不将 REF_SKIP 计入深度，sambamba 则这样做。

安装

conda install -c bioconda perbase
# OR
cargo install perbase

您也可以从 版本页面 下载二进制文件。

工具

base-depth

base-depth 工具遍历 BAM/CRAM 文件中的每个位置，并计算深度以及给定位置每种核苷酸的数目。此外，它还计算每个位置上的 Ins/Dels 数量。

输出列如下

perbase base-depth ./test/test.bam

示例输出

REF     POS     REF_BASE        DEPTH   A       C       G       T       N       INS     DEL     REF_SKIP        FAIL    NEAR_MAX_DEPTH
chr1    709636  T       16      0       0       0       16      0       0       0       0       0   false
chr1    709637  T       16      0       4       0       12      0       0       0       0       0   false
chr1    709638  A       16      16      0       0       0       0       0       0       0       0   false
chr1    709639  G       16      0       0       16      0       0       0       0       0       0   false
chr1    709640  A       16      16      0       0       0       0       0       0       0       0   false
chr1    709641  A       16      16      0       0       0       0       0       0       0       0   false
chr1    709642  G       16      0       0       16      0       0       0       0       0       0   false
chr1    709643  G       16      0       0       16      0       0       0       0       0       0   false
chr1    709644  T       16      0       0       0       16      0       0       0       0       0   false
chr1    709645  G       16      0       0       16      0       0       0       0       0       0   false

如果传递了--mate-fix标志，则每个位置将首先检查是否存在配对重叠，并选择具有最高MAPQ的配对，通过选择首先通过过滤器的配对来解决平局。被丢弃的配对不计入FAIL或DEPTH。

如果提供了--reference-fasta，则将填充REF_BASE字段。参考必须进行索引并匹配输入的BAM/CRAM头。

输出可以按以下方式压缩和索引

perbase base-depth -Z ./test/test.bam -o output.tsv.gz
tabix -S 1 -s 1 -b 2 -e 2 ./output.tsv.gz
# Query all positions overlapping region
tabix output.tsv.gz chr1:5-10

用法

Calculate the depth at each base, per-nucleotide

USAGE:
    perbase base-depth [FLAGS] [OPTIONS] <reads>

FLAGS:
    -Z, --bgzip                     
            Optionally bgzip the output

    -h, --help                      
            Prints help information

    -k, --keep-zeros                
            Keep positions even if they have 0 depth

    -m, --mate-fix                  
            Fix overlapping mates counts, see docs for full details

    -M, --skip-merging-intervals    
            Skip mergeing togther regions specified in the optional BED or BCF/VCF files.
            
            **NOTE** If this is set it could result in duplicate output entries for regions that overlap. **NOTE** This
            may cause issues with downstream tooling.
    -V, --version                   
            Prints version information

    -z, --zero-base                 
            Output positions as 0-based instead of 1-based


OPTIONS:
    -B, --bcf-file <bcf-file>
            A BCF/VCF file containing positions of interest. If specified, only bases from the given positions will be
            reported on
    -b, --bed-file <bed-file>
            A BED file containing regions of interest. If specified, only bases from the given regions will be reported
            on
    -C, --channel-size-modifier <channel-size-modifier>
            The fraction of a gigabyte to allocate per thread for message passing, can be greater than 1.0 [default:
            0.15]
    -c, --chunksize <chunksize>
            The ideal number of basepairs each worker receives. Total bp in memory at one time is (threads - 2) *
            chunksize [default: 1000000]
    -L, --compression-level <compression-level>
            The level to use for compressing output (specified by --bgzip) [default: 2]

    -T, --compression-threads <compression-threads>
            The number of threads to use for compressing output (specified by --bgzip) [default: 4]

    -F, --exclude-flags <exclude-flags>                      
            SAM flags to exclude, recommended 3848 [default: 0]

    -f, --include-flags <include-flags>                      
            SAM flags to include [default: 0]

    -D, --max-depth <max-depth>
            Set the max depth for a pileup. If a positions depth is within 1% of max-depth the `NEAR_MAX_DEPTH` output
            field will be set to true and that position should be viewed as suspect [default: 100000]
    -Q, --min-base-quality-score <min-base-quality-score>
            Minium base quality for a base to be counted toward [A, C, T, G]. If the base is less than the specified
            quality score it will instead be counted as an `N`. If nothing is set for this no cutoff will be applied
    -q, --min-mapq <min-mapq>                                
            Minimum MAPQ for a read to count toward depth [default: 0]

    -o, --output <output>                                    
            Output path, defaults to stdout

        --ref-cache-size <ref-cache-size>
            Number of Reference Sequences to hold in memory at one time. Smaller will decrease mem usage [default: 10]

    -r, --ref-fasta <ref-fasta>                              
            Indexed reference fasta, set if using CRAM

    -t, --threads <threads>                                  
            The number of threads to use [default: 32]


ARGS:
    <reads>    
            Input indexed BAM/CRAM to analyze

only-depth

only-depth工具遍历输入的BAM/CRAM文件，并计算由BED文件或BAM/CRAM头中指定的所有位置的深度。具有相同深度的相邻位置将合并在一起形成非包容范围（请参阅示例输出）。

only-depth可以运行两种不同的模式，由--fast-mode标志控制。在快速模式下运行时，覆盖区域内的深度仅由读的起始和结束位置确定，不考虑任何CIGAR相关信息。--mate-fix仍可以在此模式下使用，且配对重叠的区域不会重复计数。

没有--fast-mode标志时，每个位置的深度以类似于base-depth的方式确定，其中DEL会计入深度，但REF_SKIP不会。此外，任何失败的--exclude-flags的读将不计入深度。最后，可以应用--mate-fix以避免在配对可能重叠的区域重复计数。

关于配对修正，perbase将根据读中计数的区域进行“修正”。例如，如果您有一个从“chr1:0-1000”开始的读，CIGAR为“25M974N1M”，而配对在“chr1:45-70”处良好对齐，CIGAR为“25M”，则配对将计入“chr1:45-74”的深度。这与其他拒绝配对即使它重叠R1中不计入深度的区域的工具形成对比。

为了获得最快的输出，请使用only-depth --fast-mode。

注意，相邻位置可能不会合并，如果它们位于--chunksize边界处。如果这是一个问题，可以将--chunksize设置为所涉及的最大contig的大小。将来可能会修复此问题或提供后处理工具来修复它。对于大多数用例，这应该不会成为问题。此外，您可以将输出传递到merge-adjacent，这将修复它。例如：perbase only-depth -m file.bam | perbase merge-adjacent > out.tsv。

perbase only-depth --mate-fix --zero-base /test/test.bam的示例输出

REF     POS     END     DEPTH
chr2    0       4       1
chr2    4       9       2
chr2    9       12      3
chr2    12      14      2
chr2    14      17      3
chr2    17      19      4
chr2    19      23      5
chr2    23      34      4
chr2    34      39      3
chr2    39      49      1
chr2    49      54      2
chr2    54      64      3
chr2    64      74      4
chr2    74      79      3
chr2    79      84      2
chr2    84      89      1

如果需要类似BED的输出，可以设置--bed-format -z标志，这将写入基于0、无标题的TSV输出，其中第四列空，第五列包含深度。

用法

Calculate the only the depth at each base

USAGE:
    perbase only-depth [FLAGS] [OPTIONS] <reads>

FLAGS:
        --bed-format                
            Output BED-like output format with the depth in the 5th column. Note, `-z` can be used with this to change
            coordinates to 0-based to be more BED-like
    -Z, --bgzip                     
            Optionally bgzip the output

    -x, --fast-mode                 
            Calculate depth based only on read starts/stops, see docs for full details

    -h, --help                      
            Prints help information

    -k, --keep-zeros                
            Keep positions even if they have 0 depth

    -m, --mate-fix                  
            Fix overlapping mates counts, see docs for full details

    -n, --no-merge                  
            Skip merging adjacent bases that have the same depth

    -M, --skip-merging-intervals    
            Skip mergeing togther regions specified in the optional BED or BCF/VCF files.
            
            **NOTE** If this is set it could result in duplicate output entries for regions that overlap. **NOTE** This
            may cause issues with downstream tooling.
    -V, --version                   
            Prints version information

    -z, --zero-base                 
            Output positions as 0-based instead of 1-based


OPTIONS:
    -B, --bcf-file <bcf-file>
            A BCF/VCF file containing positions of interest. If specified, only bases from the given positions will be
            reported on. Note that it may be more efficient to calculate depth over regions if your positions are
            clustered tightly together
    -b, --bed-file <bed-file>
            A BED file containing regions of interest. If specified, only bases from the given regions will be reported
            on
    -C, --channel-size-modifier <channel-size-modifier>
            The fraction of a gigabyte to allocate per thread for message passing, can be greater than 1.0 [default:
            0.001]
    -c, --chunksize <chunksize>
            The ideal number of basepairs each worker receives. Total bp in memory at one time is (threads - 2) *
            chunksize [default: 1000000]
    -L, --compression-level <compression-level>
            The level to use for compressing output (specified by --bgzip) [default: 2]

    -T, --compression-threads <compression-threads>
            The number of threads to use for compressing output (specified by --bgzip) [default: 4]

    -F, --exclude-flags <exclude-flags>                    
            SAM flags to exclude, recommended 3848 [default: 0]

    -f, --include-flags <include-flags>                    
            SAM flags to include [default: 0]

    -q, --min-mapq <min-mapq>                              
            Minimum MAPQ for a read to count toward depth [default: 0]

    -o, --output <output>                                  
            Output path, defaults to stdout

    -r, --ref-fasta <ref-fasta>                            
            Indexed reference fasta, set if using CRAM

    -t, --threads <threads>                                
            The number of threads to use [default: 32]


ARGS:
    <reads>    
            Input indexed BAM/CRAM to analyze

合并相邻

合并相邻 是一个用于合并类似 BED 文件中重叠区域的实用程序。

它将接受一个包含四列且没有标题的文件，只要列的格式如下：

<contig>\t<start>\t<stop>\t<depth>\n

或者它可以接受具有类似标题的三列文件：

<REF|chrom>\t<POS|chromStart>\t<END|chromEnd>\t<DEPTH|COV>

END|chromEnd 列是可选的。

perbase-merge-adjacent 0.7.5-alpha.0
Seth Stadick <[email protected]>
Merge adjacent intervals that have the same depth. Input must be sorted like: `sort -k1,1 -k2,2n in.bed > in.sorted.bed`

Generally accepts any file with no header tha is <chrom>\t<start>\t<stop>\t<depth>. The <stop> is optional. See
documentation for explaination of headers that are accepted.

USAGE:
    perbase merge-adjacent [FLAGS] [OPTIONS] [in-file]

FLAGS:
    -Z, --bgzip        
            Optionally bgzip the output

    -h, --help         
            Prints help information

    -n, --no-header    
            Indicate if the input file does not have a header

    -V, --version      
            Prints version information


OPTIONS:
    -T, --compression-level <compression-level>
            The level to use for compressing output (specified by --bgzip) [default: 2]

    -T, --compression-threads <compression-threads>
            The number of threads to use for compressing output (specified by --bgzip) [default: 32]

    -o, --output <output>                              
            The output location, defaults to STDOUT


ARGS:
    <in-file>    
            Input bed-like file, defaults to STDIN

示例

perbase only-depth indexed.bam | perbase merge-adjacent > out.tsv

类似项目

• Sambamba 深度
• Samtools 深度
• Mosdepth
• Bam-readcount