macでインフォマティクス

macでインフォマティクス

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

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

 

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

brew tap brewsci/science #2/25リンク修正

 

追記

--2017--

--2018--

  

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 -s input.fq

(追記)ペアエンド。

samtools fastq input.bam -1 pair1.fq -2 pair2.fq -n
  • -1    write paired reads flagged READ1 to FILE
  • -2    write paired reads flagged READ2 to FILE
  • -0    write paired reads flagged both or neither READ1 and READ2 to
  • -s     write singleton reads to FILE [assume single-end]
  • -n    don't append /1 and /2 to the read name
  • -N    always append /1 and /2 to the read name
  • FILE

ペアエンドの順番は考慮されない。ペアエンドの順番を維持したい場合は上に載せたpicardのコマンドで出力するか、いったんfastqを出力して、BBtoolsのrepairコマンド

を使って順番を修復する。↓

・ペアリードファイルの修復(BBTools)。

repair.sh -Xmx20g in1=broken1.fq in2=broken2 out1=fixed1.fq out2=fixed2.fq outs=singletons.fq repair

 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と変換するより圧倒的に早い。例えば以下のようにして使う。

#header情報だけsamに変換
samtools view -H input.bam > header.sam
# viなどでheader.bamを編集。例えばSO:unsortedをSO:sortedに修正。
# sample名をつけ忘れていたなら,@RGのSM:部分を編集する => SM:sample1

#元のbamのヘッダーを編集したheaderで置換。
samtools reheader header.sam input.bam > reheaer.bam

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

 PCR duplicateを除く。

samtools rmdup input.bam output.bam

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

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

samtools addreplacerg 

heterozygous SNPsのコール。

samtools phase 

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

samtools sort -@ 12 -m 2G pair.bam -o pair_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で指定することで、特に大きなゲノムのソートが速くなる。上では2GB/thread指定して、スレッド数を12にしている。入力はbamとsamに対応している。samを指定すると、ソートしたbamを出力してくれる。

samからポジションソート済みbamを作成。

samtools sort -@ 12 -m 8G -O BAM pair.sam -o pair_sorted.bam
  • -O < FORMAT>   Specify output format (SAM, BAM, CRAM)

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

#samからbamに変換
samtools view -@ 10 -bS input.sam > output.bam

#bamからsamに変換
samtools view -@ 10 -h input.bam > output.sam

#bamの中身を見る。headerは見れない。
samtools view -@ 10 input.bam |less

#bamのheaderだけ見る。ソートされているか
samtools view -@ 10 -H input.bam |less

#50%ランダムサンプリング
samtools view -@ 10 -s 0.5 -b input.bam > random0.5.bam

#mマッピングされたリードだけbamとして取り出す。
samtools view -@ 10 -h -b -F 4 input.bam > mapped.bam

#chr1とchr2にマッピングされたリードだけ取り出す
samtools view -@ 10 -b input.bam chr1 chr2 > chr1-2.bam

#chr1の100-300にマッピングされたリードだけ取り出す。
samtools view -@ 10 -b input.bam chr1:100-300 > chr1_100-300.bam
  • -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]
  • -@ Number of additional threads to use [0]

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

#50%ランダムサンプリング
samtools view -s 0.5 -b input.bam > random0.5_input.bam

#mマッピングされたリードだけbamとして取り出す。samなら-bを消す。
samtools view -h -b -F 4 input.bam > mapped.bam

#chr1にマッピングされたリードだけ取り出す
samtools view -b input.bam chr1 > chr1.bam

#chr1の100-300にマッピングされたリードだけ取り出す。
samtools view -b input.bam chr1:100-300 > chr1_100-300.bam

追記 samからmappingされたリードだけbamに変換。

samtools view -@ 10 -bS -F 4 input.sam > output.bam

追記 samからmappingされたリードだけbamに変換して、さらにfastq出力(unmapを除き、mapされたリードだけfastqで回収)。

samtools view -@ 10 -bS -F 4 input.sam |samtools fastq -@ 10 - > output.fq

