samtools: インストール、使い方、オプションなど

UB3/informatics/bioinformatics/samtools

このページの最終更新日: 2022/05/18

  1. 概要: samtools とは
  2. Mac への samtools のインストール
  3. samtools の各種オプション
  4. samtools の使い方

広告

概要: samtools とは

samtools とは、マッピングの結果として得られる sam ファイルや bam ファイルをさらに解析するためのソフトである。ターミナル上で動くので、このページに書いてあることがわかりにくい場合は、まず以下のリンク先を読むことをお勧めする。

  1. ターミナルの使い方
  2. シェルスクリプト
  3. 次世代シークエンス
  4. マッピング: bowtie2, bwa

Mac への samtools のインストール

Mac へのインストールは簡単である。まず homebrew のページを参考にして homebrew をインストールし、brew install samtools で OK。X code のコマンドラインアプリケーションが必要。Mac に入っていない場合は、実行するコマンドを示したエラーメッセージが出るので、それに従う。

インストールログ

2020 年 1 月に Catalina にアップデート。X-code を再インストールし、さらに homebrew による samtools の再インストールが必要だったが、ちゃんと動いている。

さらに 2019 年 11 月に別の Mac でセットアップ。bowtie2, samtools ともに homebrew を使ってインストール可能だった。OS は Mojave。

2016 年ごろ、Mac OSX High Sierra, homebrew でインストール。


広告

samtools の各種オプション

samtools は、view, sort などのメインとなるオプションと、-h, -F などの - を含むオプションの両方を指定して使う (「オプション」という言葉が正確なのかどうかわからない)。

つまり、コマンドは以下のような形になる。

samtools view -h -b -F4 sample.bam > sample.mapped.sam


view

bam または sam ファイルを見るオプションであるが、実際は上記のように書き出しに使われることが多い。マップされたリードのみ、マップされていないリードのみを検出する -F、-f オプションは頻繁に使われる。

sort

bam または sam ファイルをソートする。samtools sort a.bam -o a_sorted.bam のようになり、-o オプションでアウトプットファイルを指定する。> a_sorted.bam でアウトプットを指定している資料もあるが、私の場合はこれでは動かなかった。

index

samtools index a.bam のようにして bam ファイルにインデックスをつける。igv ブラウザによる視覚化の前に必要である。成功すると .bam.bai という拡張子のファイルができる。

merge

複数の bam files をマージするのに使う。bam はバイナリファイルなので、cat ではマージできない。

samtools merge out.bam in-1.bam in-2.bam で in-1.bam および in-2.bam という 2 つのファイルがマージされる。この部分は *.bam でも大丈夫。

fastq

samtools fastq input.bam > output.fastq のようにして、bam から fastq ファイルを作るのに使う。詳細はページ下方の bam から fastq への変換 を参照のこと。

flagstat

samtools flagstat input.sam で、マッピング情報を表示する。fasta, fastq ファイルの情報を得る seqkit stat に似ている。

表示される情報は、例えば

  • マップされたリードの数、割合。

stat

samtools stat input.bam で、マッピング情報を表示する。SN がサマリーになる。

depth

samtools depth input.bam で、マッピングの depth を表示する。


-h

ヘッダーを維持。

-F, -f

-F でマップされた配列のみを、-f でマップされていない配列のみを扱う。view とともに書き出しで使うことが多いだろう。

-c

bam または sam ファイルのリード数をカウント。

  • samtools view -c sample.bam なら全てのリード
  • samtools view -F 0x04 -c sample.bam ならマップされたリードのみ

-b

bam ファイルでアウトプットする。

-S

input に sam ファイルを使う。ただし、新しいバージョンではこのオプションをつけなくても自動でフォーマットが認識される。

samtools の使い方

要は、上記のオプションをどう組み合わせ、どう使いこなすかである。


sam と bam の変換

sam file から bam file への変換も samtools で行うことができる。コマンドは samtools view、-S は input の sam file を指定するオプション、-b は output に bam 形式を指定するオプションで、-S -b と書いても良いが、まとめて -bS とも書ける。書き出しには > が必要 (1)。

samtools view -bS input.sam > output.bam

逆に、bam から sam に変換するときには、bam output の -b オプションは必要なく、さらに output を sam にするオプションはなくてよいので、以下のようになる。ヘッダーは維持しておいた方が良いと思うので、そのオプションのみ追加している。。

samtools view -h input.bam > output.bam



マップされた/されないリードのみを書き出す

