macでインフォマティクス

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

マッピングソフト(アライナー)のパラメータ設定

 

日付は古いが、Biostarsにアライナー比較のディスカッションがあった。

When and why is bwa aln better then bwa mem?

 

質問はbwa alnの優位点は何かということみたいだが、回答者がエボラゲノムをシミュレートして10%の高いエラー率で配列を合成し、そのゲノムにbwa back-trackとmem、bowtie2でアライメントさせてマッピング率を算出すると

  • bowtie2: 30%
  • bwa aln: 25%
  • bwa mem: 85%

だったと書いている。bwa alnとbowti2のスコアが極端に悪いのはシード領域が固定で、その領域でマッチしなければアンマップ扱いになってしまうからだと考えられる。つまりエボラゲノムの不安定性を考慮したら、bwa memのようなエラーに強い手法でないとマッピングの時点で大量のリードをロスしてしまい、変異を評価できず使い物にならなくなってしまう。

ただしこれはデフォルトの設定の話で、bowtie2は--very-sensitive-localtと

-D 20 -R 3 -N 1 -L 20

上のオプション付きで実行すれば91%のマッピング率になったらしい。-Nは最初のマッチングに使うシード領域のミスマッチ許容値。-Lはシード領域の長さ。そのほかのオプションは

-D <int>           give up extending after <int> failed extends in a row (15)

-R <int>           for reads w/ repetitive seeds, try <int> sets of seeds (2)

 

ただしBWA MEMは1-Mbpのロングリードのアライメントに対応しており、また100bp以下のリードのアライメントにも優れるとされる。リードがハイクオリティなら大抵はbwa memで問題ないとは思われる。

 

 

 

メタゲノム解析ツール

 

 

Genome sequences of rare, uncultured bacteria obtained by differential coverage binning of multiple metagenomes Albertsen et al. (2013)

 メタゲノムデータから、各生物ごとのデータを大まかに仕分け、ドラフトゲノム構築などのポストプロセスを助ける手法。難培養性でデータ比率がシーケンス比率が全体の1%程度の株のゲノムでも、ドラフトゲノムを構築できるとされる。初めて聞いた時はシンプルかつエレガントなやり方に目から鱗が落ちた。

 

 まずゲノムDNAを抽出比率が変化するような2つのやり方で抽出する。論文中では下水処理の活性化汚泥からフェノール有り/無しでDNAを抽出して、各々に別のindexをつけilluminaの150 bp x 2でシーケンスしている。そして、得られた1.7億リードと9千万リードをde novoアセンブルする。

 そこから先が面白い。全コンティグのカバレッジをフェノール有り/無しそれぞれで算出し横軸と縦軸にとってプロットする。さらに各コンティグのGC含量も計算して、各contigのプロットに色をつける。こうするとどうなるかというと、同じような色のコンティグが特定の領域に収束するようになる(下の図)。

f:id:kazumaxneo:20170523104958j:plain

公式ページから引用。

上のグラフはRのggplot2で作成されており、プロットのサイズでコンティグの長さも表現している。特定の領域に集まった同じような色のコンティグは、おそらくは同一生物か、極めて似た生物のゲノムと考えられる。理由は、抽出試薬によるDNAの抽出率は、各生物の細胞壁構造等によって固有値を持つためで、DNA抽出手法が同じでも生物ごとにDNA抽出率は変わるためである。GC含量も生物ごとに一定の値を示す。

 グラフをかき、1箇所に集まった同じ色のコンティグをラフに囲み回収する(ビニングする)。論文ではそのあたりをperlスクリプトを書いて処理している。回収後、コンティグ間にまたがってマップされるペアリードを元にコンティグの隣接関係を推測し、edge-nodeのグラフとして可視化する。ゲノムを決めたければ、PCRや追加のロングリード情報を使いフィニッシングまで持っていくこともできると思われる。

 

手順はこちらに非常に丁寧に書かれている。手法についてのスクリプトやRのコマンドについて丁寧に教えてくれるので、インフォマの勉強にもなる。ページ右端にマニュアルの各ステップページへのリンクがある。

 

 

 

 

 

 

