macでインフォマティクス

macでインフォマティクス

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

ゲノムの相同性の高い領域の網羅的な検索 MUMmer

 

MUMmerはゲノム全体を高速にアライメントするオープンソースのツールである。MUMmerは、Finisihしたゲノムだけでなくドラフトゲノムでも使用でき、容易に何百あるcontigのアライメントを行うことができる。最初の論文が発表されたのは1999年であるが(ref.1)、現在でもオープンソースのまま開発が続けられており、 間も無くMUMmer4が論文化されるとアナウンスされている。

 

MUMmerのインストールはbrewで行える。

brew install mummmer

 

 

MUMmerは、主に5つのtu-rukaranaru.mummer、nucmer、promer、rum-mummer1、そしてrun-mummer3である。 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

制作途中

 

 

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

Mummer1

1、Alignment of whole genomes

Nucleic Acids Research, 1999, Vol. 27, No. 11 2369–2376

http://mummer.sourceforge.net/MUMmer.pdf

 

Mummer2

2、Fast algorithms for large-scale genome alignment and comparison.

Delcher AL1, Phillippy A, Carlton J, Salzberg SL.

Nucleic Acids Res. 2002 Jun 1;30(11):2478-83.

http://mummer.sourceforge.net/MUMmer2.pdf

 

Mummer3

3、Versatile and open software for comparing large genomes

Stefan Kurtz*, Adam Phillippy†, Arthur L Delcher†, Michael Smoot†‡, Martin Shumway†, Corina Antonescu† and Steven L Salzberg†

Genome Biology 2004, 5:R12

http://mummer.sourceforge.net/MUMmer3.pdf

 

 

 

 

 

 

 

ゲノム比較のmurasakiと結果を表示するGMV

 

murasakiは複数ゲノムの相同性ある領域の探索を高速に行うツールで、GMVはその比較結果を見るためのビューアソフトである。領域によってカラフルな色がつくので、ゲノムリアレンジメントなどの構造変化をわかりやすく示すことができる。

 

公式サイト

murasaki genomemurasaki genome

 

murasakiのインストールはこの方のブログを参考にしてください。


下のような図を出力できる。

f:id:kazumaxneo:20170623134158j:plain

4種の非常に近縁なシアノバクテリアを比較。右端にはANI値を載せた。

 

ランさせた時の様子はこのような感じになる。

動画では、genebankファイルを2つ選択してmurasakiで相同性を調べ、出力されたフォルダをGMVで表示している。

 

チュートリアル

Yasunori Osana @ University of the Ryukyus

こちらも貼っておきます。

Yasunori Osana @ University of the Ryukyus

 

 

 

引用

Murasaki: a fast, parallelizable algorithm to find anchors from multiple genomes.

Popendorf K1, Tsuyoshi H, Osana Y, Sakakibara Y. PLoS One. 2010 Sep 24;5(9):e12651.

doi: 10.1371/journal.pone.0012651.

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

 

ゲノム比較 x 変異コール x ビューア を統合したソフト Mauve

 

mauveはよく似たゲノムのアライメントを行い、その結果を見やすいビューアで表示して比較できるソフトである。MacwindowsLinux版が用意されており、無償でダウンロードできる。

 

ダウンロードは公式サイトから行う。

the Darling lab | computational (meta)genomics

コマンドライン版もあり、そちらの方が細かくアライメントパラメータを選択できるが、ここでは使いやすいGUI版 (javaのプログラム) を説明する。

 

Align with progressive mauveを選ぶ。

f:id:kazumaxneo:20170623123720j:plain

 

FilesタブのAdd sequenceを選び比較する配列ファイルを選択。入力は、fasta、gbkに対応している。

f:id:kazumaxneo:20170623123800j:plain

 

Parametersタブでパラメータを設定。デフォルトで進める。

f:id:kazumaxneo:20170623124216j:plain

 

Scoresタブ。アライメントのスコアとペナルティ設定。デフォルトで進める。

f:id:kazumaxneo:20170623124304j:plain

最後にAlignを押しジョブを開始させる。

 

ラン中はジョブの進捗が別ウィンドウで表示される。

f:id:kazumaxneo:20170623124005j:plain

5Mbくらいのゲノム3つを比較すると、解析に10分くらいかかる。

 

