macでインフォマティクス

macでインフォマティクス

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

ゲノムなどの長い配列同士を比較し、違いをレポートする NucDiff

2018 10/13 コマンドエラー修正

2019 11/13 インストール手順訂正

2021 5/9 docker image追記

 

 

 全ゲノムシーケンシング戦略およびアセンブリアプローチの進歩により、一連の配列を互いに比較する方法が必要となっている。共通のクエスチョンは、同じリードセットの異なるアセンブリプログラムから得られたアセンブリ配列が互いにどのように異なっているか、または同じバクテリア種の異なる菌株のゲノムがどのように異なっているかである。このような分析を行うために、whole genome alignment(WGA)法が​​しばしば使用され、バイオインフォマティクスで長く研究されてきた。 WGAは、一般に、2つ以上の配列間の相同性の予測である[論文より ref.1]。 WGAは、主にゲノム間の保存された配列を同定するために使用される。ゲノム間の大規模な進化的変化を検出するゲノムの(functional)アノテーション、系統発生学的推論などをサポートする遺伝子、制御領域、非コードRNA配列、その他の機能要素[2,3]この分野は1970年代から継続的に発展しており、WGAのための多くの方法とツールが作成されている。既存の方法とツールのレビューは、[ref.1,4,5]にある。

 シーケンスセット間の相違を検出する目的で、WGA解析を実行するために使用できるツールには特定の機能が必要になる。第一に、非常に断片化されたゲノム、構造上の再構成、ゲノム配列の複製、およびリピート領域に関連することが多い様々な相違に対処できるべきである。第二に、比較分析結果は、差異の種類およびその位置に関する情報を提供すべきである。この情報は、今後の分析に適した方法で保存される必要がある。そのような比較情報は、例えばリファレンスアシストのゲノムアセンブリアセンブリエラー検出、異なるアセンブリの比較のために使用することができる。第3に、さまざまなレベルで詳細なアラインメント結果を視覚化できるようにする必要がある。グローバル規模のビジュアライゼーションは、複製、再構成、カバーされていない領域を調べるために使用できるが、ローカルスケールのビジュアライゼーションでは、置換、挿入、欠失などの小さな違いに関する情報が得られる。

 MAUVE [ref.6](紹介)、QUAST [ref.7]](紹介)、およびdnadiff [ref.8]の3つの異なるツールが現在利用できる。 MAUVEは、複数のゲノムアラインメントを行い、保存されたゲノム領域、これらの領域におけるリアレンジメントおよび逆位の正確なブレイクポイントを同定し、またヌクレオチド置換およびsmall indelsを同定する[ref.6]。また、インタラクティブな視覚化による結果の分析が可能で、情報を別々のファイルに保存する。しかし、小さな違い(置換、indels)に関する情報のみ、アクセサリプログラムを実行することなく簡単にアクセスできる。

 QUASTは、ゲノムアセンブリの品質評価のためのツールであり、リファレンスゲノム存在下でアセンブリ品質に関する様々なmetricsを出力する。QUASTは構造変化のタイプとポジション情報のみ出力する。 QUASTは、Icarusと呼ばれる付随のゲノムブラウザでの可視化を可能にしている。しかし、QUASTは検出された違いのいかなる可視化も欠いており、要約統計のみを提供している。

 Dnadiffは、MUMmer [ref.9]のNUCmerアラインメントプログラムのラッパーコマンドであり、相違を定量化し、アライメント統計と他の高レベルメトリクスを提供する[ref.8]。 QUASTと同様に、dnadiffはアセンブリの品質評価とゲノムの比較に使用できるが、検出された差を可視化することはできない。

 ここでは、シーケンス比較のためにMUMmerのNUCmer、delta-filter、show-snpsプログラムを使用するNucDiffツールを紹介する。 NUCmerは配列同士をアライメントし、アライメントした領域に関する情報を出力する。これらの領域の相対的位置の厳密な分析により、ゲノムのリアレンジメントや逆位を含む様々なタイプの差異の検出を可能にし、場合によっては、それらの領域とリピート領域との関連性も確かめる。 NucDiffは、closely relatedな2つの配列の差異を特定し、その差異をいくつかのサブタイプに分類する。両方の入力シーケンスに関する全ての差異の正確な位置は、GFF3(Generic Feature Format version 3、[ref.10])ファイルとして出力される。正確なポジションにより、視覚化とさらなる分析の両方が可能になる。したがって、NucDiffによって提供される情報は、2つの配列セットがどのように異なるかを明らかにするのに役立つ。

 NucDiffは、参照ゲノムおよびクエリーと呼ばれる2セットの配列間の様々なタイプの差異を、NUCmerによる配列のアライメントとdelta-filterおよびshow-snpsプログラムの結果を解析することによって決定する。 NUCmerはDNA配列アライメントを行い、デルタフィルタは指定された基準に従ってアライメント結果をフィルタリングする。デフォルトでNucDiffによって使用される設定では、delta-filterはクエリシーケンスの一貫性のある一番長いアラインメントを選択する。  NUCmer出力は、ソース配列、対応するリファレンスフラグメントに対するクエリフラグメントの方向、およびアライメントの類似性に関するすべてのフラグメントの正確な座標を含まれている。 show-snpsの結果には、対応リファレンスフラグメントと比較した、クエリフラグメント内のすべての挿入、欠失、および置換された塩基に関する情報が含まれる。

 

 

f:id:kazumaxneo:20180913151005j:plain

Classification of the types of differences. 論文より転載

著者らは、全ての変化は3つのグループに分けて考えることができると主張している。Nucdiffはそれぞれの違いを検出できるようパイプラインが構築されている(論文図2)。

 

インストール

ubuntu18.04のpython3.5.1::Anaconda4.1.0でテストした。

依存

NucDiff can be run on Linux and Mac OS.

  • Python 2.7
  • MUMmer v3.23
  • Biopython package

本体 GIthub

#ここでは仮想環境に入れる。python2.7の必要がある。
conda create -n nucdiff -c bioconda nucdiff -y
#bioconda (link)
mamba create -n nucdiff -c bioconda python==2.7.14 nucdiff -y
conda activate nucdiff

> nucdiff -h

 $ nucdiff -h

usage: nucdiff [-h] [--reloc_dist [int]] [--nucmer_opt [NUCMER_OPT]]

               [--filter_opt [FILTER_OPT]] [--delta_file [DELTA_FILE]]

               [--proc [int]] [--ref_name_full [{yes,no}]]

               [--query_name_full [{yes,no}]] [--vcf [{yes,no}]] [--version]

               Reference.fasta Query.fasta Output_dir Prefix

 

positional arguments:

  Reference.fasta       - Fasta file with the reference sequences

  Query.fasta           - Fasta file with the query sequences

  Output_dir            - Path to the directory where all intermediate and

                        final results will be stored

  Prefix                - Name that will be added to all generated files

                        including the ones created by NUCmer

 

optional arguments:

  -h, --help            show this help message and exit

  --reloc_dist [int]    - Minimum distance between two relocated blocks

                        [10000]

  --nucmer_opt [NUCMER_OPT]

                        - NUCmer run options. By default, NUCmer will be run

                        with its default parameters values, except the

                        --maxmatch parameter. --maxmatch is hard coded and

                        cannot be changed. To change any other parameter

                        values, type parameter names and new values inside

                        single or double quotation marks.

  --filter_opt [FILTER_OPT]

                        - Delta-filter run options. By default, it will be run

                        with -q parameter only. -q is hard coded and cannot be

                        changed. To add any other parameter values, type

                        parameter names and their values inside single or

                        double quotation marks.

  --delta_file [DELTA_FILE]

                        - Path to the already existing delta file (NUCmer

                        output file)

  --proc [int]          - Number of processes to be used [1]

  --ref_name_full [{yes,no}]

                        - Print full reference names in output files ('yes'

                        value). In case of 'no', everything after the first

                        space will be ignored. ['no']

  --query_name_full [{yes,no}]

                        - Print full query names in output files ('yes'

                        value). In case of 'no', everything after the first

                        space will be ignored.['no']

  --vcf [{yes,no}]      - Output small and medium local differences in the VCF

                        format

  --version             show program's version number and exit

