macでインフォマティクス

macでインフォマティクス

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

SVシミュレーションや、SVのマージ、レポート生成ができる SURVIVOR

structural variations (SV) evaluation tool simulator Pacbio Nanopore long read テスト失敗

2019 10/29 インストール修正

2020 1/6 追記

 

　一塩基多型（SNP）、小さな挿入 - 欠失事象（indels）、トランスポゾン挿入および大きな構造変化（SV）を含む、様々な遺伝的変化が生物種に影響し得る。欠失、重複、挿入、逆位および転座を含むSVは、タイピングするのが最も困難であり、結果として最もよく説明されていない。

　それにもかかわらず、SVは様々な生物学的過程に強い影響を与えることは明らかである。特に、コピー数変化（CNV）は、農業上重要な形質、様々なヒト疾患、微生物、植物および動物の定量的形質に影響を及ぼす（論文より ref.1,2,3,4,5)。逆位は生殖的隔離（wiki）に影響することが知られている（ref.6,7,8,9,10,11,12,13）および他の進化的過程、例えば組換え（ref.8）や種間のハイブリダイゼーション（ref.14）。

我々（著者ら）は、最近、分裂酵母Schizosaccharomyces pombeをpopulation genomicsおよびquantitative trait analysis（QTL wiki）のモデルとして開発し始めた（ref.6,7,16,17,18）。このモデル生物は、小さく、しっかりアノテーションされた一倍体ゲノム（ref.19）を持ち、遺伝子操作およびハイスループット表現型解析のための豊富なツール（ref.20）があり、ゲノム規模および遺伝子中心のデータ（ref.21,22,23）の大量のリソースを持つという利点を兼ね備えている。

　分裂酵母のこれまでの分析では、自然発生的および人工的な逆位および相互転座の両方を記述し始めている（ref.6,7,18）。 SVの証拠とそのモデル種の影響を考慮すると、我々（著者ら）はSVの体系的調査が生物学的影響の理解を促進するであろうことを認識した。ここでは、161の分裂酵母ゲノムの最新の利用可能性と、定量的形質および生殖分離に関する広範なデータを利用して、S.pombeにおけるSVの性質および効果を説明する。

著者らは、SVが様々な量的形質および本質的な（もともと備わっている:内因性の）生殖隔離に強い影響を及ぼすことを示す。それらは、形質変動の平均11％に寄与する（より豊富なSNPが平均24％寄与する）が、CNVが最大の影響をもたらす。著者らは、CNVがクローン集団内で一過性であり、しばしばSNPによってうまくタグ付けされていないことを示す。また、再構成（逆位および転座）は生殖隔離に寄与するが、CNVはそうしないことを示す。

Methodより

　著者らは、ショートリードシーケンスデータのSVを評価するためいくつかのモジュールを含むSURVIVORツールキットを開発した。第一のモジュールは、リファレンスゲノムファイル（fasta）と各SVの数とサイズの範囲（挿入、削除、複製、逆位および転座）を与えられたSVをシミュレートする。リファレンスゲノムを読み取った後、SURVIVORは、提供されたパラメータに従ってSVの位置およびサイズをランダムに選択する。続いて、SURVIVORはリファレンスゲノムをそれに応じて変更してプリントする。さらに、SURVIVORは、シミュレートされたSVの位置を報告する拡張ベッドファイルを提供する。

　第2のモジュールは、Variant Call Format（VCF）ファイルおよび任意の既知SVリストに基づいてSVコールを評価する。 SVが正しいと同定されるのは、（i）同じSVサブタイプ（例えば、欠失）である、 （ii）同じ染色体上に報告される、（iii）シミュレートされ同定されたSVの開始および停止座標が1kb以内である（ユーザ定義可能）、を全て満たす場合である。

