macでインフォマティクス

macでインフォマティクス

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

fastqのクオリティ分析を行う Quack

 

 ハイスループットDNAシーケンシングツールによって生成されたデータの品質は、そのデータが生物学的発見にどの程度役立つかを判断するために迅速に評価されなければならない。データセットのサイズがますます大きくなり、迅速な品質評価が重要になっているため、シーケンシングデータを解析するツールは、解釈しやすいグラフィックを迅速に作成する必要がある。Quackは、シーケンス技術に依存しない方法で、FASTQファイルから情報量の多いビジュアライゼーションを、他の一般的に利用可能な品質保証ツールをはるかに上回る速度で生成することで、これらの問題に対処する。

 

インストール

リリースからlinux向けバイナリをダウンロードしてテストした。

 依存

  • zlib
  • klib (pulled by the submodule update below)

Github

git clone https://github.com/IGBB/quack.git
cd quack/
make && make test

./linux.quack

$ ./linux.quack 

Usage: quack [OPTION...]

quack -- A FASTQ quality assessment tool

 

  -1, --forward file.1.fq.gz      Forward strand

  -2, --reverse file.2.fq.gz      Reverse strand

  -a, --adapters adapters.fa.gz   (Optional) Adapters file

  -n, --name NAME                 (Optional) Display in output

  -u, --unpaired unpaired.fq.gz   Data (only use with -u)

  -?, --help                      Give this help list

      --usage                     (use alone)

  -V, --version                   Print program version (use alone)

Report bugs to <thrash@igbb.msstate.edu>.

Usage: quack [OPTION...]

Try `quack --help' or `quack --usage' for more information.

 

 

実行方法

fastqを指定する。アダプターのFASTAファイルがあるなら-aで指定する。

#ペアエンド
quack -1 input_1.fq.gz -2 input_2.fq.gz -n sample_name > output.svg

#シングル
quack -u input.fq.gz > output.svg

f:id:kazumaxneo:20200915001706p:plain



一番上のパネル - 配列の各列に含まれる各ヌクレオチドの割合

f:id:kazumaxneo:20200915002100p:plain

 

中央のパネル - 配列のポジションごとのクオリティ分布を示すヒートマップと、配列全体のクオリティスコアの平均値を示すライン

f:id:kazumaxneo:20200915002529p:plain
中央両サイドのパネルは中央のパネルのクオリティ割合を示すスコア分布グラフ。グラフの左が100%、右が0%。一番下のパネルは側面のパネルはリードの長さの分布示すグラフ。

 

引用

Quack: A quality assurance tool for high throughput sequence data

Adam Thrash, Mark Arick 2nd, Daniel G Peterson

Anal Biochem. 2018 May 1;548:38-43

 

 

 

GenbankアクセッションIDからtaxonomyを返す acc2tax

 

タイトルの通りのツール。

 

インストール

macos標準のclangでビルドした。

Github

git clone https://github.com/richardmleggett/acc2tax.git
cd acc2tax/
cc -o acc2tax acc2tax.c

#パスの通ったディレクトリにコピー
cp acc2tax /usr/locasl/bin/

acc2tax -h

$ acc2tax -h

 

acc2tax v0.6

 

Bugs/comments: richard.leggett@tgac.ac.uk

 

Provide batch taxonomy information for Genbank IDs or Accessions.

 

Options:

    [-h | --help]       This help screen.

    [-a | --accession]  Query is accession IDs [default].

    [-c | --column]     1-based column number of ID in input file (default 1).

    [-d | --database]   Directory containing NCBI taxonomy files.

    [-e | --entries]    Max GI entries (default 1050000000).

    [-g | --gi]         Query is Genbank IDs.

    [-i | --input]      File of IDs (GI or Accession), one per line.

    [-k | --keep]       Copy columns from input to output file, then append taxonomy as new column.

    [-n | --nucleotide] Query IDs are nucleotide [default].

    [-o | --output]     Filename of output file.

    [-p | --protein]    Query IDs are protein.

    [-s | --strip]      Strip version from input acession IDs (ie. everything after .)

  

データベースの準備

taxdumpと.accession2taxid をダウンロードする。

ftp://ftp.ncbi.nih.gov/pub/taxonomy/accession2taxid/

f:id:kazumaxneo:20200908134727p:plain

1、/pub/taxonomy/accession2taxid/から次のファイルをダウンロードする。

nucleotide

  • nucl_gb.accession2taxid
  • nucl_wgs.accession2taxid
  • dead_nucl.accession2taxid
  • dead_wgs.accession2taxid

protein

  • prot.accession2taxid
  • dead_prot.accession2taxid 

2、1つ上の階層にはtaxdumpファイルもあるので、こちらもダウンロードする。

taxdumpを解凍したところ。解凍したフォルダのnodes.dmp、names.dmp、gi_taxid_nucl.dmp、gi_taxid_prot.dmpが必要。

f:id:kazumaxneo:20200908141051p:plain

 

3、ダウンロードしたファイルをマージしてソートする。

nucleptide ( 数十GBあるので注意)

cat nucl_gb.accession2taxid \
nucl_wgs.accession2taxid \
dead_nucl.accession2taxid \
dead_wgs.accession2taxid \
dead_nucl.accession2taxid dead_wgs.accession2taxid| sort > acc2tax_nucl_all.txt

上記のうち、nucl_gss.accession2taxidとnucl_est.accession2taxidはNucleotide database に統合が進んでいて(*1)、現在NCBIFTPサイトで利用できなくなっている。

protein

cat prot.accession2taxid dead_prot.accession2taxid | sort > acc2tax_prot_all.txt

解凍したtaxdump/に上記のファイルを移動させる。

mv acc2tax_nucl_all.txt taxdump/
mv acc2tax_prot_all.txt taxdump/

 

実行方法

 nucleotideのgenbank accesion IDテキストを指定する。

acc2tax -a genbankID.txt -c 1 -d taxdump/ -o putput

 

引用

 GitHub - richardmleggett/acc2tax: Tool for quick offline batch conversion of Genbank IDs or accessions to taxonomy strings

 

1、Upcoming Changes to EST and GSS Databases

https://ncbiinsights.ncbi.nlm.nih.gov/2018/07/30/upcoming-changes-est-gss-databases/

 

2、new NCBI taxonomy repository

https://github.com/DerrickWood/kraken2/issues/101

 

参考

NCBI taxonomy databases

https://www.uppmax.uu.se/resurser/databases/ncbi-taxonomy-databases/

(ESTとGSSレコードの廃止)

