macでインフォマティクス

macでインフォマティクス

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

BAMを感心対象のみにフィルタリングする VariantBam

 

 ゲノムシーケンスのコストが減少するにつれて、大規模なシーケンスデータセットを取り扱う際のストレージおよび計算上の負担が増大する懸念がある。ヒトゲノムの全ゲノムシーケンシングを30倍のカバレッジにすると、およそ10億リードのシーケンスが可能になり、圧縮BAMフォーマットでも100GB以上のディスクスペースが必要になる(論文より Li et al、2009)。これを軽減するために、BAMファイルのリファレンスベース圧縮を提供するCRAMフォーマットが確立され、ファイルサイズが2〜3倍減少した(Hsi-Yang Fritz et al、2011)。CRAMを用いれば全ゲノムの30カバレッジで30GBにできるが、現在のCRAMフォーマットは多くの解析ツールではサポートされていない。

 別の方法として、科学的な探求対象のシーケンス情報のみを含むトリミングされたBAMファイルを作成する方法がある。例えば、あるグループはエキソンの変異のみに関心があり、別のグループは構造変異を支持するリードにしか関心がないかもしれない。さらに、シークエンシングライブラリのクオリティによっては、低クオリティのリードの大部分を下流の分析にほとんど影響を与えずに除去できる。アライメント・タグを削除するか、カバレッジの高い領域でサブサンプリングを実行することで、さらにサイズを縮小することができる。これらのニーズをすべて満たしているリードのフィルタリングは、最も関連性の高い情報をコンパクトなファイルに保存しながらデータフットプリントを削減できる可能性がある。その後、元のBAMを低コストのアーカイブストレージに移動するか、または外部サーバ(例えばdbGaP)上に保持し、検索中にフィルタリングできる。

 リードのフィルタリングは、リードのサブセットのみで動作する解析パイプラインで直接使用することもできる。異なる分野(例えば、エピゲノム、メタゲノミクスなど)の研究では、データの異なる種類のサブセットに焦点を当てる傾向がある。シーケンスリードの効率的かつ正確な検索ができれば、そのような分析の開発を容易にするであろう。

 これらの問題に対処するため、カスタムリードセットを取得する堅牢で柔軟なメソッドを提供するBAM / CRAM / SAMフィルタリングツールVariantBamを開発した。VariantBamは階層型のフィルタリングシステムで、指定されたリードを出力に含めるかどうかを決定する。様々な分析(例えば、SNP、エキソン変異)で、異なるゲノムに異なるフィルタリングルールを適用しシングルパスで処理することができる。 VariantBamは、シーケンスディレクトリ、CIGAR解析、クオリティスコアの考慮事項、最大シーケンスカバレッジにサブサンプルする機能など、他のフィルタリングツールでは利用できない高度なフィルタリングオプションを提供する。VariantBamはC ++で書かれており、高速のC htslib(github.com/samtools/htslib)を使用して高効率でI / Oを実行する。 BAM、CRAM、SAM、および標準出力への出力がサポートされている。

 

 

インストール

cent os6でテストした。

Optional

  • samtools
  • bedtools

 

Github

https://github.com/walaj/VariantBam

git clone --recursive https://github.com/jwalabroad/VariantBam.git 
cd VariantBam
./configure
make
cd src/

 >  ./variant

 ./variant 

 

Usage: variant <input.bam> [OPTIONS] 

 

  Description: Filter a BAM/SAM/CRAM/STDIN according to hierarchical rules

 

 General options

  -h, --help                           Display this help and exit

  -v, --verbose                        Verbose output

  -t, --num-threads                    Add additional threads from pool for reading/writing. Per htslib, -t 1 adds one additional thread to main. [0]

  -x, --no-output                      Don't output reads (used for profiling with -q)

  -r, --rules                          JSON ecript for the rules.

  -k, --proc-regions-file              Samtools-style region string (e.g. 1:1,000-2,000) or BED/VCF of regions to process. -k UN iterates over unmapped-unmapped reads

  -Q, --mark-as-qc-fail                Flag reads that don't pass VariantBam with the failed QC flag, rather than deleting the read.

 Output options

  -o, --output                         Output file to write to (BAM/SAM/CRAM) file instead of stdout

  -C, --cram                           Output file should be in CRAM format

  -b, --bam                            Output should be in binary BAM format

  -T, --reference                      Path to reference. Required for reading/writing CRAM

  -s, --strip-tags                     Remove the specified tags, separated by commas. eg. -s RG,MD

  -S, --strip-all-tags                 Remove all alignment tags

  -Z, --write-trimmed                  Output the base-quality trimmed sequence rather than the original sequence. Also removes quality scores

 Filtering options

  -q, --qc-file                        Output a qc file that contains information about BAM

  -m, --max-coverage                   Maximum coverage of output file. BAM must be sorted. Negative values enforce a minimum coverage

  -p, --min-phred                      Set the minimum base quality score considered to be high-quality

 Region specifiers

  -g, --region                         Regions (e.g. myvcf.vcf or WG for whole genome) or newline seperated subsequence file.

  -G, --exclude-region                 Same as -g, but for region where satisfying a rule EXCLUDES this read.

  -l, --linked-region                  Same as -g, but turns on mate-linking

  -L, --linked-exclude-region          Same as -l, but for mate-linked region where satisfying this rule EXCLUDES this read.

  -P, --region-pad                     Apply a padding to each region supplied with the region flags (specify after region flag)

 Command line rules shortcuts (to be used without supplying a -r script)

      --min-clip                       Minimum number of quality clipped bases

      --max-nbases                     Maximum number of N bases

      --min-mapq                       Minimum mapping quality

      --min-del                        Minimum number of deleted bases

      --min-ins                        Minimum number of inserted bases

      --min-length                     Minimum read length (after base-quality trimming)

      --motif                          Motif file

  -R, --read-group                     Limit to just a single read group

  -f, --include-aln-flag               Flags to include (like samtools -f)

  -F, --exclude-aln-flag               Flags to exclude (like samtools -F)

 

 

ラン

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が出力される。

 

bamtobedに渡し各リードをbed出力する。

samtools view -h input.bam | variant - --min-mapq 10 -b | bamToBed > out.bed

 

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

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

VCF以外にBEDで領域指定することも可能 (e.g. "1:1,000,000-1,000,010")。 

VCFでコールされた部位周辺だけ出力されている。

f:id:kazumaxneo:20180606183919j:plain

リファレンスと合致するリードは必要ない情報と考えれば、たいへん効率的なbamのフィルタリング方法といえる。bam自身がシュリンクされているので、ビューアの負担も劇的に減る。

 

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

variant input.bam -L input.vcf -b -o mini.bam -v 

 

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

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

 

マッピングクオリティ(MAPQ wiki)が20以上で、10以上の大きな挿入か欠失をもつリードだけbam出力する。

variant input.bam --min-mapq 20 --min-ins 10 --min-del 10 -b -o mini.bam -v

 

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

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

 

 

Githubには多数のexampleが紹介されています。

https://github.com/walaj/VariantBam

 

引用

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.