bam file からマップされたリードのみ、またはマップされなかったリードのみを抽出するには、以下のようにする。-F4 がマップされた配列、-f4 がマップされなかった配列 (2)。

単にこれをやると、bam file のヘッダーが失われてしまい、その後の解析で missing SAM header などのエラーが出て困ることが多い。これを防ぐために、ヘッダーを維持する -h オプション を入れておくとよい。

bam から bam でも、アウトプットオプションに bam を指定する -b を入れておかないと、.bam という拡張子をもつ sam ファイルができる気がする。ファイルサイズ、text かバイナリかで見分けがつく。

samtools view -b -h -F4 sample.bam > sample.mapped.bam
samtools view -b -h -f4 sample.bam > sample.unmapped.bam

実はこれは sam file からの抽出も可能で、-S オプションと組み合わせれば良い。

samtools view -S -F4 sample.sam > sample.mapped.sam
samtools view -S -f4 sample.sam > sample.unmapped.sam

-F 4 のようにスペースを空けるように書いてあるページも多い。たぶんスペースありが正式だが、なしでも同じ結果になる。



マップされたリードの位置を書き出す

samtools view で bam file の内容を表示し、これを awk にパイプするという形をとる (8)。

samtools view inpu.bam | awk '{print $3 "\t" $4 "t\" $4+length($10)-1}' > output.txt

スペースで区切られた 3 列目はマップされた配列の名前、4 列目 ($4) はリードの始点。その次は、始点にリードの長さ (length 関数で 10 行目) を加えて 1 を引くことで、リードの終点を表している。


特定の contig にマップされたリードを書き出す

samtools view



bam から fastq への変換

bam ファイルを fastq に変換し、配列を見たり、さらに解析を続けたりしたい状況もあるだろう。samtools を使う場合は、fastq を使って以下のようにする。

samtools fastq input.bam > output.fastq

経験的には、20 GB ぐらいの bam は 80 GB ぐらいの fastq になる。samtools で行うと、この書き出しは 10 分程度。

ただしこの場合、最初の fastq が paired-end で 2 つのファイルがある場合でも、bam および変換された fastq は一つのファイルになり、paired-end の情報が失われてしまう気がする。

他の変換方法もあり (7)、paired-end の場合には bedtools または picard を使うのが良さそうだ。どちらを使っても 20 GB ぐらいの bam は 40 GB x 2 の fastq になる。1.5 時間ぐらい。

bedtools

brew install bedtools でインストールできる。single read の場合は 1 行目、paired read の場合は 2 行目。

bedtools bamtofastq -i input.bam -fq single.fastq

bedtools bamtofastq -i R1R2.bam -fq R1.fastq -fq2 R2.fastq

picard

Java ベースなのか? 変換コマンドは以下。input は bam でも SamToFastq を使う。

picard SamToFastq INPUT=input.bam F=R1.fastq F2=R2.fastq


十分に検討したわけではないが、bedtoold を使ったとき ** is marked as paired, but its mate does not occur next to it in your BAM file. Skipping. というエラーが出た。同じファイルで picard を変換するとエラーなし。しかし、生成された fastq ファイルを seqkit stat すると、含まれる配列数は完全に同一である。つまり、picard はエラーメッセージなしで配列を除いている・・ということか?


広告

References

  1. Converting SAM to BAM with samtools “view”. Link: Last access 2019/11/11.
  2. Separating mapped and unmapped reads from libraries. Link: Last access 2019/11/11.
  3. Decoding sam flags. Link: Last access 2019/11/20.
  4. SAM file format. Link: Last access 2019/11/20. SAM file の見方が載っている便利なページ。
  5. SAMtools ワンライナー覚書. Palmsonntagmorgen. Link: Last access 2020/02/02.
  6. sam/bam 関係のツールまとめ. Link: Last access 2020/12/07.
  7. フォーマット変換 bam=> Fastq アライメントされなかったリードの取り出し方など. Link: Last access 2020/12/14.
  8. Bam File: Extract Chromosome Number And Start Position Reads. Link: Last access 2021/03/08.

コメント欄

各ページのコメント欄を復活させました。スパム対策のため、以下の禁止ワードが含まれるコメントは表示されないように設定しています。レイアウトなどは引き続き改善していきます。「管理人への質問」「フォーラム」へのバナーも引き続きご利用下さい。

禁止ワード: http, the, м (ロシア語のフォントです)


このページにコメント

Name:


Comment:



これまでに投稿されたコメント

Date Name Comment