PhyloPhlAn Segata et al. (2014)

比較する生物間に保存されている全orfを使って系統樹を作成する手法。16Sのような単独のorfを使った系統樹より生物間の距離を正しく評価できるとされる。保存されている全遺伝子を使うため、わずかな遺伝子の違いでも評価できる。そのため株間の系統解析のような16Sで分解能が出ない非常に近縁な種同士の分析にも有用と思われる。

 

解析の流れはこちらに書かれている。その中のリンクからダウンロードする。

ランにはFasttreeが必要。公式ページこちらからダウンロードするかbrewでインストール。

brew install FastTree

 

アミノ酸ファイル(faa)をダウンロードして一箇所のフォルダに集め、そのフォルダ名をinputとして実行する。

./phylophlan.py -u example --nproc 12

  -u, --user_tree: Build a phylogenetic tree using user genomes only 

  --nproc N:  The number of CPUs to use for parallelizing the blasting

 

 

exampleデータを実行するとoutput/exampleにnwk形式のファイルが出力される。系統樹ソフトMEGAに読み込ませて結果を確認。

f:id:kazumaxneo:20170523104808j:plain

 

 

 

ConStrains Luo et al. (2015)

メタゲノムデータから、遺伝子間のSNPなどを指標に種内間の多様性を評価するためのツール。ある種の生物がどのくらいの多様性を持ってどのような比率で存在しているか解析することが可能。

 

gitからダウンロードする。 (Wikiリンク)

git clone https://bitbucket.org/luo-chengwei/constrains

ランにはpytohn2.7環境が必要。さらにpythonのBiopython、Numpy、Scipy、NetworkXを必要とする。そのほか、Bowtie、Metaphlan2などが必要。ないものだけ以下のコマンドでインストール。

pip install Biopython
pip install Numpy
pip install Scipy
pip install networkx
brew install Bowtie

Metaphlan2は2015年に発表された手法。ダウンロードリンクは下のMetaphlan2紹介の方に記載している。

 

テストコマンド

はじめに以下のコマンドでfastqデータとconfig_fileを自動ダウンロードする。

python utils/download_tutorial.py -o test #30分くらいかかる。

ダウンロードされたfastqはE.coliのいくつかの株の混ぜ物からARTでシミュレートされたペアリードデータらしい。データの比率は以下のように書いてある。

The two samples contain ART simulated E.coli Illumina reads. In sample_1, the composition is 5% (ETEC H10407) and 95% (HS), and in sample_2, the composition is 15.2% (IAI1), 26.3% (CFT073), and 58.5% (O157:H7).

テストではこの株の同定を試みている。

解凍する。 

tar xzvf constrains.test/fq.tar.gz

できた4つのfastqファイルをターミナルのカレントディレクトリに移動する。またはダウンロードされたtest.confを編集して、fqのパスを修正する。

ランは以下のように行う。

python ConStrains.py -c test/test.conf -o Constrains_test -t 12

 -c: コンフィグファイルの指定

-t: スレッド数

 

テストは10分くらいで終わる。Constrains_test/results/Intra_sp_rel_ab.profilesを開くと

# Species   strain_ID   masked_samples  sample_1    sample_2
Escherichia_coli    str-1   NA  0.00000000  14.88058489
Escherichia_coli    str-2   NA  3.75299123  0.00000000
Escherichia_coli    str-3   NA  96.24700877 85.11941511

となっていた。Intra_sp_rel_ab.profilesは同一種内の株の比率 (%) で(つまり合計は常に100)、Overall_rel_ab.profilesはメタゲノムのコミュニティ全体から見た株の比率を表す。

 

SRAに登録されていったbiogass producing micorobialのメタゲノムデータでも試してみた。サンプル数1は12GBのfastqが2つある。-t=16でランすると1時間程度で終了した。結果のOverall_rel_ab.profilesを見ると下のようになっていた。