追記 bamからchr1にmappingされたリードだけ抽出して、さらにfastq出力。

samtools view -@ 10 -b input.bam chr1 > chr1.bam
samtools index chr1.bam
samtools fastq -@ 10 chr1.bam > chr1_mapped.fq

追記 mapping quality (MAPQ) が1以上のリードだけ出力。

samtools view -@ 10 -bSq 1 input.sam > filtered.bam

mapping qualityについてはこちらのページが参考になる(リンク)。 

bamのindex。

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

追記 bamがおかしくなっていないか、カレントの全てのbamを調べる。

samtools quickcheck -v *.bam > bad_bams.fofn && echo 'all ok' || echo 'some files failed check, see bad_bams.fofn'

bamがおかしくなっていないか(download中に切れたことを気づかないまま進めたり、変換中におちたbamをつかったりなど)。 "||"をつかうことで、エラーがある時だけ||以降が実行されるようにしている。

追記

bamがcoordinateソート(chr順&ポジション順)されているかどうか。

samtools view -H iput.bam | grep @HD

ソートされていれば、 SO:coordinateが表示される。

sample名取得。

 samtools view -H iput.bam | grep '^@RG' | sed "s/.*SM:\([^\t]*\).*/\1/g" | uniq

追記。

たくさんのbamの並列処理。GNU parallelsを使っている。同時処理数を"-j"オプションで6つまでに制限している。samtools sortコマンド自体も4スレッドの並列処理をしている。

ls *.bam | parallel -j 6 "samtools sort -@ 4 {} > {.}_sorted.bam; samtools index {.}_sorted.bam"

splitした小さなbamなどを並列処理する時に一番高速化が期待できます。ただし、samtoolsのコマンドはI/O intensitiveなので、同時処理数が増えすぎるとかえって速度が遅くなったり、最悪、コンピュータがフリーズする恐れがあります。注意して使ってください。parallelのあとに"--dry-run"をつけると(parallel --dry-run)、コマンド内容は実行せず、実行コマンドだけ表示されます。まず--dry-runをつけて実行してください。時間短縮を考えるなら、elprep splitなどでbamを分割し、timeをつけて実行時間を測定し、最速のくみあせわを検討することをお勧めします。I/OやCPUなど環境の違いによって最速の組み合わせは変化するはずです。

追記。

minimap2でmappingし、coordinate sortしたbamを出力。スレッドは12、メモリは2G / thread。リードグループ(@RG)のサンプル名はsample1。

minimap2 -R "@RG\tID:X\tLB:Y\tSM:sample1\tPL:ILLUMINA" -t 12 -ax sr \
ref.fa pair1.fq.gz pair2.fq.gz \
|samtools sort -@ 12 -m 2G -O BAM - > sorted.bam \
&& samtools index sorted.bam

minimap2でmappingし、elprepでフィルタリングしてbam出力。elprep以外は上と同じだが、

minimap2 -ax sr -R "@RG\tID:X\tLB:Y\tSM:sample1\tPL:ILLUMINA" -t 12 \
-ax sr ref.fa pair1.fq.gz pair2.fq.gz \
|elprep filter /dev/stdin /dev/stdout \
--mark-duplicates --remove-duplicates --filter-mapping-quality 1 \
--clean-sam --nr-of-threads 12 \ --sorting-order coordinate \
|samtools sort -@ 12 -m 2G -O BAM - > sorted.bam \
&& samtools index sorted.bam

elprepの mark duplicatesとsortはin memoryで動作するため、メモリ要求量はデータサイズに直接比例します(samファイルサイズの10倍近くのメモリが必要)。メモリが足りなければ、--mark-duplicatesと--remove-duplicatesを外すか、samのsplitを行ってから順次フィルタリングするよう変更して下さい。

 

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)

クロモソームごとに分割するなら、 -referenceのフラグをつける。

bamtools split -in input.bam -reference

 input_chrromosome1.bam、input_chromosome2.bam...unmapped.bamのようにリファレンスのFASTAのヘッダー名付のbamが出力される。

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

 

 

