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

UB3/informatics/bioinformatics/samtools

このページの最終更新日: 2024/02/03

  1. 概要: samtools とは
  2. Mac への samtools のインストール
  3. samtools の各種オプション
  4. samtools の使い方
  5. 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 から 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

  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.

コメント欄

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