3ゲノムを比較した動画が以下になる。一番上がcanuでE.coliのロングリードをアセンブルして作ったcontig.fa、2つ目がE.coliのgene bankファイル、3つ目がspadesでショートリードとロングリードをhybridして作ったcontig.faである。

動画のように、気になる領域に移動してアライメントを調べられる。genebankを入れておけば、orf情報を表示させることもできて便利である。動画の後半では,SNPs、indelを出力している。

 

ゲノムサイズが100Mbを超えると重く使いづらくなるので、ゲノムが小さい生物の比較に向いている印象を受けた。ただし、mauveの論文 (ref.1) ではmouseやヒトゲノムを比較しており、使えないことはない。

アライメントできるtoolは多いが、mauveのように各OSに対応しており、アライメントして結果をGUIで確認できるソフトはなかなかないと思う。よかったら試してみてください。

 

 

引用

1、Mauve: Multiple Alignment of Conserved Genomic Sequence With Rearrangements

Aaron C.E. Darling,1,2,6 Bob Mau,2,3 Frederick R. Blattner,4,5 and Nicole T. Perna2,5

Genome Res. 2004 Jul; 14(7): 1394–1403. doi: 10.1101/gr.2289704

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

 

 

Oxford NanoporeリードのFAST5 => FASTA / FASTQ変換

 

MNIONのシーケンスデータはFAST5というフォーマットで出力される。FAST5はHDF5という時系列データ関係でよく使われる形式に乗っ取っているらしい。塩基配列になっていないバイナリーなデータのため、ビューアソフトで開いても文字化けしてしまう。

このデータから塩基配列に変換するツールがいくつか発表されている。導入して検証してみる。

  

Poretools

公式ページ https://poretools.readthedocs.io/en/latest/index.html

 

インストールはbrewでいける。

brew install poretools

 

fastqの抽出。

poretools fastq <directory> > output.fq

--type bestの指定で2Dリードのみ出力する。他にtemplateかcomplementなどがある。デフォルトではallになっている。ディレクトリはnanopore/fast5/0/にあるなら、~0/まで記載しないとエラーになる。

 

fastaの抽出。

poretools fasta <directory> > output.fq

 

シーケンスのイールド。電圧の変更でポアが急激に減っていっていないかチェックするのに使える。

poretools yield_plot --plot-type reads <directory> --saveas output.pdf

 --saveas: PDFかpngで保存。オプション無しだと画面出力される。

f:id:kazumaxneo:20170621220955j:plain

 

リード情報。

poretools stats <directory>

例えば以下のような画面が出力される。

total reads 4000

total base pairs 9561055

mean 2390.26

median 1182

min 11

max 148048

N25 10196

N50 4566

N75 1919

 

リード長のヒストグラム

poretools hist <directory> --saveas output.pdf

--num-bins、 --max-length、 --min-lengthで間隔や下限、上限を指定することもできる。

f:id:kazumaxneo:20170621222115j:plain

塩基組成。

poretools nucdist <directory>

例えば以下のような画面が出力される。

A 2242285 9561055 0.234522759256

C 2455787 9561055 0.256853140161

T 2375819 9561055 0.248489209611

G 2487164 9561055 0.260134890972

 

クオリティスコア分布。

poretools qualdist <directory>

以下のような画面が出力される。

" 1 19774 9561055 0.0020681818063

# 2 47996 9561055 0.00501994811242

$ 3 118476 9561055 0.0123915195551

.

.

.

H 39 1074 9561055 0.000112330699907

 

クオリティスコア分布のbox plot。

poretools qualpos <directory> --saveas output.pdf

f:id:kazumaxneo:20170621222741j:plain

 

フローセルの各ポアの総合パフォーマンス。

poretools occupancy <directory>

f:id:kazumaxneo:20170621223139j:plain

他にも、塩基変換前の生データやsquiggle plot を出力する機能などを持っている。公式サイトで確認してください。

 

Nanopolish

こちらでインストールは紹介している。


 

PoreSeq

Githubページ

GitHub - tszalay/poreseq: Error correction and variant calling algorithm for nanopore sequencing

githubからダウンロードして、PoreSeqのルートでsudo pip install -e .でインストールできる。

