BAMファイルの切り取り bedtools intersect ? samtools ?

bam, vcf, bed ファイルなどから特定の領域を切り出す場合、bedtoolsが便利です。

以前のブログで intersectBed を紹介しましたが、現在のバージョンでは bedtools intersect を使います。

staffblog.amelieff.jp

先日、ある遺伝子の領域だけを見ようと、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.31.0 documentation

SAMtools

timeコマンドで計算時間とメモリを調査 - アメリエフの技術ブログ

アメリエフでは、バイオデータ解析やそのシステム・インフラ環境の開発に興味のあるエンジニア・リーダー候補を募集しています。