2020 4/17 インストール追記、help更新
ゲノムシーケンスのコストが減少するにつれて、大規模なシーケンスデータセットを取り扱う際のストレージおよび計算上の負担が増大する懸念がある。ヒトゲノムの全ゲノムシーケンシングを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
https://github.com/walaj/VariantBam
git clone --recursive https://github.com/jwalabroad/VariantBam.git
cd VariantBam
./configure
make
cd src/
#bioconda (link)
conda install -c bioconda -y variantbam
> 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
-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でコールされた部位周辺だけ出力されている。
リファレンスと合致するリードは必要ない情報と考えれば、たいへん効率的な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.