　SURVIVORの 第3モジュールは、3つのSVコーラーのVCFファイルからのコールをフィルタリングして結合するために使用された。著者らの場合、これらのファイルはDELLY、LUMPY、Pindelの結果であった。このモジュールには、メソッド固有の出力フォーマットをVCFフォーマットに変換するメソッドが含まれている。 SVは、3つのVCFファイルのいずれかのみに固有であれば除外された。 2つだった場合、同じ染色体上に存在し、その開始および停止座標は1kb以内であり、それらは同じ型の場合に重複すると定義された。最後に、SURVIVORはフィルタリングされたコールを１つのVCFファイル統合した。 SURVIVORはgithub.com/fritzsedlazeck/SURVIVORから入手できる。

 

SURVIVORは以下のような機能を持つ（Githubより）。

  • Simulate SVs and evaluate existing callers
  • Merge and compare SVs within a sample and among populations/samples.
  • Convert different formats to vcf files
  • Summarize the results within vcf files or results from SURVIVOR.

そのほか、リアルデータ（PacbioやONTのデータ）からエラープロファイルを取得し、ロングリードをシミュレーションする機能をもつ。

 

wiki

https://github.com/fritzsedlazeck/SURVIVOR/wiki

 

インストール

mac os10.13、python3.4.0環境でテストした。

本体　Github

git clone https://github.com/fritzsedlazeck/SURVIVOR.git
cd SURVIVOR/Debug
make

 Anaconda環境ならcondaでも導入可能。linuxのみ

#bioconda (link)
conda install -c bioconda survivor

> ./SURVIVOR

# SURVIVOR

Program: SURVIVOR (Tools for Structural Variations in the VCF format)

Version: 1.0.6

 

Usage: SURVIVOR <command> [options]

 

Commands:

-- Simulation/ Evaluation

simSV Simulates SVs and SNPs on a reference genome.

scanreads Obtain error profiles form mapped reads for simulation.

simreads Simulates long reads (Pacio or ONT).

eval Evaluates a VCF file after SV calling over simulated data.

 

-- Comparison/filtering

merge Compare or merge VCF files to generate a consensus or multi sample VCF files.

genComp Generates a pairwise comparison matrix based on any multi sample VCF file

filter Filter a vcf file based on size and/or regions to ignore

stats Report multipe stats over a VCF file

compMUMMer Annotates a VCF file with the breakpoints found with MUMMer (Show-diff).

 

-- Conversion

bincov Bins coverage vector to a bed file to filter SVs in low MQ regions

vcftobed Converts a VCF file to a bed file

bedtovcf Converts a bed file to a VCF file 

smaptovcf Converts the smap file to a VCF file (beta version)

bedpetovcf Converts a bedpe file ot a VCF file (beta version)

hapcuttovcf Converts the Hapcut2 final file to a VCF file using the original SNP file provided to Hapcut2

convertAssemblytics Converts Assemblytics to a VCF file

 

ラン

推奨パラメータ

https://github.com/fritzsedlazeck/SURVIVOR/wiki/Methods-and-Parameter

1、configファイルの作成。

SURVIVOR simSV parameter_file

parameter_fileファイルができる。中身を確認する。

> cat parameter_file

$ cat parameter_file 

PARAMETER FILE: DO JUST MODIFY THE VALUES AND KEEP THE SPACES!

DUPLICATION_minimum_length: 100

DUPLICATION_maximum_length: 10000

DUPLICATION_number: 3

INDEL_minimum_length: 20

INDEL_maximum_length: 500

INDEL_number: 1

TRANSLOCATION_minimum_length: 1000

TRANSLOCATION_maximum_length: 3000

TRANSLOCATION_number: 2

INVERSION_minimum_length: 600

INVERSION_maximum_length: 800

INVERSION_number: 4

INV_del_minimum_length: 600

INV_del_maximum_length: 800

INV_del_number: 2

INV_dup_minimum_length: 600

INV_dup_maximum_length: 800

INV_dup_number: 2

この定義に従いSVが導入される。変更が必要ならエディタで開いて修正、保存する（順番やスペースをいじってはならない）。

 