——

 

 

ラン

クエリの配列と比較対象のリファレンス配列、出力、名前を指定して実行する。

nucdiff --vcf yes ref.fa query.fa outdir prefix

--vcf yesのフラグをつけると、SNV/indelについてvcfフォーマットでも出力される。

 

出力

resultsディレクトリに、リファレンスとクエリの各々から見たSNV、indel、構造変化それぞれの結果が出力される。

f:id:kazumaxneo:20180913151331j:plain

詳細はwikiに記載されています。

https://github.com/uio-cels/NucDiff/wiki

 

VCFやGFF3で出力されるので、そのままIGVに読み込んで可視化することもできる。

f:id:kazumaxneo:20180913170205j:plain

 

2021 5/10

 以前動作しないというお話を伺ったので、イメージ(*1)をdockerhubにpushしておきます。

docker pull kazumax/nucdiff
cd <change>/<to>/<genome_dir>/
docker run --rm -itv $PWD:/data/ -w /data kazumax/nucdiff
nucdiff -h

> nucdiff

(nucdiff) root@6695c6ada349:/data# nucdiff   

usage: nucdiff [-h] [--reloc_dist [int]] [--nucmer_opt [NUCMER_OPT]]

               [--filter_opt [FILTER_OPT]] [--delta_file [DELTA_FILE]]

               [--proc [int]] [--ref_name_full [{yes,no}]]

               [--query_name_full [{yes,no}]] [--vcf [{yes,no}]] [--version]

               Reference.fasta Query.fasta Output_dir Prefix

nucdiff: error: too few arguments

 

引用

NucDiff: in-depth characterization and annotation of differences between two sets of DNA sequences
Khelik K, Lagesen K, Sandve GK, Rognes T, Nederbragt AJ

BMC Bioinformatics. 2017 Jul 12;18(1):338. 

 

関連ツール


*1

Dockerfile

FROM kazumax/miniconda:1.0 AS build-image
RUN mamba update -n base -c defaults conda
RUN conda create -n nucdiff -c bioconda python==2.7.14 nucdiff -y && mamba clean -a -y
RUN echo "conda activate nucdiff" >> ~/.bashrc
ENV PATH=/opt/miniconda3/envs/nucdiff/bin:$PATH

Dockerfileの場所で

docker build -t kazumax/nucdiff .

ベースイメージkazumax/miniconda:1.0はこちら 

 

phylogenetic marker genesを検出し、marker genes全てを使って系統比較する自動化されたパイプライン ezTree

2019 3/9 docker pullリンク追記、インストールの流れ修正

2019 10/28誤字修正

2020 4/7 docker commnadの誤字修正

 

 メタゲノミクスおよびシングルセルゲノミクスは、様々な環境からの新規生物の発見および調査のための有望な方法として確立されている。 "microbial dark matter"という用語は、培養できない、微生物コミュニティからシーケンシングすることのみで研究される未培養生物を記述するために提案されたものであり[論文より ref.1]、新たに回収された1000以上の未培養生物のゲノムを既存の系統樹に取り込み、新しい視点で系統樹を調べる研究も報告されている[ref.2]。より多くの研究が、多種多様な微生物コミュニティから抽出された新規ゲノムの解析に焦点を当てており[ref.3-11]、これらが環境で果たす役割について我々の知識が拡張されつつある。

 微生物コミュニティを調査するための最も一般的な手法の1つは、環境からゲノム配列を直接得ることを目指すメタゲノミクスである。計算科学のビニング技術[ref.13-22]は、メタゲノムから直接、個々の生物ゲノムを抽出するために開発された。回収されたゲノムの微生物多様性を理解し、新たに同定された種をtree of lifeに配置するために、phylogenetic marker genesが使用されてきた。最も広く採用されているphylogenetic marker geneの1つである16SリボソームRNAスモールサブユニットの遺伝子は、新しく回収された生物の分類法を探り、系統樹を構築するための「ゴールドスタンダード」として確立されている[ref.23,24]。しかし、16S rRNA遺伝子の互いに非常に相同性が高い領域のために、Meta-IDBA [ref.25]、SPAdes [ref.26](紹介)、Ray Meta [ref.27]、MEGAHIT [ref.27](紹介)などのBruijnグラフベースのメタゲノムアセンブラでメタゲノムからインタクトな16S rRNA遺伝子を集めるのは依然として非常に困難な課題である[ ref.28]。その結果、メタゲノムから回収されたゲノムは、通常、16S遺伝子を欠く(または非常に短い遺伝子断片のみからなる)ため、16S配列を用いて系統樹を構築することは不可能または非常に困難である。

 系統樹上の個体間の関係を精緻化するため、全ゲノム情報を使うことが提案されている[ref.30-33]。 16S遺伝子ベースの系統樹を補うために、コンカテネート(連結)されたタンパク質のツリー(trees based on combined protein data alignments)が提案され、潜在的により堅牢で有益である[ref.34]。連結されたタンパク質ツリーを構築するためには、研究で考慮されるすべての生物において一度だけ出現する遺伝子として定義されるphylogenetic marker genesを同定する必要がある[ref.35]。この基準を満たす遺伝子は、以前の研究で示されているように、原核生物種の系統樹関係を確実に再構築するためのマーカーとして用いられている[ref.36,37]。そのようなマーカー遺伝子セットを発見するための様々な試みがなされている。例えば、Ciccarelli et alは 191種のバクテリア種において31のマーカー遺伝子を同定し、高い分解能を持つ系統樹を構築した[ref.36]。他の人々によって、異なるマーカー遺伝子セットも報告されている[ref.35,37,38]。また、checkMソフトウェア(紹介)は、系統特異的マーカー遺伝子セットを発見し、メタゲノムからビニングされた原核生物ゲノムの完全性およびコンタミ率をチェックするために使用される[ref.39]。

 メタゲノムからビニングされた個々のゲノムはめったに完全ではないので、予め定められたマーカー遺伝子セット由来の遺伝子のいくつかは、回収されたゲノムから失われている可能性がある。さらに、通常は数十または数百のゲノムが含まれるので、系統樹を構築するためのマーカー遺伝子セットを同定するために、各ゲノムの各遺伝子のコピー数を綿密にチェックする必要がある。 新たに回収された原核生物のゲノムから遺伝子を予測する努力を軽減するために、Prodigal [ref.40]やFragGeneScan [ref.41]のような信頼性の高い遺伝子予測ツールが開発されているが、ゲノムのセットの分類学的関係を推測するため、一群のゲノム中のマーカー遺伝子を自動的に同定するためのツールは依然として必要とされている。

 ここでは、原核生物ゲノムのセットからマーカー遺伝子および系統樹を推定するための計算パイプラインを紹介する。パイプラインは、新しく回収された、断片化された、または不完全なゲノムを含む一連のゲノムをとり、入力ゲノムからタンパク質コード遺伝子を予測し、すべてのゲノムによって共有されるマーカー遺伝子を同定し、マーカー遺伝子の連結タンパク質のアラインメントに加え、最尤(ML)法の系統樹を作成する。新しく回収されたいかなるクオリティのゲノムを持つユーザーも、このパイプラインで非常に簡単に系統樹を構築し、回収された種の分類を推論することができる。

 

 このパイプラインは、fasta形式の一連の原核生物ゲノム配列を取り込むように設計されている。 ゲノム配列は、完全、断片、または不完全であってもよい。 ユーザーが好む場合、ゲノム全体ではなくタンパク質配列を入力することもできる。 パイプラインのワークフローには、論文図1に示すように、ゲノムからのタンパク質コード遺伝子の予測、遺伝子への機能プロファイルの割り当て、ゲノムセットのシングルコピーマーカー遺伝子の特定、および系統樹のアライメントが含まれる。 

 

