macでインフォマティクス

macでインフォマティクス

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

Indel検出ツールまとめ

small indel structural variations (SV)

 

 

 250bp x 2のペアリードでテストした時の結果をツールごとにまとめた。

インストール編はこちら 

 

Pindel Ye et al. (2009)

split-read approaches.

　splitリード法はいくつか報告されているが、その中でおそらく最もよく知られたツール。原理はインストール編に記載した。

  • 1-数十bpのindel検出。検出率は100%近いが、insertionサイズに反比例する傾向にある。
  • 1-300 bpのdeletion検出。
  • 10 bp-5 kbのinversion検出。
  • Tandem duplicationの検出。単純なInsertionでは検出されない1-kb以上まで検出可能で、検出率は50%以上。

まとめ: 数十bpのidnel検出とinversionの検出。split-read法はエラーに弱いため、リード長が極端に短いNGS初期のサンプルや、クオリティの悪いデータではfalse callの処理で悩まされる恐れがある。それでもinversion、deplicatioponなど他の多くのツールで検出できない構造変化を予測できるため、他の手法を使ったツールと併用することで漏れを防げる可能性がある。

 

pindel -f ref.fa -T 12 -i config.txt -o output

-T: スレッド数

-f: リファレンス（fasta形式）

-i: コンフィグファイル(install編を参照)

-o: 出力ファイル名

 

GATK Haplotypecaller

split-read approaches. 

  • 1~100 bpのindel検出。検出精度はほぼ100%。
  • 1=100 bpのindel検出。リード長以上のサイズのdeletionは検出対象にならない。
  • リード長以上のサイズのInsertionも検出対象にならないが、無理やり複数のin/delとして検出することがある。

まとめ:  ハプロタイプ解析向けで1~100 bpのindel検出が可能。バクテリアの不完全indel検出などにも利用できる。long indelの扱いが不透明なため、ゲノムのリアレンジメントが想定される場合long indelやそれ以外の構造変化検出ツールとの併用が推奨される。

 

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

-T: 手法の指定(Unifiedgenotyper / haplotypecaller)

-R: リファレンス（fasta形式）

-I: ソート済みbamファイル

-o: 出力ファイル名

 

FreeBayes Garrison and Marth (2012)

gapped alignment approaches.  

　ハプロタイプ前提でindelを検出するツール。ファミリー間のハプロタイプ解析などの他、バクテリアゲノムの不完全変異の解析にも利用できると思われる。

  • 1~10 bpのindelの検出。検出精度はほぼ100%。
  • 数十bp以上のlong indelは無視される。無理やり検出されることはほぼない。
  • ハプロタイプ解析のためdefaultで高感度な条件になっており、.4%のSNV（50カバレッジで2リードのみ塩基置換）もコール対象になる。
  • 検出感度はオプション"-C"で変更可能。
  • Inversion、Tandem duplicationは10 bp前後なら検出可能。

まとめ:  SNVと1-数十bpのindelの検出。ハプリタイプの解析など置換率のわずかなSNVの高感度検出。　GATKやVarscanと異なりSNVとSVをまとめてコールするが、他のSNV検出ツールのようにlong indelをSNVとして強引にcallすることはない。

 

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

 -f: リファレンス（fasta形式）

-b: ソート済みbamファイル

-C: 検出の閾値 (defautは2で高感度)

 

Break-dancer-max Chen et al. (2009)

red-pair approaches.  

