macでインフォマティクス

読者です 読者をやめる 読者になる 読者になる

indelの検出ツール

バクテリアでパフォーマンス比較したペーパーが出ている。

実際に導入して、パフォーマンスを比較してみる。

 

 

はじめにbrew tap でsiience系のコマンドをオフィシャルコマンドとしてインストールできるようにしておく。

brew tap homebrew/science

 

ほとんどのスクリプトはインプットとしてソート済みのbamファイルを使う。BWA を使いリファレンスにアライメントしておく。BWAアルゴリズムは古いbacktrack、alnなどより、1Mbpのロングリードまで対応可能なmemが良いと思われる。下がソート済みbamとそのindexを作るまでの流れ。

bwa index -a is ref.fa 
bwa mem -R "@RG ID:X LB:Y SM:Z PL:ILLUMINA" -t 12 ref.fa R1.fq R2.fq > R1R2.sam
samtools view -@ 20 -bS R1R2.sam > R1R2.bam
samtools sort -@ 20 R1R2.bam > R1R2_sorted.bam
samtools index R1R2_sorted.bam
samtools faidx $f

 R1.fqとR2.fqがfastqファイル。ref.faはリファレンス。

リファレンスのindexも必要とするものが多い。samtools faidxでfasta.faiを作っておく。

samtools faidx $f

 これらの結果を1つのフォルダにまとめておき、各スクリプト実行時にそのフォルダを参照するようにして実行する。 

 

ショートリードから構造変化を予測する手法は大きく4つに分けることができる。

1、Alignment based 

Dindel

Stampy

SAMtools

GATK

Varscan

 

2、Split-read mapping

Pindel

SV-M

 

3、Paired-end mapping

Hydra

PEMer

Breakdancer

 

4、haplotype-based

GATK-HaplotypeCaller

Platypus

 

 

 

Breakdancer Chen et al. (2009)

discordant paired read approaches

 

macOSX環境でのインストールも可能。ペアードエンドリードの挙動から構造変化を予測する。インストールはめんどくさかった、cpanでたくさんのperlモジュールをインストールする必要があった。

 

 breakdancerはペアードエンド法で一番知られているプログラムなので、原理に触れておく。Pindelのようなsplit-mappingは、片方のリードを正確にアライメントし、もう片方のリードのアライメント状況から構造変化の種類を推測するが、breakdancerのようなread pairの手法はインサートサイズのずれからindelを予測する。つまり、ペアリードのインサートサイズが一般的な正規分布に従うと仮定し、ペア間の距離が分布の外れ値になるような部位にindelがあるだろうと考える手法となる。この時、リードの向き(正常なら--> <--)も考えれば、構造変化の種類についても情報も得ることができる。Breakdancerはこのような異常なペアリード; discordant paired read(DPR)情報に基づき構造変化を予測するツールとなる。実際にはそうゆう異常なペアが1箇所に2組以上集積している部位を信頼値付きでコールする。

 ただしこの手法には大きな問題点がある。それはペア間の距離が大きく変わっていることが前提なので、small indelは検出できないこと。さらにインサートサイズの分布によって検出部位が変化しうること、またindelの配列に関する情報が得られない上に正確なSVのbreakpointもわからないことである。

 そのため2017年現在では、このread pairの手法とsplit-mappingやlocal de novo assemblyを組み合わせて精度を上げたアルゴリズムが色々発表されている。

 

解析手順は、まずbamファイルをサンプリングしてインサートサイズなどを分析したファイル(.cfgファイル)を作り、それからコールを行う。

#cfgファイルを作成する(orted.bamはbwa memで作成した)。

perl bam2cfg.pl R1R2_sorted.bam > output.cfg

cfgファイルの中身は以下のようになっていた

 >cat output.cfg
readgroup:X platform:ILLUMINA map:R1R2_sorted.bam readlen:296.80 lib:Y num:9990 lower:48.86 upper:29146.26 mean:444.35 std:3920.08 SWnormality:minus infinity exe:samtools view

サンプリング自体は最初の数百行だけ行うらしく、すぐに終わる。ただしこの時qualityがpoorなデータだと解析できない (ランした時にpoor quality dataと出る)。 どうしてもランしたければ、オプション"-v"をdefaultの1から2-10まであげると一応ランはできる。

 

#構造変化のコール

 breakdancer-max output.cfg > output.txt

 人のchr20で試すと、解析は10分くらいかかった。500箇所ほどコールされた。数が極端に少ないのは、SNVやshort indelを検出するのが難しいと思われる。

 

 

Pindel Ye et al. (2009)

split-read approaches.

 

ショートリードのペアードエンドデータからdeletion(1bp-10kb)、medium insertion (1-20bp)を予測するツール。リードを分割してアライメントし、そのアライメントのパターンから構造変化を予測している。このようなアライメントをスプリットアライメントと呼び、BWA MEMやRNA-seqで使うtophatなども使っている方式となる。

手法の詳細は論文に書いてあるが、split-read mappingの代表的な論文なのでもう少し手順を説明する。