f:id:kazumaxneo:20180912202536p:plain

ezTreeのワークフロー。論文より転載。

 

ezTreeに関するツイート

 

インストール

mac os 10.13、Anaconda3.5.1でテストした。

依存

  • HMMER3 (Latest.  Assuming that we get version 3.1b2)
  • muscle (Latest.  Assuming that we get version 3.8.31)
  • Gblocks (Latest.  Assuming that we get version 0.91b)
  • Prodigal (Latest.  Assuming that we get version 2.63)
  • FastTree (Latest.  Assuming that we get version 2.1.9)
#Anaconda環境ならcondaで導入できる
conda install -c bioconda -y hmmer==3.1b2 muscle==3.8.31 gblocks==0.91b prodigal==2.6.3 fasttree==2.1.9

 本体 Github(マニュアルPDFもあり)

git clone https://github.com/yuwwu/ezTree.git
cd ezTree/

 > perl ezTree

$ perl ezTree 

Please input both a list file consisting of genomes and a output header.

 

ezTree - building phylogenetic trees for a set of genomes

version 0.1

Usage:

  ezTree

    -list (list file of genomes)

    -out (output header)

   (Either -list or -dir is required for running ezTree)

 

   (Other parameters)

    [-thread (thread num; default 4)]

    [-evalue (evalue for HMMER; default 1e-10)]

    [-model (JTT, WAG, or LG evolutionary models; default JTT)]

 

  Please read README file for more details.

kamisakBookpuro:ezTree kamisakakazuma$ less ezTree 

kamisakBookpuro:ezTree kamisakakazuma$ 

kamisakBookpuro:ezTree kamisakakazuma$ perl ezTree 

Please input both a list file consisting of genomes and a output header.

 

ezTree - building phylogenetic trees for a set of genomes

version 0.1

Usage:

  ezTree

    -list (list file of genomes)

    -out (output header)

   (Either -list or -dir is required for running ezTree)

 

   (Other parameters)

    [-thread (thread num; default 4)]

    [-evalue (evalue for HMMER; default 1e-10)]

    [-model (JTT, WAG, or LG evolutionary models; default JTT)]

 

  Please read README file for more details.

dockerイメージも上げておきます(version 0.1)。

docker pull kazumax/eztree

#currentと/dataをシェアしてラン
docker run -itv $PWD:/data/ kazumax/eztree
source ~/.profile

cd ~/ezTree/
perl ezTree -h

  

 

実行方法

テストデータを解析する。

1、リストファイルの準備

cd ezTree/test_example/

#ランにはfastaファイル名のリストが必要。lsで作成する。
ls *.fasta > listfile

 

2、実行

../ezTree -list listfile -out output -thread 8 -evalue 1e-5
  • -thread   (thread num; default 4)]
  • -evalue   (evalue for HMMER; default 1e-10)]

初回はPFAMのデータベースダウンロードも実行される。eztreeのルートディレクトリにdataディレクトリができる。

テストデータには23のProcaryotesゲノムが含まれており、テストジョブが完了するまでかなりの時間がかかる。実行時は9時間ほどかかった。

 

出力されるファイル

1. .aln: the concatenated alignment file of all marker proteins.
2. .nwk: the Newick tree for the genomes defined in the list file
3. .pfam: the identified single copy marker genes in terms of PFAM
4. .work directory: this is the work directory of ezTree. If one needs to re-run ezTree,
simply input the same “-out” parameter. ezTree will locate all temporary files in the 

 

output.pfamが同定されたマーカー遺伝子になる。今回は1378遺伝子見つかった。

>head  output.pfam

$ head out.pfam 

PF00238.18 Ribosomal_L14, Ribosomal

PF00398.19 RrnaAD, Ribosomal

PF00366.19 Ribosomal_S17, Ribosomal

PF00075.23 RNase_H, RNase

PF00828.18 Ribosomal_L27A, Ribosomal

PF04997.11 RNA_pol_Rpb1_1, RNA

PF02219.16 MTHFR, Methylenetetrahydrofolate

PF04551.13 GcpE, GcpE

PF14437.5 MafB19-deam, MafB19-like

PF00475.17 IGPD, Imidazoleglycerol-phosphate

 

MEGA7で.nwkを開いた。

f:id:kazumaxneo:20180913070129p:plain

 

引用

ezTree: an automated pipeline for identifying phylogenetic marker genes and inferring evolutionary relationships among uncultivated prokaryotic draft genomes
Wu YW

BMC Genomics. 2018 Jan 19;19(Suppl 1):921. 

 

PhyloPhlAn (Segata et al,. 2014) も簡単に紹介しています。


 

gANIを計算するツール ANIcalculator

 

 微生物は数と多様性の両方で生命の樹木を支配しており、その自然分類を困難かつ重要なものにしている。動物では、種は一般に交配可能な生物群と定義されるが(biological species concept)、この定義は無性生物の集合体に直接適用することはできない。結果として、微生物分類学は、生物に関する遺伝子型、表現型および化学感受性の情報を統合し、利用可能なデータのコンセンサスに基づいて微生物種を描写する多相性のアプローチ(論文よりref.2,3)を採用する。この多相性アプローチは、DNA-DNAハイブリダイゼーション(DDH)、G + C含有量の変動、選択されたDNAマーカー(16S rRNAを含む)の配列比較、脂肪酸、極性脂質、細胞壁などの特定の代謝産物の同定組成物およびエキソポリサッカライド、ならびに形態学的、生化学的および酵素学的特徴付けを含む。多相性分類学ツールボックスにおけるこれらの異なる方法が、同じ分離株の示唆された分類について完全に同意しないことは珍しいことではない。例えば、株のタイトな遺伝子型クラスター化は、それらの生化学的多様性および表現型変動性と矛盾する可能性があり、逆もまた同様である。これは、潜在的に異常で興味深い生物学的プロセスならびに実験上の誤りを示す可能性がある。しかし、特定の方法の選択および各データ型に割り当てられたウエイトは研究者の裁量に任されており、生物および研究の範囲によって異なることがあるので、分類の不一致が生じる可能性がある。高スループットシーケンシング技術の出現の前に進化したこの多相性アプローチは、これまでのところ原核生物分類を可能にした。しかし、生物の完全なゲノムは、その最終的な遺伝子シグネチャーであり、完全なゲノムは迅速で正確であり、もはや費用がかからない。全ゲノムベースのグループ化は、依然として単一の遺伝子または遺伝子の小さなサブセットから生じる生理学的変動を説明することができないかもしれないが、ゲノム距離を計算する際に考慮される情報を最大にし、よりバイアスを減らすために全ゲノム情報を活用し、新規微生物の発見とペースを同じにする。

(一部略)

  Whole-genome based Average Nucleotide Identity (gANI) は、KonstantinidisとTiedjeにより、2ゲノム間の類似性の尺度として提案された(ref.4)。ここでは、このアプローチを強化し、統合微生物ゲノム(IMG)データベース(ref.5)で公開されている86.5ミリオンのゲノムについてgANIを計算するスケーラブルな方法を開発した。著者らは、迅速ではあるが高感度な相同性検索(ref.6)であるbidirectional best hits (BBHs) でオーソロガス遺伝子のヌクレオチド同一性を計算し、平均化することによってgANIを計算する。 gANIに加えて、著者らは遺伝子含量に基づいてペアのゲノムの遺伝的関連性の相補的尺度としてオーソロガス遺伝子の割合(Alignment Fraction, AF)を考慮する。著者らは、多様なバクテリアゲノムおよびアーキアゲノムの大規模なセットについて、gANIとAFとの関係を系統的に探索し、統計的に特徴付ける。これら2つの測定値は、微生物の生物多様性のサンプリングが増加する中でゲノム間の遺伝的距離を正確に捕捉し、堅牢であることを実証している。それらは、一般に、既存の分類とよく相関し、既存の分類群に一様に適用されると、種内の遺伝的多様性のレベルを予測する一貫した種のアサインにつながる可能性のある閾値の選択を可能にする。(以下略)

 

