macでインフォマティクス

macでインフォマティクス

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

SVtools

 

 近年の全ゲノムシークエンシング(WGS)の劇的なコスト削減により、数万から数十万のディープシーケンシングされた(> 20倍)個体の包括的な形質関連の解析を行うことを目的とする大規模なヒト遺伝学研究が進行中である。その中で最も重要なものは、NHGRI’s Centers for Common Disease Genomics (CCDG) やNHLBI’s Trans-Omics for Precision Medicine (TOPMed)で、これまでに150,000を超えるディープWGSデータセットが生成されている。さらに、ゲノム生物学への洞察のためマイニングし大きなゲノム変異マップ作成を模索している進行中のゲノムアグリゲーションエフォートは、個人のゲノムおよびまれな疾患研究を解釈するのを助けることができる。さらにこれらの努力は、世界中の他の多くの人々と共に、データ中心のヒト遺伝学研究の新時代を迎えるだろう。
 WGSの重要なプロミスは、あらゆる形態のゲノム変異を評価する可能性である。しかし、多くのグループによるかなりの努力とcreativity(特に1000 Genomes Project、1KGP(Mills、et al、2011; Sudmant、et al、2015))にもかかわらず、WGSから高品質のSVマップをアセンブリすることは非常に難しく、特に何千もの個人の大規模コホートデータではそうである。初期の障害は、ブレイクポイント位置を推定しコピー数を推定するため使用されるイルミナのショートリードWGSデータではアライメント信号の完全性には根本的な制限があることであり、そのためSV検出は小規模の研究でも難しい。スプリットリード(SR)、クリップリード(CR)、リードペア(RP)、リードデプス(RD)を含むこれらのアライメント信号は、シーケンスやアライメントのアーティファクトと区別するのが難しく、統合するのが困難である。最高のパフォーマンスを発揮するツールでも、一般的に感度が低く、誤検出率(FDR)が高く、計算コストが高いという問題がある。
 第二の問題は、現在のWGSベースのマルチサンプルのバリアント検出アプローチでは、各推定バリアント部位において、各サンプルについて生の(またはほぼ生の)アラインメントデータの“joint” analysisを必要とすることである。ただし、メモリと計算の制限により、native joint callingのアルゴリズムでは数百ゲノムもあると上手く動作しない。 GATKおよびVTのようなSNV / indel検出ツールは、中間ファイル(例えば、gVCF)および並列の「scatter-gather」コンピューティングスキームを使用することによって、大きなコホートで変異検出信号を取り除き結合するための分散ワークフローを実装した。当然の目標は、SVのための同様のアプローチを開発することである。 

 しかし、SVの場合、問題は異なり、間違いなくはるかに困難である。ツールは、より高いエラー率を許容し、分類が困難な場合がある。balanced、complex and/or repetitive variantsを含むさまざまなバリアントサイズおよびアーキテクチャに対応する必要がある。単一座標にマッピングするSNVとは異なり、SVブレイクポイントは、新規のDNAジャンクションを記述する不連続で潜在的に離れたstrand指向の基準ゲノム座標(正式には“break-ends”)ペアによって定義されるという事実により、並列化スキームは複雑である。gVCFに類似した中間データ構造は、それぞれが異なる解像度およびアーチファクトモダリティを有する少なくとも4つの異なるアラインメント信号(SR、CR、RP、RD)からの情報をカプセル化しなければならないので設計が困難である。 SVブレイクポイントマッピングの分解能は一般にサンプルごとに不正確であり(平均10〜100 bp)、シーケンスやアライメントの効果はバリアントクラス、サンプル、バッ​​チによって大きく異なるため、クロスサンプルのマージ方式は位置の不確定性に対してロバストでなければならない。明らかに新しいアプローチが必要である。
 もちろん、サンプルの集まりにわたって空間的に不正確なSV /CNVコールを組み合わせるというタスクは、従来のマイクロアレイにおけるアドホックな方法を通して効果的に扱われてきた古い問題である(Conrad et al、2009; Redon et al、 2006; Wellcome Trust Case Control et al、2010)およびWGSに基づく(Mills et al、2011; Sudmant et al、2015)。しかしながら、アレイベースの方法は、balanced SV、あるいはディープWGSの分解能、複雑さおよび規模の増大に容易には及ばない。 1KGPは、複数のアルゴリズムおよびプラットフォームからの結果をマージするために巧妙なアプローチを採用した(Mills et al、2011;Sudmant et al、2015)が、記念すべきエフォートであり、その中の方法は日常的な使用には非実用的である。 GenomeSTRiPには、サンプルの集団からSVを検出するための2つの公開されたワークフローがあるが、両方とも制限がある。初期のバージョンでは欠失に焦点を合わせ、RPベースの検出とRDジェノタイピングを連続的に組み合わせている(Handsaker、et al、2011)。第2のRDベースのCNVパイプラインは、計算量が多く、低解像度(>1kb)であり、中程度のサンプルサイズ(<1000)に限定されている(Handsaker et al、2015)。著者らの知る限りでは、ここに示すように、何万ものディープWGSデータセットにおけるマルチプルアライメント信号の共同分析から高解像度SVコールセットを体系的に組み立てるための公に利用可能なツールまたは再現可能なワークフローは存在しない。

 サンプルごとのバリアント発見、解像度を意識したクロスサンプルマージ、ブレイクポイントジェノタイピング、コピー数アノテーション、バリアント分類、およびコールセット改良を組み合わせた、大規模なSVコールセット生成のためのソフトウェアツールキットと分散ワークフローsvtoolsを開発した(Preprint 図1)。 

 svtoolsはfamily解析など複数ゲノムの構造変化検出結果を読み込み、特異的なSV探索する作業を支援するために設計された一連のユーティリティ。数千から数万のゲノムにわたるspeedseq svからのコールを効率的にマージして遺伝子型を決定するように設計されている。

 

