2020 5/28 追記、10/5 タイトル修正、12/24 linuxでのインストール追記
2021 1/13 誤字修正
2024/02/21 追記,2/25追記
GenBank のような公開データベースに登録されているゲノムの数が増え続けていることから、種の遺伝子レパートリーを比較するためのツールが開発されている。このような比較には、種分化イベントの後に共通の祖先から発散したと仮定され、パラログよりも生物間で機能が保存されている可能性が高いオルソログ遺伝子の同定が含まれる（ref.1）。このため、オルソログはゲノムアノテーションや進化研究の重要な要素となっている（ref.2,3）。生命の他のどの領域よりも高速に配列決定されている細菌間では(ref.4)、直交配列を検出するための一般的なヒューリスティックな方法は、単にBLASTの相互ヒットを探すことであり(ref.5, 6)、このタスクのために様々なソフトウェアが利用可能である(ref.7)。これらのツールと増加するゲノム配列を組み合わせることで、最近の研究では、細菌のゲノムは、関心のあるグループのすべての分離株で共有される遺伝子（コアゲノム）と、菌株特異的／部分的に共有される遺伝子を含むモザイクであることを示唆する証拠がいくつか提供されている（ref.8）。コアゲノムとグループ内の残りの遺伝子の合計がパンゲノムと定義される(ref.9)。
ここでは、GNU General Public Licenseの下でリリースされているオープンソースのソフトウェアパッケージであるGET_HOMOLOGUESを紹介する。このソフトウェアは、Linux/Mac OS Xコンピュータシステム上で、系統的に異なる距離にある細菌株のパンゲノムおよび比較ゲノム解析のために特別に設計され、テストされている。このソフトウェアはいくつかの点でユニークである。ゲノムデータのダウンロード、ユーザーが選択した配列特徴の抽出、BLASTおよびHMMERジョブの実行、結果の索引付け、クラスタリング、および解析を含む、完全に自動化された高度にカスタマイズ可能な解析パイプラインを実装している。最新のマルチプロセッサアーキテクチャやコンピュータクラスタを利用して、時間のかかるBLASTやHMMERジョブを並列化することができる。Berkeley DBを使用して一時的なデータをディスクに書き込むか、双方向ベストヒット（BDBH）アルゴリズムのヒューリスティックバージョンを呼び出すことで、適度なマシン（8GB RAM未満）で大規模なデータセット（例えば、101のEscherichia coliゲノムを解析している）を扱うことができる。サポートされている配列クラスタリングアルゴリズムの組み合わせによって復元されたコンセンサスクラスタの計算を含む、遺伝子ファミリーの解析と生成を容易にするための補助スクリプトが統合されている。また、R関数を呼び出して統計解析を行い、コアプロットやパンゲノムプロットなどの結果をグラフィカルに表示するためのスクリプトも用意されている。また、多様な比較ゲノム解析も可能である。最後に、インストール作業を簡単にするためのインストールスクリプトが用意されており、また、ハンズオンチュートリアルを含む非常に詳細なマニュアルも用意されているため、使い勝手の良いパッケージとなっている。
ここでは、公開されている完全なゲノムの中からオルソログを同定するデータベースであるOMA（Orthologous Matrix）の最新版からダウンロードした50のストレプトコッカスゲノムのセットを解析することで、これらの機能の一部を紹介する（ref.10）。著者らはいくつかの理由からこの属を選んだ。この属は非常に高いレベルのゲノム可塑性を示している（ref.11）。最初のパンゲノミクス解析は、Tettelinらの先駆的な研究であるStreptococcus agalactiae (ref.12)で行われ、その後、ヒトの主要な病原体であるS. pyogenes (ref.13)やS. pneumoniae (ref.14)を含むこの属の多様な種について、非常に詳細な比較ゲノミクス研究が行われており、StreptococcusはGET_HOMOLOGUESソフトウェアの優れたテストケースとなっている。
このソフトウェアは、BLAST+（ref.16）とOrthoMCL（バージョン1.4）のコードベース（ref.17）の上に構築されており、3つの一般的なシーケンスクラスタリングアルゴリズムをサポートしている。OrthoMCL（OMCL）、COGtriangles（ref.18）、および双方向ベストヒット（BDBH）アルゴリズムの著者ら自身の実装をサポートしている（論文補足資料の図S1参照）。これらのアプローチは、戦略は異なるがBLASTの双方向ベストヒットを証拠として直交配列を呼び（ref.19）、同一ゲノム内でのベストヒットを持つ遺伝子、すなわち最近のinparalogues（ref.20）を区別している。さらに、タンパク質ドメインのPfamアノテーションを容易にするためにHMMER（http://hmmer.org）が統合されており（ref.21）、オルソログの割り当てを混乱させる可能性のある異なるドメイン構造を持つ配列を含むクラスタをフィルタリングして除外することができる。入力ファイルが GenBank 形式の場合、ヌクレオチド配列とアミノ酸配列の両方のクラスタが生成され、必要に応じて、オルソログ遺伝子に挟まれたオルソログ遺伝子間領域を抽出することもできる。また、GenBankファイル内のゲノム座標は、少なくとも1つの保存された近傍領域を有する同調性クラスタの選択に利用することができる（論文補足資料の図S2参照）。最後に、GenBankファイルには様々なDNAフィーチャーが含まれているため、例えばtRNA遺伝子に着目して検索することも可能である。この場合、デフォルトのBLASTPジョブの代わりにBLASTN検索が実行される。
メインのPerlスクリプトget_homologues.plに加えて、このソフトウェアには、その後の解析に役立ついくつかの補助スクリプトがバンドルされている。例えば、複数のアルゴリズム（COGtrianglesやOMCLなど）によって生成されたインターセクションクラスタは、compare_clusters.plで簡単に選択することができ、ユーザーはコンセンサスクラスタのみを扱うことができる。さらに、パンゲノムとコアゲノムをプロットするためのスクリプト plot_pancore_matrix.pl が提供されている。plot_pancore_matrix.pl により、Tettelinら(9)とWillenbrockら(ref.22)の指数モデルをフィッティングしてコアとパンゲノムのサイズを推定できる。本論文で図示されているように、PHYLIPパッケージ(redf.23)のPARSプログラムを用いて、パンゲノミックマトリクスを表形式とPHYLIP形式で提供しており、 parsimony criterionの下でパンゲノミックツリーを自動生成することができる。parse_pangenome_matrix.plスクリプトは、系統特異的な遺伝子ファミリーや拡張の同定に焦点を当てた比較ゲノミクスや、コアゲノム、クラウドゲノム、シェルゲノムのコンパートメントの計算やグラフ化に役立つ。コア（core）は考慮されているすべてのゲノム/taxaに含まれる遺伝子、ソフトコア（soft core）はKaasと共同研究者の研究（25）で記述されているように考慮されているゲノム/taxaの95％に含まれる遺伝子、クラウド（cloud）は少数のゲノム/Ttaxaにのみ存在する遺伝子（カットオフは、最もぽピュレーションの多い非コアクラスタクラスとそのすぐ隣のクラスとして定義されている）、シェル（shell）は残りの遺伝子で、いくつかのゲノム/taxaに存在している遺伝子と定義される。このスクリプトはまた、Snipenら(ref.26)の二項混合モデルの下でコアおよびパンゲノムサイズの推定値を計算する。ソースのゲノムデータがGenBank形式で提供されている場合、parse_pangenome_matrix.plは、その系統から選択されたリファレンスゲノムの線形化された遺伝地図上にクレード特異的な遺伝子をプロットするように動かすことができる。compare_clusters.plで計算されたパンゲノミックツリーは、parse_pangenome_matrix.plで比較するグループのメンバーを選択するのに便利である。
Manual
GET_HOMOLOGUES
インストール
リリースからmacosのビルドをダウンロードし、macos10.14でテストした。
ubuntu18でcondaを使って導入した。
本体 Github
リリースから、各プラットフォーム向けにGET_HOMOLOGUES-ESTとGET_HOMOLOGUESのバイナリがダウンロードできる。その後、データベースをダウンロードしてインストールする。
#way1 condaでインストールする
mamba create -n gethomologues python=3.8 -y
conda activate gethomologues
mamba install -c conda-forge -c bioconda get_homologues -y
#way2 mac installスクリプトを使用する場合
cd get_homologues-macosx-20200226/
./install.pl
うまく導入されたかテストする。
> ./get_homologues.pl -v
パスが見つかるかどうかチェックされる。リリースからOSXのビルドパッケージをダウンロードし、linuxでランするとblastなどのバイナリが動かないとのエラーが出る。./get_homologues.pl -vで表示されるパスにlinuxでコンパイルされた実行ファイルをコピーする。再び./get_homologues.pl -vをランしてOKになっているか、それともneed to be installedなのか確認する。全てのコマンドが動作することを確認できたら再び./install.plをランしてP-fam、SWISSPROTデータベースをダウンロードする。perlのライブラリが足りなければ、それも導入する必要がある。
> ./get_homologues.pl
$ ./get_homologues.pl
usage: ./get_homologues.pl [options]
-h this message
-v print version, credits and checks installation
-d directory with input FASTA files ( .faa / .fna ), (overrides -i,
GenBank files ( .gbk ), 1 per genome, or a subdirectory use of pre-clustered sequences
( subdir.clusters / subdir_ ) with pre-clustered sequences ignores -c, -g)
( .faa / .fna ); allows for new files to be added later;
creates output folder named 'directory_homologues'
-i input amino acid FASTA file with [taxon names] in headers, (required unless -d is set)
creates output folder named 'file_homologues'
Optional parameters:
-o only run BLAST/Pfam searches and exit (useful to pre-compute searches)
-c report genome composition analysis (follows order in -I file if enforced,
ignores -r,-t,-e)
-R set random seed for genome composition analysis (optional, requires -c, example -R 1234,
required for mixing -c with -c -a runs)
-s save memory by using BerkeleyDB; default parsing stores
sequence hits in RAM
-m runmode [local|cluster] (default local)
-n nb of threads for BLAST/HMMER/MCL in 'local' runmode (default=2)
-I file with .faa/.gbk files in -d to be included (takes all by default, requires -d)
Algorithms instead of default bidirectional best-hits (BDBH):
-G use COGtriangle algorithm (COGS, PubMed=20439257) (requires 3+ genomes|taxa)
-M use orthoMCL algorithm (OMCL, PubMed=12952885)
Options that control sequence similarity searches:
-X use diamond instead of blastp (optional, set threads with -n)
-C min %coverage in BLAST pairwise alignments (range [1-100],default=75)
-E max E-value (default=1e-05,max=0.01)
-D require equal Pfam domain composition (best with -m cluster or -n threads)
when defining similarity-based orthology
-S min %sequence identity in BLAST query/subj pairs (range [1-100],default=1 [BDBH|OMCL])
-N min BLAST neighborhood correlation PubMed=18475320 (range [0,1],default=0 [BDBH|OMCL])
-b compile core-genome with minimum BLAST searches (ignores -c [BDBH])
Options that control clustering:
-t report sequence clusters including at least t taxa (default t=numberOfTaxa,
t=0 reports all clusters [OMCL|COGS])
-a report clusters of sequence features in GenBank files (requires -d and .gbk files,
instead of default 'CDS' GenBank features example -a 'tRNA,rRNA',
NOTE: uses blastn instead of blastp,
ignores -g,-D)
-g report clusters of intergenic sequences flanked by ORFs (requires -d and .gbk files)
in addition to default 'CDS' clusters
-f filter by %length difference within clusters (range [1-100], by default sequence
length is not checked)
-r reference proteome .faa/.gbk file (by default takes file with
least sequences; with BDBH sets
first taxa to start adding genes)
-e exclude clusters with inparalogues (by default inparalogues are
included)
-x allow sequences in multiple COG clusters (by default sequences are allocated
to single clusters [COGS])
-F orthoMCL inflation value (range [1-5], default=1.5 [OMCL])
-A calculate average identity of clustered sequences, (optional, creates tab-separated matrix,
by default uses blastp results but can use blastn with -a recommended with -t 0 [OMCL|COGS])
-P calculate percentage of conserved proteins (POCP), (optional, creates tab-separated matrix,
by default uses blastp results but can use blastn with -a recommended with -t 0 [OMCL|COGS])
-z add soft-core to genome composition analysis (optional, requires -c [OMCL|COGS])
This program uses BLAST (and optionally HMMER/Pfam) to define clusters of 'orthologous'
genomic sequences and pan/core-genome gene sets. Several algorithms are available
and search parameters are customizable. It is designed to process (in a HPC computer
cluster) files contained in a directory (-d), so that new .faa/.gbk files can be added
while conserving previous BLAST results. In general the program will try to re-use
previous results when run with the same input directory.
> ./compare_clusters.pl
$ ./compare_clusters.pl
[options]:
-h this message
-d comma-separated names of cluster directories (min 1 required, example -d dir1,dir2)
-o output directory (required, intersection cluster files are copied there)
-n use nucleotide sequence .fna clusters (optional, uses .faa protein sequences by default)
-r take first cluster dir as reference set, which might contain (optional, by default cluster dirs are expected
a single representative sequence per cluster to be derived from the same taxa; overrides -t,-I)
-s use only clusters with syntenic genes (optional, parses neighbours in FASTA headers)
-t use only clusters with single-copy orthologues from taxa >= t (optional, default takes all intersection clusters; example -t 10)
-I produce clusters with single-copy seqs from ALL taxa in file (optional, example -I include_list; overrides -t)
-m produce intersection pangenome matrices (optional, ideally expects cluster directories generated
with get_homologues.pl -t 0)
-x produce cluster report in OrthoXML format (optional)
-T produce parsimony-based pangenomic tree (optional, requires -m)
> ./plot_matrix_heatmap.sh
$ ./plot_matrix_heatmap.sh
# ERROR: no input tab file defined!
USAGE synopsis for: [plot_matrix_heatmap.sh v.v1.0.4_31Jan18]:
plot_matrix_heatmap.sh <string (tab-delimited matrix file)> [-t|-m|-v|-o|-p|-H|-W|-C|-r|-k|-N|-M]
REQUIRED
-i <string> *Avg_identity.tab file
OPTIONAL:
* filter out excessive redundancy in the tab-delimited ANI matrix file
-c <float> (maximum) identity cut-off value (e.g. 97.3) [def: 100]
-d <int> maximum number of decimals in matrix display [0-2; def: 0]
* tweak the graphical output:
-a <integer> angle to rotate node angles [def 45]
-t <string> text for plot title [def input_tab_file_name]
-m <integer> margins_horizontal [def 18]
-v <integer> margins_vertical [def 18]
-o <string> output file format [def svg]
-p <integer> points for plotting device [def 15]
-H <integer> output device height [def 10]
-W <integer> output device width [def 15]
-C <flag> do not reorder clusters [def reorder clusters and plot dendrogram]
-r <flag> remove column names and cell contents [def names and cell contents are printed]
-k <string> text for scale X-axis [def "Value"]
-b <integer> right border (margin) in dendrograms [def 10]
-X <float> character expansion factor [def 1.0]
* Goodness of clustering
-K <integer> max. number of clusters [def: ; NOTE: 2 <= K <= n-1 ]
RUN NJ ANALYSIS using ANI matrix (average identity matrix) generated by get_homologues.pl -A -t 0 -M|G
-N <flag> will compute a distance matrix (ANdist) from the input similarity matrix
and use the former to construct a NJ tree and save it as newick string,
as well as a complete linkage dendrogram for convenient cluster delimitations
at ANdist cutoffs of 6,5 and 4, which correspond to ANI values of 94%, 95% and 96%, respectively
Do not use with a binary presence-absence matrix
SUBSET MATRIX WITH REGULAR EXPRESSIONS
-x <string> regex, like: 'species1|species2' [def ]
MORE DETAILED HELP:
-M <flag> prints gplot installation instructions and further usage information
EXAMPLE:
plot_matrix_heatmap.sh -i Avg_identity.tab -c 99.5 -d 1 -t "Genus X ANIb (OMCL all clusters)" -N -o pdf -m 22 -v 22 -p 20 -H 20 -W 30 -x 'sp1|sp2' -a 45 -b 12 -X 0.9 -K 20
#------------------------------------------------------------------------------------------------------------------
AIM: Plot ordered heatmaps with row and col. dendrogram, from squared numeric (distance or presence-absence) matrix,
like the *Avg_identity.tab matrices produced by get_homologues.pl with the -A flag
OUTPUT: svg|pdf files
DEPENDENCIES:
R packages ape and gplots. Run plot_matrix_heatmap.sh -M for installation instructions.
NOTES:
1. To improve the display, shorten the genome/taxon names as much as possible!
2. Use -m and -v to control horizontal and vertial margins to the plot
to fit the taxon labels, and -X character_expansion_factor to control font size
#------------------------------------------------------------------------------------------------------------------
# Run check_dependencies() ... looks good: R is installed.
>./plot_pancore_matrix.pl
$ ./plot_pancore_matrix.pl
[options]:
-h this message
-i input .tab data file (required, generated by get_homologues.pl)
-f type of fit [pan|core_Tettelin|core_Willenbrock|core_both]
(default: core_Tettelin, PubMed:16172379,18088402)
-F font scale for fitted formulas (optional, default=0.8)
-a save snapshots for animations in dir (optional, example: -a animation)
(2020/12/24)
debianでdocker imageを挙げておきます。すぐに使えます。
docker pull kazumax/get_homologues:latest
cd /home/get_homologues-x86_64-20200226/
#run check
perl get_homologues.pl -v
condaで導入した方がよいと思います。
実行方法
解析するゲノムのGenbankファイル（*１）またはアミノ酸のfastaファイルを集め、そのディレクトリを指定する。"-d"で指定するgenbankファイルは、gzip圧縮した状態で提供してもよい。
#default settings
perl get_homologues.pl -d genbank_dir/ -n 12
#orthoMCL使用。inparalogue含めず、10genome中９ゲノムヒット
perl get_homologues.pl -d genbank_dir/ -M -e -t 9 -n 12
- -n number of threads for BLAST/HMMER/MCL in 'local' runmode (default=2)
- -G use COGtriangle algorithm (COGS, PubMed=20439257) (requires 3+ genomes|taxa)
- -M use orthoMCL algorithm (OMCL, PubMed=12952885) Options that control sequence similarity searches:
- -X use diamond instead of blastp (optional, set threads with -n)
- -e exclude clusters with inparalogues (by default inparalogues are included)
- -t report sequence clusters including at least t taxa (default t=numberOfTaxa,
t=0 reports all clusters [OMCL|COGS])
テストデータを用いた解析の流れ
１、データのダウンロードとget_homologues.plのラン。
#1 download
wget https://github.com/eead-csic-compbio/get_homologues/releases/download/v3.3.2/get_homologues-x86_64-20200226.tgz
tar -xvf get_homologues-x86_64-20200226.tgz && rm get_homologues-x86_64-20200226.tgz
cd get_homologues-x86_64-20200226/
#2-1 run get_homologues(BDBH algorithm)
#4genomeのアミノ酸faaファイルが収められているsample_buch_fastaを対象にget_homologuesをラン
./get_homologues.pl -n 12 -d sample_buch_fasta .
#run後、カレントにsample_buch_fasta_homologues/ができる。
#2-2 run get_homologues(COG triangles algorithm)
./get_homologues.pl -n 12 -d sample_buch_fasta -G
#2-3 run get_homologues(OMCL algorithm)
./get_homologues.pl -n 12 -d sample_buch_fasta -M
#3 venn diagram(#4のため-Tと-mも指定)
./compare_clusters.pl -o sample_intersection -T -m -d sample_buch_fasta_homologues_BDBH/BuchaphCc_f0_alltaxa_algBDBH_e0_,sample_buch_fasta_homologues_COG/BuchaphCc_f0_alltaxa_algCOG_e0_,sample_buch_fasta_homologues_OMCL/BuchaphCc_f0_alltaxa_algOMCL_e0
#3の出力
２、コア遺伝子マトリクスファイルの出力
compare_clusters.plをサンプルと"-m -T"のみ指定して実行すると、閾値を満たしたコア遺伝子のgene presence/absenceを示したmatrxixファイルとその系統樹ファイルが得られる。
#get_homologues.plの出力ディレクトリにできるサブディレクトリ（日付で始まる）を指定する
perl compare_clusters.pl -o sample_intersection -T -m -d GET-HOMOLOGUES_output/xxx_algBDBH_e0_/
- -m produce intersection pangenome matrices (optional, ideally expects cluster directories generated with get_homologues.pl -t 0)
- -T produce parsimony-based pangenomic tree (optional, requires -m)
-
-d comma-separated names of cluster directories (min 1 required, example -d dir1,dir2)
-
-o output directory (required, intersection cluster files are copied there)
#4 heatmap (Rのライブラリが必要、ないなら導入）
install.packages("gplots")
install.packages("dendextend")
install.packages("factoextra")
install.packages("ape")
#run
./plot_matrix_heatmap.sh -i sample_intersection/pangenome_matrix_t0.tab -o pdf \
-r -H 8 -W 14 -m 28 -t "sample pangenome (clusters=180)" -k "genes per cluster"
#4の出力(別のサンプルデータ)
#5 cloud, shell、 core遺伝子の計算
./parse_pangenome_matrix.pl -m sample_intersection/pangenome_matrix_t0.tab -s
#5の出力
#6
./get_homologues.pl -d sample_buch_fasta -c -n 12
#7
./plot_pancore_matrix.pl -i sample_buch_fasta_homologues/core_genome_algBDBH.tab -f core_Tettelin
- -c report genome composition analysis (follows order in -I file if enforced,
ignores -r,-t,-e)
#ANI推定
./get_homologues.pl -d sample_plasmids_gbk -a 'CDS' -A -t 0 -M
cd sample_plasmids_gbk_homologues
../plot_matrix_heatmap.sh -i EcoliST131plasmidpKC394_f0_0taxa_CDS_algOMCL_e0_Avg_identity.tab \
-d 2 -t "CDS ANI matrix with 2 decimals"
- -A calculate average identity of clustered sequences, (optional, creates tab-separated matrix, by default uses blastp results but can use blastn with -a recommended with -t 0 [OMCL|COGS])
- -a report clusters of sequence features in GenBank files (requires -d and .gbk files, instead of default 'CDS' GenBank features example -a 'tRNA,rRNA', NOTE: uses blastn instead of blastp, ignores -g,-D)
追記
GET HOMOLOGUESの実践的な使い方についての論文です。分類が難しいLactobacillus属（乳酸菌）の菌株をコア遺伝子の存在パターンに基づいて分類を試みたという論文です。コア遺伝子探索とクラスタリングにGET HOMOLOGUESが使用されています。
https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-018-4601-5
引用
GET_HOMOLOGUES, a Versatile Software Package for Scalable and Robust Microbial Pangenome Analysis
Bruno Contreras-Moreira, Pablo Vinuesa
Appl Environ Microbiol. 2013 Dec; 79(24): 7696–7701
*１
DFASTでアノテーション付けすると上手く認識したが、PROKKAでアノテーション付けしたGenBnakファイルを使うとエラーになった。prodigalを使って得たgenbankファイルも認識しなかった（proteinのfastaは問題なく認識した）。
#自分用メモ
#pan-genome解析でコア遺伝子を探索。50genomeとしたら50で全ゲノムに存在。refも別途指定しているなら(refもgenbank-dirに含めないと解析されないので注意)
./get_homologues.pl -d genbank-dir/ -M -t 50 -r ref.gbk -e -n 20
- -r reference proteome .faa/.gbk file (by default takes file with least sequences; with BDBH sets first taxa to start adding genes)
-
-t report sequence clusters including at least t taxa (default t=numberOfTaxa, t=0 reports all clusters [OMCL|COGS])
