macでインフォマティクス

macでインフォマティクス

NGS関連のインフォマティクス情報についてまとめています。

sam/bamファイルを変換、編集したり分析するためのツール

 

samとbamのハンドリングに関するツールを紹介する。

似たような機能と名前を持つツールが多くて混乱するが、データ処理時に使えると便利なものもある。代表的な機能に限って紹介するので、確認してみて下さい。

 

brewでツールをインストールするので、はじめにbrew tap でsiience系のコマンドをオフィシャルコマンドとしてインストールできるようにしておく。

brew tap homebrew/science

 

9/1追記  samblaster。

12/16 追記

 

picard-tools

brewでインストールできる。

brew install picard-tools

BamIndexStats リファレンスにアラインされたリードの数。

picard BamIndexStats I=output_sorted.bam > bam_index_stats.tsv

I=でbamを指定。

>のリダイレクトで出力ファイル名を指定。

クロモソームごとのアライメントされたリードの数を出力できる。

Lambda_NEB length= 48502 Aligned= 4613 Unaligned= 0

 

 CollectAlignmentSummaryMetrics アライメントされたリードの平均サイズ。

picard CollectAlignmentSummaryMetrics I=R1R2_sorted.bam OUTPUT=aln_sum_metric.tsv

 

CollectInsertSizeMetrics インサートサイズの分析。

Picard CollectInsertSizeMetrics INPUT=R1R2_sorted.bam OUTPUT=insert_size_metric.tsv H=insert_size_metric.pdf

Rで描画された図が出力される。

 

f:id:kazumaxneo:20170723191440j:plain

QualityScoreDistribution クオリティスコアの分析。

Picard QualityScoreDistribution  INPUT=R1R2_sorted.bam OUTPUT=quality_score_dist.tsv CHART_OUTPUT=quality_score_dist.pdf

Rで描画された図が出力される。

f:id:kazumaxneo:20170723191557j:plain

gatkのBQSR補正したbamを使った場合、補正前のスコアがOriginal quality scoreで表示され、補正前後のクオリティの変化を比較できる。

 

SamToFastq fastqの抽出

Picard SamToFastq  INPUT=R1R2_sorted.bam F=R1.fq F2=R2.fq

ペアリードをFとF2で指定する。シングルならF=で指定。 ペアリードの片方だけアンマップのペアや、セカンドアライメントのペア、ペアでないリードだけfastqに返すようなオプションもある。https://gsl.hudsonalpha.org/information/software/bam2fastqにはunmapのリードを返せるとあるが、そのようなオプションはないように見える。

 

samtools

アライメントレポートを出す。

samtools flagstat R1R2_sorted.bam

出力結果の例