# Species strain_ID masked_samples ERR843249
Lactobacillus_curvatus NA,insufficient 1 0.28228
Enterococcus_cecorum NA,insufficient 1 2.77777
Coprococcus_catus NA,insufficient 1 0.11833
Weissella_koreensis NA,insufficient 1 0.37525
Prevotella_copri NA,insufficient 1 0.22341
Methanoculleus_bourgensis NA,insufficient 1 25.03193
Aminobacterium_colombiense NA,insufficient 1 1.07588
Streptococcus_macedonicus NA,insufficient 1 0.02517
Lactobacillus_johnsonii NA,insufficient 1 0.11143
Thermoanaerobacter_thermohydrosulfuricus NA,insufficient 1 5.57953
Streptococcus_lutetiensis NA,insufficient 1 0.6998
Lactobacillus_reuteri NA,insufficient 1 0.45737
Streptococcus_infantarius NA,insufficient 1 8.61296
Leuconostoc_citreum NA,insufficient 1 0.15594
Lactobacillus_amylovorus NA,insufficient 1 3.39372
Agrobacterium_sp_H13_3 NA,insufficient 1 34.99582
Aerococcus_viridans NA,insufficient 1 0.49274
Methanoculleus_marisnigri NA,insufficient 1 0.24412
Staphylococcus_lentus NA,insufficient 1 0.12412

 シーケンス率の異なるかなり多様な菌が混じっていることがわかる。

 

  

 

 

MetaPhlAn2 Luo et al. (2015)

 

 解説ページ。

Metaphlanはメタゲノムデータから種のレベルで分類する手法。clade特異的なマーカー遺伝子を元に仕分けを行なっており、カバレッジを計算することでデータ中の比率も自動計算されるようになっている。Metaphlan2ではこれをさらに高精度化して、strainの分類まで検出することが可能になっている。7500種以上の生物から100万近いマーカー遺伝子を準備し、ウィルスも検出できるようになったらしい。 

こちらのリンクからダウンロードする(大量のマーカ遺伝子を含むため1GB以上のサイズになっている)。

 

 

 

 

Metacoder Foster et al. (2017)

 Rのパッケージとして提供されている。

インストールは、Rを起動して

install.packages("metacoder")

 

 

 

 

 

 

------------------------------------------引用------------------------------------------

Genome sequences of rare, uncultured bacteria obtained by differential coverage binning of multiple metagenomes

Genome sequences of rare, uncultured bacteria obtained by differential coverage binning of multiple metagenomes. - PubMed - NCBI

 

ConStrains identifies microbial strains in metagenomic datasets.

ConStrains identifies microbial strains in metagenomic datasets. - PubMed - NCBI

 

PhyloPhlAn is a new method for improved phylogenetic and taxonomic placement of microbes

Nicola Segata, Daniela Börnigen, Xochitl C. Morgan, and Curtis Huttenhower. Nature Communications 4, 2013

https://www.ncbi.nlm.nih.gov/pubmed/23942190

  

MetaPhlAn2 for enhanced metagenomic taxonomic profiling

Duy Tin Truong, Eric A Franzosa, Timothy L Tickle, Matthias Scholz, George Weingart, Edoardo Pasolli, Adrian Tett, Curtis Huttenhower & Nicola Segata. Nature Methods 12, 902–903 (2015)

MetaPhlAn2 for enhanced metagenomic taxonomic profiling. - PubMed - NCBI

 

Metacoder: An R package for visualization and manipulation of community taxonomic diversity data

Metacoder: An R package for visualization and manipulation of community taxonomic diversity data

 

 

 

 

Indel検出ツールまとめ

 

 

 

 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つの部位の変異は1つと限定せずにコールする。ファミリー間のハプロタイプ解析などの他、バクテリアゲノムの不完全変異の解析にも利用できると思われる。

  • 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%。ただしインサートサイズとそのばらつきにより検出可能サイズは変動する。
  • 構造変化の種類は予測されるが、配列情報は得られない。
  • 1つの構造変化を複数変異として検出することが多い(排他性は考慮されていないと思われる)。

 