メタゲノムのビニングされた真核生物由来コンティグの品質を調べる EukCC

 

 微生物のDNAは日常的に抽出され、配列決定され、ゲノムにアセンブリされている。回収されたゲノムの品質を推定することは、不完全なゲノムや汚染されたゲノムが公表されるのを防ぐために非常に重要である。シングルコピーマーカー遺伝子(SCMG)は、新たにアセンブルされたゲノムの品質を推定するために日常的に使用されている。これらの遺伝子はゲノム内に一度しか存在しないと予想されるため、ドラフトゲノム内で発見されたSCMGの数と予想されるマーカー遺伝子の数を比較することで、完全性の推定が可能となり、マーカー遺伝子の追加コピーの存在は汚染の指標として使用することができる。このアプローチは原核生物でも真核生物でも広く受け入れられている[1,2,3,4]。原核生物のゲノムについては、CheckM [4]が完全性と汚染度を推定するために最も広く使用されているツールだが、他のアプローチも使用されており、BUSCO [2,5]やanvi'o [6]によって原核生物のSCMGのセットも提供されている。CheckMは、ゲノムのクレードを識別するために普遍的なSCMGの初期セットを使用し、その後、品質を推定するためにクレード固有のセットを使用する。

 同様に、SCMGは真核生物の単離ゲノムの評価にも使用されている。CEGMA [1]は、6つのモデル生物から同定された240のユニバーサルシングルコピーマーカー遺伝子を用いてゲノムの完全性を評価している。BUSCOはCEGMAに取って代わり、真核生物の単一普遍マーカー遺伝子セットに加えて、複数の真核生物および原核生物のクレードのマーカー遺伝子セットをキュレーションしたことが大きな進歩であった。しかし、BUSCO(バージョン3.1)は真核生物、原生生物、植物、真菌の完全性を推定するためのセットを提供しているが、ゲノムの品質を評価する際にどのセットが最も適しているかは利用者の選択に委ねられている。これらのより普遍的なアプローチの他にも、特定の真核生物のクレードに焦点を当てたアプローチも存在する。FGMPは、SCMGと真菌ゲノム内の高度に保存された領域の両方を利用して、真菌単離ゲノムの品質を推定している[7]。原生生物については、anvi'oは、真核生物のゲノムの品質を推定するために、BUSCOの「protist」セットからのプロファイルの数を減らしたものを提供している(未発表)。

 過去20年の間に、海や土壌、ヒトの腸のような宿主関連サイトなどの環境(マイクロバイオームとして知られている)に存在する微小な生物の理解が飛躍的に進歩した。このような知見のほとんどは、メタバーコーディング(マーカー遺伝子の増幅)やメタゲノミクス(ショットガンシーケンス)などの手法を用いて、微生物の集合的な遺伝物質に最新のDNAシーケンス技術を適用したことに由来している。このような配列データの解析から、全微生物の最大99%が培養されていないと考えられている[8]。

 これまでのところ、圧倒的な数のメタバーコーディングおよびメタゲノミクス研究は、サンプル内に存在する細菌に焦点を当ててきた。しかし、ウイルスや真核生物もまた、その数と機能の両面で、生態系の重要なメンバーである[9,10,11,12]。実際、単細胞原生生物と真菌は、世界の微生物バイオマスの約17%を占めると推定されている。微生物の真核生物バイオマスの中では、原生生物として知られる遺伝的に多様な単細胞生物が25%もの割合を占めている[13]。

 利用可能な完全なゲノム配列や完全に近いゲノム配列の数が増加しているにもかかわらず、微生物真核生物の解析を含むメタバーコーディングやメタゲノムアプローチは、原生生物の真の多様性が、現在のゲノム参照データベース(RefSeqやENAなど)に反映されているものよりもはるかに大きいことを実証してきた。例えば、メタバーコーディングシークエンシングに基づく最近の推定では、海洋だけでも15万種の真核生物が存在することが示唆されている[14]が、GenBankにエントリーがある代表種はわずか4551種である(15.Nov 2019)。したがって、マイクロバイオームの機能的役割を完全に理解するためには、これらのまだ解明されていない生物が何であるか、そしてそれらが実行している可能性のある機能的役割を知る必要がある。

 現在、マイクロバイオームの機能を理解するための最良のアプローチの一つは、ショットガンリード(通常は200〜500bpの長さ)を集めて、より長いコンティグ(通常は2000〜50万bpの範囲)を得ることである。これらのコンティグは、完全なタンパク質へのアクセスを提供し、そのタンパク質は周囲の遺伝子の文脈の中で解釈される可能性がある。ここ数年、この種の解析を拡張して、メタゲノムアセンブルゲノム(MAG)と呼ばれる仮説ゲノムを回収することが一般的になってきた。MAGは、単一の生物に由来すると考えられるコンティグをグループ化することによって生成されるが、これはビニングと呼ばれるプロセスである[15]。しかし、ビニング後であっても、生物学的な理由(微生物の豊富さなど)、実験的な理由(シーケンシングの深さなど)、技術的な理由(アルゴリズムなど)の組み合わせにより、MAGの完全性にはばらつきがあり、一般的には断片化されている。さらに、コンティグをビニングするために使用される計算手法は、異なる生物に由来するコンティグの区別に失敗することがあり、キメラゲノム(コンタミネーション)につながることがある。これは、MAGを分析する際に考慮しなければならない問題である[16]。そのため、高品質のMAGの同定と選択を可能にするためには、ゲノムの品質評価がプロセスの重要な部分を占めている。

 ここでは、異なる真核生物のクレードにまたがる現在のアプローチの性能を調査し、真核生物のMAGの品質評価にこのツールを適用することを特に視野に入れ、完成度とコンタミネーションの観点から真核生物のゲノム品質を推定するための教師なし手法であるEukCCについて説明する。

 BUSCOでは、特定のクレードのSCMGセットを作成しているため、より正確な品質推定が可能である。この観察を基に、著者らは原生生物や真菌類の広い範囲をカバーする複数のSCMGセットを定義することを目指した。マーカー遺伝子ライブラリの主な用途は、特徴の乏しいゲノムへの応用であり、そのような場合、de novo予測によって遺伝子が同定される可能性が高いと予想している。そこで、メタゲノミクスのユースケースをより正確に表現するために、GeneMark-ESで生成された遺伝子予測を使用して、すべての真核生物のRefSeq種(両生類や維管束植物に属さない)に(再)アノテーションを行うことを選択した。De novoアプローチで予測可能なマーカー遺伝子の選択にも有利である。GeneMark-ESを選択して適用したのは、このツールが大規模な種の範囲でうまく機能し、一般的にRefSeqのアノテーションベンチマークに近い性能を発揮することを以前に実証したからである。さらに、UniProtKBにリファレンスゲノムとして使用されている真核生物の全種を追加した。結果として得られたタンパク質は、hmmer (version 3.2) [22]を用いて、PANTHER 14.1 [21]のファミリーレベルプロファイルHMMを用いてアノテーションを行った。PANTHERを選択したのは、テスト済みのデータベースの中で、分析したタンパク質のカバレッジが最大であることが示されていること [23] と、PANTHERのプロファイルHMMが構成グローバルドメインではなく、完全長のタンパク質ファミリーをモデル化しているからである。

 パラログの分離を高め、共通ドメインによる局所的な一致を最小限に抑えるために、プロファイル固有のビットスコアしきい値を定義することを目的とした。これを達成するために、分類学的にバランスのとれた種のセットに依存し、各プロファイルについて、シングルコピーマッチの数が最も多いビットスコアしきい値を特定した(「方法」のセクションを参照)。

 その後、クレード固有のSCMGを定義するために、まず、55の広く存在するSCMGを用いて、与えられたゲノムの参照ツリーを構築した(ここから「参照セット」と呼ぶ、「方法」の項を参照)。ツリーの各クレードにおいて、率が98%以上のSCMGをチェックした。次に、前述の閾値に一致するクラッド内に20以上のPANTHERファミリーが見つかった場合には、マーカー遺伝子のセットを定義した。この方法を用いて、木全体で477個のSCMGセットを定義することができた。BUSCOやCEGMAとは対照的に、PANTHERモデルに基づいて真核生物キングダム全体に適用可能なSCMGを特定することはできなかったが、代わりに多くのサブクレードに適用可能なセットを発見した。これは特異性を高めるためには望ましいことだが、明らかな欠点は、どのセットを使用するのが最も適切であるかを知ることである。そこで著者らは、最適なSCMGを選択し、それを用いてゲノムの品質を推定するソフトウェアパッケージであるEukCCを開発した。
 分類学的系統が不明な新規ゲノムに対して最も特異的なSCMGのセットを選択するために、EukCCは、55の広く存在するSCMGの参照セットを使用して、de novoで予測されたタンパク質をアノテーションすることにより、初期の分類学的分類を行う。次に、Pplacer [24]を適用して、参照ツリー内の各マッチを系統的に文脈化する。ツリー内の各配置を追跡すると、EukCCは、SCMGセットがデータベースに定義されている最低共通祖先(LCA)ノードを決定する。

 予想されるように、pplacerは多くの場合、すべての配列を参照ツリーの単純で狭い領域に配置するが、時折、整合性のない、遠くに関連するクレード内に配置されることがある。このような場合、SCMGの単一のセットがすべての場所を網羅しているとは限らない。このようなケースを克服するために、配置の最大の割合をカプセル化したSCMGセットが配置される。このプロセスは、不正確な配置または一貫性のない配置のために外れる配置が発生する場合を克服するが、このアプローチは、新規ゲノムからの参照SCMGとの一致がツリー内に確実に配置できない場合、不正確なSCMGセットを選択する可能性がある。これを制御するために、EukCCは常にセット内のプロファイルの数を報告し、配置位置をプロットするオプションを提供している(図5S)。したがって、参照SCMGの一部しか網羅していないセットが選択された状況では、このMAGのより詳細な分析を行うことができるので、これを行うべきである。(以下省略)

 

Documentation

EukCC — EukCC documentation

 

インストール

ubuntu18.04LTSでcondaの仮想環境を作ってテストした。

依存

Github

#1 ここでは3.7の仮想環境を作成
conda create -n eukcc -y python=3.7
conda activate eukcc
conda install -c bioconda -c conda-forge eukcc -y

#GeneMarkはライセンスの関係でHPより直接ダウンロードする
#2 genemakr ESを入れる準備
apt install -y cpanminus make gcc dialog
#3 perl libray
cpanm inc::Module::Install::DSL Hash::Merge MCE::Mutex FindBin Test::Pod Logger::Simple Parallel::ForkManager.pm YAML
#cpanmで入らない場合 => *1

#4 genemakr ESのダウンロード => link

#5 プログラム本体とキーをダウンロードして、
cd GeneMark_dir
#6 Install the key
cp gm_key ~/.gm_key
#7 test run
bash check_install.bash
#8 パスを通す
export PATH=$PWD:$PATH

eukcc -h

$ eukcc -h

