macでインフォマティクス

macでインフォマティクス

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

複数のSV検出結果を統合し、精度の高いSVコールを行う MetaSV

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/

f:id:kazumaxneo:20180731173653j:plain

 

MetaSV: An Accurate and Integrative Structural-Variant Caller for Next Generation Sequencing

f:id:kazumaxneo:20180731174637j:plain

 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実行前

f:id:kazumaxneo:20180731145324j:plain

実行後。out/にvariants.vcfができている。

f:id:kazumaxneo:20180731145348j:plain

variants.vcfの中身。

f:id:kazumaxneo:20180731145507j:plain

動作している。

 

ラン

いくつかのSVコーラーは以前まとめて紹介してます(リンク)。ここでは各ツールの SVコールコマンドと、それを入力としてmetaSVを実行する手順を確認する。

 

1、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を除き、リアライメントまで実行したいときは、こちらを参照。

 

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

 

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

 

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

 

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

 

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

 

7、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実行する。

方法1

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ファイルが出力される。

f:id:kazumaxneo:20180803163629j:plain

 

方法2

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については、バージョンが異なると内容が変わり、余計な混乱を招くので載せてません(GATK3と最新のGATK4の違いとか)。

MetaSV 

 

より新しい手法も発表されてきています。また紹介します。

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")。