mac OSXで動作はするが、内部でLASTとsamtoolsを動かしているらしく、samtoolsのsort付近でエラーになりファイルが消える。現在のところ未修正。

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Oxford Nanoporeリードのアセンブリ パフォーマンス比較

 

S. cerevisiaeとE.coli K12をilluminaとnanoporeで読んだ論文のデータ が公開されている。

http://schatzlab.cshl.edu/data/nanocorr/

 

このデータを使い、2017年6月現在のナノポアリードのパフォーマンスを調べてみる。

上記URLからMinionとilluminaのシーケンスデータ (fastq) をダウンロードした。 

 

 

E.coliのリード数などは以下の通り。

・Oxford nanopore Minion: 59008 reads、total_length 240Mb、average 4068 bp、max 43000bp、min 120bp

・illumina Miseq (301x2): 22Mx2 reads、total_length 13.4Gb、max 301 bp

ナノポアのリードは平均4000bpを超えているが、リードのトータルサイズは2.4億bpしかないため、5Mbゲノムとした時のカバレッジは48xにしかならない。一方でMiseqのトータルサイズは301bp計算で67億bpを超えている!。カバレッジを計算すると2680xで、かなりディープに読まれたデータだ。gzから解凍したファイルサイズはそれぞれ15GB近くある。

ナノポアのlong readのサイズ分布をヒストグラムで調べる。

f:id:kazumaxneo:20170622150029j:plain

横軸がリードの長さ(bp)で、縦軸がその長さのリードが出てくる頻度である。(右端は単位を5000bpに変えていることに注意してください)。一番多いのは1000bp程度のリードであるが、綺麗な二項分布にはなっていない。6000-7000付近にも弱いピークがあるように見える。

 

S. cerevisiaeのリード数は以下の通り。

・Oxford nanopore Minion: 88218 reads、total_length 526Mb、average 5969 bp、max 72879bp、min 103bp

・illumina Miseq (301x2): 25Mx2 reads、total_length 15Gb、max 301 bp

S. cerevisiaeのゲノムサイズを大雑把に12.5Mbとした時のナノポアのカバレッジは44x、Miseqのカバレッジは1300xとなる。

ナノポアのlong readのサイズ分布をヒストグラムで調べる。

f:id:kazumaxneo:20170622150005j:plain

横軸がリードの長さ(bp)で、縦軸がその長さのリードが出てくる頻度である。(同じく右端は単位を5000bpに変えていることに注意してください)。一番多いのは1000bp程度のリードであるが、7000付近にも強いピークがある。

推測だが、7000bp付近のリードが本来のナノポアMinionのパフォーマンスで、左端のピークはDNA抽出時に物理的に切れてしまったDNAではないだろうか?(海外研究者のナノポアのワークショップの投稿で、λ HInd III cutをナノポアで読むと、等カバレッジにならずサイズが長くなるほどカバレッジが下がという話がある)。

 

 

アセンブル

以下はE.coliの例であるが、 S. cerevisiaeもゲノムサイズ以外は同じパラメータで実行した。 

1、canuを使いONTリードだけでアセンブル。

canu -p test -d long_read genomeSize=4.8m -nanopore-raw oxt.fq

E.coli=> 3 contig

S. cerevisiae=> 27 contig

contigの長さをグラフにすると以下のようになる。

f:id:kazumaxneo:20170623112456j:plain

グラフの左からcontigの長さ順に示す。縦軸がcontigの長さ(bp)である。

 

S.cerevisiaeは17クロモソーム持つので、クロモソームのサイズとcontigの長さを比較した。下の表は左側がS. cerevisiaeのゲノムをサイズ順に並べたもので、右側がナノポアリードをアセンブルして作ったcontigをサイズ順に並べたものである。

chr Size (bp)   contig Size (bp)
IV 1530000   tig00000001  1546985
VII 1090000   tig00000003  1091026
XV 1090000   tig00000007  1086779
XII 1080000   tig00000006  939389
XVI 948066   tig00000013  805656
XIII 924431   tig00000022  737318
II 813184   tig00013668  678294
XIV 784333   tig00000024  644519
X 745751   tig00000053  573971
XI 666816   tig00013665  519780
V 576874   tig00000181  487116
VIII 562643   tig00000055  462121
IX 439888   tig00000058  429563
III 316620   tig00000072  271056
VI 270161   tig00000185  259055
I 230218   tig00000066  239610
MT 85779   tig00000183  237080
      tig00000050  190760
      tig00000075  188868
      tig00000035  187180
      tig00000184  162871
      tig00000086  155346
      tig00000042  140091
      tig00000036  110210
      tig00013670  86062
      tig00000009  52963
      tig00000046  25807