1、ペアリードの片側だけマップされるペアを探索する。この片側のリードのアライメントが肝なので、ユニークでかつミスマッチがないパーフェクトマッチのアライメントが行われる。

2、相方リードがマップされないのは、この相方リードが構造変化部位を横断するようにマップされているからと考える。この相方リードのアライメント部位を探すため、アライメントできた方のリード(アンカーと呼んでいる)からインサートサイズ平均の2倍の範囲内でマッピングされなかったリードの方がマッチする部位を探す。ただし相方リードはそのままではアライメントできないので、リードを5'側と3'側、真ん中の3ピースに分けてそれぞれでアライメントを実行する(slit alignmnet)。

3、5'側と3'側のアライメント部位にずれがあれば、即ちそのずれのサイズがdeletionのサイズになる(下図a)。5'側と3'側のアライメント部位にズレはないが、真ん中のピースが浮いた状態でマッピングされていれば、その真ん中のピースがinsert配列そのもののはずである。

4、2つ以上のペアで同じ情報が検出された部位をコールする。

f:id:kazumaxneo:20170327095850j:plain

片方のリードを正確にアライメントし、その相方の位置も限定することで偽陽性を抑え且つ解析をスピードアップしている。原理上リードの長さによって検出可能なindelの長さの上限が決まる。初期の30-40bpのリードでは20bpのinsertionが上限ということらしい。

 

以下の図は解析の流れ。

f:id:kazumaxneo:20170509104834j:plain

 

http://gmt.genome.wustl.edu/packages/pindel/install.html にインストールの仕方が書いてある。macOSX環境ではインストールできなかったのでcentOSにインストールした。

git clone git://github.com/genome/pindel.git
cd pindel
./INSTALL /path/to/samtools_FOLDER/

 

デモランは以下の通りに行う。

cd demo
./pindel -f simulated_reference.fa -i simulated_config.txt -o output

ランに必要なものはリファンレス配列(.fasta)とそのindexファイル(.fasta.fai)。またそのリファレンスにmapして作ったbamファイルとそのindexファイル(.bam.bai)。このあたりは他のツールと変わりない。

pindelは情報をファイルに記載したものを作って、それを元にランするタイプのスクリポトになる。サンプルファイルの場合これに相当するのはsimulated_config.txtファイルとなる。ファイル内には以下のことが書かれている。

simulated_sample_1.bam  250     SAMPLE1

simulated_sample_2.bam  250     SAMPLE2

simulated_sample_3.bam  250     SAMPLE3

データが1つなのでこれを以下のように変更した。

R1R2_sorted.bam  600     SAMPLE1

この600はインサートサイズ。ランするには

./pindel -f input.fasta -T 12 -i config.txt -o output

-T: thread

-f: ref.fa

-i: config_file

-o: output

終わると構造変化の種類ごとにアライメントファイルが出力される。結果を見るにはpindel2vcfで変換する。pindel2vcfはビルドしてでできたpindel/に入っている。

pindel2vcf -p pindel_output.txt_D -r ref.fa -R ref -d 20170509 -v pindel_output_D.vcf -e 5
pindel2vcf -p pindel_output.txt_INT_final -r ref.fa -R ref -d 20170509 -v pindel_INT_final.vcf -e 5
pindel2vcf -p pindel_output.txt_INV -r ref.fa -R ref -d 20170509 -v pindel_output_INV.vcf -e 5
pindel2vcf -p pindel_output.txt_SI -r ref.fa -R ref -d 20170509 -v pindel_output_SI.vcf -e 5
pindel2vcf -p pindel_output.txt_TD -r ref.fa -R ref -d 20170509 -v pindel_output_TD.vcf -e 5
pindel2vcf -p pindel_output.txt_BP -r ref.fa -R ref -d 20170509 -v pindel_output_BP.vcf -e 5

-r/--reference The name of the file containing the reference genome: required parameter -R/--reference_name The name and version of the reference genome: required parameter -d/--reference_date The date of the version of the reference genome used: required parameter -p/--pindel_output The name of the pindel output file containing the SVs  

-e/--min_supporting_reads  The minimum number of supporting reads required for an event to be reported (default 1)

出力のvcfフォーマットはvcf higher versionに準拠しているらしい。

-eは必須パラメータではないが、デフォルトの1では高感度すぎて明確な変異なしのシミュレートデータ(50カバレッジ)でも1000近くの偽陽性コールが検出される。indelソフトのパフォーマンス比較の論文で-e 5になっていたのでここではそれに習った。ただし、目的やリードのカバレッジによってチューニングが必要と思われる。

 

 

 

LUMPY layer et al. (2014)

split-read + discordant paired read + CNV based approaches. 

 解析の流れはこのURLに丁寧に書いてある。

 

macOSX環境でのインストールも可能。homebrewでインストールできる。

brew install homebrew/science/lumpy-sv

