bam, vcf, bed ファイルなどから特定の領域を切り出す場合、bedtoolsが便利です。
以前のブログで intersectBed を紹介しましたが、現在のバージョンでは bedtools intersect を使います。
先日、ある遺伝子の領域だけを見ようと、bamファイルを切り取っていたのですが
40GBもあるbamだけあって、なかなか終わらない...（最終的に、25分かかった）
のですが、bamを扱うソフトウェア、 samtools を使うと、高速に切り取ることができます！
※ 特定の少数の領域だけを切り取る場合です。
多数の領域や、重なり方があれこれ...などはbedtoolsが多機能で便利です。
samtools view でBAMの領域切り取り
切り取りたいbamは、ソートされてインデックスファイルが作成されている必要があります。
# sort samtools sort -O BAM -o sorted.bam unsorted.bam # index samtools index sorted.bam
特定の領域の切り取りは、以下のコマンドで実行します。
# samtools view <option> <入力bam> <領域> samtools view -o result.bam -O BAM sorted.bam chr1:12345678-23456789
領域の指定方法は、染色体:開始-終了 ですが、
染色体番号のみでも、開始・終了位置にはカンマが入っていてもOKです。
領域は、スペース区切りで複数並べることもできます。(染色体順ではなく、コマンドで指定した順番に出力されます)
samtools view -o result.bam -O BAM sorted.bam chr1:12345678-23456789 chrX chr2:1,000,000-2,000,000
ただし、手打ちやコピー&ペーストではミスもしやすいので
領域が多数ある場合はおとなしくbedtools を使うといいと思います。
viewは速いよ
samtools vs bedtools で 時間、使用メモリを測ってみました。
上のコマンドを使うと、なんと！
# samtools view $ /usr/bin/time -f "Memory %M KB Time %E" \ > samtools view -o result.bam -O BAM sorted.bam chr1:8,888,000-8,888,888 Memory 17756 KB Time 0:00.16 # bedtools intersect $ /usr/bin/time -f "Memory %M KB Time %E" \ > bedtools intersect -a sorted.bam -b test.bed > result.bam Memory 4492 KB Time 25:04.56
bedtools で25分かかったのに対し、samtools では1秒足らずで計算が終わりました。
ただしメモリ使用はbedtools の方が小さかったです。
更に！
samtools は view を含め多くのコマンドでマルチスレッドで計算することができます。
-@ ／ --threads で使用スレッド数を増やせばもっと高速で計算できます。
samtools view -o result.bam -O BAM -@ 3 sorted.bam chr1:8,888,000-8,888,888 # → 1 + 3で4スレッド使用
ぜひお役立てください。
浮いた時間でブログが書けますよ♪
追記：samtoolsで切り取る領域
お客様よりご質問いただきました。
samtools の仕様で、指定領域にオーバーラップするリードが 出力bam に残されるため、
"指定領域にマップ+少しはみ出る部分もある"リードがあった場合に、指定領域よりも長い領域で抽出されます。
余談。 samtools の隠しコマンド?
指定する領域の書式（chr1:12345678 のような）を調べるため
samtools view のヘルプメッセージを見てみると、[region] と書いてあるだけです。
$ samtools view
Usage: samtools view [options] <in.bam>|<in.sam>|<in.cram> [region ...]
Options:
-b output BAM
-C output CRAM (requires -T)
-1 use fast BAM compression (implies -b)
-u uncompressed BAM output (implies -b)
-h include header in SAM output
-H print SAM header only (no alignments)
-c print only the count of matching records
-o FILE output file name [stdout]
-U FILE output reads not selected by filters to FILE [null]
-t FILE FILE listing reference names and lengths (see long help) [null]
-L FILE only include reads overlapping this BED FILE [null]
-r STR only include reads in read group STR [null]
-R FILE only include reads with read group listed in FILE [null]
-q INT only include reads with mapping quality >= INT [0]
-l STR only include reads in library STR [null]
-m INT only include reads with number of CIGAR operations consuming
query sequence >= INT [0]
-f INT only include reads with all of the FLAGs in INT present [0]
-F INT only include reads with none of the FLAGS in INT present [0]
-G INT only EXCLUDE reads with all of the FLAGs in INT present [0]
-s FLOAT subsample reads (given INT.FRAC option value, 0.FRAC is the
fraction of templates/read pairs to keep; INT part sets seed)
-M use the multi-region iterator (increases the speed, removes
duplicates and outputs the reads as they are ordered in the file)
-x STR read tag to strip (repeatable) [null]
-B collapse the backward CIGAR operation
-? print long help, including note about region specification
-S ignored (input format is auto-detected)
--input-fmt-option OPT[=VAL]
Specify a single input file format option in the form
of OPTION or OPTION=VALUE
-O, --output-fmt FORMAT[,OPT[=VAL]]...
Specify output format (SAM, BAM, CRAM)
--output-fmt-option OPT[=VAL]
Specify a single output file format option in the form
of OPTION or OPTION=VALUE
-T, --reference FILE
Reference sequence FASTA FILE [null]
-@, --threads INT
Number of additional threads to use [0]
はっ( ﾟ ω ﾟ ) ! !
" -? print long help, including note about region specification "
領域指定についての説明など、長いヘルプを見るには -? をつける
隠れ機能がありました ¯\_(ツ)_/¯
$ samtools view -?
...
...
Notes:
1. This command now auto-detects the input format (BAM/CRAM/SAM).
Further control over the CRAM format can be specified by using the
--output-fmt-option, e.g. to specify the number of sequences per slice
and to use avoid reference based compression:
samtools view -C --output-fmt-option seqs_per_slice=5000 \
--output-fmt-option no_ref -o out.cram in.bam
Options can also be specified as a comma separated list within the
--output-fmt value too. For example this is equivalent to the above
samtools view --output-fmt cram,seqs_per_slice=5000,no_ref \
-o out.cram in.bam
2. The file supplied with `-t' is SPACE/TAB delimited with the first
two fields of each line consisting of the reference name and the
corresponding sequence length. The `.fai' file generated by
`samtools faidx' is suitable for use as this file. This may be an
empty file if reads are unaligned.
3. SAM->BAM conversion: samtools view -bT ref.fa in.sam.gz
4. BAM->SAM conversion: samtools view -h in.bam
5. A region should be presented in one of the following formats:
`chr1', `chr2:1,000' and `chr3:1000-2,000'. When a region is
specified, the input alignment file must be a sorted and indexed
alignment (BAM/CRAM) file.
6. Option `-u' is preferred over `-b' when the output is piped to
another samtools command.
参考
bedtools: a powerful toolset for genome arithmetic — bedtools 2.30.0 documentation
timeコマンドで計算時間とメモリを調査 - アメリエフの技術ブログ
アメリエフでは、バイオデータ解析やそのシステム・インフラ環境の開発に興味のあるエンジニア・リーダー候補を募集しています。
www.wantedly.com