おそらく多くのクロモソームは1つにアセンブルされている。quastでエラーチェックも行なったが、リファレンスゲノムと比べて明確な構造変化は起こっていない。つまり、あとはPCRを総当たりすれば、ゲノムの構造は決まってしまう。かなり衝撃的な結果である。(ただしエラーチェックは入念に行う必要がある!)

E.coli K12の方は、3つのcontigにアセンブルされていた(E.coli K12 MG1655のゲノムは4.64Mbのクロモソームのみを持つ。プラスミドはない)。

 

  

2、spadesを使いilluminaのペアリードをアセンブル。

spades

E.coli=> 763 contig

S. cerevisiae=> 1389 contig

f:id:kazumaxneo:20170623111700j:plain

 

lengthのtop100に絞って表示する。縦軸がcontigの長さ(bp)である。

 

 

3、spadesを使いONTリードとilluminaリード両方でhybirdアセンブル。

spades

E.coli=>  602 contig

S. cerevisiae=>  1145 contig

E.coliのcontigは一番長いもので4.64Mbで、残りは2kb以下のcontigだった。グラフにすると下のようになる。

f:id:kazumaxneo:20170623112142j:plain

lengthのtop30に絞って表示する。縦軸がcontigの長さ(bp)である。

 

 

 

 

 

 

 

引用

 

ショートリードのアダプタートリミングツール

 

これまで様々なアダプタートリミングツールが報告されてきた。OMIC toolsで検索すると、2017年6月で35件ヒットする(OMIC toolリンク)。ただしFastQC、cutadapt、Fastx-toolkitなどよく使われるツールは限定されている。Trim Galore!はFastQCとcutadaptを内部で動かし、fastqから自動でアダプター配列を認識してトリムするツールで、Biostarのhit数を見る限り使っている人も多そうである。インストールして動作を見ていく。

 

Trim Galore!

公式ページリンク Babraham Bioinformatics - Trim Galore!

Trim Galore!本体はperlスクリプトで、内部でFastQCとcutadaptを動かしてトリムを行う。そのため動作にはFastQCとcutadaptが必要である。この2つを先にインストールしておく。

 

インストー

FastQCとcutadaptのインストールでbrewで行える。

brew install cutadapt
brew install FastQC

 

FastQCとcutadaptのインストールが終わったら、公式サイトからTrim Galore!をダウンロードして解凍する。解凍したディレクトリのトップ階層にあるtrim_galoreを/usr/local/bin/に移動する。

mv TrimGalore-0.4.3/trim_galore /usr/local/bin

/usr/local/binにはパスが通っているので、 これでパスは通る。

 

ラン

trim_galore <fastq>

でランできる。

ペアリードなら両方のデータを指定し、--pairedのオプションをつけて実行する。

trim_galore --paired R1.fq R2.fq

出力は入力ファイル名_val_1.fq、入力ファイル名_val_2.fqとなる。

 

自動で認識させないでアダプター配列を明示するには-a <seq>をつける。配列を明示させるオプションは以下のようなものがある。

  • --illumina Adapter sequence to be trimmed is the first 13bp of the Illumina universal adapter 'AGATCGGAAGAGC' instead of the default auto-detection of adapter sequence.
  • --nextera Adapter sequence to be trimmed is the first 12bp of the Nextera adapter 'CTGTCTCTTATA' instead of the default auto-detection of adapter sequence.
  • --small_rna Adapter sequence to be trimmed is the first 12bp of the Illumina Small RNA 3' Adapter 'TGGAATTCTCGG' instead of the default auto-detection of adapter sequence. 

他にも様々なオプションがある。オプションを見るには

trim_galore --help

 

Trim Galore!のデフォルト動作はphread scoreをphred33と認識して行うが、古いIllumina 1.5 encodingのシーケンスデータを使う場合、--phred64をつけて明示する必要がある。

 

 