3153604 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
18472 + 0 mapped (0.59% : N/A)
3153604 + 0 paired in sequencing
1576802 + 0 read1
1576802 + 0 read2
5708 + 0 properly paired (0.18% : N/A)
16576 + 0 with itself and mate mapped
1896 + 0 singletons (0.06% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

properly pairedは、同じクロモソームに適切な向きと位置でアライメントされたペアリードの数になる。詳細に興味がある人はsamtoolsのマニュアルを参考にしてください。

http://www.htslib.org/doc/samtools.html

レポート結果が正しいかどうかは簡単なカウントで判断できる。

samtools view R1R2_sorted.bam |wc #1行目と同じ数値になる。

   

アライメントされたリードのデプスを出す。

samtools depth R1R2_sorted.bam

sam/bamのマージ。

samtools merge 

2つ以上のbam、samを1つにまとめる。splitコマンドでクロモソームごとに分割し、別のクラスターでランして、後で統合するような用途にも使えそうである。

 

sam/bamの分割。

samtools split 

 

bam/samからfastqを抽出。

samtools fastq input.bam > input.fq

 

bam/samからfasta抽出。

samtools fasta input.bam > input.fa

 

bamのヘッダーを変更する。

samtools reheader header.sam input.bam > output.bam

input.bamのheader.samのヘッダー情報で置換する。bam=>sam=>bamと変換するより圧倒的に高速である。例えば以下のようにして使う。

samtools view -H input.bam > header.sam #header情報だけsamに変換
# viなどでheader.bamを編集。例えばSO:unsortedをSO:sortedに修正。
samtools reheader header.sam input.bam > reheaer.bam #元のbamのヘッダーを編集したheaderで置換。

samtools view -H reheaer.bam |head #headerが交換されたことを確認。

 

 

PCR duplicateを除く。

samtools rmdup input.bam output.bam

ペアリードのみに使える。ペアリードがFRの向きで同じクロモソームの同じ位置にアライメントされているリードが2つ以上あると適用される。duplicate判定されるとマッピングクオリティがベストのものだけ残し、あとのペアは排除される。ペアの片方しかアライメントされていないものには適用されない。duplicateがあるかどうかは、fastqcでfastqを分析すれば調べられる。

 

リードグループ情報を追加したり、消したりする。

samtools addreplacerg 

 

heterozygous SNPsのコール。

samtools phase 

 

リードをポジション順にソート。

samtools sort -@ 12 -m 8000000000 R1R2.bam -o R1R2_sorted.bam
  • -@ Set number of sorting and compression threads [1]
  • -m Set maximum memory per thread; suffix K/M/G recognized [768M]
  • -n Sort by read name

デフォルトではメモリーを768Mまでしか使わない。メモリーがもっとあるなら、-mで指定することで、特に大きなゲノムのソートが速くなる。上では8GB指定して、スレッド数を12にしている。入力はbamとsamに対応している。samを指定すると、ソートしたbamを出力してくれる。

 

bam <=> sam変換やフィルタリング。

samtools view -bS input.sam > input.bam #samからbamに変換
samtools view -h input.bam > input.sam #bamからsamに変換
samtools view input.bam |less #bamの中身を見る。headerは見れない。
samtools view -H input.bam |less #bamのheaderだけ見る。ソートされているか

samtools view -s 0.5 -b input.bam > random0.5_input.bam #50%ランダムサンプリング
samtools view -h -b -F 4 input.bam > mapped.bam #mマッピングされたリードだけbamとして取り出す。
samtools view -b input.bam chr1 > chr1.bam #chr1にマッピングされたリードだけ取り出す
samtools view -b input.bam chr1:100-300 > chr1_100-300.bam #chr1の100-300にマッピングされたリードだけ取り出す。
  • -b output BAM
  • -C output CRAM requires -T)
  • -h include header in SAM output
  • -S ignored (input format is auto-detected)
  • -T Reference sequence FASTA FILE [null]
  • -s FLOAT integer part sets seed of random number generator [0]; rest sets fraction of templates to subsample [no subsampling]
  • -F INT   only include reads with none of the bits set in INT set in FLAG [0]

他にも以下のような使い方がある。

samtools view -s 0.5 -b input.bam > random0.5_input.bam #50%ランダムサンプリング
samtools view -h -b -F 4 input.bam > mapped.bam #mマッピングされたリードだけbamとして取り出す。samなら-bを消す。
samtools view -b input.bam chr1 > chr1.bam #chr1にマッピングされたリードだけ取り出す
samtools view -b input.bam chr1:100-300 > chr1_100-300.bam #chr1の100-300にマッピングされたリードだけ取り出す。

 

 

bamのindex。

samtools index sorted.bam
  • -b Generate BAI-format index for BAM files [default]
  • -c Generate CSI-format index for BAM files

 

bamtools 

bamファイルを使い色々なことができるツール。samtoolsの機能と重複するものが多いが、よりデータ解析よりに作られている。

brewでインストールできる。

brew install bamtools

 

BamTools convert bamtoolsの一番基本的な使い方はbamからのformat変換である。

bamtools convert -format <format> -in input.bam -out output

formatの部分に出力するフォーマットを記載する。pileup、bedのほかfastaやfastqに直接変換することも可能になっている(bedtoolsのbamtofastqと同じ)。ただしbamtoolsはシングルスレッドで動くので、samの出力はマルチスレッド化したsamtools viewが速いと思われる。

 

 BamTools stats bamの情報を出力できる。

bamtools stats -in input.bam

bamtools statsの出力

**********************************************

Stats for BAM file(s): 

**********************************************

 

Total reads:       525231

Mapped reads:      524170 (99.798%)

Forward strand:    262861 (50.0467%)

Reverse strand:    262370 (49.9533%)

Failed QC:         0 (0%)

Duplicates:        0 (0%)

Paired-end reads:  525231 (100%)

'Proper-pairs':    517998 (98.6229%)

Both pairs mapped: 523577 (99.6851%)

Read 1:            262603

Read 2:            262628

Singletons:        593 (0.112903%)

 

samtoolsのflagstatの出力

525231 + 0 in total (QC-passed reads + QC-failed reads)

