macでインフォマティクス

macでインフォマティクス

HTS (NGS) 関連のインフォマティクス情報についてまとめています。

高効率なカバレッジ計算ツール BamToCov

2022/02/25 論文引用

 

 多くのゲノミクスアプリケーションでは、リファレンスのヌクレオチドカバレッジを計算したり、リファレンス領域に何本のリードがマッピングされているかをカウントしたりする必要がある。本発表では、BamToCovを紹介する。このツールは、メモリ効率の良いアルゴリズムに依存し、カスタムパイプラインに柔軟に統合できるように設計された、迅速かつ柔軟なカバレッジ計算のためのツールスイートである。このツール群は、ソートされたBAMファイルやCRAMファイルを処理し、様々なフィルタリングアプローチを用いてカバレッジ情報を抽出することができる。

 BamToCovツールは、既存のツールとは異なり、最小限のメモリで、ワークフローに容易に統合でき、ストランドに特化したカバレッジ解析ができるように開発されている。独自のカバレッジ計算アルゴリズムにより、ロングリードのアラインメント解析に最適になっている。プログラムとそのドキュメントは、https://github.com/telatin/bamtocov で自由に利用することができる。

 

 アライメントファイル(BAM形式)からカバレッジ情報を抽出するツールは、すでにいくつか存在する。Samtools (Li et al., 2009), Bedtools (Quinlan, 2014), Sambamba (Tarasov et al., 2015) 。そして新しく、より機能豊富なMosdepth (Pedersen and Quinlan, 2018b) とMegaDepth (Wilks et al., 2021)がある。既存のツールの共通の限界は、mate-pairs ライブラリを使用してアセンブリの完全性を決定する際に重要な物理的カバレッジを計算できないことである。また、鎖ごとのカバレッジを分離することができない。ある位置がフォワードリードのみ、あるいはリバースリードのみによってカバーされている場合、それはおそらくミスアライメントに起因する。これらの制限を解決するために、Covtobed (Birolo and Telatin, 2020)を開発した。これは、コンピュータプログラミングのUNIX哲学に触発され、入出力ストリームをサポートする単一タスクに焦点を当てたシンプルかつ効率的なC++プログラムである。ここでは、Nim言語で記述されたBamToCovプログラムとその補助ユーティリティを紹介する。これは、入力ストリームを読み込む機能を維持しながら、インターバルターゲット、新しい出力フォーマット、カバレッジ統計、複数のBAMファイルをサポートする新しい機能を備えたCovtobedのコアアルゴリズムを用いてカバレッジ計算を行い、全体的にパフォーマンスの向上(すなわち、より小さなメモリフットプリントと最大3倍の速度向上)を達成している。

 

Documentation

https://telatin.github.io/bamtocov/

Wig format

https://telatin.github.io/bamtocov/notes/wig.html

 

 

ペアエンドライブラリを使用する場合、物理カバレッジも計算することができる。物理カバレッジとは、ロングインサートのペアエンド(メイトペア)で-->... <--のようにリードはカバーしていない領域(...)が発生するが、この領域もカウント対象としたカバレッジの事(プレプリント図1))

 

特徴

  • UNIX哲学の入力ストリームをサポートに対応しており、bamのインデックスは必要ない
  • ストランドバイアスをチェックするために、ストランドごとのカバレッジを計算可能
  • 少ないメモリ使用量
  • CRAMファイルをネイティブにサポート
  • ロングリードのアラインメントにも対応(Table.1)
  • 高速(MegaDepthに次ぐ速度)

 

インストール

Github

mamba create -n bamtocov
conda activate bamtocov
mamba install -y -c bioconda bamtocov

> bamtocov -h

BamToCov 2.3.0

 

  Usage: bamtocov [options] [<BAM>]...

 

Arguments:                                                                                                                                                 

  <BAM>         the alignment file for which to calculate depth (default: STDIN)

 

Core options:

  -p, --physical               Calculate physical coverage

  -s, --stranded               Report coverage separate by strand

  -q, --quantize <breaks>      Comma separated list of breaks for quantized output

  -w, --wig <SPAN>             Output in WIG format (using fixed <SPAN>), 0 will print in BED format [default: 0]

  --op <func>                  How to summarize coverage for each WIG span (mean/min/max) [default: max]

  -o, --report <TXT>           Output coverage report

  --skip-output                Do not output per-base coverage

  --report-low <min>           Report coverage for bases with coverage < min [default: 0]

 