　ペアリード間の距離の異常や、向きの異常（普通は--> <--）を足がかりにしてlarge indelなどの構造変化を捉える手法。1つの部位に同じような種類の異常ペア(discordant read pair）が2つ以上存在するるとコールするアルゴリズムになっている。実際にはbamファイルを最初の数百行サンプリングしてインサートサイズを調べ、それを元に解析が行われる。よってインサートサイズのばらつき具合(SD) によってfalse positiveの数が変化する可能性がある。

 　コール精度はサイズが大きくなると上がるため、100 bp以上のindelなら高確率でコールされる。一方1-10 bpのindelは全く検出できない。 100 bpくらいのindelでもコールできないこともしばしばある。

  • 100 bp以上のInsertion、Deletion、Inversion、Tandem duplicationの検出。検出精度はほぼ100%。ただしインサートサイズとそのばらつきにより検出可能サイズは変動する。
  • 構造変化の種類は予測されるが、配列情報は得られない。
  • １つの構造変化を複数変異として検出することが多い（排他性は考慮されていないと思われる）。

 

まとめ:  100-200 bp以上の大きな構造変化の検出。inversionの検出精度も高い。配列情報も得られないため、主力で使うよりは他のツールの追試に使う方が現実的と思われる。

 

perl bam2cfg.pl R1R2_sorted.bam > output.cfg #パラメータファイルの作成
breakdancer-max output.cfg > output.txt #SVの検出

 

VarScan Koboldt et al. (2009)

gapped alignment approaches. 

 

  • SNV検出モードとSV検出モードがある。
  • 複数のアライナーをサポートしている（BLAT, Newbler, cross_match, Bowtie and Novoalign）。
  • SV検出モードでは1-の10bp程度のindelを検出可能。精度は他のプログラムよりやや劣る。
  • １つのin/delを複数の変異としてコールすることがある。
  • SNV検出モードを走らせると、 in/del部位のsoft clipされたリードをSNVとして解釈しコールしてしまうことが多い。
  • Inversion、Tandem duplicationは10 bp前後なら検出可能。
  • javaで書かれたVerscan2が存在する。

まとめ: 偽陽性/偽陰性コールが他のプログタムより多い傾向がある。モデルの見直しで精度が改善されたVerscan2が存在する。

 

 

VarScan pileup2indel mpileup > output.txt

 

Platypus Rimmer et al. (2014)

local realignment & local assembly based approaches. 

多型を検出する手法として構築された。他の手法ではlong indelをSNVやsmall indelとして強引にコールすることがあるが、Platypusはローカルリアライメントとアセンブリを行い変異の正確なbreakpointを評価できるようにチューニングされている。indel検出パフォーマンス比較の論文でも高評価を受けている。

  • SNVと1-数 bpのindelやInversion、duplicationの検出。検出精度はほぼ100%。
  • assemblyを行なっているが構造変化検出はできない。例えば250 bpx2データでも100 bpの構造変化は全く検出されない。
  • 変異が疑われる領域（active region）のみをターゲットにするため、解析速度は非常に高速（100Mbゲノムで１分超）。
  • 感度が非常に高く、デフォルトではサポートするリードが２で検出される。

まとめ: SNVと1-数十bpのindelの検出。ローカルリアライメントとアセンブリを行うことで変異のbreakpointを正確に検出することができる。計算速度も速いため、small indel解析の主力ツールとしての利用を検討してもよい。解決できなかったlong indelを無理やりコールすることもないので、大きな構造変化を捉えるツールと併用しやすい。

ダウンロードリンク 

http://www.well.ox.ac.uk/platypus

以下のコマンドでビルドする

cd Platypus_X.X.X 
./buildPlatypus.sh
python Platypus.py callVariants --help #テストラン

 

python Platypus.py callVariants --refFile=ref.fa --bamFiles=R1R2_sorted.bam --nCPU 12 --minReads 4

--refFile= :  リファレンス（fasta形式）

bamFiles= : ソート済みbamファイル

--nCPU=: CPU数 (default 1)

--minReads: サポートするリード数。(default 2)

 

LUMPY layer et al. (2014)

split-read approaches  + read-pair + read-depth

　複数手法を使ってindelを検出するプログラムはLUMPY以前にも発表されていたが、いずれも1つの手法で変異を絞り込み、それから他の手法で制度を上げていくものだった。このやり方では検出される変異が初めの手法に依存する度合いが強くなる。LUMPYはsplit-read、read-pair、read-depthの3つの手法を同時に評価して、各種法間の整合性の高さを元にindelなどの構造変化を評価する方法論となる。

　1000 Genome projectなど複数のサンプルがある場合、サンプル間の比較を行いindel検出精度を高めることも可能となっている。その場合、1つのアルゴリズムだけを走らせることになるらしい（論文Fig. 1B）。

  • 100 bp以上のInsertion、Deletion、Inversion、Tandem duplicationが主な検出対象となる。ただし１つの構造変化を複数のindelとしてコールすることが多い。
  • 1-10 bpのsmall indelは検出されない。

まとめ: long indelなら検出されるようだが、一つの部位を複数コールしたりすることがあるためやや使いにくい。手法としては面白いが、現状では実用性はやや低い印象。

 

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

 

 

Breseq Deatherage et al. (2014)

split-read approaches.

Breseqはhaploidのバクテリアゲノムのindelをgene locusまで含めて正しく検出するために開発された。これまでのindel検出ツールは偽陽性のコールがまだまだ多く、indel breakpointも正しく評価できないものが多かった。またユニークにマップされたリードのみ解析するツールが多く、repeat 領域のindelを評価できない問題もあった。Breseqはこれらの問題に取り組み、repeat領域も検出対象とした。検出はsplit-read法を用いており、read-depthなどを調べ統計処理を行うことでin/delを判別してコールする。論文ではE.coliの変異を感度良く検出できたことが記載されている。

1-数十 bpのindel検出。検出精度は100%近い。100 bpくらいのindelでも7割くらいは検出される。それ以上長いindelは検出できず無視される。

  • 10 bp程度のTandem duplication検出。
  • long Inversion検出。ただしInversion末端のBeakpointを2つのindelとして検出する。
  • ランはfastqとリファレンスを指定するだけでワンライナーで実行でき、また出力はdefaultでhtml形式になっており、ビューアなしに変異のリードの分布を確認できる。
  • 極力false callが出ないよう工夫されている。
  • シングルリードでもランは可能。
  • プログラムはilluminaのほか、454にも対応している。

まとめ: 変異部位のリードのアライメント状況、deletion部位のカバレッジなどもhtmlファイルから全てアクセスすることができるため、素早く変異の状況を確認できる。論文の手法としてよくここまで作り込んだと思わせる完成度で、NGSに不慣れな人に特にお勧めできる。

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

-r; リファレンス（fasta形式、.gbk (genbank形式)、GFF3に対応)。