２、SVを含むゲノムfastaの作成。先ほど作ったconfigファイルを指定する。SNPs頻度(0-1の間）は0.1にした。

SURVIVOR simSV ref.fa parameter_file 0.1 0 simulated

上のコマンドは以下の順番に記載している。

1: Reference fasta file

2: Parameter file

3: SNP mutations frequency (0-1)

4: 0= simulated reads; 1= real reads

5: output prefix

4番目の数値(0 or 1)はシミュレーションのモードを表す。0はシミュレーションリードを発生させ、それをリファレンス（WTゲノム）にマッピングするモード。1はリアルシーケンスデータをSVを導入したリファレンス（MTゲノム）にマッピングするモード。いずれにしてもSVは検出されるが、挿入と欠失のコールなどは逆転する。

 

human chr20を使いテストしたが、simSVのランが2日経っても終わらなかった。CPU usageは100%のままになっているので、何か計算していると思われるが、logやエラーを立てるフラグがないためどうなっているかわからない。コマンドを勘違いしているのかもしれない。

 

３、コールされたSVの評価。

SURVIVOR eval caller.vcf simulated.bed 10 eval_res

 ２のランが終われば追記します。

 

 

他にも以下のようなコマンドがある。

-- Comparison/filtering --

  • merge   Compare or merge VCF files to generate a consensus or multi sample vcf files.
  • filter   Filter a vcf file based on size and/or regions to ignore
  • stats   Report multipe stats over a VCF file
  • compMUMMer   Annotates a VCF file with the breakpoints found with MUMMer (Show-diff).

-- Conversion --

bincov   Bins coverage vector to a bed file to filter SVs in low MQ regions

  • vcftobed   Converts a VCF file to a bed file
  • bedtovcf   Converts a bed file to a VCF file 
  • smaptovcf   Converts the smap file to a VCF file (beta version)
  • bedpetovcf   Converts a bedpe file ot a VCF file (beta version)
  • hapcuttovcf   Converts the Hapcut2 final file to a VCF file using the original SNP file provided to Hapcut2
  • convertAssemblytics   Converts Assemblytics to a VCF file

 

ここではリアルデータのエラープロファイルを取り込んだシミュレーションロングリードを発生させる流れを記載する。

1、リアルデータを読み込みエラープロファイルファイルを作成。

SURVIVOR scanreadsにsamtools viewをパイプしてリードを読み込む。ターゲットは最小1000bpとする（ロングリードの場合、短かすぎるリードには異なるbiasが発生しているかもしれない。）。時間節約のため、viewに"-s 0.1"をつけて例えば10%だけサンプリングする。

samtools view -s 0.1 your_file.bam | ./SURVIVOR scanreads 1000 error_file

error.txtができる。ちなみにONT由来のbamでは動作したが、ショートリードのbamを使うと、パラメータを変えてもsegmentation faultを起こした。

 

２、エラープロファイルに従いリードを発生させる。

指定するパラメータはカバレッジだけになる。

SURVIVOR simreads ref.fa error_file 100 output.fa

ロングリードのFASTAファイルが出力される（fastqは出力不可）。

 

追記

SV callのマージ と視覚化（チュートリアルに載っています）

ls *vcf > sample_files
#1kb以内　のコールを統合、two SV caller以上がサポート、10bp以下除去
SURVIVOR merge sample_files 1000 2 10 0 0 10 merged.vcf

#視覚化
SURVIVOR genComp merged.vcf merged.mat.txt
#Rにて
t=read.table("merged.mat.txt",header=F)
dst <- data.matrix(t(t))
image(1:nrow(dst), 1:ncol(dst),log(dst), col=rev(heat.colors(10)),axes=F, xlab="", ylab="")
grid(col='black',nx=nrow(dst),ny=ncol(dst),lty=1)

f:id:kazumaxneo:20200106212437p:plain



 

引用

Transient structural variations have strong effects on quantitative traits and reproductive isolation in fission yeast.
Jeffares, Daniel C; Jolly, Clemency; Hoti, Mimoza; Speed, Doug; Shaw, Liam; Rallis, Charalampos; Balloux, Francois; Dessimoz, Christophe; Bähler, Jürg; Sedlazeck, Fritz J.
Nature communications, Vol. 8, 14061, 24.01.2017, p. 1-11. DOI:10.1038/NCOMMS14061