まとめ:  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を検出可能。精度は他のプログラムよりやや劣る。
  • 1つの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検出パフォーマンス比較の論文でも高評価を受けている。

  • 変異はハプロタイプを前提に評価される。
  • 1-数十 bpのInsertion、Deletion、Inversion、Tandem duplication検出。検出精度はほぼ100%。
  • assemblyを行なっているが、long indelは検出対象外。例えば250 bpx2データでも100 bpの構造変化は全く検出できない。
  • 変異が疑われる領域(active region)のみをターゲットにするため、解析速度は非常に高速。

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

 

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

--refFile= :  リファレンス(fasta形式)

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

 

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が主な検出対象となる。ただし1つの構造変化を複数の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 (genebank形式)、GFF3に対応)。

-j: スレッド数

 

Scanindel Yang et al. (2015)

 gapped alignment + split reads + de novo assembly)

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

 検出の流れは、1)BWA-MEMを使ったギャップを考慮したアライメント。2A)BLATによる高精度なアライメント。2B)2Aと同時に、soft-clippedされたリードとマッピングされなかったリードを使ったde novo assemblyを行いcontigを作成する。

3)1のアライメントではsmall indelが、2の高精度なアライメントでは1で見つかりにくいmedium indelが、3のcontigとリファレンスの比較では1、2で見つからないリード長以上のindelが見つかる可能性がある。それらの候補部位をハプロタイプ検出のFreebayesでコールする。

 1のアライメントでsoft-clippedされたリード領域というのはおそらくsmall indelを横断しており、indelをリード内に完全に含んでいると考えられる。そこで、その領域を対象に2でBLATアライメントを行うというわけである。BLATはBLASTと似た手法だが、NGSのアライナーより高精度なアライメントが可能となっている。同時にInchwormアセンブラを使ってアセンブルを行なう(InchwormはTrinitiyのアセンブルでも使われている)。その結果から様々なindelを検出している。論文ではハイクオリティなリードのみ使うことが強調されている。おそらくそれがfalse callを減らすキーなのだと思われる。

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

まとめ: 100 bp以上のlndelも高確率でコールされる。ただし長さは間違っている場合が多い。しかしながら他のindelプログラムがlong indelをSNVやsmall indelとしてコールしてしまう傾向が強いのと比べると、ScanIndelは間違っているなりにlong indelはlong indelとしてコールするあたり改善が見える。ただしリード長より長いindelが近接して2つ以上存在すると、2つのindelを1つとして捉えるなど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データでアセンブルがおよそ1日、コールが1時間半と実用的な時間でランできるらしい(16 CPU、peak memory 85 GB)。

 解析の流れは、1)BFCでエラーコレクション、2)ropeBWT2でBWT constraction、3)Fermiでアセンブリ、4)BWA-MEMでマッピング、5)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を検出する。解析手順は

1)ギャップを考慮したアライメントでsmall indelを検出、2)clipped readを検出しmedium indelを検出、3)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を探し、それがクラスターを作っている部位をまとめる。2)1の部位からsoft-clippedされたリードの情報を抽出しsplit-alignmentを行う。3)アライメント結果から構造変化を分析してコールする。

 

  • 検出可能な構造変化は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

 

 

 

 

 

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

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

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

 

 

 

 

 

 

条件

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

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

 

fastq配列のシミュレーション

ARTを使ったやり方を記載する。

https://www.niehs.nih.gov/research/resources/software/biostatistics/art/

 

ARTはIllumina's Solexa, Roche's 454 and Applied Biosystems' SOLiDのリードジェネレートができる。illuminaはMiseq v3の250bpをサポートしている。gccのライブラリをインストールしておけばbinary版をダウンロードして叩くだけで動く。

ARTはこちらからダウンロード。

から一番新しいmacOSX版(ソースと書いてない方がbinary版)

 

テストスクリプト