0 + 0 secondary

1683 + 0 supplementary

0 + 0 duplicates

524170 + 0 mapped (99.80% : N/A)

523548 + 0 paired in sequencing

261774 + 0 read1

261774 + 0 read2

517042 + 0 properly paired (98.76% : N/A)

521908 + 0 with itself and mate mapped

579 + 0 singletons (0.11% : N/A)

846 + 0 with mate mapped to a different chr

791 + 0 with mate mapped to a different chr (mapQ>=5)

大きな差はないが、個人的にはbamtools statの方が分かりやすい表記だと思う。

 

BamTools count アライメントされたリード数のカウント。

bamtools count -in input.bam

 

BamTools sort bamをname順に並べ替える。

bamtools sort -in input.bam -out input_sort.bam -byname

 リード名順に並べ替えられる。samtools view -Hでヘッダーを確認すると

$ samtools view -H input_sort.bam |head 1

@HD VN:1.3 SO:queryname

になっており、名前順にソートされていることを確認できる。(-HはSAMheaderのみprintするsamtools viewのオプション)。samtools sort -nと同じと考えられる。

 

BamTools index bamのindexを作る。

bamtools index -in input_sort.bam

samtools indexと同じである。~.baiが出力される。

 

BamTools coverage アライメントされたリード数のカウント。

bamtools convert -in input.bam -out coverage.txt

version2.4.1ではなぜかsortしてないというエラーメッセージが出て検証できていない(bamtools sortとindexをかけていてもエラーになってしまう)。

 

 BamTools filter bamのフィルタリング。

bamtools filter -in input.bam -out output.bam -length 50

 様々なパラメータでbamをフィルタリングできる。上記では50bpのリードのみ出力される。ただし~以上や~以下、また複数指定もできない(正直、使い勝手はイマイチである)。以下のようなフィルターがある。

  • -alignmentFlag  keep reads with this *exact* alignment flag (for more detailed queries, see below)
  • -insertSize   keep reads with insert size that matches pattern -length <int> keep reads with length that matches pattern
  • -mapQuality   keep reads with map quality that matches pattern -name <string> keep reads with name that matches pattern
  • -queryBases  keep reads with motif that matches pattern
  • -tag   keep reads with this key=>value pair

 

 BamTools header ヘッダー情報を出力。

bamtools header -in input.bam

bamのheaderを出力する。samtools view -Hと同じである。

 

 BamTools merge 複数のbamを1つにマージする。3つ以上も可。

bamtools merge -in input1.bam -in input2.bam -out merge.bam

 

 BamTools random bamを指定数ランダムサンプリングして出力。

bamtools random -in input1.bam -out output.bam -n 10000

この例だと10000リードランダム出力する。bamtools formatやsamtools view -hでsamに戻し、リード数をカウントしてチェックする (wc -l input.sam)。

 

 BamTools split ユーザー定義の情報に従いbamを分割する。

bamtools split -in input.bam -mapped

上記ではリファレンスにマッピングされたリード(> input.MAPPED.bam)、マッピングされていないリード(> input.UNMAPPED.bam)に分けて出力される。利用できる定義は

  • -mapped split mapped/unmapped alignments
  • -paired split single-end/paired-end alignments
  • -reference split alignments by reference
  • -tag <tag name> splits alignments based on all values of TAG encountered (i.e. -tag RG creates a BAM file for each read group in original BAM file)

 

他にもいくつかの機能を持つ。bamtools <command> -hで各コマンドの全機能は確認できる。

  

BAMStats

bamのアライメント情報をGUIで確認できるツール。内部でsamtools statコマンドを利用してるらしい。

http://bamstats.sourceforge.net

公式サイトのリンクからダウンロードする。フォルダを解凍すると、以下のファイルが出てくる。

f:id:kazumaxneo:20170620004702j:plain

java -Xmx6g -jar BAMStats-GUI-1.25.jarコマンドラインで叩くか、このファイルをダブルクリックするとアプリが起動する。

GUI版では下のような情報を見ることができる。 

このように、平均や中央値でなく、分布でリードのアライメント状況を見ることができる。 bam/samを追加して複数データを同時比較することも可能である。

 

BAMStatsはコマンドライン版もある。ランするには以下のように打つ。

java -Xmx4g -jar BAMStats-1.25.jar -i input.bam
  •  -iでbam/samを指定する。

 

 

 

samstat

公式サイト

SAMStat

