2020 4/19 help更新, コマンド追記
2021 5/14 インストール追記, help更新
カバレッジデプスの測定は、コピー数変動（CNV）の検出などのゲノム解析にとってクリティカルで、例えばcn.mops（Klambauer et al、2012）、Indexcov（Pedersen et al、2017）などのどのゲノム領域のカバレッジが低すぎるか高すぎるか（Li、2014）調べる方法がある。カバレッジプロファイルの適用範囲を考えると、ゲノム全体のカバレッジを計算する既存のツールがいくつかある。 Samtools depth（Li et al、2009）は、ベースごとのカバレッジを出力する。 BEDTools genomecov（Quinlan and Hall、2010; Quinlan、2014）は、リージョンごとまたはベースごとのカバレッジを出力する。 Sambamba（Tarasov et al、2015）（紹介）はまた、ベースごとおよびウィンドウごとのデプスを提供する。効率的なカバレッジ計算の必要性は、全ゲノム配列の数およびシーケンスデプスと共に増加し、既存の方法は、30倍のカバレッジを有する典型的なヒトゲノムについて約1時間またはそれ以上の計算を必要とする。ここでは、既存のメソッドより高速で、追加のユーティリティーを持ったmosdepthを示す。
Mosdepthは、nimプログラミング言語（https://nim-lang.org）を経由しHTSLib（http://www.htslib.org/）を使用している。入力したBAMまたはCRAMファイルがポジションによってソートされることが期待される。 mosdepthは、リードごとに各ヌクレオチドを追跡する「パイルアップ」エンジンを使用するsamtoolsとは対照的に、アライメントのチャンク（塊）のみを追跡する。アラインメントの各チャンクの開始位置と終了位置（削除またはその他のイベントによって分割された場合、各アライメントは複数のチャンクを持つ可能性がある）は、サイズが染色体の長さである配列（32ビット整数）で追跡される。リファレンスゲノムへのアラインメントの各チャンクについて、mosdepthはその染色体位置（論文 図1）に対応する配列内のインデックスの値の開始とインクリメントをインクリメントする。ペアエンドシーケンシングフラグメントの両端に重複するアライメントがある場合（論文図1、ブラックアライメント）、ダブルカバレッジカバレッジを回避する。カバレッジアレイがBAMまたはCRAMファイル内のすべてのアラインメントの開始と終了を追跡すると、特定の位置の深さは、それに先行するすべての位置の累積合計として計算される（同様のアルゴリズムがBEDToolsで使用されている）。
染色体に沿って、カバレッジは、配列の各要素までの複合開始および終了カウントを累積合計で置き換えることによって、適切な位置で計算される。完了すると、領域のカバレッジは、配列の要素の最初から最後までの平均値になる。これにより、数百万の小領域であっても、極めて迅速にカバレッジを計算することができる。このセットアップは、ゲノムのカバレッジ分布の迅速な計算、すなわち、ゲノムまたは所定の領域にわたる所定数のリードによってカバーされる塩基の数にも影響を受けやすい。分散計算には、各カバレッジ値の発生をカウントする配列を介した余分な繰り返しが必要となる。 mosdepth法はヒトゲノムの249メガベースのchr1に対してより多くのメモリを必要とし、約1GBのメモリが必要だが、その数はカバレッジのデプスやアラインメントされたリード数に依存しない。その柔軟性にもかかわらず、mosdepthは使いやすく、理解しやすくなっている（例えば、論文のsupplementary documentを参照）。
インストール
cent os6でテストした。後日、ubuntu18.04LTSにcondaを使って導入した（pyhthon3.7）。
依存
- htslib version 1.4 or later
- PCRE
PCREが古いと動作しない。その場合はアップグレードする。
cd ~/src/
wget ftp://ftp.csx.cam.ac.uk/pub/software/programming/pcre/pcre-8.41.tar.gz
tar zxvf pcre-8.41.tar.gz
cd pcre-8.41/
./configure
make
#bioconda (link) (linux)
mamba create -n mosdepth -y
conda activate mosdepth
mamba install -c bioconda -y mosdepth
#docker
docker pull quay.io/biocontainers/mosdepth:0.2.4--he527e40_0
> mosdepth -h
# mosdepth -h
mosdepth 0.3.1
Usage: mosdepth [options] <prefix> <BAM-or-CRAM>
Arguments:
<prefix> outputs: `{prefix}.mosdepth.dist.txt`
`{prefix}.mosdepth.summary.txt`
`{prefix}.per-base.bed.gz` (unless -n/--no-per-base is specified)
`{prefix}.regions.bed.gz` (if --by is specified)
`{prefix}.quantized.bed.gz` (if --quantize is specified)
`{prefix}.thresholds.bed.gz` (if --thresholds is specified)
<BAM-or-CRAM> the alignment file for which to calculate depth.
Common Options:
-t --threads <threads> number of BAM decompression threads [default: 0]
-c --chrom <chrom> chromosome to restrict depth calculation.
-b --by <bed|window> optional BED file or (integer) window-sizes.
-n --no-per-base dont output per-base depth. skipping this output will speed execution
substantially. prefer quantized or thresholded values if possible.
-f --fasta <fasta> fasta file for use with CRAM files [default: ].
Other options:
-F --flag <FLAG> exclude reads with any of the bits in FLAG set [default: 1796]
-i --include-flag <FLAG> only include reads with any of the bits in FLAG set. default is unset. [default: 0]
-x --fast-mode dont look at internal cigar operations or correct mate overlaps (recommended for most use-cases).
-q --quantize <segments> write quantized output see docs for description.
-Q --mapq <mapq> mapping quality threshold. reads with a quality less than this value are ignored [default: 0]
-T --thresholds <thresholds> for each interval in --by, write number of bases covered by at
least threshold bases. Specify multiple integer values separated
by ','.
-m --use-median output median of each region (in --by) instead of mean.
-R --read-groups <string> only calculate depth for these comma-separated read groups IDs.
-h --help show help
ラン
キャプチャされた領域のbedファイルはイルミナからダウンロードした（本来あり得ないが、このテストではキットの種類は気にしない）。
exomeデータ（ERR035486）を使い、ヒトゲノムアセンブリhg19に対してマッピングする。minimap2を使う（mminimap2紹介）。
minimap2 -t 10 -ax sr ref.fa pair1.fq pair2.fq | samtools sort -@ 10 -O BAM -o exome.bam - && samtools index exome.bam
mosdepthを使ってchr1のデプスを計算する。
mosdepth --by capture.bed output exome.bam -Q 0 -t 10 -c chr1
- -Q mapping quality threshold [default: 0]
- -t number of BAM decompression threads [default: 0]
- -c chromosome to restrict depth calculation.
ポジションごとのカバレッジと、bedの領域ごとのカバレッジがgz出力される。-nをつけると ポジションごとのカバレッジファイルは出力されなくなる。
ファイル内容を確認。
> head output.regions.bed
$ head output.regions.bed
chr1 13402 13639 CEX-chr1-13403-13639 149.31
chr1 69088 70010 CEX-chr1-69089-70010 824.77
chr1 324438 325605 CEX-chr1-324439-325605 0.79
chr1 664484 665108 CEX-chr1-664485-665108 1.34
chr1 721404 721912 CEX-chr1-721405-721912 0.21
chr1 762079 762571 CEX-chr1-762080-762571 2.47
chr1 861319 861395 CEX-chr1-861320-861395 15.22
chr1 865532 865718 CEX-chr1-865533-865718 7.06
chr1 866416 866471 CEX-chr1-866417-866471 46.58
chr1 871149 871278 CEX-chr1-871150-871278 29.43
> head output.per-base.bed
$ head output.per-base.bed
chr1 0 10004 0
chr1 10004 10007 1
chr1 10007 10008 0
chr1 10008 10029 1
chr1 10029 10051 2
chr1 10051 10052 1
chr1 10052 10054 2
chr1 10054 10055 3
chr1 10055 10061 4
chr1 10061 10073 5
output.mosdepth.global.dist.txtには、カバレッジの度数分布が出力される。
追記
multiqcで統合レポートを作成する。全スレッド使用。MQ≥1。
mosdepth --by gencode.v33.annotation.bed sample1.mosdepth sample1.bam -Q 1
mosdepth --by gencode.v33.annotation.bed sample2.mosdepth sample2.bam -Q 1
mosdepth --by gencode.v33.annotation.bed sample3.mosdepth sample3.bam -Q 1
multiqc .
引用
Mosdepth: quick coverage calculation for genomes and exomes.
Pedersen BS, Quinlan AR
Bioinformatics. 2018 Mar 1;34(5):867-868.