bamutils

keepbest(リンク best mappingを抽出。

bamutils keepbest input.bam output.bam 

export(リンク  指定した情報を全リードから出力

bamutils export -name input.bam > output #リード名を出力

bamutils export -ref input.bam > output #mappingされたクロモソーム名を出力

expressed(リンク マッピングされた領域をBED形式で出力。(samtools indexでbam.baiファイルを準備しておく必要がある)

bamutils expressed input.bam > output.bed

filter(リンク 指定した条件を満たすリードだけbamで出力。

#chr1にマッピングされており、リード長が100以上でペアエンドが-> <-の向きに両方マップングされているリードだけbamで出力。
bamutils filter input.bam output.bam -minlen 100 -properpair -includeref chr1

merge(リンク bamをマージ(best matchを選択してくる)。

bamutils merge output.bam input1.bam input2.bam

pair(リンク別々にマッピングされているペアエンドについて、指定条件内でペアを再度調べる。

bamutils pair output.bam input1.bam input2.bam

removevlipping(リンクhard-clipされた領域を捨てる。

bamutils removeclipping input.bam output.bam

renamepair(リンク一部アライナーはペアエンドのリード名に/1や/2をつけてbsmを作り、これが下流の解析でエラーになる。これを消す(ペアエンドは同じリード名になる)。

bamutils renamepair input.bam output.bam

tag(リンクtagをbamに追加する。

bamutils tag -xs input.bam output.bam

check(リンク破損のチェック。

bamutils check input.bam 

cleancigar(リンクCIGAR情報の修復。

bamutils cleancigar  input.bam

stats(リンクbamを分析。

bamutils stats input.bam

$ bamutils stats input.bam

Calculating Read stats...

Done! (0:36)                                                              

out.bam

Reads: 1083286

Mapped: 1083286

Unmapped: 0

 

Flag distribution

[0x001] Multiple fragments       1083286 100.00%

[0x002] All fragments aligned    888871 82.05%

[0x008] Next unmapped            58635 5.41%

[0x010] Reverse complimented     542299 50.06%

[0x020] Next reverse complimented 539537 49.81%

[0x040] First fragment           1083286 100.00%

[0x800] Supplementary            1 0.00%

 

Template length: 568.90 +/- 824.04

 

Reference counts count

1 1083286

 

 

 追記

elprep(紹介)

filter 以下の作業をワンライナーで行う。

  • アンマップリードを完全に除く(--filter-unmapped-reads-strict)。
  • duplicationのタグを付け(--mark-duplicates)duplicationリードを捨てる(--remove-duplicates)。
  • mapping qualityがゼロのリード(bwaならmultiple mappingはゼロ)を捨てる(--filter-mapping-quality
  • samのcleaning(--clean-sam)。
elprep filter input.sam output.sam --mark-duplicates --remove-duplicates --filter-mapping-quality 0 --clean-sam

split   インプット (sam/bam/cram)をクロモソームごとに分割する。unmapリードは別出力する。

elprep split input.sam output/ --output-prefix "split-sam" --output-type sam --nr-of-threads 8 --single-end 

 ソートやindexファイルがなくても動作する。シングルエンドのinputは--single-endをつける。デフォルトでは入力をペアエンドとして扱う。

merge  splitで分割した複数のインプット (sam/bam/cram)をマージする。

elprep merge input/ output.sam --nr-of-threads 8 

シングルエンドのinputは--single-endをつける。デフォルトでは入力をペアエンドとして扱う。bamとcram出力にも対応している。

 

VariantBam(紹介)

chr1の1-1,100,000にマッピングされたリードのみbam出力する。

variant input.bam -g 'chr1:1-1,100,000' --min-mapq 10 -b -o mini.bam -v

sam出力なら-bを外す。mini.bamが出力される。

VCFでコールされたポジション周囲100-bpのリードのみbam出力する。

variant input.bam -l input.vcf -P 100 -b -o mini.bam -v 

sub-samplingしてmax-coverageを100に減らす(マッパーによっては、リピート領域にMAPQ=0のリードが多数マッピングされる)。

variant input.bam -m 100 -o mini.bam -v -b

clippingされたハイクオリティなリードのみbam出力する。

variant input.bam --min-phred 4 --min-clip 5 -b -o mini.bam -v

 

goleft(紹介)

covstats   bamをサンプリングしてカバレッジとインサートサイズをレポートする。

goleft covstats pair.bam

カバレッジ、インサートサイズとそのSDなどが出力される。 

f:id:kazumaxneo:20180905091943j:plain

insert sizeは-> <-の内側サイズが計算される。250x2でinner inset sizeが100なら、outer inset sizeは600。

 

depth  一定のウィンドウサイズでbamのカバレッジを計算する。 

goleft depth --reference ref.fa --prefix output input.bam 

indexcov   bam.baiからのカバレッジの超高速な推定。30サンプルを30秒で解析可能。複数のbamをディレクトリに準備してランする。

goleft indexcov --directory output/ bam_data/*.bam

関係ないファイルがあるとエラーになるのでbamとbaiだけ集める。htmlで結果は出力される。

 

biobambam2(紹介)

1、duplicateにタグをつけ、coordinateソートも行って出力する。

bamsormadup < input.bam level=1 threads=8 inputformat=bam SO=coordinate outputformat=bam indexfilename=output.bam.bai > output.bam 

処理の大半が並列化されており、I/Oがリミットになりうるため、I/Oの高速なSSDなどのストレージ使用が推奨されている。threadsフラグがなければ利用できる全ての論理コアを使う。

2、ペアエンドを照合し、secondary hitは除く。samtools collateと同様の機能(samtools1.8 manual)。 level=0だと非圧縮になる。

bamcollate2 <input.bam collate=1 level=1 exclude=SECONDARY > out.bam 

3、bamを再圧縮する。

bamrecompress <input.bam level=9 numthreads=8 > output.bam

 4、ポジションソートしてbam出力する。

./bamsort SO=coordinate < input.bam inputthreads=8 level=1 outputformat=bam > out.bam 

 5、bamからペアエンドfastqを出力。

bamtofastq < out.bam F=1.fq F2=2.fq

 

Sambamba(紹介)

view bamの基本情報

sambamba view --reference-info input.bam
  •  -I, --reference-info  output to stdout only reference names and lengths in JSON 

view samの基本情報

sambamba view -S --reference-info input.bam
  • -S specify that input is in SAM format

bamのリファレンス"chr19"に適切にアライメントされるリードだけsam出力。(カウントなら-cをつけ、-o以下を消す)

sambamba view -F "proper_pair" input.bam chr19 -t 8 -f sam -o output.sam
    • -o specify output filename

sortリンク ソート。

指定がなければindexが自動でつくcoordinate-sorted。

sambamba sort input.bam -o sorted.bam -t 20 -p -l 6
  • -t use specified number of threads
  • -p show progressbar in STDERR
  • -l level of compression for sorted BAM, from 0 to 9(指定がなければ"6"くらいのサイズになる)

名前でソート。(indexは付けられない)

sambamba sort input.bam -o sorted.bam -t 20 -p -l 5 

flagstatリンク bamの情報表示。

bamの情報表示。

sambamba flagstat input.bam -t 8 -p

sambamba flagstat input.bam

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

0 + 0 secondary

538 + 0 supplementary

0 + 0 duplicates

456058 + 0 mapped (89.55%:N/A)

508742 + 0 paired in sequencing

254371 + 0 read1

254371 + 0 read2

451458 + 0 properly paired (88.74%:N/A)

454670 + 0 with itself and mate mapped

850 + 0 singletons (0.17%:N/A)

0 + 0 with mate mapped to a different chr

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

表示された情報については上のflagstatのリンク先参照。bamのBGZFフォーマットしか対応していないのでsamはviewで変換してから使う。

indexリンク bamのindexファイル作成

sambamba index input.bam -t 12 -p

mergeリンク bamのマージ。

圧縮率最大の設定で、12スレッド使い3つのbamをマージ。

sambamba merge merge.bam input1.bam input2.bam nput3.bam -t 12 -p -l 9
  • -l Specify compression level of output file as a number from 0 to 9

mpileupリンク samtoolsのmpileupのparallelバージョン。

20スレッド使い、pileupを実行。samtoolsの条件は--samtoolsをつけて書く。

sambamba mpileup input.bam -t 12 --samtools samtools -gu -S -D -d 1000 -L 1000 -m 3 -F 0.0002

 

Samblaster(紹介)

 bwa memのアライメントの過程でduplicationにタグをつける。

bwa index -a is input.fa
bwa mem -t 12 -R "@RG\tID:X\tLB:Y\tSM:Z\tPL:ILLUMINA" input.fa *.fastq | samblaster |samtools view -@ 12 -Sb - |samtools sort -@ 12 - > samp.sorted.bam

discordant-readとsplit-readは別ファイルに出力。

bwa mem -t 12 -R "@RG\tID:X\tLB:Y\tSM:Z\tPL:ILLUMINA" input.fa *.fastq | samblaster -e -d samp.disc.sam -s samp.split.sam > output.sam
  • -e Exclude reads marked as duplicates from discordant, splitter, and/or unmapped file.
  • -d FILE Output discordant read pairs to this file. [no discordant file output]
  • -s FILE Output split reads to this file abiding by paramaters below. [no splitter file output]

  duplication-readは全出力から除く。

bwa mem -t 12 -R "@RG\tID:X\tLB:Y\tSM:Z\tPL:ILLUMINA" input.fa *.fastq | samblaster -r -e -d samp.disc.sam -s samp.split.sam > output.sam
  •  -r Remove duplicates reads from all output files.  

注意;bwa memで-M(mark shorter split hits as secondary)をつけている時は、samblasterにも-Mをつけてランを行う

Alfred(紹介)

リファレンスのFASTAと.bamを指定する。

alfred qc -r ref.fa -o qc.tsv.gz input.bam

#結果の可視化
#なければ本体をcloneして取ってくる
git clone https://github.com/tobiasrausch/alfred.git
#可視化
Rscript alfred/R/stats.R qc.tsv.gz

出力例(direct link)。

 

fastqp(紹介)

fastqp input.bam

複数のmetricsで分析が行われる。下はその一部。

f:id:kazumaxneo:20180622090028p:plain

 

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:20170620183455j:plain

 

sam/bamの分析

samstat input.bam

 ~.htmlが出力される。

f:id:kazumaxneo:20170620183839j:plain

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

 

bcftools

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

brewインストールできる。

brew install bcftools

時間があればまとめます。 

 

 

 

 

 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が高速で使い易いと思います。ランダムサインプリング、逆相補鎖に変換、ソート、重複の除去、指定した範囲だけ抽出など、なんでもできます。

以下にまとめています。

bamからのfastqの抽出はtophat (ver2.1)でもできるようです。

bam2fastx -A -o output.fastq input.bam

 

 

引用

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

 

VariantBam: filtering and profiling of next-generational sequencing data using region-specific rules.

Wala J, Zhang CZ, Meyerson M, Beroukhim R

Bioinformatics. 2016 Jul 1;32(13):2029-31.

 

biobambam: tools for read pair collation based algorithms on BAM files

German Tischler, Steven Leonard

Source Code Biol Med. 2014; 9: 13.

 

elPrep: High-Performance Preparation of Sequence Alignment/Map Files for Variant Calling.

Herzeel C, Costanza P, Decap D, Fostier J, Reumers J.

PLoS One. 2015 Jul 16;10(7):e0132868.

 

 

parallel


samtools

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

 

bamのsplit

https://www.biostars.org/p/46327/

 

Biostars

https://www.biostars.org/p/247903/

 

Tip and tricks for BAM files

https://github.com/IARCbioinfo/BAM-tricks