fastq、sam、bamなどを分析できる。インストールは

tar -zxcf samstat.tgz 
cd samstat
./configure
make clean
make

ビルドできたらsrc/のsamstatにパスを通す。

 

fastqの分析

samstat input.fq

 ~.htmlが出力される。

f:id:kazumaxneo:20170620183443j:plain

f:id:kazumaxneo:20170620183448j:plain

f:id:kazumaxneo:20170620183452j:plain

 

 

 

f:id:kazumaxneo:20170620183455j:plain

 

sam/bamの分析

samstat input.bam

 ~.htmlが出力される。

f:id:kazumaxneo:20170620183839j:plain

f:id:kazumaxneo:20170620183844j:plain

f:id:kazumaxneo:20170620183848j:plain

f:id:kazumaxneo:20170620183852j:plain

f:id:kazumaxneo:20170620183853j:plain

 

ショートリード向けに設計されているので、ロングリードを使うとエラーになる。

 

 

 samblasterを使えば、samファイルから異常なアライメントをされているリードを取り除いたり、duplicationを除くことができます。samblasterは別に紹介しています。


 

 

 bedファイルを出力するにはbedtoolsを使う。

 bamをbedに変換する(bed)。 

bedtools bamtobed -i reads.bam > reads.bed

 

bedtoolsについては以下にまとめています。


 

bam/samが増えてきていちいちコマンドを実行するのが面倒になったら、スクリプトを組む事を考える。ただし簡単な作業ならparallelを使い並列化した方が簡単で高速な場合も多い。

例えば特定のディレクトリにbamが複数あり、それぞれについてpicardを使いマッピングスコアを分析する。

find *bam | parallel 'Picard QualityScoreDistribution INPUT={} OUTPUT={}.quality_score_dist.tsv CHART_OUTPUT={}.quality_score_dist.pdf'

bamそれぞれについてpdfファイルとtsvファイルが出力される。

findでbamファイルを全て拾い、parallelコマンドに渡している。parallelコマンドが' 'で囲まれた部位を並列実行する。ただし負荷はファイル分増加するので気をつけること。

 

 

fasta、fastqなどを扱うなら、seqkitが高速で使い易いと思います。ランダムサインプリング、逆相補鎖に変換、ソート、重複の除去、指定した範囲だけ抽出など、なんでもできます。

以下にまとめています。


 

追記

SAMtoolsより高速なsambambaというツールがあります。SAMtoolsより全体的に速くなっておりお勧めです。


 

 

 

 追記

bcftools

非常にたくさんのコマンドがある。よく使うだろうコマンドを紹介する。 

brewインストールできる。

brew install bcftools

 

BamUtil

Github

https://github.com/statgen/bamUtil

 

 作成中。 別にまとめてリンクします。

  

 

引用

BCFtools/RoH: a hidden Markov model approach for detecting autozygosity from next-generation sequencing data.

Narasimhan V1, Danecek P1, Scally A2, Xue Y1, Tyler-Smith C1, Durbin R1.

Bioinformatics. 2016 Jun 1;32(11):1749-51. doi: 10.1093/bioinformatics/btw044. Epub 2016 Jan 30. 

https://www.ncbi.nlm.nih.gov/pubmed/26826718

 

The Sequence Alignment/Map format and SAMtools.

Li H1, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R; 1000 Genome Project Data Processing Subgroup.

Bioinformatics. 2009 Aug 15;25(16):2078-9. doi: 10.1093/bioinformatics/btp352. Epub 2009 Jun 8.

https://www.ncbi.nlm.nih.gov/pubmed/19505943

 

SAMStat: monitoring biases in next generation sequencing data.

Lassmann T1, Hayashizaki Y, Daub CO.

Bioinformatics. 2011 Jan 1;27(1):130-1. doi: 10.1093/bioinformatics/btq614. Epub 2010 Nov 18.

https://www.ncbi.nlm.nih.gov/pubmed/21088025

 

BamTools: a C++ API and toolkit for analyzing and managing BAM files. 

Barnett DW1, Garrison EK, Quinlan AR, Strömberg MP, Marth GT.

 Bioinformatics. 2011 Jun 15;27(12):1691-2. doi: 10.1093/bioinformatics/btr174. Epub 2011 Apr 14.

BamTools: a C++ API and toolkit for analyzing and managing BAM files. - PubMed - NCBI

 

parallel


samtools

http://davetang.org/wiki/tiki-index.php?page=SAMTools