./art_illumina -ss MSv3 -sam -i ref.fa -p -l 250 -f 100 -m 600 -s 10 -o 250bp_paired

-l: リード生成元のfasta

-l: リード長(最大250)

-m:インサートサイズ(リード込みのサイズ)

-f: リードデプス(100なら-lで指定したゲノムのx100カバレッジ

合成速度はそれなりに速い。バクテリアサイズなら数分で終わる。ペアードエンドのfastqの他にそのalnファイルとsamファイルもできる。indel解析のシミュレーションなどする場合に便利。illuminaのデータをシミュレートすると、5'末のクオリティが非常に低いことを確認できる。

 

 

パラメータを貼っておく。

===== PARAMETERS =====

 

  -1   --qprof1   the first-read quality profile

  -2   --qprof2   the second-read quality profile

  -amp --amplicon amplicon sequencing simulation

  -c   --rcount   number of reads/read pairs to be generated per sequence/amplicon (not be used together with -f/--fcov)

  -d   --id       the prefix identification tag for read ID

  -ef  --errfree  indicate to generate the zero sequencing errors SAM file as well the regular one

                  NOTE: the reads in the zero-error SAM file have the same alignment positions

                  as those in the regular SAM file, but have no sequencing errors

  -f   --fcov     the fold of read coverage to be simulated or number of reads/read pairs generated for each amplicon

  -h   --help     print out usage information

  -i   --in       the filename of input DNA/RNA reference

  -ir  --insRate  the first-read insertion rate (default: 0.00009)

  -ir2 --insRate2 the second-read insertion rate (default: 0.00015)

  -dr  --delRate  the first-read deletion rate (default:  0.00011)

  -dr2 --delRate2 the second-read deletion rate (default: 0.00023)

  -k   --maxIndel the maximum total number of insertion and deletion per read (default: up to read length)

  -l   --len      the length of reads to be simulated

  -m   --mflen    the mean size of DNA/RNA fragments for paired-end simulations

  -mp  --matepair indicate a mate-pair read simulation

  -M  --cigarM    indicate to use CIGAR 'M' instead of '=/X' for alignment match/mismatch

  -nf  --maskN    the cutoff frequency of 'N' in a window size of the read length for masking genomic regions

                  NOTE: default: '-nf 1' to mask all regions with 'N'. Use '-nf 0' to turn off masking

  -na  --noALN    do not output ALN alignment file

  -o   --out      the prefix of output filename

  -p   --paired   indicate a paired-end read simulation or to generate reads from both ends of amplicons

                  NOTE: art will automatically switch to a mate-pair simulation if the given mean fragment size >= 2000

  -q   --quiet    turn off end of run summary

  -qL  --minQ     the minimum base quality score

  -qU  --maxQ     the maxiumum base quality score

  -qs  --qShift   the amount to shift every first-read quality score by

  -qs2 --qShift2  the amount to shift every second-read quality score by

                  NOTE: For -qs/-qs2 option, a positive number will shift up quality scores (the max is 93)

                  that reduce substitution sequencing errors and a negative number will shift down

                  quality scores that increase sequencing errors. If shifting scores by x, the error

                  rate will be 1/(10^(x/10)) of the default profile.

  -rs  --rndSeed  the seed for random number generator (default: system time in second)

                  NOTE: using a fixed seed to generate two identical datasets from different runs

  -s   --sdev     the standard deviation of DNA/RNA fragment size for paired-end simulations.

  -sam --samout   indicate to generate SAM alignment file

  -sp  --sepProf  indicate to use separate quality profiles for different bases (ATGC)

  -ss  --seqSys   The name of Illumina sequencing system of the built-in profile used for simulation

       NOTE: sequencing system ID names are:

            GA1 - GenomeAnalyzer I (36bp,44bp), GA2 - GenomeAnalyzer II (50bp, 75bp)

           HS10 - HiSeq 1000 (100bp),          HS20 - HiSeq 2000 (100bp),      HS25 - HiSeq 2500 (125bp, 150bp)

           HSXn - HiSeqX PCR free (150bp),     HSXt - HiSeqX TruSeq (150bp),   MinS - MiniSeq TruSeq (50bp)

           MSv1 - MiSeq v1 (250bp),            MSv3 - MiSeq v3 (250bp),        NS50 - NextSeq500 v2 (75bp)

BLAST2GOでアノテーションをつける

basic版をここからダウンロード。インストールが終わったら立ち上げる。

統合TVの説明を参考にした。

リンク

 

 

BLAST2GOを起動したら、上のメニューからSTART -> Load Sequences -> Broseを選択し読み込むFASTAファイルを選択 -> Openボタンで読み込み。

 

基本的に左のアイコンから順に実行していく。解析が終わるとその色がつく。例えばblastはオレンジアイコンだが、blastがおわった配列はオレンジになってかつTagの列が紫になる。

 

 

 左端から3番目のBLAST

上のメニューからBLAST -> Run blast -> NCBI blast -> blastxを選択。データベースは"nr"。メールアドレスを入れて、そのほかはデフォルトのまま実行。

1配列で10分前後かかる。数千配列あるなら10日くらいかかるかもしれない。ちなみに

 local データベースを作成して高速化も可能らしい。

 ただし、nrのデータベースでもかなり容量がある。2016年頭にダウンロードしたときは解凍後200GBくらいにはなった。ダウンロードにも相当な時間がかかったし、保存場所も確保しないといけない。ローカルサーバーがあれば問題ないが、そうでなく時折しかblast解析しないのであればネットワーク版で問題ないと思う(ローカル版のデータベースの差分更新も既知の方法はありそう)。

 

 左端から4番目のINTERPRO SCAN

BLAST解析と並列化することが可能。役割はこちらのスレッドが分かりやすい。

やってることはモチーフ検索。blastでhypothetical proteinしかヒットしなくてもInterpro scanならみつかるかもしれない。統合TVの説明。

 

 上のメニューからInterpro -> Run interpro -> メールアドレスを記入してそのほかはデフォルトのまま実行。blastより短い傾向にあるがblastx同様1配列で数分かかるので、やはりこれも気長に終わるのを待つ。

 

 

  左端から5番目のmapping

blast情報に従ってgo termを関連付けする作業。

上のメニューからmapping-> Run mapping ->Runボタンをクリック。これだけ。

 

 

  左端から6番目のアノテーション

Go termは複数つくので、適切なものを選ぶ作業。

上のメニューからannotation-> Run annotation-> デフォルトのままRunボタンをクリック。

 

 

  左端から7番目のStatics

様々な情報を見ることができる。

 

 

エンリッチメント解析

特定のサンプルのみ変動している遺伝子群をGo termでまとめることができる。

メニューバーのAnalysis -> Enrichment Analysisを選択。

 

 

 

 

 

 

 

 

 

Waifuで拡大時の画質劣化を抑える

 

 

sudo apt update

sudo apt install cmake git libgtk2.0-dev pkg-config libavcodec-dev libavformat-dev libswscale-dev

sudo apt install python-dev python-numpy libtbb2 libtbb-dev libjpeg-dev libpng-dev libtiff-dev libjasper-dev libdc1394-22-dev

sudo apt install cmake freeglut3-dev libglew-dev git

git clone https://github.com/glfw/glfw && mkdir build && cd build && cmake ../glfw && make -j4 && sudo make install

wget https://github.com/opencv/opencv/archive/3.2.0.zip

unzip 3.2.0.zip

cd opencv-3.2.0/

cmake .

make

sudo make install

cd

git clone https://github.com/khws4vl/waifu2x-converter-glsl.git

cd waifu2x-converter-glsl/

CLC genomics workbench (7.0)でRNA seq解析

CLCで行う前提でワークフローを書いてみる。

 

アダプター配列の確認

アダプター配列はシーケンス->画像からのシグナル取得 -> ベースコールファイル(.bcl) -> FASTQ変換の過程で自動除去されるようだが、インサートが短いペアリードファイルなどでは3'側にアダプターが存在している可能性がある。例えば以下の

https://www.researchgate.net/post/PCR_primers_trimming_for_MiSeq

でそのことが論じられている。アダプター残存による影響は解析によって変わってくる。例えばCLCの資料の

https://jp.illumina.com/content/dam/illumina-marketing/apac/japan/documents/pdf/2015_techsupport_session10.pdf

には、BWA MEMでマッピングした時には影響がないが(たぶんクリッピングされるから)、同じBWA でもショートリード向けのBWA Backtrackでは影響が多いと書かれている。アセンブルする場合も影響は大きいだろう。

 

 

アダプターが残存しているかCLCを使って調べてみる。

アダプター配列が存在しているか確認するには(オフィシャルサポートリンク)

 

 1、アダプターリストの作成

 2、fastqをダウンサンプリング

 3、モチーフ検索ツールでダウンサンプリングfastqファイル中のモチーフを検索

 4、アダプター配列が見つかることが確認できたら、アダプタートリミングツールでアダプターを除去

 

1、アダプターリストの作成

CLC7を起動

File -> New -> Trim motif list を選択。

ウィンドウが出てきたら下のAdd Rowボタンを押してアダプター配列情報を記入していく。Mismatchコストは初期値のままにしてFinish。パラメータの詳細はこちらを参考

f:id:kazumaxneo:20170418140342p:plain

 

アダプター配列分これを繰り返す。

最後にFile -> Save で保存。

 

 

2、fastqの一部データのサンプリング

・モチーフ検索ツールは1000リードまでしか対応していないみたいなので、FASTQのリードを1000だけランダムサンプリングする。手順は

・前もって読み込んだfastqファイルを右クリック -> Toolbox -> NGS core tool -> Sample reads

リード数を1000に指定してRun

10秒くらいでサンプリングは終わる。~ sampledというファイルが作られているはず。

 

 

3、モチーフ検索ツールでサンプリングしたfastqファイル中のモチーフを検索

・~ sampledファイルを右クリック -> Toolbox -> Classical Sequence Analysis -> General Sequence Analysis -> Motif search を選択。

・Motif Search TypeはMotif Listを選択。

・Motif listは先ほどで作成したリストを選択。

・Include negative strandにチェック

f:id:kazumaxneo:20170418141636p:plain

 

・次のページでCreate TablesとAdd anotations to sequencesにチェック

・resultはsaveを選択しfinish。

 ・解析が終わると、モチーフが見つかったリードを記載したリストができる。

 

 

4、アダプター配列が見つかるか確認

・~ sampledファイルの中身をダブルクリックで開き、右側のリストからAnnotation typeを選択する。分かりにくいが写真の青色のMotifというアイコンの右上にAnnotation typeというタブがある。

f:id:kazumaxneo:20170418142416p:plain

 

・青色のMotifボックスにチェックをつけるとモチーフが見つかった部位が矢印で表示される。今回は人工的に下の画像の矢印部位が100%マッチするようにモチーフ配列を設計した。

f:id:kazumaxneo:20170418144401p:plain

 

同じことをRでやるなら

library(qrqc)

fastq <- readSeqFile("input.fastq")

basePlot(fastq)

f:id:kazumaxneo:20170418154709j:plain

先頭だけATGC頻度が極端なのでアダプターと考えられる。

 

kmerKLPlot(fastq)  #6-mer時の塩基出現頻度

f:id:kazumaxneo:20170418155116j:plain

高頻出の配列がアダプターということになる。

 

k-mer長を変えるなら

fq <- readSeqFile("1.fq", kmer=TRUE, k=15) #k=15の時 

kmerKLPlot(fq) #グラフ描画

 

 

 

アダプター配列のトリミング

上の画像の部位のトリミングを行ってトリミングツールの挙動を確かめてみる。以下の手順で進める。

File -> New -> Trim Adapter list を選択。

・アダプター配列を順番に記載してsave。リストファイルを作る。(上のモチーフ検索リスト作成と同じ要領)

・NGS core tools -> Trim sequences を選択。

・いまアダプタートリミングだけ行いたいので、クオリティトリミングや5' or 3'末端トリミングのチェックを全て外す。先ほど作ったリストファイルを読み込ませ、Finishでラン開始。

 

終わったら、FASTQファイルのtrimmedというのができるので、開いてみる。

f:id:kazumaxneo:20170418145842p:plain

先ほどの青色矢印部分の配列が除去されている。

アダプターの上流側に配列がある場合、それもまとめて除去されている。5'側がまとめて除去されるので、index配列が分かればその上流も含めて全て除去できるのは都合がいい。ただし、配列の途中にアダプターが出てくる場合、強制的に5'側がトリムされる仕様は不都合が出てくるかもしれない。例えばインサートが短くてペアリードが調整したライブラリの末端までシーケンスされていた場合(特にsmall RNAシーケンスなど)、3'側にアダプター配列が出てくる可能性がある。その場合、下流側の配列を除く必要があり、強制的に5'側が削除されるのはまずい。

 

 

CLC以外のコマンドツールではどうか?

1、cutadaptの場合

 

@DRR057320.1 MG00HS11:593

TTAGGTTTCTACAAAATGAAGATTTCGAAAGTTTATCAAAACAAAGAATCT

+

@@@DDA=DHFFBDA:EFHIII@FHHHEG?+??CFEDEDHD:FHCG;DCGIF

上のリードの太字部分がトリム対象なら、cutadapt -gコマンドで

cutadapt -g TACAAAATGAAGATTTCG 1.fq > 1_trimmed.fq #5'側トリミング

 

トリム後、アダプターの5'側も除去される。

@DRR057320.1 MG00HS11:593

AAAGTTTATCAAAACAAAGAATCT

+

G?+??CFEDEDHD:FHCG;DCGIF

こうなる。3'側をトリムするなら-gを-aに変えて

cutadapt -a TACAAAATGAAGATTTCG 1.fq > 1_trimmed.fq #3'側トリミング

リードを見ると

@DRR057320.1 MG00HS11:593

TTAGGTTTC

+

@@@DDA=DH

今度はアダプターとその下流側が除去されている。

 

 

 

 

 

 

 

 

トリミング

詳細はこちらを参照した。 

CLC genomics workbench のトリムは累積のクオリティスコアが一定以上の領域が続いた場合に取り除くらしい。つまり1塩基だけクオリティスコアが悪くても必ずしもトリミングはされないアルゴリズムを用いている。

#アダプター配列の除去

toolboxNGS core tools → Trim sequences を選択。

 

de novo アセンブリ

trinityの結果と比較する。

Trinity --seqType fq --left S1_L001_R1.fastq --right S1_L001_R2.fastq --max_memory 10G --CPU 18

 

>行をカウントしてコンティグ数を見積もる。

Trinity 1656 contigs

CLC genomics workbench 7 1956 contigs

rnaspades 1952 contigs

 

数についてはTrinityが一番少なく、rnaspadesが一番多い。

長さの分布やlongest conntigを見るとどうか?

 

trinity

f:id:kazumaxneo:20170415141857p:plain

TOP10の長さ

f:id:kazumaxneo:20170415142259p:plain

 

 

CLC genomics workbench 7(default condition)

f:id:kazumaxneo:20170415142525p:plain

TOP10の長さ

f:id:kazumaxneo:20170415142214p:plain

 

 

 

 

rnaspades

rnaspades.py -t 30 -k auto -o data -1 1_R1.fq -2 1_R2.fq )

TOP10の長さ

Name length
NODE1 19799
NODE2 18969
NODE3 17335
NODE4 15793
NODE5 15760
NODE6 13466
NODE7 13232
NODE8 12802
NODE9 12596
NODE10 12065

 

 

top10の長さはrnaspadesとtrinityが優れている感じ。今回のデータはMiseq 301bpのペアデータである。他のデータでも同様の傾向が出るのかはわからない。