f:id:kazumaxneo:20190407210008p:plain

 Figure 1. The svtools pipeline.  Preprintより転載。

 

チュートリアル

https://github.com/hall-lab/svtools/blob/master/Tutorial.md

 

twitter

 

インストール

mac os 10.12のanaconda2-4.2.0でテストした。

依存

  • A Linux-like environment with bash, awk, and sort
  • A working installation of cnvnator-multi (a component of speedseq sv)
  • Python2.7Numpy
  • Scipy
  • Pandas
  • Statsmodels
  • Pysam (≥0.8.1)

本体 Github

pip、またはcondaで依存も含め導入できる。

#Anaconda環境ならcondaで導入できる。
conda install -c bioconda svtools

#Anaconda環境ではないならpipを使う。
pip install svtools

$ svtools -h

usage: svtools [-h] [--version] [--support] subcommand ...

 

Comprehensive utilities to explore structural variation in genomes

 

optional arguments:

  -h, --help    show this help message and exit

  --version     show program's version number and exit

  --support     information on obtaining support

 

  subcommand    description

    lsort       sort N LUMPY VCF files into a single file

    lmerge      merge LUMPY calls inside a single file from svtools lsort

    vcfpaste    paste VCFs from multiple samples

    copynumber  add copynumber information using cnvnator-multi

    genotype    compute genotype of structural variants based on breakpoint

                depth

    afreq       add allele frequency information to a VCF file

    bedpetobed12

                convert a BEDPE file to BED12 format for viewing in IGV or the

                UCSC browser

    bedpetovcf  convert a BEDPE file to VCF

    vcftobedpe  convert a VCF file to a BEDPE file

    vcfsort     sort a VCF file

    bedpesort   sort a BEDPE file

    prune       cluster and prune a BEDPE file by position based on allele

                frequency

    varlookup   look for variants common between two BEDPE files

    classify    reclassify DEL and DUP based on read depth information

 

 

テストラン(large data-set)

1、データの準備

NA12878の家系からNA12877、NA12878、NA12889の bamをダウンロードする。3つで300 GB以上ある(*1)(このデータセットになる。fastqもある)。

wget ftp://ftp.sra.ebi.ac.uk/vol1/run/ERR194/ERR194146/NA12877_S1.bam
wget ftp://ftp.sra.ebi.ac.uk/vol1/ERA172/ERA172924/bam/NA12878_S1.bam
wget ftp://ftp.sra.ebi.ac.uk/vol1/ERA172/ERA172924/bam/NA12879_S1.bam

 リファレンスのGRCh37

wget ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/reference/human_g1k_v37.fasta.gz

 

2、speedseqのrealignでbamを変換し、svコマンドでSV検出。

#realign
speedseq realign -o NA12877 human_g1k_v37.fasta NA12877_S1.bam
speedseq realign -o NA12878 human_g1k_v37.fasta NA12878_S1.bam
speedseq realign -o NA12879 human_g1k_v37.fasta NA12879_S1.bam

#sv call
speedseq sv -o NA12877 -R human_g1k_v37.fasta.gz -B NA12877_realign.bam -D NA12877.discordants.bam -S NA12877.splitters.bam -v -d -P -g -k
speedseq sv -o NA12878 -R human_g1k_v37.fasta.gz -B NA12878_realign.bam -D NA12878.discordants.bam -S NA12878.splitters.bam -v -d -P -g -k
speedseq sv -o NA12879 -R human_g1k_v37.fasta.gz -B NA12879_realign.bam -D NA12879.discordants.bam -S NA12879.splitters.bam -v -d -P -g -k

 

3、svtools lsortを使ってvcfを結合しソートする。

svtools lsort NA12877.sv.vcf.gz NA12878.sv.vcf.gz NA12879.sv.vcf.gz | bgzip -c > sorted.vcf.gz

 

 

より大規模なデータセットの例や、他のコマンドの使い方はgithubを参照してください。

引用

svtools: population-scale analysis of structural variation
David E. Larson, Haley J. Abel, Colby Chiang, Abhijit Badve, Indraniel Das, James M. Eldred, Ryan M. Layer, Ira M. Hall

bioRxiv preprint first posted online Dec. 13, 2018

 

関連


*1

テスト時はNA12879がダウンロードできなかったので、NA12877とNA12878だけで実行した。