JGIのANI HP

https://ani.jgi-psf.org/html/home.php?

 

インストール

ローカルで使うには、実行ファイル(linux向けにコンパイルされている)をJGIよりダウンロードする。

https://ani.jgi-psf.org/html/download.php?

(リンク先の一番下のacceptをクリック)

> ANIcalculator --help

$ ANIcalculator --help

********************************************************************************************************************

This tool will calculate the bidirectional average nucleotide identity (gANI) between two genomes.

Required input is the full path to the fna file (nucleotide sequence of genes in fasta format) of each query genome.

Either the rRNA and tRNA genes can be excluded, or provided in a list with the -ignoreList option.

********************************************************************************************************************

Usage:ANIcalculator

   -genome1fna <fna file of the first query genome> *REQUIRED*

   -genome2fna <fna file of the second query genome> *REQUIRED*

   -outfile <the output file> OR -stdout <output to screen> *Default: ANIcalculator.out*

   -outdir <output directory> *Default: Current directory*

   -ignoreList <file containing list of genes to ignore (Should include TRNA and RNA genes)>

   -logfile <log file> *Default: ANIcalculator.log*

   -help <prints this page>

 

Results are shown in tab-delimited format with following headers:

   Genome1 <name of FNA file of genome1>

   Genome2 <name of FNA file of genome2>

   ANI(1->2) <Average nucleotide idenitity of the first genome to the second>

   ANI(2->1) <Average nucleotide idenitity of the second genome to the first>

   AF(1->2) <Alignment Fraction of the first genome to the second>

   AF(2->1) <Alignment Fraction of the second genome to the first>

 

 

 ローカルでのラン

ANIcalculator -genome1fna input1.fasta -genome2fna input2.fasta -outdir outdir

 出力。

> outdir/ani.blast.dir/cat input1.fasta.input2.fasta.blout

# cat input1.fasta.input2.fasta.blout 

NC_000911.1 AE005174.2 66.38 2484 835 16 2368750 2716600 35208 2230375 1.1e+03 27248.00

NC_000911.1 AE005174.2 49.47 3307 1671 11 644387 1351186 436561 4738477 1.5e+03 7825.00

NC_000911.1 AE005174.2 40.02 2431 1458 13 3367689 3566614 3495590 4924147 3.1e+02 2266.00

NC_000911.1 AE005174.2 61.02 1411 550 2 3129308 3431892 794424 2112869 4.3e+02 861.00

NC_000911.1 AE005174.2 65.86 1775 606 8 3349356 3351139 136681 138470 1.6e+03 563.00

NC_000911.1 AE005174.2 74.11 1016 263 5 2448950 2449972 4226221 4227238 1.3e+03 490.00

NC_000911.1 AE005174.2 71.15 780 225 9 3327399 3328183 4902776 4903565 8.6e+02 330.00

NC_000911.1 AE005174.2 71.43 658 188 0 1667696 1668354 4778438 4779095 7.7e+02 282.00

NC_000911.1 AE005174.2 63.48 660 241 1 3550315 3550974 3614286 3614948 5.4e+02 178.00

NC_000911.1 AE005174.2 67.34 395 129 1 1520691 1521088 1794691 1795085 4.1e+02 137.00

NC_000911.1 AE005174.2 65.18 336 117 0 2152291 2152626 3477274 3477609 2.9e+02 102.00

NC_000911.1 AE005174.2 65.45 275 95 0 1073491 1073765 4329764 4330038 2.6e+02 85.00

NC_000911.1 AE005174.2 62.96 324 120 3 990415 990744 4135533 4135859 2.6e+02 84.00

NC_000911.1 AE005174.2 63.55 299 109 0 808917 809215 755514 755812 2.5e+02 81.00

NC_000911.1 AE005174.2 65.86 249 85 0 1026588 1026837 2198237 2198485 2.4e+02 79.00

NC_000911.1 AE005174.2 60.37 323 128 1 2887886 2888209 4117005 4117330 2.2e+02 67.00

NC_000911.1 AE005174.2 64.19 215 77 0 3505100 3505315 2231875 2232089 1.7e+02 61.00

NC_000911.1 AE005174.2 79.79 94 19 0 2598253 2598347 4208942 4209035 1.5e+02 56.00

NC_000911.1 AE005174.2 85.71 70 10 1 633308 633377 4845001 4845073 1.4e+02 50.00

NC_000911.1 AE005174.2 92.98 57 4 0 320590 320647 4517677 4517733 1.4e+02 49.00

NC_000911.1 AE005174.2 80.25 81 16 0 2087795 2087875 776366 776446 1.4e+02 49.00

NC_000911.1 AE005174.2 72.38 105 29 0 139872 139977 4352722 4352826 1.3e+02 47.00

NC_000911.1 AE005174.2 73.61 72 19 0 2708477 2708548 2059415 2059486 1e+02 34.00

NC_000911.1 AE005174.2 86.96 46 6 0 3224161 3224207 4225817 4225862 96 34.00

NC_000911.1 AE005174.2 72.46 69 19 0 1082309 1082377 777015 777083 90 31.00

NC_000911.1 AE005174.2 79.17 48 10 0 2411957 2412005 4153597 4153644 79 28.00

NC_000911.1 AE005174.2 80.00 45 9 0 2791704 2791748 3737114 3737158 73 27.00

NC_000911.1 AE005174.2 68.33 60 19 0 1082308 1082368 2421426 2421485 68 22.00

Results are shown in tab-delimited format with following headers:
Genome1 <name of FNA file of genome1>
Genome2 <name of FNA file of genome2>
ANI(1->2) <Average nucleotide idenitity of the first genome to the second>
ANI(2->1) <Average nucleotide idenitity of the second genome to the first>
AF(1->2) <Alignment Fraction of the first genome to the second>
AF(2->1) <Alignment Fraction of the second genome to the first>
 

 

オンラインで計算するにはgANIにアクセスする。

https://ani.jgi-psf.org/html/calc.php?

f:id:kazumaxneo:20180910205045p:plain

 出力。

f:id:kazumaxneo:20180910214058p:plain

グループ

JGI MiSI - Clusters

 

グループをネットワークで可視化したもの

https://ani.jgi-psf.org/html/clusters.php?page=cliqueGroups

f:id:kazumaxneo:20180911235033p:plain

 

引用
Microbial species delineation using whole genome sequences
Varghese NJ, Mukherjee S, Ivanova N, Konstantinidis KT, Mavrommatis K, Kyrpides NC, Pati A

Nucleic Acids Res. 2015 Aug 18;43(14):6761-71. 

 