usage: eukcc [-h] --db DB [--outdir OUTDIR] [--config CONFIG] [--ncores int]

             [--ncorespplacer int] [--hmm HMM] [--training] [--proteins]

             [--bed file.bed] [--force] [--keeptemp] [--fplace] [--noglob]

             [--quiet] [--debug] [--HPA] [--nPlacements n] [--minGenomes n]

             [--fullineage] [--minPlacementLikelyhood float] [--mindist n]

             [--touch] [--gmes] [--pygmes] [--diamond DIAMOND] [--plot] [-v]

             fasta

 

Evaluate completeness and contamination of a MAG. Args that start with '--'

(eg. --db) can also be set in a config file (specified via --config). Config

file syntax allows: key=value, flag=true, stuff=[a,b,c] (for details, see

syntax at https://goo.gl/R74nmi). If an arg is specified in more than one

place, then commandline values override config file values which override

defaults.

 

positional arguments:

  fasta                 Run script on this bin (fasta file)

 

optional arguments:

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

  --db DB               Path to EukCC DB

  --outdir OUTDIR, -o OUTDIR

                        Location for the output. Names will be prefixed using

                        the bin filenames

  --config CONFIG, -c CONFIG

                        Config file to define parameters, YAML

  --ncores int, -n int  set number of cores for GeneMark-ES, pplacer and Hmmer

  --ncorespplacer int   Pplacer requires a lot of memory. If you want you can

                        set less cores for pplacer, which improves memory

                        consumption significantly

  --hmm HMM             run hmmer on all these HMMs instead

  --training            Run EukCC in training mode (needed to create a new

                        release of the DB)

  --proteins            Input fasta is proteins

  --bed file.bed, -b file.bed

                        You can pass a bedfile of the protein location to omit

                        fragmented proteins being detected twice

  --force, -f           Force rerun of computation even if output is newer

                        than input. Don't resume previous run.

  --keeptemp            Keep all temporary files, by default EukCC will remove

                        some temp files

  --fplace, -p          Force rerun of placement and subsequent steps

  --noglob, -g          Do not expand paths using glob

  --quiet, -q           Silcence most output

  --debug, -d           Debug and thus ignore safety

  --HPA                 Set placement method to HPA

  --nPlacements n       Set number of proteins to support location in tree

                        (default: 2)

  --minGenomes n        Minimal number of genomes to support a set (default:

                        3)

  --fullineage          Output full lineage for MAGs

  --minPlacementLikelyhood float

                        minimal pplacer likelyhood (default: 0.4)

  --mindist n           Distance to collapse hits (default: 2000)

  --touch               Do not run, but touch all output files

  --gmes                only run GeneMark-ES

  --pygmes              Use pygmes, will improve eukccs capability of running

                        on highly fragmented bins but will take longer

  --diamond DIAMOND     required to use pygmes option

  --plot                produce plots

  -v, --version         show program's version number and exit

 

データベースの準備

mkdir Eukcc
cd Eukcc/
wget http://ftp.ebi.ac.uk/pub/databases/metagenomics/eukcc/eukcc_db_v1.1.tar.gz
tar -xzvf eukcc_db_v1.1.tar.gz
mv eukcc_db_20191023_1 eukccdb
#remove tar ball
rm eukcc_db_v1.1.tar.gz

eukcc実行時にはこのeukccdbディレクトリを指定する。 

 

テストラン

22 Contig配列からなる小さなメタゲノムアセンブリを使用する。

wget -O testgenome.fa.gz ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/002/251/995/GCF_002251995.1_ASM225199v2/GCF_002251995.1_ASM225199v2_genomic.fna.gz
gunzip testgenome.fa.gz
eukcc --db eukccdb \
--ncores 4 \
--ncorespplacer 1 \
--outdir eukcc_testgenome \
testgenome.fa
  • --db   Path to EukCC DB
  • --outdir   Location for the output. Names will be prefixed using the bin filenames
  • --ncores    set number of cores for GeneMark-ES, pplacer and Hmmer
  • --ncorespplacer    Pplacer requires a lot of memory. If you want you can set less cores for pplacer, which improves memory consumption significantly

テストデータのランでエラーになる。修正できたら追記します。

 

実際にビニングから始める場合のチュートリアルも用意されています。上のDocumentationリンクを確認して下さい。

引用

Estimating the quality of eukaryotic genomes recovered from metagenomic analysis with EukCC

Paul Saary, Alex L. Mitchell & Robert D. Finn
Genome Biology volume 21, Article number: 244 (2020)


関連


 

*1

cpanmは標準では~/perl5/lib/perl5/に入れるが、@INCに~/perl5/lib/perl5/が含まれていないならライブラリを認識できない。対策は単純に

export PERL5LIB=~/perl5/lib/perl5/

すればよい。

 

 

バクテリアゲノムの比較ゲノム解析を行うwebサービス EDGAR

2022/03/28 URL 更新, POCP matrix追記

 

Bergey's Manual of Systematics of Archaea and Bacteriaより

 次世代シークエンシングアプローチの展開により、完全にシークエンシングされたゲノムの数は急速に増加している。その結果、単一ゲノムだけでなく、関連する大規模なゲノム群を比較的に解析することが可能になった。また、タイプストレインのゲノムの全ゲノムシークエンシングは、より高分解能の系統分類・分類学的分類を得るための大きな可能性を秘めている。EDGARプラットフォームは、過去9年間で比較ゲノムの分野で最も確立されたソフトウェアツールの一つとなった。この間、ソフトウェアは継続的に改良され、多くの新しい解析機能が追加されてきた。近年では、コアゲノムに基づく系統図・分類学的解析にEDGARを利用することが、このソフトウェアの主な応用分野となっている。原核生物種のすべてのタイプストレインのゲノム配列を生成することに焦点を当てることで、基本的な16S rRNA遺伝子配列の系統図を大幅に拡張し、より高分解能のコアゲノムベースの分類学に拡張することができ、ラボでの集中的なDNA-DNAハイブリダイゼーション(DDH)を、DDHと同様に種の境界を反映するゲノム配列ベースの指標に置き換えることができる。

 EDGARのウェブベースのユーザーインターフェースは、新種の提案に必要な系統学的な種間・種内分類学的解析に必要なすべてのツールを提供している。EDGARは、隣接結合法と最尤法を用いてコアゲノムベースの系統樹を計算し、アミノ酸同一性(AI)と平均ヌクレオチド同一性(ANI)マトリックスを計算する。さらに、ベン図、シンテニープロット、直交遺伝子のゲノム近傍の比較表示などの便利な可視化機能を提供する。また、統計解析、レプリコングループ化オプション、メタ遺伝子セットの第2レベル解析など、様々な新機能が追加された。進化の関係を迅速に調査することができ、近親ゲノムの差異遺伝子の内容についての新たな生物学的洞察を得るプロセスを簡略化することができる。

 また、EDGARでは、比較ゲノムや系統図の結果を事前に計算したプロジェクトを公開データベースとして提供している。このプラットフォームでは、8,079の完全ゲノムからなる322の属ベースの公開データベースを提供している。ここでは、これらの属ベースのプロジェクトの他に、ファミリーレベルでクラスタリングされた226の新しい公開プロジェクトを紹介する。これらの新規公開プロジェクトは、さらに4,400のゲノムを含んでいる。EDGARは学術利用は無料で、ドイツのバイオインフォマティクス・インフラストラクチャ・ネットワーク(de.NBI)のサービスとして提供されている。EDGARは公開ウェブサーバー http://edgar.computational.bio から利用できる。

 

EAGAR2の論文より

 微生物のゲノム配列の利用可能性が急速に高まっていることから、密接に関連したゲノムの比較に基づく機能解析を支援するバイオインフォマティクスソフトウェアツールの需要が高まっている。遺伝子レベルでの比較アプローチを利用することで、研究対象となる生物の共通の特徴を表すコア遺伝子の洞察を得ることができる。また、逆にシングルトン遺伝子を同定することで、個々のゲノムの特異的な性質を明らかにすることができる。EDGARプラットフォームは、最初の出版以来、比較ゲノムの分野で最も確立されたソフトウェアツールの一つとなっている。ここ数年、EDGARは継続的に改良され、多くの新しい解析機能が追加されてきた。新バージョンのEDGAR 2.0では、遺伝子オルソロジー推定アプローチが新たに設計され、完全に再実装された。他の新機能の中でも、EDGAR 2.0では、AAI(Average Amino Acid Identity)やANI(Average Nucleotide Identity)行列のような系統解析機能の拡張、ゲノムセットサイズの統計、インタラクティブシンテニープロットやベン図のような最新の可視化機能を提供している。これにより、微生物ゲノム間の進化関係を迅速かつユーザーフレンドリーに調査することが可能となり、微生物ゲノムに含まれる遺伝子の違いについての新たな生物学的洞察を得るプロセスを簡素化することができる。すべての機能は、ウェブベースのプラットフォームに依存しないユーザーインターフェースを介して研究者に提供されており、事前に計算されたデータセットを簡単に閲覧することができる。ウェブサーバーは http://edgar.computational.bio からアクセスできる。

 

help

https://www.uni-giessen.de/fbz/fb08/Inst/bioinformatik/software/EDGAR/Features

 

webサービス

https://edgar3.computational.bio.uni-giessen.de/cgi-bin/edgar_login.cgiにアクセスする。

f:id:kazumaxneo:20200911012700p:plain

 

ここではサルモネラを選択。

f:id:kazumaxneo:20200911014238p:plain

 

まずはコアゲノムを探索してみる。パンゲノムも同じ感じに使えるが、コアゲノムとは異なり、パンゲノム解析では、特定のゲノムにしか見つからない遺伝子もプリントされる。

左のメニューからコアゲノムを選択。

f:id:kazumaxneo:20200911014317p:plain

SELECT ALLを選択するかクリックしてゲノムを選ぶALLの表記はchrとplasmid全てを含む。参照のゲノムは一度選択後に右クリックすることで選べる。

f:id:kazumaxneo:20200911194949p:plain

CP015724を選んだ。リファレンスのゲノムはより黒い色に変わる。

 

SHOW CORE GENOMEボタンを押すと解析がスタートする。 ゲノムが多すぎると計算が終わらないことがあるので注意。

 

出力は遺伝子のリストとなる。1行に1セットのオーソロガス遺伝子が表示される。参照のゲノムは結果は1列目に表示される。パンゲノム解析では、特定の遺伝子のオルソログが他のゲノムに見当たらない場合はセルに"-"が表示される。

f:id:kazumaxneo:20200911221302p:plain

 

各行の最後の2列(下図)には、その行の遺伝子のマルチプルアラインメントと、翻訳開始部位の前の400bpの上流領域のマルチプルアラインメントを作成するためのリンクになっている。DNAまたはアミノ酸配列どちらでも作成できる。

f:id:kazumaxneo:20200911223357p:plain

 

遺伝子はFASTAファイル(DNAまたはタンパク質配列)や座 標タグと機能記述をTABで区切ったフラットファイルとして保存できる。

f:id:kazumaxneo:20200911223712p:plain

 

 

他の解析も、同じようにゲノムを選んでスタートさせる。

f:id:kazumaxneo:20200911224320p:plain

 

結果だけ順番に見ていく。

Upset plot

f:id:kazumaxneo:20200911200959p:plain

 

Venn diagram

f:id:kazumaxneo:20200911201154p:plain

GET LISTから結果のリストをダウンロードできる。

ゲノムの数が多すぎると視覚化結果はよくわからなくなるので注意。その場合は複雑な部分集合を調べることができるCalculate genesetsを使う。

 

Circular Plot

上では3ゲノム選択して描画した。

f:id:kazumaxneo:20200911201643p:plain

f:id:kazumaxneo:20200911201647p:plain

Core development plot

より多い数のゲノムに対して「本当の」コアゲノムのサイズやシングルトン数を推定する。
ある特定のゲノムセットについて計算されたコアゲノムは、常にそのゲノムセットの状況のスナップショットに過ぎない。種の "本当の "コアゲノムを知るためには、無限のサイズのゲノムセットについてコア遺伝子の数を外挿する。これは、Tettelinら(2005)による近似的なアプローチによって行われる。

f:id:kazumaxneo:20200911201804p:plain

Singleton development plot

f:id:kazumaxneo:20200911202319p:plain

 

Pan development plot

f:id:kazumaxneo:20200911202423p:plain

 

Pan vs. Core development plot

f:id:kazumaxneo:20200911202218p:plain

 

Synteny Plots

下図のように、3つ位以上のゲノムでも実行可能。

f:id:kazumaxneo:20200911230029p:plain

 

Synteny Matrix

f:id:kazumaxneo:20200911230145p:plain

 

Create AAI matrix

f:id:kazumaxneo:20200911231622p:plain

f:id:kazumaxneo:20200911231635p:plain

 

Create ANI matrix

f:id:kazumaxneo:20200911231836p:plain

 

POCP (percentage of conserved genes) matrix

 

Phylogenetic trees

EDGARの主な目的は複数のゲノム間のオルソログの推定であるため、そのオルソログ情報から系統知識を推論することは容易にできる。このパイプラインでは、選択されたゲノムのコアゲノムを計算する。全ゲノム中に見出されたオロソロガスな遺伝子のセットは、マルチプルアラインメントツールMUSCLE(Edgar, 2004)を用いて個別にアラインメントされている。アライメントは、数十万残基になる巨大なマルチプルアラインメントに連結される。このアラインメントから距離行列を計算し、最終的にこの距離行列に基づいてNeighbor-Joining法を用いて系統樹が構築されている。後者の二つの方法は、Felsenstein (1995) による PHYLIP の実装で使用されている。Neighbor-Joining法が選択されているのは、非常に計算効率の良いヒューリスティックなアプローチであるため、コアゲノムに基づく系統樹構築の結果として得られる大規模なデータセットに適しているからである。

f:id:kazumaxneo:20200911232028p:plain

 

補足

マニュアルより

EDGARは、より高いレベルの比較のために "メタコンティグ "やレプリコングループを定義する。基本的にEDGARは単一のコンティグや完全なゲノムの比較を提供しているが、より複雑な比較が必要な場合(例えばある生物のすべてのプラスミドと別の生物のすべてのプラスミドの遺伝子内容を比較したいなど)、生物とコンティグの間の抽象化レベルが必要となる。このような場合、ユーザーが選択したコンティグのグループを作成して使用することができる。

引用

EDGAR: A Versatile Tool for Phylogenomics

Jochen Blom Stefanie P. Glaeser Tobias Juhre Julian Kreis Patrick H.G. Hanel Jascha G. Schrader Peter Kämpfer Alexander Goesmann

Bergey's Manual of Systematics of Archaea and Bacteria, First published: 06 June 2019

https://doi.org/10.1002/9781118960608.bm00038


EDGAR 2.0: an enhanced software platform for comparative gene content analyses
Jochen Blom, Julian Kreis, Sebastian Spänig, Tobias Juhre, Claire Bertelli, Corinna Ernst, Alexander Goesmann

Nucleic Acids Res. 2016 Jul 8; 44(Web Server issue): W22–W28


EDGAR: a software framework for the comparative analysis of prokaryotic genomes

Jochen Blom, Stefan P Albaum, Daniel Doppmeier, Alfred Pühler, Frank-Jörg Vorhölter, Martha Zakrzewski, Alexander Goesmann

BMC Bioinformatics. 2009 May 20;10:154

 

EDGAR3.0


 

 

Average amino acid identity(AAI)を計算する CompareM

2020 10/5 追記

2020 10/9 修正

 

CompareMは大規模な比較ゲノム解析をサポートするソフトウェアツールキットである。ゲノムのセット(アミノ酸の同一性など)と個々のゲノム(コドン使用率など)の統計を提供する。数千ゲノムへのスケーラビリティを可能にするために、並列化された実装が提供されている。一般的なワークフローは、ユーザーが簡単に導入できるように単一のメソッドとして提供され、経験豊富なユーザーが特定の機能を利用できるようにより詳細なインターフェースが提供されている。CompareMはオープンソースであり、GNU General Public License (Version 3)の下でリリースされている。

 CompareMには2つの共通のワークフローがある。1つ目はゲノム間のペアワイズAAIを計算するためのもので(--aai_wf)、2つ目はリファレンスゲノムのセットに対してAAIの値を計算してクエリゲノムを分類するものである(--classify_wf)。これらのワークフローはCompareMが提供するより詳細なコマンドのいくつかを呼び出している。具体的には、aai_wf は call_genes, similarity, and aai コマンドを実行し、classify_wf は call_genes, similarity, and classify コマンドを実行する。CompareM は Prodigal を使用してゲノムのセットそれぞれの遺伝子をコールするが、ユーザーが提供したタンパク質配列セットを使用することもできる。

 

ユーザーガイド

https://github.com/dparks1134/CompareM/blob/master/users_guide.pdf

 

インストール

yamlファイルが3.6以上になっているので、conda mambaでpython3.6の仮想環境を作ってテストした( macos10.14で試したが、エラーメッセージが出ないまま正しい数値が出力されない危険なエラーが起きた。ubuntu18.04で動作確認した)。

依存

  • CompareM makes use of the numpy, scipy, matplotlib, and biolib python packages
  • prodigal >= 2.6.2
  • diamond >= 0.9.0

Github

#bioconda (link)
mamba create -n Comparem -y python=3.6
conda activate Comparem
mamba install -c bioconda comparem -y

comparem

$ comparem 

 

                ...::: CompareM v0.1.1 :::...

 

    Common workflows:

     aai_wf      -> Calculate AAI between all pairs of genomes

                    (runs call_genes => similarity => aai)

     classify_wf -> Identify similar genomes based on AAI values

                    (runs call_genes => similarity => classify)

                    

    Gene prediction:

     call_genes -> Identify genes within genomes

                     

    Gene homology and genome similarity:

     similarity -> Perform reciprocal sequence similarity search between proteins

     aai        -> Calculate AAI between all pairs of genomes

     classify   -> Identify similar genomes based on AAI value

 

    Usage profiles:

     aa_usage    -> Calculate amino acid usage within each genome

     codon_usage -> Calculate codon usage within each genome

     kmer_usage  -> Calculate kmer usage within each genome

     stop_usage  -> Calculate stop codon usage within each genome

 

    Lateral gene transfer:

     lgt_di    -> Calculate dinuceotide (3rd,1st) usage of genes to identify putative LGT events

     lgt_codon -> Calculate codon usage of genes to identify putative LGT events

 

    Visualization and exploration:

     diss      -> Calculate the dissimilarity between usage profiles

     hclust    -> Perform hierarchical clustering

 

  Use: comparem <command> -h for command specific help.

 

  Feature requests or bug reports can be sent to Donovan Parks (donovan.parks@gmail.com)

    or posted on GitHub (https://github.com/dparks1134/comparem).

    

comparem aai_wf -h

$ comparem aai_wf -h

usage: comparem aai_wf [-h] [-e EVALUE] [-p PER_IDENTITY] [-a PER_ALN_LEN]

                       [-x FILE_EXT] [--proteins] [--force_table FORCE_TABLE]

                       [--blastp] [--sensitive] [--keep_headers] [--keep_rbhs]

                       [--tmp_dir TMP_DIR] [-c CPUS] [--silent]

                       input_files output_dir

 

Calculate AAI between all pairs of genomes

 

positional arguments:

  input_files           genome files

  output_dir            output directory

 

optional arguments:

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

  -e, --evalue EVALUE   e-value cutoff for identifying initial blast hits

                        (default: 0.001)

  -p, --per_identity PER_IDENTITY

                        percent identity for defining homology (default: 30.0)

  -a, --per_aln_len PER_ALN_LEN

                        percent alignment length of query sequence for

                        defining homology (default: 70.0)

  -x, --file_ext FILE_EXT

                        extension of files to process (default: fna)

  --proteins            indicates the input files contain protein sequences

  --force_table FORCE_TABLE

                        force use of specific translation table

  --blastp              use blastp instead of diamond

  --sensitive           use sensitive mode of DIAMOND

  --keep_headers        indicates FASTA headers already have the format

                        <genome_id>~<gene_id>

  --keep_rbhs           create file with reciprocal best hits

  --tmp_dir TMP_DIR     specify alternative directory for temporary files

                        (default:

                        /var/folders/v0/rdv6fcx11fz1dy2jjvthmx1r0000gn/T)

  -c, --cpus CPUS       number of CPUs to use (default: 1)

  --silent              suppress output

comparem classify_wf -h

$ comparem classify_wf -h

husage: comparem classify_wf [-h] [-k NUM_TOP_TARGETS] [-t TAXONOMY_FILE]

                            [-e EVALUE] [-p PER_IDENTITY] [-a PER_ALN_LEN]

                            [-x FILE_EXT] [--proteins]

                            [--force_table FORCE_TABLE] [--blastp]

                            [--sensitive] [--keep_headers] [--keep_rbhs]

                            [--tmp_dir TMP_DIR] [-c CPUS] [--silent]

                            query_files target_files output_dir

 

Identify similar genomes based on AAI value.

 

positional arguments:

  query_files           query genome files

  target_files          target genome files

  output_dir            output directory

 

optional arguments:

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

  -k, --num_top_targets NUM_TOP_TARGETS

                        number of top scoring target genomes to report per

                        query genome (default: 1)

  -t, --taxonomy_file TAXONOMY_FILE

                        file indicating taxonomic identification of all target

                        genomes

  -e, --evalue EVALUE   e-value cutoff for identifying initial blast hits

                        (default: 0.001)

  -p, --per_identity PER_IDENTITY

                        percent identity for defining homology (default: 30.0)

  -a, --per_aln_len PER_ALN_LEN

                        percent alignment length of query sequence for

                        defining homology (default: 70.0)

  -x, --file_ext FILE_EXT

                        extension of files to process (default: fna)

  --proteins            indicates the input files contain protein sequences

  --force_table FORCE_TABLE

                        force use of specific translation table

  --blastp              use blastp instead of diamond

  --sensitive           use sensitive mode of DIAMOND

  --keep_headers        indicates FASTA headers already have the format

                        <genome_id>~<gene_id>

  --keep_rbhs           create file with reciprocal best hits

  --tmp_dir TMP_DIR     specify alternative directory for temporary files

                        (default:

                        /var/folders/v0/rdv6fcx11fz1dy2jjvthmx1r0000gn/T)

  -c, --cpus CPUS       number of CPUs to use (default: 1)

  --silent              suppress output

comparem call_genes -h

$ comparem call_genes -h

usage: comparem call_genes [-h] [-x FILE_EXT] [--force_table FORCE_TABLE]

                           [-c CPUS] [--silent]

                           input_genomes output_dir

 

Identify genes within genomes.

 

positional arguments:

  input_genomes         genome files to process

  output_dir            output directory

 

optional arguments:

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

  -x, --file_ext FILE_EXT

                        extension of files to process (default: fna)

  --force_table FORCE_TABLE

                        force use of specific translation table

  -c, --cpus CPUS       number of CPUs to use (default: 1)

  --silent              suppress output

comparem similarity -h

$ comparem similarity -h

usage: comparem similarity [-h] [-e EVALUE] [-p PER_IDENTITY] [-a PER_ALN_LEN]

                           [-x FILE_EXT] [--blastp] [--sensitive]

                           [--keep_headers] [--tmp_dir TMP_DIR] [-c CPUS]

                           [--silent]

                           query_proteins target_proteins output_dir

 

Perform sequence similarity search between genes.

 

positional arguments:

  query_proteins        query protein files to process

  target_proteins       target protein files to process

  output_dir            output directory

 

optional arguments:

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

  -e, --evalue EVALUE   maximum e-value for reporting an alignments (default:

                        0.001)

  -p, --per_identity PER_IDENTITY

                        minimum percent identity for reporting an alignment

                        (default: 30.0)

  -a, --per_aln_len PER_ALN_LEN

                        minimum percent coverage of query sequence for

                        reporting an alignment (default: 70.0)

  -x, --file_ext FILE_EXT

                        extension of files to process (default: faa)

  --blastp              use Blastp-fast instead of DIAMOND

  --sensitive           use sensitive mode of DIAMOND

  --keep_headers        indicates FASTA headers already have the format

                        <genome_id>~<gene_id>

  --tmp_dir TMP_DIR     specify alternative directory for temporary files

                        (default:

                        /var/folders/v0/rdv6fcx11fz1dy2jjvthmx1r0000gn/T)

  -c, --cpus CPUS       number of CPUs to use (default: 1)

  --silent              suppress output

comparem aai -h

$ comparem aai -h

usage: comparem aai [-h] [-e EVALUE] [-p PER_IDENTITY] [-a PER_ALN_LEN]

                    [--keep_rbhs] [-c CPUS] [--silent]

                    query_gene_file sorted_hit_table output_dir

 

Calculate the AAI between all genome pairs.

 

positional arguments:

  query_gene_file       file with all query genes

  sorted_hit_table      sorted file indicating genes passing sequence

                        similarity criteria

  output_dir            output directory

 

optional arguments:

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

  -e, --evalue EVALUE   maximum e-value for reporting an alignments (default:

                        0.001)

  -p, --per_identity PER_IDENTITY

                        minimum percent identity for reporting an alignment

                        (default: 30.0)

  -a, --per_aln_len PER_ALN_LEN

                        minimum percent coverage of query sequence for

                        reporting an alignment (default: 70.0)

  --keep_rbhs           create file with reciprocal best hits

  -c, --cpus CPUS       number of CPUs to use (default: 1)

  --silent              suppress output

comparem classify -h

$ comparem classify -h

usage: comparem classify [-h] [-k NUM_TOP_TARGETS] [-t TAXONOMY_FILE]

                         [-e EVALUE] [-p PER_IDENTITY] [-a PER_ALN_LEN]

                         [-x FILE_EXT] [--keep_rbhs] [-c CPUS] [--silent]

                         query_gene_file target_gene_file sorted_hit_table

                         output_dir

 

Identify similar genomes based on AAI value.

 

positional arguments:

  query_gene_file       file with all query genes

  target_gene_file      file with all target genes

  sorted_hit_table      sorted file indicating genes passing sequence

                        similarity criteria

  output_dir            output directory

 

optional arguments:

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

  -k, --num_top_targets NUM_TOP_TARGETS

                        number of top scoring target genomes to report per

                        query genome (default: 1)

  -t, --taxonomy_file TAXONOMY_FILE

                        file indicating taxonomic identification of all target

                        genomes

  -e, --evalue EVALUE   e-value cutoff for identifying initial blast hits

                        (default: 0.001)

  -p, --per_identity PER_IDENTITY

                        percent identity for defining homology (default: 30.0)

  -a, --per_aln_len PER_ALN_LEN

                        percent alignment length of query sequence for

                        defining homology (default: 70.0)

  -x, --file_ext FILE_EXT

                        extension of files to process (default: fna)

  --keep_rbhs           create file with reciprocal best hits

  -c, --cpus CPUS       number of CPUs to use (default: 1)

  --silent              suppress output

comparem aa_usage –h

$ comparem aa_usage –h

usage: comparem aa_usage [-h] [--counts] [-x FILE_EXT] [-c CPUS] [--silent]

                         protein_gene_files output_file

comparem aa_usage: error: the following arguments are required: output_file

comparem codon_usage -h

$ comparem codon_usage -h

usage: comparem codon_usage [-h] [--counts] [-x FILE_EXT] [--keep_ambiguous]

                            [-c CPUS] [--silent]

                            nucleotide_gene_files output_file

 

Calculate codon usage within each genome.

 

positional arguments:

  nucleotide_gene_files

                        input files with genes in nucleotide space

  output_file           output file indicating codon usage of each genome

 

optional arguments:

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

  --counts              output raw counts instead of frequencies

  -x, --file_ext FILE_EXT

                        extension of files to process (default: fna)

  --keep_ambiguous      keep codons with ambiguous bases

  -c, --cpus CPUS       number of CPUs to use (default: 1)

  --silent              suppress output

comparem kmer_usage -h

$ comparem kmer_usage -h

usage: comparem kmer_usage [-h] [--counts] [-k K] [-x FILE_EXT] [-c CPUS]

                           [--silent]

                           genome_files output_file

 

Calculate kmer usage within each genome.

 

positional arguments:

  genome_files          input files with genomes in nucleotide space

  output_file           output file indicating kmer usage of each genome

 

optional arguments:

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

  --counts              output raw counts instead of frequencies

  -k K                  length of kmers, e.g., 4 -> tetranucleotides (default:

                        4)

  -x, --file_ext FILE_EXT

                        extension of files to process (default: fna)

  -c, --cpus CPUS       number of CPUs to use (default: 1)

  --silent              suppress output

comparem stop_usage -h

$ comparem stop_usage -h

usage: comparem stop_usage [-h] [--counts] [--mean_gene_length] [-x FILE_EXT]

                           [-c CPUS] [--silent]

                           nucleotide_gene_files output_file

 

Calculate stop codon usage within each genome.

 

positional arguments:

  nucleotide_gene_files

                        input files with genes in nucleotide space

  output_file           output file indicating stop codon usage of each genome

 

optional arguments:

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

  --counts              output raw counts instead of frequencies

  --mean_gene_length    report mean gene length for genes with each stop codon

  -x, --file_ext FILE_EXT

                        extension of files to process (default: fna)

  -c, --cpus CPUS       number of CPUs to use (default: 1)

  --silent              suppress output

comparem lgt_di -h

$ comparem lgt_di -h

usage: comparem lgt_di [-h] [-x FILE_EXT] [--crit_value CRIT_VALUE] [-c CPUS]

                       [--silent]

                       nucleotide_gene_files output_dir

 

Calculate dinuceotide (3rd,1st) usage of genes to identify putative LGT

events.

 

positional arguments:

  nucleotide_gene_files

                        input files with genes in nucleotide space

  output_dir            output directory to write dinucleotide usage for each

                        gene in each genome

 

optional arguments:

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

  -x, --file_ext FILE_EXT

                        extension of files to process (default: fna)

  --crit_value CRIT_VALUE

                        critical value for defining deviant genes (default:

                        0.001)

  -c, --cpus CPUS       number of CPUs to use (default: 1)

  --silent              suppress output

comparem lgt_codon -h

$ comparem lgt_codon -h

usage: comparem lgt_codon [-h] [-x FILE_EXT] [-c CPUS] [--silent]

                          nucleotide_gene_files output_dir

 

Calculate codon usage of genes to identify putative LGT events.

 

positional arguments:

  nucleotide_gene_files

                        input files with genes in nucleotide space

  output_dir            output directory to write dinucleotide usage for each

                        gene in each genome

 

optional arguments:

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

  -x, --file_ext FILE_EXT

                        extension of files to process (default: fna)

  -c, --cpus CPUS       number of CPUs to use (default: 1)

  --silent              suppress output

comparem diss -h

$ comparem diss -h

usage: comparem diss [-h]

                     [--metric {euclidean,minkowski,cityblock,seuclidean,sqeuclidean,cosine,correlation,hamming,jaccard,chebyshev,canberra,braycurtis,mahalanobis,yule,matching,dice,kulsinski,rogerstanimoto,russellrao,sokalmichener,sokalsneath,wminkowski}]

                     [--full_matrix] [--silent]

                     profile_file output_file

 

Calculate the dissimilarity between usage profiles.

 

positional arguments:

  profile_file          file with usage profile for each genome

  output_file           output file with pairwise dissimilarity between

                        genomes

 

optional arguments:

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

  --metric {euclidean,minkowski,cityblock,seuclidean,sqeuclidean,cosine,correlation,hamming,jaccard,chebyshev,canberra,braycurtis,mahalanobis,yule,matching,dice,kulsinski,rogerstanimoto,russellrao,sokalmichener,sokalsneath,wminkowski}

                        distance metric to use (default: euclidean)

  --full_matrix         output full dissimilarity matrix

  --silent              suppress output

comparem hclust -h

$ comparem hclust -h

usage: comparem hclust [-h]

                       [--method {single,complete,average,weighted,centroid,median,ward}]

                       [--similarity] [--max_sim_value MAX_SIM_VALUE]

                       [--name_col1 NAME_COL1] [--name_col2 NAME_COL2]

                       [--value_col VALUE_COL] [--silent]

                       pairwise_value_file output_tree

 

Perform hierarchical clustering.

 

positional arguments:

  pairwise_value_file   file with pairwise similarity or dissimilarity values

                        between genomes

  output_tree           name for output hierarchical cluster tree

 

optional arguments:

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

  --method {single,complete,average,weighted,centroid,median,ward}

                        clustering method to use. (default: average)

  --similarity          indicates file contain similarity values

  --max_sim_value MAX_SIM_VALUE

                        maximum similarity value (default: 100)

  --name_col1 NAME_COL1

                        index of first column with genome names (default: 0)

  --name_col2 NAME_COL2

                        index of second column with genome names (default: 1)

  --value_col VALUE_COL

                        index of column with similarity or dissimilarity

                        values (default: 2)

  --silent              suppress output

 

 

実行方法

ランするにはfasta形式のゲノムファイルが必要。

f:id:kazumaxneo:20200910194851p:plain

 

1、comparem aai_wf  ペアワイズAAIの計算

比較するゲノムのディレクトリ、出力ディレクトリを指定する。 

comparem aai_wf --cpus 12 -x fna input_dir/ outdir

#または翻訳済みのproteinディレクトリを指定
comparem aai_wf --cpus 12 --proteins -x faa input_dir/ outdir
  • -e    e-value cutoff for identifying initial blast hits (default: 0.001)
  • -p    percent identity for defining homology (default: 30.0)
  • -a    percent alignment length of query sequence for defining homology (default: 70.0)
  • --sensitive   use sensitive mode of DIAMOND
  • -x    extension of files to process (default: fna)
  • -blastp   use blastp instead of diamond
  • --proteins    indicates the input files contain protein sequences

出力

f:id:kazumaxneo:20200910195147p:plain

aai/aai_summary.tsv

f:id:kazumaxneo:20200910195841p:plain

出力はall versus allのマトリックスファイルではなく、ペアワイズの比較結果のテーブルになる。使用したデータでは、上の表のようにgenome3と他のゲノムとの比較(2−5行目)、次にゲノム5と他のゲノムとの比較(6−8行目)、次にゲノム1(9行目)、最後にゲノム4(10−11行目)という順番でプリントされている。ゲノム5から見たゲノム3との比較結果は2行目にあるが、ゲノム3から見たゲノム5の比較結果はない。レシプロカルヒットの平均からAAIを出しているなら不要になる。

 

*2

 

2、comparem classify_wf AAIの結果から近いゲノムを推定。

  クエリのゲノムと比較するゲノムのディレクトリ、出力ディレクトリを指定する。 

comparem classify_wf -x fna --cpus 12 query_genome.fna input_dir outdir

テスト時はエラーを起こした。

 

補足

--keep_headers をつけるとAAIだけでなく全比較結果が残る。aai_summary.tsvは巨大なテキストファイルになるので、ラン後にユーザーがパースする必要がある。

 

現在のところ、CompareMの実装にはエラーが出る可能性のある部分がある。詳細はGithubで確認してください。

引用

GitHub - dparks1134/CompareM: A toolbox for comparative genomics.

 

参考

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

 

 

*1 エラーが起きるときはファイル名が問題の可能性がある。シンプルなファイル名にリネームする。

 

*2

Teaching-Dundee-BS32010/06-RBBH.ipynb at master · widdowquinn/Teaching-Dundee-BS32010 · GitHub

 

 

入力プロテオームから類似したタンパク質のデータベースを自動検索し、プロテオームから近い種を調べる AAI-profiler

 

 全ゲノムショットガンシーケンスは、分類学的分類の再評価を推進し、シングルセルゲノミクスの出現は生物多様性に関する知識を大きく広げている(1)。これらすべての応用分野において、分類学的分類に関するオリジナルの文献を検索するよりも、配列データを直接比較する方が、分類学的・系統的関係の概要を迅速に把握することができる。残念ながら、配列データベースのメタデータは、古い同義語を使用していたり、完全に誤って分類されていたりすることがある。多くの推論手法は、配列の種の割り当てが正しいと仮定して、配列のツリーと種のツリー(分類学)との整合性をテストするため、正しいメタデータは重要である。このようなアプリケーションには、種分化、遺伝子複製イベント(2)および水平方向の遺伝子移動イベント(3)を識別するためのツリーの再構成、メタゲノミクスにおける分類学的プロファイリングのためのlowest common ancestor(LCA)アプローチ(4)、およびLCA分類群の配列クラスタへの割り当て(5)などがある。

 近年、Pairwise overall genomic relatedness indices(OGRI(6))は、種の発見と同定において人気を博している。OGRIは、16S rRNA遺伝子配列の分解能が限られていることや、Multilocus Sequence Analysis (MLSA)におけるデータの欠落など、遺伝子ベースの計算検定の限界をいくつか克服したものである(7)。ゲノム全体の関連性の指標としては、Karlinゲノムシグネチャー、平均ヌクレオチド同一性(ANI)、平均アミノ酸同一性(AAI)、スーパーツリー、in silico Genome-to-Genome Distance Hybridization(GGDH)などがある(8)。特に、従来の多相性分類からゲノム微生物分類への転換が期待されており、時間のかかる実験室での検査による表現型の特徴付けを避けることができる(8-11)。

 ここでは、クエリプロテオームとUniprotデータベース内の全標的種との間のAIを計算する、ユーザーフレンドリーなウェブサーバーであるAI-profilerを紹介する。既存のツール(12-15)では、2つの種を比較したり、あらかじめ定義された少数の種のセットを比較したり、メタゲノミクスを目的としているのとは対照的に、AAI-profilerは1つの種のプロテオームを入力とし、類似したタンパク質を持つ種のタンパク質データベースを自動的に検索する。AAI-profilerは、AAI距離検索を含む一連のメタゲノムおよびゲノム比較分析を実装しているMiGAのウェブサイト(http://microbial-genomes.org/)に似ている。AAI-profilerはBLASTの代わりにSANSparallel(16)を使用しているため、応答時間が速くなっているが、最も重要な違いは、AAI-profilerで検索されたUniprotデータベースは、MiGAで検索された原核生物のリファレンスゲノムコレクション(NCBI RefSeqで1927件、NCBI Prokで11 566件)に比べて、より多くの種の表現(809 540件のユニークラベル)を持っていることである。

 AAI-profilerは、高速相同性検索ツールであるSANSparallel(16)を利用している。相同性検索では、配列同一性が約50%までの隣接する細菌属や哺乳類のファミリーを識別するのに十分なレベルまで検出する。検索対象の種は、そのプロテオーム(FASTA形式のタンパク質配列)で表される。ヌクレオチド配列ではなくアミノ酸配列を比較することで、AAI-profilerを真核生物に適用することができる。真核生物のゲノムは細菌のゲノムの数百倍から数千倍の長さがあるが、遺伝子密度が低いため、真核生物のプロテオームは細菌の10倍程度の大きさしかない。例えば、大腸菌ゲノムは約5Mb、フェレットゲノムは約2.4Gbで、それぞれ約5000と48000のタンパク質をコードしている(NCBIゲノム、https://www.ncbi.nlm.nih.gov/genome)。SANSparallelを用いて、AAI-profilerは細菌のプロテオームを数分、真核生物のプロテオームを1時間で処理することができる。

 AAI-profilerの結果は2次元散布図で示される。横軸はAAIで、いくつかのソースが種の境界を∼95% AAIとしている(11)。姉妹種と姉妹属は連続して低いAAI値を持っている(論文図1A)。縦軸はカバレッジ、すなわち一致したタンパク質ペアの割合である。ゲノムが完全に配列決定されている種は、サンプル数が少ない種よりも、異なるタンパク質ファミリーが異なる速度で進化するため、よりロバストな AAI 推定値を与える。散布図のデータポイントは、分類学的なグループ分けに応じて色分けされており、予想される単系統のパターンの例外を視覚的に見分けることが容易になる。一つは、分類群は単系統であると予想され、したがって、分類群内の距離は、異なる分類群からの種間の距離よりも小さいはずである。したがって、期待されるパターンは、クエリプロテオームの種と属が、最も高いAAI値で一様な色のクラスタで表示される(図1B。*1 )。しかし、多くの場合、この「自己」クラスターの中には、ラベル付けされていない種が散在している。クエリプロテオームのカバレッジが高い場合(「一致率」が高い場合)、これらの例外は、誤って分類されたサンプルや誤ってラベル付けされたサンプルに起因している。また、カバレッジが低い場合は、コンタミネーションや遺伝子の水平伝播の可能性があることを示唆している。

 

 クエリープロテオーム(FASTA形式のタンパク質配列)について、一方向性と双方向性のAAIプロファイルを計算する。SANSparallelを使用して、Uniprotから相同タンパク質を検索する。Species情報は、UniprotヘッダーのOSタグから取得している。分類学メタデータはDictServerから取得している。各クエリタンパク質について、最も高いビットスコアを持つデータベースの全種との一致を保持する。蛋白質のデータベース化を行い、その蛋白質とデータベース種の標的蛋白質との many-to-oneのマッピングを行う。多重度は、ターゲット種で一致するクエリタンパク質の数を、異なるターゲットタンパク質の数で割ったものと定義している。例えば、あるタンパク質ファミリーがクエリ種の中で多数のパラログにまで拡大している場合、多重度は1よりも大きくなる。この効果は、転移可能なエレメントによってコードされた(疑似)タンパク質が問い合わせプロテオームに含まれている場合に顕著になることがある。双方向性AAIプロファイルはone-to-oneマッピングに基づいており、どちらかの配列に高いスコアの一致がある場合には、クエリタンパク質とデータベースタンパク質の一致を除外している。双方向ヒットの多重度は定義上1である。種ごとの一致数は、幅1%の配列同一性ビンで集計される。配列同一性はSANSparallelによって返されたアラインメントの位置ごとに計算される。AAIは,クエリプロテオームとデータベースの種の間でマッチしたすべてのペアの配列同一性の平均で,各クエリタンパク質は1の重みを持つ。ノーマッチのクエリタンパク質(SANSparallelによって報告)は0の重みを持つ。

 

 

tutorial PDF

http://ekhidna2.biocenter.helsinki.fi/AAI/AAI.pdf

example result

http://ekhidna2.biocenter.helsinki.fi/AAI/

 

AAI-profilerのダウンロード(ダイレクトリンク)

http://ekhidna2.biocenter.helsinki.fi/AAI/#download

 

webサービス

http://ekhidna2.biocenter.helsinki.fi/AAI/ にアクセスする。

f:id:kazumaxneo:20200910124057p:plain

クエリのプロテオーム配列(multi fasta)を入力する。

f:id:kazumaxneo:20200910133105p:plain

 タイトルとメールアドレスを記載してsubmitをクリック。

f:id:kazumaxneo:20200910133216p:plain

クエリのサイズとサーバーの混雑度によってランタイムは変わる。バクテリアのプロテオーム配列(4000配列ほど)を使ってテストした時は数分で計算は終わった。

 

出力

リンクが表示される。

f:id:kazumaxneo:20200910134419p:plain

example出力

出力図の横軸はクエリとデータベース種間のAAIを示している。平均は、SANSparallelが一致を報告しているクエリ蛋白質について、データベースの種ごとの最良の一致を計算している。縦軸は,その種で一致するクエリタンパク質の割合を示している。プロットは、属(細菌)または目(真核生物)によって色分けされている。真核生物の種はダイヤモンド、バクテリアは丸、古細菌は十字、それ以外のもの(ウイルス、メタゲノム、未分類サンプル)は四角で表示される。

f:id:kazumaxneo:20200910134521p:plain

 

スコア高かった特定の種が自動で選ばれ、クエリとターゲットの1対1で全タンパク質のAAI比較結果が示される。下では6つの図があるが、もっと多い場合がある。分布が狭く急激なピークが観察されるのは、登録されているタンパク質が少ないためである。完全に配列決定されているプロテオーム同士を比較する場合、ブロードでなだらかな分布が描かれる。

一番上は自分自身との比較(上の図で一番右上の位置)なので、下のグラフになる(データベースにクエリ配列があれば)。

f:id:kazumaxneo:20200910140517p:plain

 

その下がスコアが高かったいくつかの種との比較になる。

f:id:kazumaxneo:20200910134556p:plain

f:id:kazumaxneo:20200910135104p:plain

データベースへのタンパク質登録数がいずれも少ないのか、分布形状はシャープ。他のexampleの方がわかりやすいかもしれない(リンク)。

 

Taxonomic profile excluding top species

f:id:kazumaxneo:20200910134650p:plain

 

Taxonomic profile including top species

f:id:kazumaxneo:20200910134713p:plain

 

Onesided AAI profiles

excelで開いた。

f:id:kazumaxneo:20200910134901p:plain

Bidirectional AAI profiles

excelで開いた。

f:id:kazumaxneo:20200910135008p:plain

 

 

引用

AAI-profiler: fast proteome-wide exploratory analysis reveals taxonomic identity, misclassification and contamination
Alan J Medlar, Petri Törönen, Liisa Holm
Nucleic Acids Research, Volume 46, Issue W1, 2 July 2018, Pages W479–W485

 

関連


 

*1 

ここでのクエリはFrankliniella occidentalisという農業上の害虫の共生菌であるBFo1である。もう1つは、地理的に隔離された F. occidentalislから見つかった共生菌BFo2でこちらの結果も比較されている。詳細は論文のRESULTSのセクション参照。

ところで、図1Bでの横軸のAAIでのベストマッチ(一番右下のプロット)は、AAIがほぼ1のErwinia aphidicolaだが、図のプロットでは縦軸がほぼゼロのところに位置している。これは、E. aphidicolaはuniprotにタンパク質が数十個しか登録されていないためと述べられている。

ラージデータセットのコアゲノムを高速に構築する CoreCruncher

 

 コアゲノムとは、原核生物のある集団や種のすべての、あるいはほぼすべての系統が共有する遺伝子の集合を意味する。コアゲノムを推定することは多くのゲノム解析に不可欠だが、ほとんどの手法はすべてのゲノムのペアを比較することに依存している。ここでは、数百から数千のゲノムのコアゲノムをロバストかつ迅速に構築するプログラムであるCoreCruncherを紹介する。CoreCruncherはすべてのペアワイズゲノム比較は計算せず、同一性スコアの分布に基づいたヒューリスティックな手法を用いて、配列をorthologsまたは paralogs/xenologsに分類する。現在の方法よりも高速だが、他のツールよりも保守的で、paralogsやxenologsの存在に敏感ではないことが分かった。CoreCruncher は https://github.com/lbobay/CoreCruncher で自由に利用できる。CoreCruncherはPython 3.7で書かれているが、Python 2.7でもそのまま実行できる。一部のオプションでは、プログラム muscle または mafft が必要である。

 

インストール

condaでpython3.7の仮想環境を作ってテストした(macos10.14)。

依存

  • Usearch (32 bits) or BLAST (tested successfully with BLAST 2.7.1+)
  • python library numpy
  • OPTIONAL: can align core genomes and provide sequence concatenate if alignment program is available Alignment program muscle or mafft must be executable and in /usr/local/bin

usearchを使うなら実行ファイル'usearch61' として、パスの通った /usr/local/bin/などに置く。

mamba create -n CoreCruncher -y python=3.7
conda activate CoreCruncher
mamba install numpy -y
#mafft
mamba install -c bioconda mafft -y
#muscle
mamba install -c bioconda muscle -y

Github

git clone https://github.com/lbobay/CoreCruncher.git
cd CoreCruncher/

python corecruncher_master.py -h

$ python corecruncher_master.py -h

USAGE: python corecruncher_master.py  -in input_folder   -out  output_folder   OPTIONAL: -freq  frequency_across_genomes  -ref ref_genome -prog usearch/blast -id unique/combined   -ext .fa/.fasta/.prot/.faa   -list path_to_genome_list   -score identify_score   -length sequence_length_conservation  -restart yes/no    -align muscle/mafft  -stringent yes/no  -batches 1

Use -h to see the different options

 

Location= 

 

USAGE:

-in      input folder

-out     output folder

 

OPTIONS: 

-freq     Minimum frequency of the gene across genomes to be considered core (default= 90%, an ortholog is considered a core gene even if it is missing in 10% of the set of genomes)

-score    Identity score used by usearch or blast to define orthologs in % (Default= 90)

-length   Minimum sequence length conservation used by to define orthologs (default= 80%)

-prog     Program to use to compare sequences: usearch or blast (default= usearch)

-ref      Reference genome (default: first genome in folder will be used as reference). If you want to specify the reference genome to use, specify the name of the file in the folder (e.g. -ref genome1.prot)

-id       Type of gene IDs in output files. Choose 'unique' if the same gene IDs are not found in different genomes or 'combined' to combine genome ID & gene ID (default= 'combined').  

-ext     File extensions .fa/.fasta/.prot/.faa (default: will try to find it automatically)

-list     Path to a file containing the list of genomes to analyze (default: none, all the genomes in the folder will be analyzed by default)

-restart  Restart analysis from scratch: yes or no (default= no). If yes is chosen, the program will erase the usearch output files and relaunch usearch or blast

-align    Align core gene sequences with specified program (muscle or mafft) and merge all the core genes into a single concatenate. Example=  -align musclev0.0.0 or -align mafft)

-stringent   Define a stringent core genome: yes or no (default= no). By default, core genes with paralogs will be conserved in the core genome and the paralogous sequences will be removed. If stringent is chosen, the core gene will be entirely removed from the core genome)

-batches     Number of batches of genomes. You can divide the analysis in multiple batches of genomes when analyzing large datasets and if your computer can't process all the genomes at once (default= 1, all genomes are anlyzed together)

 

 

テストラン

example/を使う。17  ファイル(protein multi-fasta)ある。

f:id:kazumaxneo:20200909162650p:plain

Serratia_ureilytica_Lr5-4.prot

f:id:kazumaxneo:20200909162717p:plain

ディレクトリを指定して実行する。”-align”をつけるとマルチプルシーケンスアラインメントと連結タンパク質配列作成まで行う。

python corecruncher_master.py -in example -out out_folder -prog usearch -align mafft
  • -freq    Minimum frequency of the gene across genomes to be considered core (default= 90%, an ortholog is considered a core gene even if it is missing in 10% of the set of genomes)
  • -score   Identity score used by usearch or blast to define orthologs in % (Default= 90)
  • -length    Minimum sequence length conservation used by to define orthologs (default= 80%)
  • -prog    Program to use to compare sequences: usearch or blast (default= usearch)
  • -ref    Reference genome (default: first genome in folder will be used as reference). If you want to specify the reference genome to use, specify the name of the file in the folder (e.g. -ref genome1.prot)
  • -align    Align core gene sequences with specified program (muscle or mafft) and merge all the core genes into a single concatenate. Example= -align musclev0.0.0 or -align mafft)

alignフラグなしだと計算は数十秒で終わった。

出力f:id:kazumaxneo:20200909162227p:plain

CC/

f:id:kazumaxneo:20200909162521p:plain

1つ開いてみる。

Serratia_plymuthica_4Rx13.prot-Serratia_sp_AS12.prot

f:id:kazumaxneo:20200909162557p:plain

families_core.txt

f:id:kazumaxneo:20200909162327p:plain

core/

f:id:kazumaxneo:20200909162401p:plain

1つ開いてみる。

fam617.prot

f:id:kazumaxneo:20200909162428p:plain

 

summary.txt

f:id:kazumaxneo:20200909162244p:plain

 

alignフラグありだと最後にconcat.protができる。

 

引用

CoreCruncher: fast and robust construction of core genomes in large prokaryotic datasets
Connor D Harris, Ellis L Torrance, Kasie Raymann, Louis-Marie Bobay
Molecular Biology and Evolution, Published: 04 September 2020 

 

関連