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

このページの最終更新日: 2024/09/30

概要: samtools とは
Mac への samtools のインストール
samtools の各種オプション
samtools の使い方
samtools のエラー

概要: samtools とは

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

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 から 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 はエラーメッセージなしで配列を除いている・・ということか？

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

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

マップされた/されなかったリードの数を書き出す

リード数のカウントには -c オプションを使用する。

samtools view -F 0x04 -c sample.bam

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

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 にマップされたリードを書き出す

-b は bam output のオプション。-b なしで sam でやっても動きそうだが、実際には試していない。

contig 名は、index を作った fasta ファイルの配列名であり、アセンブリーがよくできているなら染色体名になっていることもある。

# "samtools view -b 対象bamファイル対象contig" の形になる。
samtools view -b all.bam contig1 > contig1.bam

samtools のエラー

[E::sam_hrecs_update_hashes] Duplicate entry "文字列" in sam header
samtools view: failed to add PG line to the header

このページにあるように、sam file のヘッダーに重複があるというエラーであることはわかる。

このエラーが出ると、samtools view による bam への変換、samtools sort なども動かず、八方塞がりになる。sam ファイルはテキストファイルなので、重複を探して削除することはできるが、そもそもなぜ重複が生じるのかよくわからない。

fastp 済みのあるペアエンドのデータを bowtie2 でマップしたときにこのエラーが出た。BWA でもマップしてみたが、結果は同じだった。

ペアエンドなのが原因かと思い、一つのファイルを seqkit で　10 個に分割し (これは処理時間を短くするため)、それをマップしてみた。

References

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

コメント欄

サーバー移転のため、コメント欄は一時閉鎖中です。サイドバーから「管理人への質問」へどうぞ。