ANItoolsをwebで実行できるANItools web

 

 単離されたバクテリア株の迅速かつ正確な分類は、医療微生物学、特に全国的または全世界的な広がりの脅威を伴う感染症発症中の最も重要な課題である(論文より ref.1)。しかし、現在の分類方法はすべて、表現型の類似性や化学的性質に基づく方法のみならず、フラグメントのヌクレオチド配列(16Sおよびマルチシーケンス配列タイピング[MLST])に基づく現代の遺伝的方法も解像度の点で欠点を有する(ref.3-5)。closely relatedな種(97%以上の類似性)を区別するためには、16S rRNAの分子構造は保存され過ぎている(ref.6-8)。さらに、原核生物の早期分類は、表現型の類似性および化学的特性のみに基づいており、温度およびpHなどの環境因子によってある程度影響され、偏りの可能性がある(ref.4)。 16S rRNAおよびMLST法を用いた分類はまた、1つ以上のシーケンシングエラーによってバイアスがかかる可能性がある。

 最近では、任意の2つの株の間で共有される全ての配列のペアワイズ比較から計算された平均ヌクレオチド同一性(ANI)が、細菌種の定義および分類のための新しい測定基準として提案されている。 2005年には、Pro. Konstantinidisは最初に70の関連した種を評価し、2種間の共通遺伝子のANIが株間の遺伝的関連性を比較するための堅牢な手段であることを見出した。 〜94%のANI値は、現在の種定義(ref.9-11)の従来の70%のDNA-DNA再会合(ハイブリダイゼーション)基準に対応することが示された。 2012年にはさらに、ChanがテストケースとしてAcinetobacter属の38株を用いて、ANI結果がコアゲノム系統学および伝統的なアプローチと合致し、既存の分類学(ref.12)と適合することを証明した。著者らの以前の研究(ref.2)では、1226菌株における2つのゲノム比較の正確なANI値を計算し、列挙した。これは、ANIに基づく種の分類がNCBIバクテリア分類学と優れた一致を示すことを示している。この研究は、ANIがバクテリア分類学のために有用であることを証明し、既存および新規のバクテリア種の定義のための強力な候補方法を表している(ref.2)。

 他の方法と比較して、2つの株間の全ゲノム比較に基づくANI分析は、より高い解像度を有し、配列選択およびエラーによって引き起こされるバイアスを回避することができる。 2つのclosely relatedなバクテリア種であっても、ゲノムレベルでのDNA発散に基づいて区別することができ、シーケンシングリードのカバレッジの助けを借りて、1つまたは少数のシーケンシングエラーを容易に調整することができる(ref.2)。 ANItoolsに加えて、ANI 値計算には他のプログラムも用意されている (JSpecies) (ref.6,13)。しかし、以前のバージョンのANItoolsを使用するには、パラメータ調整に加えて、パーソナルコンピュータ上のいくつかの追加プログラム(BLASTやHmmerなど)のインストールと同様にインストールが必要である。さらに、何千ものバクテリアゲノムのために現在ANIデータベースを一般に公開することはできない。 JSpeciesWS(ref.13)はいくつかのバクテリア株間のANI値を計算するためのWebバージョンもサポートしているが、系統数制限(最大15ゲノム)は属または種レベルでANIマトリックスを得る可能性を妨げ、また、同じ属または種の系統間の関係を図式的に示すこともできない。 ANItoolsの他に、JSpecies(ref.6,13)という別のANI値計算プログラムがある。しかし、両方のツールは期待どおり完璧ではない。 ANTIoolsは、BLASTやHmmerなどの特定のAdd-Inプログラムも、その間に常に関連するパラメータの設定や調整作業を伴い、パーソナルコンピュータまたはサーバーにローカルにインストールする必要がある。それは、ITやバイオインフォマティクスについての背景がないユーザーにとってあまり親切ではないことを意味する。 JSpeciesWS(ref.13)に関して、最初のツールはANIを計算するためにオンラインで使用でき、パラメータのセットアップや調整作業を一切必要としない。しかしキャパシティ制限(最大15)のために、ユーザは属または種レベルでANIマトリックスを得る機会がなく、その間に多すぎる菌株を分析する必要がある。さらに、同じ属または種の系統間の関係を図式的に示す系統樹の結果は存在しない。

 そのため、上記の結論に沿って現在のツールのすべての欠点を解消するために、最終的に、ANItools 2.0(http://ani.mypathogen.cn/)のWebバージョンをプログラミングした。 ANItools Webバージョンは、ユーザーがオンラインでANI値を直接取得するのに役立ち、以前のLinuxバージョンと比較して調べられたゲノムの数を増加させる。このデータベースに含まれている (2773 strains, 1487 species and 668 genera)任意の2つの系統のANI値が比較できる 。 ANItools webは細菌株の関係を定義するのに有用であり、ゲノムデータを使用するバクテリア種の分類と同定に役立つ。現在利用可能なソフトウェアと比較して、ANItools Webはユーザーの関与を最小限に抑えている。必要なのはゲノム配列のアップロードと属データの選択のみである。入力ゲノムとデータベースのゲノム配列間の比較を自動的に実行し、ANI計算結果のグラフレポートを生成する。

 

ヘルプ

http://ani.mypathogen.cn/docs/ANItools_Help.pdf

 

使い方

ANItools webにアクセスする。

http://ani.mypathogen.cn

f:id:kazumaxneo:20180910200035p:plain

 

 

調べたいバクテリアゲノムのfastaをアップロードする。例えばE.coli O157-H7のゲノムfastaをアップロードしてみる。

f:id:kazumaxneo:20180910200141p:plain

Browseからローカルのファイルを選択し、Uploadをクリック。

 

データベースの比較したいバクテリアの属などを選択する。 

f:id:kazumaxneo:20180910200412p:plain

 

準備ができたらランする。

f:id:kazumaxneo:20180910200545p:plain

 

しばらく待つとジョブが開始される。

f:id:kazumaxneo:20180910200617p:plain

 

ANI計算結果

f:id:kazumaxneo:20180910213630p:plain

f:id:kazumaxneo:20180910213627p:plain

 

引用
ANItools web: a web tool for fast genome comparison within multiple bacterial strains
Na Han, Yujun Qiang, Wen Zhang

Database (Oxford). 2016; 2016: baw084.

 

Whole-genome sequence comparison as a method for improving bacterial species definition
Zhang W, Du P, Zheng H, Yu W, Wan L, Chen C. 

J Gen Appl Microbiol. 2014;60(2):75-8.

 

参考資料

http://www.wfcc.info/iccc12/presentations/kkonstantinidis.pdf

Pan-genome解析をwebで実行できる PanWeb

 

 次世代シーケンシング(NGS)プラットフォームは、DNAシーケンシングの大きな進歩をもたらした。これは主に、イールドの向上と精度の向上、およびコストの大幅な削減によるものである[論文より ref.1,2]。 NGS技術のために、オンラインゲノムデータベース(https://gold.jgi.doe.gov/)などの公開データベースにdepositされたcomplete genomesの数が指数関数的に増加している。利用可能なゲノム、特に原核生物のゲノムが多数存在することでpan-genome解析などの比較分析が可能になった[ref.3,4]。

 生物の遺伝的レパートリーの比較は、生物工学、生物医学および環境に関する遺伝子発見を助けることができる[ref.5]。このアプローチは、進化的イベントの発生をチャート化し、系統発生的関係を確立することができる[ref.6]。比較ゲノミクスを用いて、特定の種に内在するいくつかのゲノム特性を調べることが可能である[ref.7]。Pan-genomeアプローチにより同じ種または属のいくつかの生物のゲノム間の類似性を比較して、病原性生物の病原性メカニズムを解明することができる[ref.8]。

 さらに、比較ゲノミクスは、異なるライフスタイルを有する微生物に対して用い、intracellular pathogensがしばしばreductive evolutionおよびgene lossを受ける遺伝子レパートリーおよびゲノムサイズを相関させるために使用することができる[ref.9]。

 真核生物では、比較ゲノミクスを使用して、異なる種類の疾患(心臓血管、視覚、聴覚、内分泌および骨疾患)に関連するヒト遺伝子のホモログをショウジョウバエ(Drosophila melanogaster)などのモデル生物において同定し、遺伝子治療のテストに使うことができる[ref.10,5]。

 pan-genomics解析を実行するために利用可能なソフトウェアツールの中で、Panseq(Pan-Genome Sequence Analysis Program)[ref.11]は、ユーザ定義パラメータに基づいてゲノム配列の集合からコア領域およびアクセサリ領域を決定したり、SNP(single-nucleotide polymorphisms)を検出することができる。しかしPanseqはpan-genomicプロファイル解析やfunctional enrichment解析を行うことができない[ref.11,12]。 PGAT(Prokaryotic Genome Analysis Tool)[ref.13]は、多くの微生物ゲノム間の遺伝的な内容を比較するために設計されたウェブデータベースである。しかしPGATは現在、そのデータベース内の種については限られた結果しか提供しておらず、ユーザー提供のゲノムデータ解析もサポートしていない[ref.13,12]。最後に、PGAP(Pan-Genome Analysis Pipeline)(紹介)[ref.12]は、いくつかの機能モデルを統合し、病原性メカニズムを発見し、流行を防除するためのバクテリアの進化の歴史を研究するために使用することができる。しかし、PGAPはPerlスクリプト言語を使用して開発されており、パイプラインをインストール、設定、パラメータ化する必要があり、計算処理の経験がないユーザは使用することが難しくなる。さらに、その出力ファイルは複雑で解釈が難しい[ref.14]。この研究では、分析を簡素化し、結果の解釈に役立つ新しいグラフィックスを組み込むPGAPのためのグラフィカルインターフェイスを提供するWebプラットフォームの開発について説明する。

Implementation

PanWebはプログラミング言語PHPhttp://php.net/)とR(https://www.r-project.org/)を使用して実装された。PGAPの出力結果を処理し、グラフを生成してデータ解釈を容易にする。さらに、サイト構造をマークするためのHTML(Hypertext Markup Language)5、サイトのスタイルと外観を定義するCSSバージョン3、JavaScriptプログラミング言語など、他のテクノロジを使用した。

Pan-genome解析

 PGAPではアクセサリゲノム、セントラルゲノムおよび種特異的領域を同定し、さらにオーソロガスおよびパラロガス遺伝子を同定することができる。機能遺伝子のクラスター分析は、複数のゲノム間のオルソログおよびパラログの検索における不可欠なステップの1つである。この目的のために、PGAPは2つの方法、すなわちMP(MultiParanoid)とGF(gene family)を使用する[ref.12 pubmed]。MP法は、2つのアルゴリズム、すなわち、InParanoidとMultiParanoid [ref.15]を使用する。InParanoidアルゴリズムは、まず各株ペア間のオーソログ検索を行い、各株のパラログは、BLASTを用いてゲノム中のホモログを検索することによって検討される。MultiParanoidアルゴリズムは、複数の系統の遺伝子クラスター探索を行う[ref.12]。GF法は、タンパク質配列のmixuture中でBLASTALLを使用し、クラスタリングプロセスはMCLアルゴリズムによって実行される[ref.16,17]。BLASTの最小スコア値とe valueは、両方の方法(GFとMP)で40と10^-8である。 Blast検索に費やされる時間を短縮するために、PGAPスクリプトはDiamond  [ref.18]を実行するように変更された。これはBlastより高速で、PGAPパイプラインでは主要な低速ステップの1つだった。

パイプライン

PanWebはEMBL形式で各ゲノムのアノテーションファイルを入力として受け取り、WebアプリケーションはNUCおよびPEPファイルを自動的に生成する。NUCおよびPEPファイルにはそれぞれヌクレオチド形式とタンパク質形式のコードDNA配列(CDS)が含まれている。各FUNCファイルも生成される。これらのファイルは、その後PGAP [ref.12]によって処理され流。分析のバイアスを避けるために、すべてのゲノムに対して同じアノテーションツールを使用することによって、ゲノムの予測およびアノテーションを標準化することが重要になる。 PGAPを実行すると、PanWebは結果を処理し、パイプラインによって生成された出力ファイルから得られたグラフを表示する。

 

Introduction

http://www.computationalbiology.ufpa.br/panweb/instructions.html

 

使い方

 テストデータを解析してみる。 

http://www.computationalbiology.ufpa.br/panweb/analysistestfiles.php

f:id:kazumaxneo:20180909111627p:plain

 E-value、Identity、Coverage、Analysisタイプなどを指定する。最後にメールアドレスを記載したらAnalyseをクリックする。Exampleはすぐ終わるが、普通はジョブ完了のメールが届くまで数分~数十分かかる。

 

Boxplot

f:id:kazumaxneo:20180909111918p:plain

Pie

4ゲノム全てに見つかるオーソロガスな遺伝子が86.6%を占めた。

f:id:kazumaxneo:20180909111949p:plain

Barplot

f:id:kazumaxneo:20180909111956p:plain

 

 

 

手持ちのデータを 解析するには、genbank、embl形式のファイルをアップロードする(*1)。またはすでに登録されているバクテリアのゲノムデータを比較する。E.coliのK-12系統4株と0-157 H7を1株比較してみる。Analysisのタブに戻り、左下から生物種を選択する。

f:id:kazumaxneo:20180909131507p:plain

繰り返すと複数登録できる。

f:id:kazumaxneo:20180909124225p:plain

パラメータをセットして分析を開始する。

 

ユニークなORFの数。アウトグループとして加えたNC_011353が圧倒的に多かった。

f:id:kazumaxneo:20180909123822p:plain

シェアされている遺伝子数の表

f:id:kazumaxneo:20180909123854p:plain

 

引用

PanWeb: A web interface for pan-genomic analysis.
Pantoja Y, Pinheiro K, Veras A, Araújo F, Lopes de Sousa A, Guimarães LC, Silva A, Ramos RTJ

PLoS One. 2017 May 24;12(5):e0178154

 

*1

NCBI Refseqなどからgenbankをダウンロードして使ってもいい。ただしテスト時はエラーを起こした。genbankの自動ダウンロードはこの辺りの質問を参照。

https://www.biostars.org/p/223758/

 

PGAP



パンゲノム解析ツール PGAP

2019 7/6 誤字修正

 

 DNAシーケンシング技術の急速な発展に伴い、「Ten Thousand Microbial Genomes Project」や「NIH Human Microbiome Project(HMP)」(Peterson et al、2009)など多くの大規模な微生物ゲノムプロジェクトが処理されている。バクテリア全ゲノムシーケンスの蓄積はまた、生物学者に以前よりも大きなスケールで進化仮説を探索し、試験する機会を与える。 2005年、Tettelinらは新しい概念「pan-genome」を導入した(Tettelin et al、2005)。その後しばらくして、pan-genomeは、肺炎連鎖球菌(Hiller et al、2007)、インフルエンザ菌(Haemophilus influenzae)(Hogg et al、2007)、大腸菌(Escherichia coli)(Rasko et al、 2008年)などがの解析の洞察に使われ始めた。進化への洞察に加え、pan-genomeはレジオネラニューモフィラいくつかの病原体の株特異的ビルレンス因子を検出するために広く使用されている(D'Auria et al、2010)。コアゲノムの機能が変化する遺伝子をスキャンし、流行病の病原体を調べることで(Bayjanov et al、2010; Holt et al、2008)、逆ワクチン薬理ゲノミクス(Serruto et al、 2009)(wiki)により病原性細菌のワクチンを開発することの助けにもなる。

 1つのバクテリア集団についてpan-genome解析を可能な限り容易にするために、バクテリアパンゲノム解析のための高効率ツールを開発する大きな必要性がある。pan-genome解析のために使えるツールは、これまでPanseq(Laing et al、2010)およびPGAT(Brittnacher et al。、2011)だけだった(論文執筆時点)。 Panseqは、ゲノム配列の中の「コア」領域と「アクセサリー」領域を抽出し、コア領域間でSNPを検出することに優れている。しかし、それは与えられた株のpan-genomeプロファイルを提示し、進化の歴史を複数の物質で追跡し、機能的遺伝子の変異および function enrichmentを指摘する能力には不十分である。PGATは、ウェブベースのデータベースとして、オルソログのアサインメント、遺伝子内容の照会、配列の多型および代謝経路情報を統合している。しかしながら、これまでのところ、限られた種の分析結果のみをデータベースに提供し、ユーザからのゲノムデータを分析することはできない。著者らは、パンゲノム解析パイプライン(PGAP)と呼ばれる新しいスタンドアローンプログラムを開発した。PGAPは複数の機能モデルを統合し、バクテリアの進化の歴史を研究し、病原性のメカニズムを発見し、流行を予防し、防除するために使用できる。

 

 

f:id:kazumaxneo:20180908181921j:plain

Supplementaryより。

 

f:id:kazumaxneo:20180908181909j:plain

同名のツール NCBI prokaryotic genome annotation pipeline (PGAP) とは異なります。

 

 

インストール

依存

  • BLAST (2.2.12 or higher).
  • mafft
  • dnaml, dnadist, neighbor, seqboot, consense in PHYLIP (version 3.69)
  • mcl

SourceForge

https://sourceforge.net/projects/pgap/

ここではdockerイメージを使う。 

docker pull kastman/pgap

 > docker run --rm kastman/pgap

$ docker run --rm kastman/pgap

 

====== Pan-Genome Analysis Pipeline (PGAP) ======

                   Version 1.2.1

 

Usage:   perl PGAP.pl [Options]

 

Options: 

  --strains    String    Input strains nicknames, and join them with '+', for example: A+B+C

  --input      String    Input data directory 

  --output     String    Result output directory

 

  --cluster              Run homologous gene clustering

  --pangenome            Run pan-genome analysis

  --variation            Run homologous clusters variation analysis

  --evolution            Run evolution analysis

  --function             Run Function analysis

 

  --method     String    GF for GeneFamily method,  and MP for MultiParanoid method

                           for GF: fast, but not very accurate

                               evalue, score, indentity, coverage are employed

                           for MP: slow, but more accurate

                               score, coverage, local, global are employed

  --thread     Int       Number of processors to use in blastall. [default:1]

  --score      Int       Minimum score in blastall. [default:40]

  --evalue     Decimal   Maximal E-value in blastall. [default:1e-10]

  --coverage   Decimal   Minimum alignment coverage for two homologous proteins. [default:0.5]

  --local      Decimal   Minimum local alignment overlap in MP method. [default:0.25]

  --global     Decimal   Minimum global alignment overlap in MP method. [default:0.5]

  --identity   Decimal   Minimum alignment indentity for two homologous proteins. [default:0.5]

  --bootstrap  Int       Bootstrap times for phylogenetics tree. [default:1]

 

  --h or help            Display this message

 

 

ラン

解析するには、対象生物種それぞれについて、3種類のファイルが 必要(.pep for protein sequences, .nuc for nucleotide sequences and .function for annotation file. For each strain, these three type files should have the same prefix, )。,nucと.pepはそれぞれ遺伝子の塩基配列アミノ酸配列、.functionは1情報につき1行3列で情報を示したフォーマット。以下のように記載する(マニュアルより)。

f:id:kazumaxneo:20180908193758j:plain

1列目は.nucと.pepのヘッダー名と同じにする。2列目はCOGのID(無いなら"-")。

 

ジョブの開始

seq_dir=/path/to/input
out_dir=/path/to/output

docker run --rm\
-v "${seq_dir}":/input -v "${out_dir}":/output -w /pgap \
kastman/pgap perl ./PGAP.pl --strains $strains \
--input /input --output /output \
--cluster --pangenome --variation --evolution \
--function --method MP --thread 1

kastmanさんのnoteより: * the -w /pgap default working directory - this is important for PGAP.pl to correctly find sub-modules in the /pgap directory.

出力の詳細はPDFマニュアルを確認してください。

 

引用
PGAP: pan-genomes analysis pipeline
Zhao Y, Wu J, Yang J, Sun S, Xiao J, Yu J

Bioinformatics. 2012 Feb 1;28(3):416-8. 

Procaryotesの自動アセンブリパイプライン Mypro

Pos 

 最近の全ゲノムシーケンシング(WGS)技術のコストの低下は、様々な原核生物のシーケンシングの増加をもたらした。典型的なゲノミクスプロジェクトでは、データマイニングの前にシーケンシングリードを処理する必要がある(Hasman et al、2014; Rhoads et al、2014)。このプロセスには、品質管理チェックや前処理対策、新規アセンブリおよび/またはリファレンスベースアセンブリ、マニュアル改良の有無にかかわらず自動アノテーション、ギャップ充填などの方法によるゲノム品質の向上が含まれる。ゲノムアセンブリアノテーションおよび(Koren et al、2014; Magoc et al、2013; Seemann、2014; Swain et al、2012)を達成するために、現在、様々なソフトウェアが利用可能である。しかし、適切なバイオインフォマティクスの訓練と支援がない科学者や小規模な研究所の能力は、有意義な結果を得るためにゲノムインフォマティクスツールの実行を制限するかもしれない(Nocq et al、2013)。

 ここでは、最小限のプログラミング技術しか必要としない原核生物ゲノム用のユーザーフレンドリーなゲノミクスソフトウェアパイプラインであるMyProを紹介する。パイプラインは、品質管理と前処理ツール、デノボアセンブラ、contigインテグレーションソフトウェア、アノテーションとリファレンスベースアセンブリのツールで構成されている。

 使いやすさを向上させるため、MyProは自動でデノボアセンブリアセンブリ統合(genome assembly reconciliation toolsを使う)、ゲノムアノテーションを行うよう設計されている。これは、3つの主要なモジュールで構成されている:アセンブリ、統合、アノテーション(論文 図1)。 VelvetOptimiser 2.2.5(Zerbino and Birney、2008)、Edena V3.131028(Hernandez et al、2008)、Abyss 1.5.2(Simpson et al、2009)、SOAPdenovo 2.01(Luo et al 、2012)およびSPAdes 3.1.1(Bankevich et al、2012)はBioLinux 8(Field et al、2006)(pubmed)にインストールおよび設定されており、MyProのAssembleモジュール内で複数のアセンブリを作成できる。

 5つのアセンブラのほとんどは、de-Bruijnグラフベースのアルゴリズムを使用している。これらのアルゴリズムは、k-merパラメータに大きく依存している。 Fastq形式のシーケンシングデータを入力した後、VelvetOptimiserを使用して最初に0.6から0.9までのk-merのパラメータスイープを実行する。最適なk-merは、最も高いN50コンティグの長さに対して自動的に選択され、このk-merはAbyssおよびSOAPdenovoアセンブリに使用される。ペアエンドリードの場合、VelvetOptimiserはインサートサイズとその標準偏差を推定し、後で使用するためにMyProに保存する。入力リード長に応じて、SPAdesアセンブリは、リード長150 bpでk-mer長21, 33, 55が設定され、150〜250bpのリードに対しては21, 33, 55, 77のk-mer長が設定される。250bp以上のリード長では21,33,55,77,99および127のk-mer長が設定される。 Edenaはオーバラップレイアウトに基づいたアセンブリを行うため、k-merの代わりに最小オーバラップサイズ(m)を必要とする。 mの値はリード長の半分に設定され、リード長まで徐々に増加する。最適なmは、最も高いN50コンティグの長さに対して自動的に選択され、対応するアセンブリはEdenaアセンブリとして保持される。

 複数のアセンブリが取得されると、必要に応じて事前計算されたリード情報(例えば、ペアエンドリードの最小および最大インサートサイズ)を用い、SOAP2(Li et al、2009)によって入力リードがマッピングされる。アライメント率はアセンブリを評価するために使用される。高品質のコンティグを提供し、可能性のある汚染配列を除去するために、100アライメント未満のコンティグは除去され、コンティグ数、最も長いコンティグ長、N50および全ゲノムサイズなどの基本的なアセンブリ統計がレポートされる。 MyPro(Integrate.py)は、CISA 1.3(Lin and Liao、2013)を使用してAssembleモジュールによって生成されたアセンブリを統合することにより優れたアセンブリを提供する。自動化されたパイプラインでは、アセンブリ統計に基づくワーストのアセンブリは、Integrateモジュールを実行する前に削除される。ユーザーは5つのアセンブリの中から選択することも、他のアセンブリを統合に追加することもできる。

 アノテーションモジュールでは、迅速なアノテーションを実装するProkka(シーマン2014年)のコアデータベース制限を、高品質のリファレンスゲノムデータベースと最新のSwiss- Protデータベース(Tatusova et al、2013)(注: 論文当時のデータ)を使用して増強した。 MyPro(Annotate.py)はCISAアセンブリに自動的にアノテーションを付けるが、ユーザーはフォルダに配置された他のアセンブリアノテーションを付けることもできる。

 3つの主要なモジュールに加えて、MyProは、前処理、exploration(以後、探索)、およびpost-assemblyのための機能を提供する。特定のゲノムサイズの前処理では、MyPro(Preprocess.py)はクオリティトリミングを実行し、入力リードを所望のカバレッジデプス(デフォルト100×)にサブサンプルする。アセンブリの探索のために、MyProはリードのアライメントをBAM形式で出力し、Tablet(Milne et al、2012)によるアセンブリの外観検査を可能にする。closely relatedなリファレンスゲノムが利用可能である場合、r2cat(Husemann and Stoye、2010)を使い、関連したリファレンスを用いコンティグを順序付け、不一致コンティグを見つける。MyPro(Postassemble.py)は、コンティグとリードの両方をセットを使用してポストアセンブルを行う。リファレンスガイドのオーダーで重複するコンティグ(15 bp以上)はマージされる。マージされたコンティグ上のSOAP2によってマッピングされないリードは、Edenaにより様々なmを使用してアセンブルされたローカルアセンブリは、2つのコンティグをブリッジするために使用される。

  MyProを検証するために、著者らは、GC-コンテンツが両極端な3つの菌種:GC-poor Staphylococcus aureus (33% GC)、GC-rich Rhodobacter sphaeroides 2.4.1 (69% GC) 、 Escherichia coli MG1655 (GC 51%)でゲノムアセンブリを行った。公開されているシーケンシングリードをダウンロードし、MyProによって分析した(詳細は補足データを参照)。 QUAST 2.3(Gurevich et al、2013)によって評価されたアセンブリ結果は、MyProがアセンブリ精度とアノテーション精度の連続性に関して高品質の注釈付きアセンブリを生成したことを示している。 QUASTはまた、最近シーケンシングされた8つの口内連鎖球菌株についてMyProによるアセンブリを評価するためにも使用された(論文 表1)。これらの株については、データの質を維持するために500bp未満のコンティグは捨てられた。MyPro(AutoRun.py)(表2)によってデノボアセンブリされた10の口内ストレプトコッカス株についても評価された。利用可能なリファレンスゲノムを有する口内連鎖球菌株およびリファレンスゲノムがない株に関して、MyProはデノボアセンブリツールおよびコンティグインテグレーションソフトウェア(CISA)より一貫して優れたN50またはコンティグ数を示した。 (以下略)

 

Supplementary Documentリンク

Myproに関するツイート。

 

 

インストール

SourceFourge(マニュアルPDFあり)

sb2nhri - Browse /MyPro at SourceForge.net

 

Mypro.ovaファイルダウンロードリンク

http://sb.nhri.org.tw/MyPro/index.html

> Preprocess.py 

manager@bl8vbox[manager] Preprocess.py                                                                          [10:02AM]

Usage:

Preprocess.py -read1 [/path/read1.fastq] -read2 [/path/read2.fastq] -g [genomesize] -d [depth] -l [limit] -r [reduce] -min [min length]

-g genomesize. [default:9000000]

-d Int. depth. [default:100]

-l Float. cutoff value for quality trimming. [default:0.01]

-r Int. step size for reducing target length. [default: 10] 

-min Int.  Discard reads below length.(defult 50)

-h Help.

 > AutoRun.py    

manager@bl8vbox[manager] AutoRun.py                                                                                                   [10:20AM]

Usage:

AutoRun.py [projectName] -read1 [/path/read1.fastq] -read2 [/path/read2.fastq] [options]

-p Str. parameters for annotate.py

-thr Int. To remove ctg with lower than thr depth.

-h Help.

Postassemble.py 

manager@bl8vbox[manager] Postassemble.py                              [10:21AM]

Usage:

Postassemble.py -o [filename] -u [filename] -read1 [/path/read1.fastq] -read2 [/path/read2.fastq] -m mean -s std

-o ordered.fa.

-u unmatched.fa

-m mean of insert size. default:200

-s std of insert size. default:50

-h Help.

 

ラン

Virtual PCovaファイルをインポートする。起動前に、設定からCPUやメモリなどをチューニングしておく。ここではメモリ8GB、CPU数は写真では12となっているが4とした。

f:id:kazumaxneo:20180907194304j:plain

立ち上げ、端末エミュレータを起動する。

f:id:kazumaxneo:20180907194031j:plain

DockにはTabletとfastqcのショートカットも準備されている。

 

端末にコマンドを打ち込むことで実行する。

1、データをエミュレータにコピーする。

 

2、端末を立ち上げ、作業ディレクトリを作成する。

#作業ディレクトリtestを~(home/)に作成
mkdir ~/test

#testディレクトリに移動
cd ~/test

 

3、前処理。ペアエンドのファイルpair1.fq とpair2.fq のクオリティトリミングを行う。

#ゲノムサイズは5Mとする
Preprocess.py -read1 pair1.fq -read2 pair2.fq -g 5000000

ジョブが終わるとtrimmed_pair1.fqとtrimmed_pair2.fqができる。

 

4、アセンブリアセンブリ結果のInteggration

AutoRun.py test1 -read1 pair1.fq -read2 pair2.fq

 アセンブリしたfastaと似たリファレンスを加え、dot plot解析を行う。 OSをエミュレーションしてランしている上に重たいパイプラインを走らせるので、ジョブが終わるまでにはかなりの時間がかかる(*1)。

ランが終わると複数のサブディレクトリができる。

f:id:kazumaxneo:20180910110514j:plain

 

--Optional step--

5、r2catを立ち上げ、同一、またはよく似たリファレンスのゲノムと比較する。

#ゲノムサイズは5Mとする
java jar ~/Desktop/r2cat.jar

 

結果をソートする。

f:id:kazumaxneo:20180910113941j:plain

unmatchとsortしたcontigをそれぞれ出力する。

f:id:kazumaxneo:20180910114231j:plain

 

6、ポストアセンブリ

Postassemble.py -o cisa.ordered.fa -u unmatched.fa  -read1 pair1.fq -read2 pair2.fq

ジョブが終わると、ランディレクトリにBridged.xxx,faファイルが出力される。

 

引用

MyPro: A seamless pipeline for automated prokaryotic genome assembly and annotation
Yu-Chieh Liao, Hsin-Hung Lin, Amarpreet Sabharwal, Elaine M. Haase, Frank A. Scannapiecob

J Microbiol Methods. 2015 Jun; 113: 72–74

 

*1

mac proで4コア割り当てランしたところ、4M 50xリードの解析に2日ほどかかった。