8/19 pindel、lumpy、cnvkit、breakdanerコマンドミス修正
SVは、ゲノムの多様性およびゲノムの障害に寄与することに関与している（Stankiewicz and Lupski、2010）。したがって、SVの検出には相当量の作業が行われている。一般に、SVを検出するためのツールは、リードのアライメントからの以下の信号のうちの1つまたは複数を使用する：split-read [分割アライメント。例: Pindel（Ye et、2009）]、read-pair[異常なペアエンドアライメント。例: BreakDancer（Chen et al、2009）]、カバレッジデプス[異常なカバレージ。例: CNVnator（Abyzov et al、2011）]、ジャンクションマッピング[既知のSVブレークポイントへのアライメント。例: BreakSeq2（Abyzov et al、2015; Lam et al、2010）]またはブレークポイント周辺のアセンブリ[MindTheGap（Rizk et al、2014）]。しかし、SVタイプとサイズのニッチがあり、すべてのタイプのSVを包括的に検出するシグナルはない。これは、SV検出を改善するための複数の方法を統合するツールの開発を必要とする。
以前の研究（Lam et al、2012; Mills et al、2011）は、複数のツールと方法によって作成されたバリアントコールが一般により正確であることを示している。このため、複数の方法を採用するためのツールが開発されている。 DELLY（Rausch et al、2012）、LUMPY（Layer etal、2014）、SVMerge（Wong etal、2010）などの複数のツールの結果をマージするものなどがある。ただし、LUMPYとDELLYは挿入を検出することができず、SVMergeはコールのマージの際に個々のツールのSV解像度を無視する。したがって、この研究は、異なるタイプとサイズのSVを高い精度と解像度で検出するための既存のSVマージングツールの限界に取り組もうとしている。
MetaSVは、複数のSV検出方法とツールを使用して、信頼性が高く正確なSVコールセットを検出する。 MetaSVの新しさは、以下の重要なアイデアの組み合わせにある。複数の独立した方法の組み合わせで報告されたコールは一般に品質が良く、動的プログラミングによるローカルアセンブリを使用してSVブレークポイントを絞り込むことができる。
MetaSVは、以下のステップ（論文 図1）で進行する。
Intra-tool merging：同じツールで生成された潜在的なduplicate callsがここでマージされる。十分なオーバーラップがある場合、2つのコールはマージされる。
Inter-tool merging：複数のツールで生成されたコールがマージされる。複数のツールに共通するコールのブレークポイントを決定する際には、精度が高いことが知られているメソッドが優先される（e.g., スプリット・リード情報の方がリード・ペア情報より優先）。このメソッド認識マージは、報告されたSVのブレークポイントが可能な限り正確であることを確実にするために重要である。
Local assembly：追加の証拠を収集してSVの配列を決定するため、ツールによって報告されたSV領域に対してローカルアセンブリが実行される。アセンブリのため独自のペア情報を使用する機能があるため、SPAdes（Bankevich etal、2012）アセンブラが使用される。
Breakpoint resolution：アセンブリされたSV配列は、ダイナミックプログラミングを使用してブレークポイントを検出または調整するためリファレンスとアライメントされる（Abyzov and Gerstein、2011）。
Genotyping：SVブレークポイント周辺のカバレッジがSVのzygosityを決定するために使用される。最終的な出力は、各SVのgenotypeを持つVCFファイルとして生成される。
Annotation：MetaSVは入力と出力をVCFで標準化している。各SVには、個々のツールによって行われた対応するコーラーを示し、その信頼レベルを分類するために注釈が付けられる。複数のツールで検出されたSVは、高い信頼性があるとみなされる。
MetaSV HP
http://bioinform.github.io/metasv/
MetaSV: An Accurate and Integrative Structural-Variant Caller for Next Generation Sequencing
MetaSVパイプライン。論文より転載。
MetaSVに関するツイート。
インストール
mac os 10.12のAnaconda環境でテストしたが、アセンブリをonにするとエラーを起こした（spades 3.6.2と3.12でテストした）。
依存
- pybedtools: Note that bedtools (version 2.21 or greater) has to be installed separately in order for pybedtools to work. Please use pybedtools version 0.6.9.
- pysam-0.7.7: MetaSV has only been tested with version 0.7.7 of pysam.
- pyvcf: Please use pyvcf version 0.6.7.
In addition, paths to the following tools must be provided as MetaSV arguments:
- SPAdes: Used for assembly around the SV breakpoints
- AGE: Used for determining the SV breakpoints using assembled contigs. Please compile AGE without OpenMP support using make OMP=no.
spades（bioconda）、age（bioconda）もcondaで導入可能（"spades.py"と "age_align"でヘルプが出るなら、whichでパス確認）。
spadesは"conda search spades"でインストール可能なバージョンを調べ、3.6.2を導入した(*5)。
conda search spades
conda install -c bioconda spades==3.6.2
上記に加え、metaSVに入力するSVをコールするツールも必要になる。 ヘルプには以下のSV検出ツール（SV caller）のオプションが記載されている。
- Pindel
- breakdancer (breakdancer-max)
- breakseq (2も含む?)
- cnvnator
- gatk
- manta
- lumpy
- cnvkit
- wham
SV検出に使うcallerだけ導入すればよい。インストール方法については、各コマンド実行部分に簡単な説明を残した。
本体 Github
condaで導入できる（リンク）。
conda install -c bioconda metasv
#permittionエラーが出たので、てテスト時はsudoで0.4.0を入れた (" metasv==0.4.0")
> run_metasv.py
$ run_metasv.py -h
usage: run_metasv.py [-h] --sample Sample
[--pindel_vcf pindel_vcf [pindel_vcf ...]]
[--pindel_native File list [File list ...]]
[--breakdancer_vcf breakdancer_vcf [breakdancer_vcf ...]]
[--breakdancer_native File list [File list ...]]
[--breakseq_vcf breakseq_vcf [breakseq_vcf ...]]
[--breakseq_native breakseq_native [breakseq_native ...]]
[--cnvnator_vcf cnvnator_vcf [cnvnator_vcf ...]]
[--cnvnator_native File list [File list ...]]
[--gatk_vcf file [file ...]]
[--manta_vcf MANTA_VCF [MANTA_VCF ...]]
[--lumpy_vcf LUMPY_VCF [LUMPY_VCF ...]]
[--cnvkit_vcf CNVKIT_VCF [CNVKIT_VCF ...]]
[--wham_vcf WHAM_VCF [WHAM_VCF ...]]
[--mean_read_length MEAN_READ_LENGTH] --reference
reference [--chromosomes CHROMOSOMES [CHROMOSOMES ...]]
[--gaps gaps] [--filter_gaps] [--keep_standard_contigs]
[--bams BAMS [BAMS ...]] [--isize_mean ISIZE_MEAN]
[--isize_sd ISIZE_SD] [--wiggle WIGGLE]
[--inswiggle INSWIGGLE] [--minsvlen MINSVLEN]
[--maxsvlen MAXSVLEN] [--overlap_ratio OVERLAP_RATIO]
[--min_avg_base_qual MIN_AVG_BASE_QUAL]
[--min_mapq MIN_MAPQ] [--min_soft_clip MIN_SOFT_CLIP]
[--max_nm MAX_NM] [--min_matches MIN_MATCHES]
[--min_support_ins MIN_SUPPORT_INS]
[--min_support_frac_ins MIN_SUPPORT_FRAC_INS]
[--max_ins_intervals MAX_INS_INTERVALS]
[--mean_read_coverage MEAN_READ_COVERAGE]
[--min_ins_cov_frac MIN_INS_COV_FRAC]
[--max_ins_cov_frac MAX_INS_COV_FRAC]
[--sc_other_scale SC_OTHER_SCALE] [--spades SPADES]
[--spades_options SPADES_OPTIONS]
[--spades_timeout SPADES_TIMEOUT] [--disable_assembly]
[--svs_to_assemble {DUP,INV,DEL,INS} [{DUP,INV,DEL,INS} ...]]
[--svs_to_softclip {DUP,INV,DEL,INS} [{DUP,INV,DEL,INS} ...]]
[--extraction_max_read_pairs EXTRACTION_MAX_READ_PAIRS]
[--spades_max_interval_size SPADES_MAX_INTERVAL_SIZE]
[--assembly_max_tools ASSEMBLY_MAX_TOOLS]
[--assembly_pad ASSEMBLY_PAD] [--stop_spades_on_fail]
[--age AGE] [--age_timeout AGE_TIMEOUT]
[--min_inv_subalign_len MIN_INV_SUBALIGN_LEN]
[--min_del_subalign_len MIN_DEL_SUBALIGN_LEN]
[--age_window AGE_WINDOW] [--boost_sc]
[--gt_window GT_WINDOW] [--gt_normal_frac GT_NORMAL_FRAC]
[--svs_to_report {INV,CTX,INS,DEL,ITX,DUP} [{INV,CTX,INS,DEL,ITX,DUP} ...]]
[--enable_per_tool_output] [--workdir WORKDIR]
[--num_threads NUM_THREADS] --outdir OUTDIR [--version]
Merge SVs from multiple tools for accurate SV calling
optional arguments:
-h, --help show this help message and exit
Input data options:
--sample Sample Sample name (default: None)
--pindel_vcf pindel_vcf [pindel_vcf ...]
VCF file or dir for Pindel VCFs (default: )
--pindel_native File list [File list ...]
Pindel native files (default: )
--breakdancer_vcf breakdancer_vcf [breakdancer_vcf ...]
VCF file or dir for BreakDancer VCFs (default: )
--breakdancer_native File list [File list ...]
BreakDancer native files (default: )
--breakseq_vcf breakseq_vcf [breakseq_vcf ...]
VCF file or dir for BreakSeq VCFs (default: )
--breakseq_native breakseq_native [breakseq_native ...]
BreakSeq native GFF files (default: )
--cnvnator_vcf cnvnator_vcf [cnvnator_vcf ...]
VCF file or dir for CNVnator VCFs (default: )
--cnvnator_native File list [File list ...]
CNVnator native files (default: )
--gatk_vcf file [file ...]
VCF file or dir for gatk VCFs (default: )
--manta_vcf MANTA_VCF [MANTA_VCF ...]
VCF file or dir for Manta VCFs (default: )
--lumpy_vcf LUMPY_VCF [LUMPY_VCF ...]
VCF file or dir for Lumpy VCFs (default: )
--cnvkit_vcf CNVKIT_VCF [CNVKIT_VCF ...]
VCF file or dir for CNVkit VCFs (default: )
--wham_vcf WHAM_VCF [WHAM_VCF ...]
VCF file or dir for WHAM VCFs (default: )
--mean_read_length MEAN_READ_LENGTH
Mean read length (default: 100)
Reference options:
--reference reference
Reference file (default: None)
--chromosomes CHROMOSOMES [CHROMOSOMES ...]
Chromosome list to process. If unspecified, then all
chromosomes will be considered. (default: )
--gaps gaps Gap bed file (default: None)
--filter_gaps Filter out gaps (default: False)
--keep_standard_contigs
Keep only the major contigs + MT (default: False)
Input BAM options:
--bams BAMS [BAMS ...]
BAMs (default: [])
--isize_mean ISIZE_MEAN
Insert size mean (default: 350.0)
--isize_sd ISIZE_SD Insert size standard deviation (default: 50.0)
Tool output merging options:
--wiggle WIGGLE Wiggle for interval overlap (default: 100)
--inswiggle INSWIGGLE
Wiggle for insertions, overides wiggle (default: 100)
--minsvlen MINSVLEN Minimum length acceptable to be an SV (default: 50)
--maxsvlen MAXSVLEN Maximum length SV to report (default: 1000000)
--overlap_ratio OVERLAP_RATIO
Reciprocal overlap ratio (default: 0.5)
Soft-clip detection options:
--min_avg_base_qual MIN_AVG_BASE_QUAL
Minimum average base quality (default: 20)
--min_mapq MIN_MAPQ Minimum MAPQ (default: 5)
--min_soft_clip MIN_SOFT_CLIP
Minimum soft-clip (default: 20)
--max_nm MAX_NM Maximum number of edits (default: 10)
--min_matches MIN_MATCHES
Mininum number of matches (default: 50)
--min_support_ins MIN_SUPPORT_INS
Minimum read support for calling insertions using
soft-clips (including neighbors) (default: 15)
--min_support_frac_ins MIN_SUPPORT_FRAC_INS
Minimum fraction of reads supporting insertion using
soft-clips (including neighbors) (default: 0.05)
--max_ins_intervals MAX_INS_INTERVALS
Maximum number of insertion intervals to generate
(default: 10000)
--mean_read_coverage MEAN_READ_COVERAGE
Mean read coverage (default: 50)
--min_ins_cov_frac MIN_INS_COV_FRAC
Minimum read coverage around the insertion breakpoint.
(default: 0.5)
--max_ins_cov_frac MAX_INS_COV_FRAC
Maximum read coverage around the insertion breakpoint.
(default: 1.5)
--sc_other_scale SC_OTHER_SCALE
Control degree of incorporation of breakpoints from
other methods. (default: 5)
Assembly options:
--spades SPADES Path to SPAdes executable (default: None)
--spades_options SPADES_OPTIONS
Options for SPAdes (default: )
--spades_timeout SPADES_TIMEOUT
Maximum time (in seconds) for running SPAdes on an
interval (default: 300)
--disable_assembly Disable assembly (default: False)
--svs_to_assemble {DUP,INV,DEL,INS} [{DUP,INV,DEL,INS} ...]
SVs to assemble (default: ['INS', 'INV', 'DUP'])
--svs_to_softclip {DUP,INV,DEL,INS} [{DUP,INV,DEL,INS} ...]
SVs to soft-clip (default: ['INS', 'INV', 'DUP'])
--extraction_max_read_pairs EXTRACTION_MAX_READ_PAIRS
Maximum number of pairs to extract for assembly
(default: 10000)
--spades_max_interval_size SPADES_MAX_INTERVAL_SIZE
Maximum SV length for assembly (default: 50000)
--assembly_max_tools ASSEMBLY_MAX_TOOLS
Skip assembly if more than this many tools support a
call (default 1) (default: 1)
--assembly_pad ASSEMBLY_PAD
Padding base pairs to use for assembling breakpoint
with Spades and AGE (default: 500)
--stop_spades_on_fail
Abort on SPAdes failure (default: False)
--age AGE Path to AGE executable (default: None)
--age_timeout AGE_TIMEOUT
Maximum time (in seconds) for running AGE on an
assembled contig (default: 300)
--min_inv_subalign_len MIN_INV_SUBALIGN_LEN
Minimum length of inversion sub-alginment (default:
50)
--min_del_subalign_len MIN_DEL_SUBALIGN_LEN
Minimum length of deletion sub-alginment (default: 50)
--age_window AGE_WINDOW
Window size for AGE to merge nearby breakpoints
(default: 20)
--boost_sc Use soft-clips for improving breakpoint detection
(default: False)
Genotyping options:
--gt_window GT_WINDOW
Window for genotyping (default: 100)
--gt_normal_frac GT_NORMAL_FRAC
Min. fraction of reads supporting reference for
genotyping (default: 0.05)
Output options:
--svs_to_report {INV,CTX,INS,DEL,ITX,DUP} [{INV,CTX,INS,DEL,ITX,DUP} ...]
SV types to report (default: set(['INV', 'CTX', 'INS',
'DEL', 'ITX', 'DUP']))
--enable_per_tool_output
Enable output of merged SVs for individual tools
(default: False)
Running environment options:
--workdir WORKDIR Scratch directory for working (default: work)
--num_threads NUM_THREADS
Number of threads to use (default: 1)
--outdir OUTDIR Output directory (default: None)
Other options:
--version show program's version number and exit
——
テストラン
git clone https://github.com/bioinform/metasv.git
cd metasv/test/
./test_run.sh
test_run.sh実行前
実行後。out/にvariants.vcfができている。
variants.vcfの中身。
動作している。
ラン
いくつかのSVコーラーは以前まとめて紹介してます（リンク）。ここでは各ツールの SVコールコマンドと、それを入力としてmetaSVを実行する手順を確認する。
１、mapping
ペアエンドのfastqからbamを作成する。リファレンスはhuman_g1k_v37_decoy.faとする。minimap2でマッピングし（*1）、samtoolsにパイプしてソートし、BAM出力（*2）（*6, *7）。
#normal fastq
minimap2 -ax sr -R "@RG ID:X LB:Y SM:NA12877 PL:ILLUMINA"
-t 16 human_g1k_v37_decoy.fa pair1.fg pair2.fq
|samtools sort -@ 16 -m 4g -O BAM - > sample.bam
samtools index sample.bam
PCR duplicationを除き、リアライメントまで実行したいときは、こちらを参照。
２、Breakdancer-maxでSVコール（リンク）。dokcer imageを使う。
#コンテナ内で作業する
docker pull molecular/breakdancer #pull
docker run --rm -itv /Volumes/10TB/test/:/data/ molecular/breakdancer
> bam2cfg.pl sample.bam > breakdancer.cfg
> breakdancer-max breakdancer.cfg > breakdancer_output
３、PindelでSVコール（リンク）。dokcer imageを使う。
#コンテナ内で作業する
docker pull opengenomics/pindel
docker run -m "50g" --rm -itv /Volumes/10TB/test/:/data/ opengenomics/pindel
> #configファイル作成
> printf "sample.bam_350_sample1" > config
> sed -e 's/_/ /g' config > pindel_config && rm config #tabを入れたかったので少し遠回りした
> #pindel実行
> pindel -f human_g1k_v37_decoy.fa -T 10 -i pindel_config -o pindel
> cat pindel_SI pindel_D pindel_LI pindel_TD > pindel_merge
> pindel2vcf -R human_g1k_v37_decoy.fa -r human_g1k_v37_decoy.fa -p pindel_merge -v pindel.vcf -d 20180801 -e 4
> exit #コンテナを抜ける
vcffilter -f " SVTYPE = INS | SVTYPE = DEL " pindel.vcf > pindel.indel.vcf
４、MantaでSVコール（リンク）。dokcer imageを使う。
docker pull eitanbanks/manta-1.0.3
#configファイル作成
docker run --rm -itv /Volumes/10TB/test/:/data eitanbanks/manta-1.0.3
/bin/manta/bin/configManta.py
--normalBam=/data/sample.bam
--referenceFasta=/data/human_g1k_v37_decoy.fa
--runDir=/data/output
#configファイルを実行
docker run --rm -itv /Volumes/10TB/test/:/data eitanbanks/manta-1.0.3
data/output/runWorkflow.py -m local -j 8 --memGb=20
５、LumpyでSVコール（リンク）。Lumpyはbrewで導入した。
#normal bam
samtools view -@ 12 -b -F 1294 sample.bam > sample.discordants.unsorted.bam
samtools view -@ 12 -h sample.bam | extractSplitReads_BwaMem -i stdin | samtools view -Sb - > sample.splitters.unsorted.bam
samtools sort -@ 12 sample.discordants.unsorted.bam > sample.discordants.bam
samtools sort -@ 12 sample.splitters.unsorted.bam > sample.splitters.bam
#lumpyを実行
lumpyexpress -B sample.bam -S sample.splitters.bam -D sample.discordants.bam -o lumpy.vcf
６、CnvkitでSVコール（リンク）。Cnvkitはcondaで導入した。dockerイメージ使うなら"docker pull etal/cnvkit"。
cnvkit.py batch sample.bam --normal \
--targets my_baits.bed --fasta human_g1k_v37_decoy.fa \
--output-reference sample.cnn --output-dir cnvkit/
#.cnsからvcfを出力
cnvkit.py export vcf cnvkit/sample.cns -i "SampleID" -o sample.cnv.vcf
７、CnvnatorでSVコール（リンク）。dokcer imageを使う。メモリ要求量が大きいようなので、メモリlimitは100gとした（*4）。
#オフィシャルらしきものもあるが、timothyjamesbeckerさんのをpullした（リンク）。
docker pull timothyjamesbecker/cnvnator
docker run -m "100g" --rm -itv /Volumes/10TB/test/:/data/ timothyjamesbecker/cnvnator
#作成中
SVコールが終わったら、MetaSV実行する。
方法１
Only merging output of 5 SV detectors without further sof-clips based analysis or local assembly.
SV callerの結果のマージだけ実行。soft clipやローカルアセンブリを介した解析は行わない。
run_metasv.py --reference human_g1k_v37_decoy.fa
--lumpy_vcf lumpy_tumor_normal.vcf
--cnvkit_vcf Sample.cnv.vcf
--breakdancer_native breakdancer_output
--manta_vcf manta.vcf
--pindel_vcf pindel.indel.vcf
--outdir metasv_out --sample HG005 --disable_assembly --filter_gaps --keep_standard_contigs --num_threads 12 --workdir work
出力ディレクトリ。統合されフィルタリングされたVCFファイルが出力される。
方法２
Complete run of MetaSV using all 4 SV detectors, soft-clips based analysis to enhance insertion detection, and local assembly to improve breakpoint resolution.
全解析。
run_metasv.py --reference human_g1k_v37_decoy.fa
--lumpy_vcf lumpy_tumor_normal.vcf
--cnvkit_vcf Sample.cnv.vcf
--breakdancer_native breakdancer_output
--manta_vcf manta.vcf
--pindel_vcf pindel.indel.vcf
--spades SPAdes/spades.py --age AGE/age_align
--min_support_ins 2 --max_ins_intervals 500000 --isize_mean 500
--isize_sd 150 --num_threads 12 --workdir work --outdir out
--sample HG005 --bam alignments.bam
mac osでは、アセンブリしたcontigのindex作成あたりでエラーを起こした。やはりlinux機でテストした方が良かったかもしれない。
他の方法は著者らが準備した説明ページで確認して下さい。BreakSeq2のコマンドはわからなかったので載せてません。GATKについては、バージョンが異なると内容が変わり、余計な混乱を招くので載せてません（GATK３と最新のGATK４の違いとか）。
より新しい手法も発表されてきています。また紹介します。
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5859555/
引用
MetaSV: an accurate and integrative structural-variant caller for next generation sequencing
Marghoob Mohiyuddin, John C. Mu, Jian Li, Narges Bani Asadi, Mark B. Gerstein,Alexej Abyzov, Wing H. Wong, and Hugo Y.K. Lam
Bioinformatics. 2015 Aug 15; 31(16): 2741–2744.
*1 elprepでの簡単フィルタイングを考えたが、elprepの一部タスクは全データをメモリ内に収納して動作するためメモリ要求量が高く、chrごとにsplitしてもメモリが足りないのでやめた（ヒト50xWGSでも、おそらくは384GBくらいメモリが確保されていれば、splitで動く）
elprep split
↓
elprep filter --mark-duplicates --remove-duplicates --filter-mapping-quality 1 --clean-sam
*2 mac osでランしたかったので、処理速度の低下には目をつぶり、dokcer imagesを使っている。
*3 使用スレッド数を増やす場合、利用可能なメモリー量などに応じて慎重に行って下さい。
*4 実メモリに応じて変えて下さい。
*5 新しくても動くかもしれない。未検討。
*6 メモリ追加。指定値xスレッド数、だけメモリを食います。
*7 sort時にtoo many filesのエラーが出たら、このページを参考にopen filesの上限を上げて下さい（現在の値を確認するには"ulimit -n"）。