ソフトクリップのマッピングファイルから構造変化を予測する手法。マッピングソフトは1Mbpまで対応したbwa memなどのマッパーを使う。デプスのないデータなどでもsensitivityが高いらしい。また複数サンプルで差分をとって精度を上げるやり方も搭載している。

f:id:kazumaxneo:20170327095430j:plain

テストランはtestディレクトリにあるtest.shを実行する。

./test.sh

問題なければNA12878、honeybee、riceのindel予測結果がvcf形式で出力される。

 

 解析時は下のような流れでコールする。

bwa index -a bwtsw ref.fa
bwa mem -R "@RG ID:X LB:Y SM:Z PL:ILLUMINA" ref.fa R1.fastq R2.fastq |samtools view -S -b - > sample.bam
samtools view -@ 20 -b -F 1294 sample.bam > sample.discordants.unsorted.bam
samtools view -@ 20 -h sample.bam | ./scripts/extractSplitReads_BwaMem -i stdin | samtools view -Sb - > sample.splitters.unsorted.bam
samtools sort -@ 20 sample.discordants.unsorted.bam > sample.discordants.bam
samtools sort -@ 20 sample.splitters.unsorted.bam > sample.splitters.bam
lumpyexpress -B sample.bam -S sample.splitters.bam -D sample.discordants.bam -o sample.vcf

解析できるまで手間取った。最後の肝心のlumpyexpressがモジュールを読み込めずエラーになる。最終的に強引にrootで進めることでなんとかランできた。

samtoolsのソートでエラーになる人は、>(リダイレクト)とリダイレクトの後ろ側の.bamを決してランする。これはsamtoolsのバージョンによっては>と.bamが入らないバージョンがあるため。

 

 

 

Hydra-sv

discordant paired read approaches. 

 

マニュアル: Google Code Archive - Long-term storage for Google Code Project Hosting.

 