-j: スレッド数

 

Scanindel Yang et al. (2015)

 gapped alignment + split reads + de novo assembly)

　複数の手法を合わせた手法。以前から似た手法は報告されていたが、Scanindelは手法の実行順序を考えサイズの様々なindelをできる限りなく漏れなく、誤りなく検出するようにパイプラインが工夫されている。

　検出の流れは、１）BWA-MEMを使ったギャップを考慮したアライメント。２A）BLATによる高精度なアライメント。２B）2Aと同時に、soft-clippedされたリードとマッピングされなかったリードを使ったde novo assemblyを行いcontigを作成する。３）１のアライメントではsmall indelが、２の高精度なアライメントでは１で見つかりにくいmedium indelが、３のcontigとリファレンスの比較では１、２で見つからないリード長以上のindelが見つかる可能性がある。それらの候補部位をハプロタイプ検出のFreebayesでコールする。

　１のアライメントでsoft-clippedされたリード領域というのはおそらくsmall indelを横断しており、indelをリード内に完全に含んでいる。そこで、クリップされた領域を対象に２でBLATアライメントを行うというわけである。BLATはBLASTと似た手法で、NGSのアライナーより高精度なアライメントを行なっている。並行してInchwormアセンブラを使ってアセンブルを行なう（InchwormはTrinitiyのアセンブルでも使われている）。論文ではハイクオリティなリードのみ使うことが強調されており、これがfalse callを減らすキーなのだと思われる。

 

  • 1-100 bpのindel検出。検出精度はほぼ100 %で、1 kbのindelもおよそ60 %は検出される。
  • １つのlarge indelを複数回検出することがある。
  • Inversionは完全に無視される。
  • 10 bp程度のTandem duplication検出。

まとめ: 100 bp以上のlndelも高確率でコールされる。ただし長さは間違っている場合が多い。しかしながら他のindelプログラムがlong indelをSNVやsmall indelとしてコールしてしまう傾向が強いのと比べると、ScanIndelは間違っているなりにlong indelはlong indelとしてコールするあたり改善が見える。ただしリード長より長いindelが近接して２つ以上存在すると、2つのindelを１つとして捉えるなどlong indelの解釈がややおかしい事もある。プログラムもスクリプトとしてこなれていない感じがあり、他の論文ではまだ改善の余地があると指摘されている。 

python ScanIndel.py -i sample.txt -p config.txt