Nextraで調整されMiseqでランされたペアードエンドデータでテストすると、合計323万bp(データの4%)がトリムされた。

  >>> Now performing quality (cutoff 20) and adapter trimming in a single pass for the adapter sequence: 'AGATCGGAAGAGC' from file T1second_R2.fastq <<< 

This is cutadapt 1.8.3 with Python 2.7.10

Command line parameters: -f fastq -e 0.1 -q 20 -O 1 -a AGATCGGAAGAGC T1second_R2.fastq

Trimming 1 adapter with at most 10.0% errors in single-end mode ...

Finished in 5.59 s (21 us/read; 2.81 M reads/minute).

 

=== Summary ===

 

Total reads processed:                 261,774

Reads with adapters:                    69,845 (26.7%)

Reads written (passing filters):       261,774 (100.0%)

 

Total basepairs processed:    77,820,790 bp

Quality-trimmed:               3,237,687 bp (4.2%)

Total written (filtered):     74,481,043 bp (95.7%)

 

=== Adapter 1 ===

 

Sequence: AGATCGGAAGAGC; Type: regular 3'; Length: 13; Trimmed: 69845 times.

 

No. of allowed errors:

0-9 bp: 0; 10-13 bp: 1

 

Bases preceding removed adapters:

  A: 38.0%

  C: 25.5%

  G: 15.9%

  T: 20.6%

  none/other: 0.0%

 

Overview of removed sequences

length count expect max.err error counts

1 57208 65443.5 0 57208

2 9072 16360.9 0 9072

3 2290 4090.2 0 2290

4 536 1022.6 0 536

5 110 255.6 0 110

6 40 63.9 0 40

7 117 16.0 0 117

8 21 4.0 0 21

9 39 1.0 0 31 8

10 47 0.2 1 29 18

 

 

 

 

 

 

 

 

 

 

 

 

 

 

ナノポアのアセンブルデータのキュレーション及び変異の検出 nanopolish

 

ナノポアリードでアセンブルしたcontigのエラー修復と、変異のコールを行うツール。

インストールから動作まで見ていく。

 

ダウンロード

Github HP

GitHub - jts/nanopolish: Signal-level algorithms for MinION data

 

インストー

mac ではビルドできないとの情報があるので、cent OSにインストールした (dockerを使い

macでランする方法もオーサーたちは紹介してます)。

 

依存するもの

解凍したディレクトリのルートでmakeするだけでbiopython以外の依存するパッケージは自動導入してくれる。

make 

brewで導入することもできる。macでエラーが出る人は、biostarの関連スレッドを参照。

 

 代表的な機能を見ていく。

 

Minionから出力されたFAST5をfastaやfastqに変換する。

nanopolish extract: extract reads in FASTA or FASTQ format from a directory of FAST5 files

 

 

 

ドラフトゲノムをキュレートする。

nanopolish variants --consensus: calculate an improved consensus sequence for a draft genome assembly

ミスマッチが修復されるだけでなく、contigが伸びたりスキャホールドが長くなることもある。ドラフトゲノムをpolishしてクオリティを上げるワークフロ-は以下のようになる。

1、FAST5からfasta(fastq)を抽出

2、bwa memなどを使い、nanoporeのリードをドラフトゲノム(.fa) にアライメント

3、nanopolishでゲノムを50kbのセグメントに分け、コンセンサス配列を分析。

4、nanopolishで各セグメントをマージし、polishしたゲノム (.fa) を出力。

 

 

実際にテストしてみる。

S. cerevisiaeとE.coli K12をilluminaとnanoporeで読んだ論文のシーケンスデータ が公開されている。

http://schatzlab.cshl.edu/data/nanocorr/

49GBある生データをダウンロードし、変換から行ってみる。

 

 

 

 

作成途中

 

 

 

 

他にも以下のような機能がある。

 

メチル化サイトをコールする。

nanopolish call-methylation: predict genomic bases that may be methylated

 

SNVとindelをコールする。

nanopolish variants: detect SNPs and indels with respect to a reference genome

 

 

 

 

2, Oxford Nanopore sequencing, hybrid error correction, and de novo assembly of a eukaryotic genome

Sara Goodwin,1 James Gurtowski,1 Scott Ethe-Sayers, Panchajanya Deshpande, Michael C. Schatz, and W. Richard McCombie

Genome Res. 2015 Nov; 25(11): 1750–1756. doi: 10.1101/gr.191395.115

 

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