下のような流れでコールする。(https://code.google.com/archive/p/hydra-sv/wikis/TypicalWorkflow.wiki

bwa index -p $f -a is $f
bwa aln -t 12 $f $p > R1.sai
bwa aln -t 12 $f $q > R2.sai
bwa sampe $f R1.sai R2.sai $p $q |samtools view -Sb - > sample.tier1.queryorder.bam
samtools view -uF 2 sample.tier1.queryorder.bam |bamToFastq -bam stdin -fq1 R1_disc1.fq -fq2 R2_disc1.fq
novoindex ssuis.nix $f #Novoalignのindex
novoalign -d ssuis.nix -f R1_disc1.fq R2_disc1.fq -i 500 50 -r Random -o SAM |samtools view -Sb - > sample.tier2.queryorder.bam
samtools view -uF 2 sample.tier2.queryorder.bam |bamToFastq -bam stdin -fq1 R1_disc2.fq -fq2 R2_disc2.fq
novoalign -d ssuis.nix -f R1_disc2.fq R2_disc2.fq -i 500 50 -r Ex 1100 -t 300 -o SAM |samtools view -Sb - > sample.tier3.queryorder.bam 
bamToBed -i sample.tier3.queryorder.bam -tag NM |pairDiscordants.py -i stdin -m hydra -z 800 > sample.disc.bedpe
dedupDiscordants.py -i sample.disc.bedpe -s 3 > sample.disc.deduped.bedpe
hydra -in sample.disc.deduped.bedpe -out sample.breaks -mld 500 -mno 1500

 

ただし結果はひどくしょぼい。進め方かパラメータを間違ってる可能性がある。

 

 

 

 

 

DELLY2 Rausch et al. (2012)

split-read +discordant paired read approaches. 

 

Delly2 can discover,deletions, tandem duplications, inversions and translocations at single-nucleotide resolution in short-read massively parallel sequencing data.

deletionなどを予測する。コピーナンバーが変わるリピートのdeletionや、depthが少ないデータにも強いらしい。 コントロールのデータ(somatic sampleなど)と解析データ(tumor sampleなど)両方を使った差分解析を特徴とするが、単体のデータでも解析は可能。insertionはターゲットではないので注意する。こちらからダウンロードできる。

 

macOSX環境ではインストールできなかったのでcentOSサーバーにインストールした。

基本的にはgitからダウンロードしてmakeするだけ。

git clone --recursive https://github.com/dellytools/delly.git 
cd delly/
make all

 

ランは

./delly call -t DEL -x human.hg19.excl.tsv -o delly.bcf -g ref.fa R1R2_sorted.bam

-xはなくても問題ないが、コールを除外する領域を設定できる(例えばヒトのセントロメアなど)。

コントロールと比較するなら

./delly call -t DEL -x human.hg19.excl.tsv -o delly.bcf -g ref.fa R1R2_sorted.bam Control_R1R2_sorted.bam

Contorl~bamは変異がないと期待される体細胞シーケンスデータ。tumor~bamが変異の可能性がある腫瘍組織のシーケンスデータ。

 

ランが終わるとbinaryのbcfファイルが出力される。結果を見るにはbedtoolsで変換する必要がある。

./bcftools/bcftools view delly.bcf > delly.vcf

この時フィルターコマンドは実行することも可能。

 

5つのコールタイプがあるので逐一実行する。

./delly call -t DEL -o del.bcf -g ref.fa R1R2_sorted.bam &
./delly call -t DUP -o dup.bcf -g ref.fa R1R2_sorted.bam &
./delly call -t INV -o inv.bcf -g ref.fa R1R2_sorted.bam &
./delly call -t TRA -o tra.bcf -g ref.fa R1R2_sorted.bam &
./delly call -t INS -o ins.bcf -g ref.fa R1R2_sorted.bam &
./bcftools/bcftools view del.bcf > deletions.vcf 
./bcftools/bcftools view dup.bcf > duplications.vcf
./bcftools/bcftools view inv.bcf > inversions.vcf
./bcftools/bcftools view tra.bcf > translocations.vcf
./bcftools/bcftools view ins.bcf > insertions.vcf

 

 

Platypus Rimmer et al. (2014)

local assembly based approaches. 

nature geneticsに載った手法。原理の詳細はこちら http://nextgenseek.com/2014/07/platypus-a-new-variant-caller-that-integrates-mapping-assembly-and-haplotype-based-approaches/

 

候補部位をローカルリアライメントしてプロタイプを調べ、さらにlocal de novo assemblyも行いSNV、indel、complex polumorphismsを検出する手法。パフォーマンス比較の論文でも精度が高いと評価されている。

ふつうショートリードのアライメントは各々の塩基が一致するかに焦点を当てて行われる。この方式は計算リソースの面から考えると大量のデータを処理するには向いているが、indel周辺はミスマッチが多発してしまうことになる。そこでindelを考量して計算リソースを潤沢に使った丁寧なギャップアライメントを行う機能がいくつかの手法に実装されている(GATKなど)。ただしギャップアライメントでうまく補正できるのはせいぜい数塩基で、long indelが起きていた場合修正しきれない。Platypusはそのような補正が出来ない領域だけを対象にde novoでアセンブルを行い、long indelの検出を試みる手法である。対象領域を限定することで、解析スピードを大きく縮小していることも特徴の1つで、例えば数メガのバクテリアの解析だと早ければ1分以内でランは終わる。

 

macOSでも動く。インストールにはgccとcythonが必要。

brew install gcc #ない人だけ 
pip install cython #python2環境なので、pip3を使ってはダメ

samtoolsとhtslibも必要。

brew install homebrew/science/htslib #ない人だけ
brew install homebrew/science/samtools #ない人だけ
pip install cython #python2環境なのでpip3を使ってはダメ

GitHub - andyrimmer/Platypus: Platypus Variant Caller

からPlatypusをzipをダウンロードして解凍。中に入って

make

これで完了。

ランに必要なリファレンスのfaiファイルを作成。

samtools faidx ref.fa #faiファイル作成

 

最後にラン。ソートしたbamファイルを使う。

python bin/Platypus.py callVariants --refFile=ref.fa --bamFiles=R1R2_sorted.bam

パスを通しておけばPlatypus.pyだけで動く。解析速度は一番速いかもしれない。3Mゲノムだと10秒以内にrunは終わった。またtrue indelのコール率も高かった。 

 

追記;人のデータでも試してみた。hg19のchr20をリファレンスにしてテストしたところ(gz圧縮シーケンスデータは2GBx2あったが)、10分以内に解析は終了した。bam作るまでの方が時間はかかる。Platypusでのcall部位は75000くらいあった。桁違いに多い...

 

 

fermikit Li (2015)

assembly based approaches. 

 

一度アセンブルを行なって、それからアライメントをとるプログラム。そう聞くと重そうだが、ヒトゲノムデータでもピークメモリー使用量は90GBくらいで、アセンブルも24時間くらいで終わると書いてある。またその先のアライメントは1時間半くらいで終わるとのこと。

macOSでも動く。インストールは

git clone --recursive https://github.com/lh3/fermikit.git
cd fermikit
make

ビルドしてできたfermi.kitフォルダの中に実行に必要なものが入っている。アライメントにbwaを使っている。

実行にはリファレンスのbwa indexが必要なので最初に作っておく。 

bwa index ref.fa

アセンブル条件ファイルの作成

fermi2.pl unitig -s 3m  -t 12 -l 301 -p prefix R1.fastq > prefix.mak

-s: 推定ゲノムサイズ(3mは3 Mbp)

-l: リード長 (Miseq 301bp)

-p 出力されるprefix名

ペアーリードには対応していないみたいで、ペアリード指定の方法が見つからない。

 アセンブル

make -f prefix.mak

structural variationsのコール

fermi.kit/run-calling -t 12 ref.fa prefix.mag.gz | sh

prefix.mag.gzとprefix.flt.vcf.gzというファイルができる。READMEには

`prefix.flt.vcf.gz` for filtered SNPs and short INDELs and `prefix.sv.vcf.gz` for long deletions,

に書いてある通り、prefix.flt.vcf.gzの中にSV予測結果、prefix.sv.vcf.gzの中にlong del予測結果が保存されている。playtus同様、解析時間は非常に短い。playtusやfermikitはPEM法などで予測が難しいsmall indelもコールしてくれるのが1つの利点と言える。

 

 

 

FreeBayes Garrison and Marth (2012)

gapped alignment approaches.  

 

macOSでも動く。小さなindel専用の検出ツール。ベイズ的アプローチでlndelを検出する。 倍数体ゲノムの変異検出にも対応している。ショートリードデータからhaplotypeの変異を検出するために開発された。

brewでインストール可能。

brew install FreeBayes
freebayes -f ref.fa -b R1R2_sorted.bam > out.vcf 

感度は-Cで調整できる。デフォルトの"-C 2"では感度が高すぎて、例えば50リードシーケンスリードがある部位で2リードだけッミスマッチがあると検出される。腫瘍の変異などわずかな変異を検出するには都合良いが、この条件では偽陽性も大量に出てくる。クローナルなサンプルなどで偽陽性を減らすには-Cを5ぐらいまであげるといい。 ただしはプロタイプの解析など感度が必要な解析では注意して進める。

freebayes -C 5 -f ref.fa -b R1R2_sorted.bam > out.vcf 

人のchr20で検証すると、390000万箇所コールされた。解析には30分程度しかかからなかった。早い方だと思われる。

 

 

 

GATK

split-read approaches. 

多分一番有名な有名なツール。 macOSでも動く。https://software.broadinstitute.org/gatk/documentation/topic?name=methods に詳細が記載されている。

 

UnifiedGenotyperのラン

gatk -T UnifiedGenotyper -R ref.fa -I  R1R2_sorted.bam -o UnifiedGenotyper.txt 

 HaplotypeCallerのラン

gatk -T HaplotypeCaller -R ref.fa -I R1R2_sorted.bam -o HaplotypeCaller.txt

 

 

SvABA Wala et al. (2017)

local assembly approaches. 

 

gitのリンク。

macでもインストール可能とあったが、ビルド中にエラーが出たのでcentOSにインストールした。

git clone --recursive https://github.com/walaj/svaba 
cd svaba
./configure
make
make install

ビルドが終わると、src/svava/に本体がある。

 ランは

svaba run -p 12 -t R1R2_sorted.bam -G re.fa

-t: bamファイル。

-G: リファレンス。

~indel.vcfが出力ファイル。 

-p: thread数

 

bedファイルを指定することで解析領域を限定することも可能。マニュアルには300bpまでのinsertionをコールできるとあるが、検証したところ確かに300bpまでの長さのInsertionを検出できていた。

 

 

 

VarScan Koboldt et al. (2009)

gapped alignment approaches. 

 

BLAT, Newbler, cross_match, Bowtie and Novoalignなどの多くのアライナーをサポートしている。

 

macOSでも動く。brewでインストール可能。

brew install VarScan

Varscanはpileupデータを使うため、初めにソート済みbamファイルをpileupする必要がある。

samtools mpileup -f re.fa R1R2_sorted.bam > mpileup

SNVをコール

VarScan pileup2snp mpileup > output.txt

indelをコール

VarScan pileup2indel mpileup > output.txt

オプションとして以下のようなものがある。

--min-coverage Minimum read depth at a position to make a call [8]

--min-reads2 Minimum supporting reads at a position to call variants [2]

--min-avg-qual Minimum base quality at a position to count a read [15]

--min-var-freq Minimum variant allele frequency threshold [0.01]

--min-freq-for-hom Minimum frequency to call homozygote [0.75]

--p-value Default p-value threshold for calling variants [99e-02]

--variants Report only variant (SNP/indel) positions [0]

 

pileupしたファイルは大きいので、パイプを使い1ライナースクリプトで進めるやり方がHPに紹介されている。例えばindelをコールするなら

samtools mpileup -f ref.fa R1R2_sorted.bam | VarScan pileup2indel > output.txt 

 

 

Breseq Deatherage et al. (2014)

split-read approaches.

 マニュアルリンク

 

macOSでも動く。brewでインストール可能。

brew install Breseq

indelコールは

breseq -r ref.fa -j 12 R1.fq R2.fq

-r; リファレンスを指定。fast.gbk (genebank形式)、GFF3に対応。

-j: スレッド数

-0: アウトプット

遺伝子のアノテーションがついたgff3かgenebankファイルがあるなら、検出結果を遺伝子の部位情報付きで出してくれるので非常に便利。genbankファイルは、NCBIftpサイトか、もしくは該当ページの情報を配列付きで表示させて全テキストコピーしてもいい。gff3ファイルはUCSCなどからダウンロードできる。genebankリファンレスでのランは

breseq -r ref.gbk -j 12 R1.fq R2.fq

感度調整のオプションはないが、-pをつけると多型検出ができる。helpのメッセージを載せておく。

  -p,--polymorphism-prediction     The sample is not clonal. Predict polymorphic (mixed) mutations. Setting this flag changes from CONSENSUS MODE (the default) to POLYMORPHISM MODE

 

 

解析が終わるとhtmlのまとめファイルができるのが特徴。内部でRなどの描画ライブラリを動かしているぽい。/outputの中のsummary.htmlを開くと、

f:id:kazumaxneo:20170512131212j:plain

上のようにアミノ酸置換まで含めて報告されている。 

変異部位をクリックするとアライメント結果も表示される。すごい。

f:id:kazumaxneo:20170511134951j:plain

 

coverage

f:id:kazumaxneo:20170426164310j:plain

deletion予測時は下のようにリードのデプスを見ることもできる。f:id:kazumaxneo:20170513170011j:plain

 

 

 indelの抜け漏れはこの中で一番少なかった。インフォマティクスのツールに不慣れなユーザーにも使いやすいと思われる。ただし、動作が複雑なためか他のソフトの安定性かわからないが、時折解析中にエラーを起こす。そうゆう時は大体パラレルにBreseqを動作させていた時だったりする。複数走らせないほうがいいみたい。

 

 

Scalpel Narzisi et al. (2014)

assembly based approaches. 

このスクリプトは人をターゲットにしている。腫瘍サンプル特異的な変異の差分検出もできる。オプションにもdad(父)とかmam(母)とかある。人のサンプルを使うことがあれば試してみたい。一応ランの流れを書いておく。

 

解析にはbedフォーマットのポジション情報が必要(リンク)。 

bcftoolsを使ってbamファイルから変換することができる(SEQanswers)。bamファイルはあらかじめソートしてindexをつけておく。

bamToBed -i R1R2_sorted.bam > R1R2_sorted.bed 

Runにはデータベースファイルが必要。

scalpel-discovery --single --bam R1R2_sorted.bam --bed R1R2_sorted.bed --ref ref.fa

上のコマンドを実行するとoutdir/フォルダに色々なファイルができる。その中にvariants.db.dirというファイルができている。これがデータベースファイル。データベースファイルができないエラーもあるらしい

indelの解析

scalpel-export --single --db outdir/variants.db.dir --bed R1R2_sorted.bed --ref ref.fa 

人データでテストしたところ、エラーなく最後まで走ったが、コールがゼロになる。

現在のところ未解決

 

 

Breakseek Zhao et al. (2015)

breakpoint + discordant paired read approaches. 

ダウンロードリンク

 

pythonで動くため、ソースコードのビルドは必要ないがpython の2.7が必要。

またランにはsamファイルが必要。推奨される作り方がマニュアルに書いてあるのでそれに従う。bwa backtrackを使う。

bwa index -a bwtsw ref.fa
bwa aln ref.fa read1.fq > read1.sai 
bwa aln ref.fa read2.fq > read2.sai
bwa sampe ref.fa read1.sai read2.sai read1.fq read2.fq > aln_pe.sam

テスト用のsamとfaが用意されている。下記URLからダウンロードする。

https://sourceforge.net/projects/breakseek/files/toy_testset/

以下はそのファイルで進行している。 

 

bwa alnでsamを作ったら、sam_prep.pyでクロモソームごとに分割する。

python sam_prep.py -f test.sam

終わるとfasta内のクロモソームごとにディレクトリが作られ、中にsamができる。このsamファイルを使ってランする。

python breakseek.py -f testS.sam -m 600 -s 100 -c 50 -q 250 -n 12 -r test.fa

-m: インサートサイズ

-s: インサートサイズのSD。正規分布ならインサートサイズの10-20%くらいらしい。

-c: リードデプス

-q: リード長

以下のメッセージ

gathering statistics on the dataset...

が出てから、解析途中は何もlogが出てこない。初めはランに失敗したと思っていたが、強制終了せず気長に待ったところ、1時間近くかかった。3Mのバクテリアゲノムを使っただけでこれは長い...。プログラムの最適化の度合いが低いのかもしれない。

 

 samファイルはbwa alnでなくてもOK。例えばbwaの別のアルゴリズムbwa memを使うなら

bwa index -a is test.fa
bwa mem -R "@RG ID:X LB:Y SM:Z PL:ILLUMINA" -t 12 test.fa R1.fq R2.fq > R1R2.sam 

is: bwaのindexアルゴリズム。"bwasw"か"is"。

-t: スレッド数

R1.fq: fastqファイル

R2.fq: fastqファイル

 

 

 

Scanindel Yang et al. (2015)

gapped alignment + split reads and + assembly based approaches. 

performance比較の論文で高評価されているツール。

 

動作に必要なものがgitのダウンロードページに記載されている。

Softwares and Python packages:

 

上記の大半はbrewとpipコマンドですぐ導入できたが、Inchwormはmacで動作しなかった。最終的にubuntu 14.04にインストールした。まずbrewで導入できるコマンドでないものがあれば入れておく。

brew listでインストール済みリストのチェック。

brew list

 下のコマンドのうちlistででないものがあれば導入。

brew install bwa 
brew install samtools
brew install bedtools
brew install blat
brew install freebayes

pysam、PyVCF、Biopython、SciPy、NumPyもpythonのpipコマンドで導入できる。PyVCFはpythonのSNV検出スクリプト

インストールされているpythonライブラリの確認 

pip freeze 

ないものがあれば導入。

pip install pysam 
pip install PyVCF
pip install Biopython
pip install --user numpy scipy matplotlib ipython jupyter pandas sympy nose

Inchworm assemblerはTrinitiyの中のアセンブラで知ってる人も多いと思う。TrinityのソースコードHPからダウンロードしてインストールする。 上記ファイルをzipでダウンロードして解凍。中にあるInchwormディレクトリに移動。

unzip trinityrnaseq-Trinity-v2.3.2.zip
cd trinityrnaseq-Trinity-v2.3.2/Inchworm/

マニュアルに従って Inchwormをビルド。

./configure --prefix=`pwd` 
make
make install

 bin/内に実行ファイル"inchworm"ができているはず。これにパスを通すかusr/bin内にコピー。

これでようやくランができるようになる。

 

ランは

python ScanIndel.py -i sample.txt -p config.txt [options]

ScanIndelは内部でblatを動かし、自動で上のコマンドを走らせポートを開けてリファレンスとのマッチを行う。そのため、blat検索時のポート番号の指定などのオプションを備えている。

 

 

 

 

-------参考-------

scanindelはリードとリファンレスの相同部位検索にblatを用いている。blatはblastに類似した相同部位検索プログラムで、indexファイルの作成がいらないなどの特徴をもつ。詳細はこちら

 blatスタンドアローンとwebサーバー版(UCSC)がある。最初にbinaryファイルを作るところはどちらも同じ。

faToTwoBit ref.fa ref.2bit

そのあとサーバー版ではgfServerコマンドを使いポートを開ける。

gfServer -trans -mask start localhost 17775 ref.2bit & 
gfServer start localhost 17774 ref.2bit &

 localhost自分で指定する名前。

 17774は自分で指定するポート番号。

 

 

 ----エラーが出たとき-----

ただしたくさんのスクリプトを統合しているためか、エラーが出てきたときはlogファイル(gfserver.temp.log)を開きどこで止まったか確認して進める。今回はbwa memのプロセスでエラーが出た。

これはsamtools sortコマンドがバージョンにより

samtools sort input.bam input_sorted

でsorted.bamができるbwaのバージョンと

samtools sort input.bam > input_sorted.bam

にしないとsortred.bamができないバージョンがあるかららしい。この件はgithubのScanindelのところでもコメントで表示されている。logにbwa memのソートでエラーが出ていた場合、pythonスクリプトを以下のように修正。

まずScanIndel.pyをテキストエディタで開く。

94行目の

"bwa.bam"+output+".bwa.sorted"

"bwa.bam >"+output+".bwa.sorted.bam"

に変更。さらに241行目の

'temp.bam' + output + '.temp.sorted'

'temp.bam >' + output + '.temp.sorted.bam'

に変更すれば動く。

 

 

ScanIndelのランに失敗すると、gfServerが停止せずポートが開きっぱなしになる。 このままだと、次回ランでエラーになる。ポート番号を--gfServer_portオプションで指定するか、ポートを閉じておく。ポートを閉じるには、まずポート番号をlsofコマンドで確認し、

lsof -i

表示されたプロセス番号が2587だとしたら

kill 2587

 これで停止できる。その他、途中で止まるとログファイル(gfserver.temp.log)も残っている。これも次回ランでエラーを起こすぽいので、エラー内容を確認したら消しておく。 その他、シミュレーションででsmall indelを検出しようとすると、エラーログなしに解析が止まることもたくさんあった。その場合、large indelを少し加えると、ランは最後まで行くようになった。そのとき、small indelも正確に検出されていたので、何かしら小さなバグだと思われる。ランが途中で止まる人は大きなindelを擬似的に加えてやり直してみると良い(あまり懸命じゃないが)。

 

その他、シミューレーションデータで検証時、短いindelを検出しようとすると途中で止まることがあった。原因は不明だがlong indelを混ぜ込むとランは正常に終了したので、対処療法的にバグ回避はできたが原因は不明。

 

 

 

PRISM Jiang et al. (2012)

discordant paired read + split-read approaches. 

 

 インストールは

tar xvzf PRISM.x86_64.tar.gz
cd PRISM
file bin/smapper
export PRISM_PATH=$PWD

PRISMは大きく分けて4つのステップからなる。それらのtoolはtoolkitディレクトリ中に保存されている。ただし、実際のランはrun_PRISM.shを走らせるだけでよい。

 

テストランは

./run_PRISM.sh -m500 -e30 -i /home/user/PRISM_1_1_6/sample/chr100.sam -r /home/user/PRISM_1_1_6/sample/reference/chr100.fa

samファイルとfaファイルはフルパスで指定する。

ランが成功すると、カレントディレクトリにPRISM_outputができ、中にdiscordantなリードペアのクラスターがsam形式で出力される。これが構造変化の種類ごとに分けて出力されるのでかなりたくさんのsamができる(詳細はREADMEに書いてある)。肝心のindelはsplit_all.sam_ns_rmmul_cigar_sorted_svに保存されている。

 

パスを通してからのランは下のように行う。 

run_PRISM.sh -m600 -e100 -i /home/user/input.sam -r /home/user/ref.fa -I PRISM_working_directory -O output

 -Iと-Oでそれぞれ作業ディレクトリと出力ディレクトリを明示する。インサートサイズとSDを大きく間違うと途中で止まることがあるみたいので注意する。

 

 

 ざっくり総評

 上に紹介したソフトの中では、Fermikit、ScanIndelが最も検出率が高かった。コール内容も正確で、導入もunixが使えれば難しくもない。Fermikitについては、default条件で他のソフトが検出できなかった不完全な(その領域の50%くらいのリードに変異あり)large indel変異はコールできた点も注目できる。2017現在バクテリアのショートリードでのindel検出ではこの2つが最もオススメできる。

 Scanindelはlong indelの検出率も高く配列までコールする点も長所になるが、偽陽性もそれなりに多く、またソフトの柔軟性が低く、エラーが起こりやすい点はマイナス。例えばログファイルが残っているだけで次はエラーになる。

 3番目はBreseq。変異それぞれについてhtmlファイルを開ければリードの張り付きからアミノ酸置換の結果まで表示してくれる。よくここまで作り込んだなと感心するレベルで、Breseqは100bpくらいまでのsmall indel とlong deletionを高精度に検出可能で、varscanやpinldeのような漏れはほとんど出ない。他のソフトと比較してF-measure値(参照)は高いんじゃないだろうか。インフォマの知識のない人でもまとめられるレベルになっている。ただしBreseqではlong insertion、inversionなどの構造変化や不完全indelを検出するには不十分で、改善の余地もある。

 

 

 

 

 

 

 

 

 

 

 

 

 

 

引用

-------------------------------------------------------------------------------------------------------------------------

BreakDancer – Identification of Genomic Structural Variation from Paired-End Read Mapping

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4138716/

 

LUMPY: a probabilistic framework for structural variant discovery

https://genomebiology.biomedcentral.com/articles/10.1186/gb-2014-15-6-r84

 

Pindel: a pattern growth approach to detect break points of large deletions and medium sized insertions from paired-end short reads

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2781750/

 

DELLY2: structural variant discovery by integrated paired-end and split-read analysis

https://academic.oup.com/bioinformatics/article/28/18/i333/245403/DELLY-structural-variant-discovery-by-integrated

 

PRISM: PRISM: Pair-read informed split-read mapping for base-pair level detection of insertion, deletion and structural variants

https://academic.oup.com/bioinformatics/article/28/20/2576/202193/PRISM-Pair-read-informed-split-read-mapping-for

 

Platypus: Platypus: A New Variant Caller that Integrates Mapping, Assembly and Haplotype-based approaches

http://www.nature.com/ng/journal/v46/n8/full/ng.3036.html

 

FermiKit: FermiKit: assembly-based variant calling for Illumina resequencing data

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4757955/

 

FreebayesHaplotype-based variant detection from short-read sequencing

https://arxiv.org/pdf/1207.3907.pdf

 

SvAVA: SvABA: Genome-wide detection of structural variants and indels by local assembly

http://biorxiv.org/content/biorxiv/early/2017/02/01/105080.full.pdf

 

VarScan: VarScan: variant detection in massively parallel sequencing of individual and pooled samples

VarScan: variant detection in massively parallel sequencing of individual and pooled samples | Bioinformatics | Oxford Academic

 

Scalpel: Accurate de novo and transmitted indel detection in exome-capture data using microassembly

http://www.nature.com/nmeth/journal/v11/n10/full/nmeth.3069.html

 

Breakseek: BreakSeek: a breakpoint-based algorithm for full spectral range INDEL detection

https://academic.oup.com/nar/article/43/14/6701/2902746/BreakSeek-a-breakpoint-based-algorithm-for-full

 

INDELseek:  INDELseek: detection of complex insertions and deletions from next-generation sequencing data

https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-016-3449-9

 

Breseq: Identifying structural variation in haploid microbial genomes from short-read resequencing data using breseq

https://bmcgenomics.biomedcentral.com/articles/10.1186/1471-2164-15-1039

 

ScanIndel: ScanIndel: a hybrid framework for indel detection via gapped alignment, split reads and de novo assembly

https://genomemedicine.biomedcentral.com/articles/10.1186/s13073-015-0251-2

 

PRISM: PRISM: Pair-read informed split-read mapping for base-pair level detection of insertion, deletion and structural variants

https://genomemedicine.biomedcentral.com/articles/10.1186/s13073-015-0251-2