macでインフォマティクス

macでインフォマティクス

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

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

 

とにかく色々なツールがあり混乱するが、データ処理するため、使えると便利なものもある。代表的なツールの機能を紹介する。

 

bcftools

 

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

brewインストールできる。

brew install bcftools

 

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で描画された図が出力される。

 

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

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

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

f:id:kazumaxneo:20170617001516j:plain

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

 

 

 

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 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

input.bamのheader.samのヘッダー情報で置換する。bam=>sam=>bamと変換するより圧倒的に高速である。

 

PCR duplicateを除く。

samtools rmdup input.bam output.bam

ペアリードのみに使える。ペアリードがFRの向きで同じクロモソームの同じ位置にアライメントされていると、マッピングクオリティがベストのものだけ残し、あとのペアを排除する。ペアの片方しかアライメントされていないものには適用されない。

 

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

samtools addreplacerg 

 

heterozygous SNPsのコール。

samtools phase 

 

samtools depad 

 

 

 

 

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で各コマンドの全機能は確認できる。

 

 

 

bedtools

作成中。 

 

 

 

 

 

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

どうやらショートリード専用らしい。ロングリードを使うとエラーになる。

 

 

 

引用

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