Target files:

  -r, --regions <bed>          Target file in BED or GFF3/GTF format (detected with the extension)

  -t, --gff-type <feat>        GFF feature type to parse [default: CDS]

  -i, --gff-id <ID>            GFF identifier [default: ID]

  --gff-separator <sep>        GFF attributes separator [default: ;]

  --gff                        Force GFF input (otherwise assumed by extension .gff)

 

BAM reading options:

  -T, --threads <threads>      BAM decompression threads [default: 0]

  -F, --flag <FLAG>            Exclude reads with any of the bits in FLAG set [default: 1796]

  -Q, --mapq <mapq>            Mapping quality threshold [default: 0]

 

Other options:

  --debug                      Enable diagnostics

  -h, --help                   Show help

> bamtocounts -h

$ bamtocounts -h

BamToCounts 2.3.0

 

  Usage: bamtocounts [options] <Target> <BAM-or-CRAM>...

 

Arguments:                                                                                                                                                 

 

  <Target>       the BED (or GFF) file containing regions in which to count reads

  <BAM-or-CRAM>  the alignment file for which to calculate depth

 

Options:

 

  -T, --threads <threads>      BAM decompression threads [default: 0]

  -r, --fasta <fasta>          FASTA file for use with CRAM files [default: ].

  -F, --flag <FLAG>            Exclude reads with any of the bits in FLAG set [default: 1796]

  -Q, --mapq <mapq>            Mapping quality threshold [default: 0]

  -g, --gff                    Force GFF for input (otherwise autodetected by .gff extension)

  -t, --type <feat>            GFF feature type to parse [default: CDS]

  -i, --id <ID>                GFF identifier [default: ID]

  -n, --rpkm                   Add a RPKM column

  -l, --norm-len               Add a counts/length column (after RPKM when both used)

  --header                     Print header

  --debug                      Enable diagnostics    

  -h, --help                   Show help

>  covtotarget -h

covToTarget 2.3.0

 

  Usage: covtotarget [options] <Target> [<covtobed-output>]

 

Arguments:                                                                                                                                                 

 

  <Target>           the BED (or GFF) file containing regions in which to count reads

  <covtobed-output>  covtobed output, or STDIN if not provided

 

Options:

 

  -g, --gff                    Force GFF for input (otherwise autodetected by .gff extension)

  -t, --type <feat>            GFF feature type to parse [default: CDS]

  -i, --id <ID>                GFF identifier [default: ID]

  -l, --norm-len               Normalize by gene length

  -b, --bed-output             Output format is BED-like (default is feature_name [tab] counts)

  -h, --help                   Show help

 

 

実行方法

bamtocov - BAMファイルを解析してBED形式のカバレッジファイルを出力

bamを指定する。

インデックスがなくてもソートされたBAMファイルを読み込むことができる。
bamtocov input.bam > coverage.bed

#物理カバレッジ
bamtocov -p input.bam > coverage.bed
  •  -p   Calculate physical coverage

wigファイルフォーマットで出力する。

bamtocov --wig 200 input.bam > coverage.wig
  • -w   Output in WIG format (using fixed <SPAN>), 0 will print in BED format [default: 0]

strandedを指定すると、forwardとreverseそれぞれのカバレッジが計算される。

bamtocov --stranded input.bam > coverage.bed
  • -s    Report coverage separate by strand

5列のBEDライクなファイルが出力される。4列目がforward strand coverage、5列目がreverse strand coverage。

 

BamToCounts - BAMファイル中のターゲットのリード数をカウント

ターゲット領域のbedを指定する。

BamToCounts target.bed input.bam  > coverage.txt

 

BamCountsRefs - 複数のBAMファイル(同じ参照配列を持つ)からカウントテーブルを出力

wigファイル出力

bamcountrefs --tag "Chr1" input1.bam input2.bam   

 

covToTarget - ターゲット(BED または GFF3 フォーマットのアノテーションファイル)と covtobed 1.0 の出力を基に、フィーチャーごとのカバレッジレポートを作成

covtobed input/mini.bam | covtotarget input/mini.bed > output/counts.tsv

 

シミュレートされたbamを生成するコマンドも用意されています。ドキュメントを確認して下さい。

引用

BamToCov: an efficient toolkit for sequence coverage calculations
Giovanni Birolo,  Andrea Telatin

bioRxiv, Posted November 17, 2021

 

2022/02/24

BamToCov, an efficient toolkit for sequence coverage calculations 
Giovanni Birolo,  Andrea Telatin  Author Notes
Bioinformatics, Published: 23 February 2022