-i; リード情報 (fastqのパス、出力ファイル名）

-p: コンフィグファイル(リファンレスとそのindexのパス、blatのリファレンスのパス)

 

fermikit Li (2015)

assembly based approaches. 

　アセンブリベースの手法は他にも存在するが、これまではsplit-readやdiscordant read-pair情報を元に場所を限定してからlocal assemblyを実行する方法が中心だった。それに対してfermikitはde novo でフルにリードをアセンブリし、contigをリファレンスにアライメントして構造変化を検出する。この方法は重くなるのが難点だが、論文によるとFermikitは30 x humanデータでアセンブルがおよそ１日、コールが１時間半と実用的な時間でランできるらしい（16 CPU、peak memory 85 GB）。

　解析の流れは、１）BFCでエラーコレクション、２）ropeBWT2でBWT constraction、３）Fermiでアセンブリ、４）BWA-MEMでマッピング、５）HTSboxでコール、の流れとなる。車輪の再発明を避け、先行研究のツールを有効活用している。

  • 1-5000 bpのInsertion、Deletion、Inversion、Tandem duplication検出。検出精度は100 %近い。
  • リードより長いindelについては配列は出力されず長さだけレポートされる。

まとめ: アセンブリがうまくいけば、small indelだけでなく1 kb以上の構造変化もほぼ100 %検出可能である。しかしアセンブルの成功度合いに手法は依存しているため、50 bp前後のリードではわずかな繰り返し配列部分でも解析は厳しいと思われる。またコピーナンバーバリアントの変異についてはどうゆう挙動になるかはわからないし、laptopなどのメモリが少ない環境では解析できないかもしれない。とはいえ方法論としては理想的なものと思われる。 

fermi2.pl unitig -s 3m  -t 12 -l 301 -p prefix R1.fastq > prefix.mak
make -f prefix.mak 
fermi.kit/run-calling -t 12 ref.fa prefix.mag.gz | sh

-t: スレッド数

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

-l: リード長

-p: 出力prefix名

-2: 2pass エラーコレクション

 

SvABA Wala et al. (2017)

local assembly approaches. 

　2015年発表のFermikitとは異なり、local de novo assemblyを行なってindelを検出する。解析手順は

１）ギャップを考慮したアライメントでsmall indelを検出、２）clipped readを検出しmedium indelを検出、３）discordant read pairを検出してlocal de novo assemblyを実施してindelを検出、の流れになっている。indel以外に、稀にだが1000bpくらいのゲノムのリアレンジメント（転移）を捉える事もあるらしい。

 

  • 1~300 bp以下のinsertion検出。特に20-300 bpのmedium indelの検出感度が優れている。
  • 1 bp以上のdeletion検出。1 kbp以上のdeletionもほぼ100%検出できる。
  • inversionも検出するが、inversionのbreakpointを2つの変異として検出する。
  • Tandem duplicationの検出。単純なInsertionでは検出されない1 kb以上まで検出可能で、検出率は100%近い。

 まとめ:

 

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

-p: スレッド数

-t: bamファイル。

-G: リファレンス（fasta形式）

 

PRISM Jiang et al. (2012)

read-pair + split-read approaches. 

　手法としてはPindelにに似ているが、PRISMはread-pairの情報をより積極的に使う。具体的にはPRISMは以下の流れでSV検出を行っている。1）discordant read pairを探し、それがクラスターを作っている部位をまとめる。２）１の部位からsoft-clippedされたリードの情報を抽出しsplit-alignmentを行う。３）アライメント結果から構造変化を分析してコールする。

 

  • 検出可能な構造変化はPindelに似ており、small insertion。small-large deletion。inversionを見つけることができる。

 まとめ:

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

 

 

 

 

 

　どの検出法も万能ではないため、あらゆる構造変化を想定するとこれというプログラムを１つだけ上げるのは難しい。それを踏まえて使いやすい/使いにくいソフトを選別するなら、まず当然だが擬陰性(false negative)が多いソフトは使いづらい。2つ目に、構造変化を無理やり塩基置換として検出（false positive）してしまうようなスクリプトも使いづらい。リアルデータではどれがfalse positiveでどれがTrue callなのかわからないので、基本的に検出された以上全ての部位を分析する必要が出てくるからである。

　それらのことを前提に優劣をつけるなら、シミュレーションデータではアセンブリベースの手法が最も精度が高い結果となった。特にde novo assemblyベースのFermikitの検出率は高く、シミュレーションデータでの1-5000bpのindelの検出率は99 %だった。さらにTandem duplicationについても正確なサイズとブレイクポイントをダブりなく予測していた。またPlatypusはレアなsmall indelの検出感度が高く、small indelをほぼ100%捉えていた。

　local de novo assemblyのSvAVA、Scanindelも良好な結果を出した。scanindelは漏れがあったものの、5 kbのindelを60 %ほど検出し、かつ配列まで正確に予測していた。Indelのサイズが大きくなると、リードをが限定するlocal de novo assemblyより全データを使ったフルのde novo assemblyの優位が出ていると推測される。ただしFermikitはリードより長いindelは予測できても配列はコールしないため、medium Insertionに高感度で配列も予測するSvAVAを並行して使うと良いかもしれない。

　Indel検出にはアセンブリベースの手法が有効だが、long Inversionの検出精度はイマイチだった。例えばFermikit、FermikitはInversionの末端を2つのindelとして検出し、ScanindelはInversion変異を完全に無視していた。それに対してsplit-read法のPindelのInversion検出精度は非常に高く、配列まで正確にInversionを予測していた。Pindelはリードの個数を指定して検出できるので、高感度な解析を行いやすいと思われる。

やや古い手法であるがBreakdancerもInversion変異を正しく検出していた。Breakdancerでは配列に関する情報は得られないが、もっと複雑な構造変化を捉えられる可能性もある。 

 

 まとまりが悪いので表にします。

 

 

 

 

条件

リード長250bpのペアリード

インサートサイズは600、インサートサイズのSDは10に設定。