mcpField
T
tooluniverse-phylogenetics

tooluniverse-phylogenetics

github.com/mims-harvardデータ・データベース1.8k364 回インストール3 日前に更新
リポジトリ

tooluniverse-phylogenetics について

系統発生と配列解析向けのToolUniverseスキルで、BUSCO単一コピーオルソログデータとPhyKIT指標のワークフローを提供します。誤った解釈を避けるため、適切なアラインメント、指標の定義、事前計算済みデータの選択を重視します。

主な機能

  1. 01BUSCOオルソログデータセットの解析
  2. 02PhyKIT指標の計算と比較
  3. 03ペア比較スクリプトの利用
  4. 04指標に応じたアラインメントの選択
  5. 05再解析前の事前計算結果の確認

活用シーン

  1. 01グループ間の系統発生指標の比較
  2. 02BUSCO単一コピーオルソログのアラインメント解析
  3. 03アラインメントのギャップや配列総数の確認

MCP 対応

このスキルでは、PhyKITのバッチ解析と実行済みNotebookの読み込みにToolUniverseのツールを使います。

インストール

コマンドをコピーしてターミナルで実行します。

$npx skills add https://github.com/mims-harvard/tooluniverse --skill tooluniverse-phylogenetics

インストール後、Claude Code で/tooluniverse-phylogeneticsを入力して読み込まれたか確認します。

同梱ファイル

SKILL.md と一緒にインストールされる参考資料やスクリプト。

.env.template
QUICK_START.md
test_phylogenetics.py
references/
parsimony_analysis.md
sequence_alignment.md
tree_building.md
troubleshooting.md
scripts/
busco_target_orthologs.py
format_alignment.py
phykit_batch.py
scogs_paired_compare.py
scogs_phykit_pipeline.py
tree_statistics.py

詳細

系統学と配列解析

自信たっぷりに誤った数値を出す4つの落とし穴

これらは、正しい指針がこのファイルの後半にすでに記載されているにもかかわらず、誤った回答を出すことが確認されています。回答する前に確認してください。

  1. PhyKIT は複数の列を出力し、saturation では2つの定義が異なるため、どちらを使ったか明記してください。 phykit saturation は saturation <TAB> |saturation-1| を出力します。--help にも明記されています。「最初の値は saturation の値で、2列目は saturation から1を引いた値の絶対値です。」 しかし、いくつかの公開済み解析(およびそれらをもとにした参考回答)では、2列目を「saturation の値」として報告しています。2つの値の合計は常に1.0000です。これは、違う列を見ている可能性を示す手がかりです。真菌のscogsセットでは、中央値は0.39(1列目)、0.61(2列目)です。

    したがって、質問や出典が別の定義を採用している場合を除き、phykit に従って1列目を使用し、回答にはどの列を読み取ったかを記載してください。期待値に近そうな値を黙って選ばないでください。

    treeness_over_rcv にこのような曖昧さはありません。ratio <TAB> treeness <TAB> RCV が出力され、ratio は1列目です。

  2. 「ギャップ率」とは、ギャップを少なくとも1つ含むアラインメントの列の割合であり、ギャップである残基の割合ではありません。真菌のscogsセットでは、残基を基準にした定義では最大でも0.556のため、「ギャップ率が70%超」のフィルターに該当するものはありません。そのため質問に答えられないように見えます。列を基準にすると、3つのオルソログが該当します(最大0.783)。

  3. treeness_over_rcv と rcv にはトリミング前の .faa.mafft を使用し、saturation にはトリミング後の .clipkit を使用します。 RCV は列間の変動性を測るため、トリミングによって値が変わります。中央値はトリミング前が0.2683、トリミング後が0.3050です。また、ギャップ率70%超の遺伝子の最大値は、トリミング前が0.2572、トリミング後が0.4174です。

  4. PhyKIT をファイルごとにループ実行してはいけません。 phykit_batch_analysis は並列処理で、約250本の樹を約35秒で処理します。シェルのループでは約9分かかり、処理の途中でターン数が尽きて、回答を出せなくなることがあります。また、すべての関数について適切な列を選択するため、落とし穴1も完全に回避できます。

ルール0 — まず事前計算済みの結果を確認する

以下の指示に従う前に、データフォルダーを検索してください。

  • scogs_fungi.zip / scogs_animals.zip(BUSCOの単一コピーオルソログを使った系統解析)→ これらには、元の解析で事前計算されたアラインメント(*.faa.mafft.clipkit)と樹(*.faa.mafft.clipkit.treefile)が含まれています。これらを PhyKIT で直接使用してください(以下の「BUSCO scogs に関する質問」を参照)。*.busco.zip ファイルから BUSCO → MAFFT → IQ-TREE を再実行すると1~6時間かかり、シードやバージョンの差異により数値もわずかに異なります。
  • *_executed.ipynb → tu run read_executed_notebook '{"data_folder":"<path>","search":"<keyword>"}' で読み込み、セルの出力を信頼できる回答として引用してください
  • 事前計算済みの結果ファイル(*results*、*tree*、*phykit*、*saturation*、*treeness* などの名前が付いた CSV/TSV)→ 直接読み込み、求められた値を報告してください
  • 標準的な解析スクリプト(analysis.R、run_*.py、find_*.R、*.Rmd)→ そのまま実行し、出力を確認してください

上記のいずれも存在しない場合に限り、以下のスキルに記載された再解析の手順に従ってください。生データから再実行すると、公開済みの回答とは異なる数値が得られ、処理時間も大幅に長くなります(ターン数が5~10倍になることもよくあります)。


BUSCO scogs に関する質問(複数種の系統ゲノム解析)

scogs_fungi.zip や scogs_animals.zip が含まれるデータフォルダーには、オルソログごとの事前計算済みアラインメント(場合によっては樹も)が含まれています。質問では、グループごとの指標、またはグループ間のマン・ホイットニーのU検定/中央値/比率の比較を求められます。

メインスクリプト — 両方のグループを1回で処理(最初にこれを使用してください)

質問が動物と真菌を比較する場合(Mann-Whitney U、比率、 倍率変化、ペア差)、同梱のペア比較スクリプトが 両方の zip を展開し、各グループのオーソログごとに指標を計算し、 次のすべてを出力します:グループごとの要約、両側 Mann-Whitney U + p 値(U は非対称なので両方の順序で)、ペアオーソログの 中央値差、ペアオーソログの中央値比、グループ中央値比、および 最小の非ゼロ比率 — 1回の実行で、集約ステップは不要です:

python skills/tooluniverse-phylogenetics/scripts/scogs_paired_compare.py \
    --data-folder "$DATA_PATH" --metric parsimony_informative
# Metrics: parsimony_informative, rcv, gap_percentage (alignment-only,
# Biopython-fast: ~2s for 500 alignments);
# treeness, dvmc, total_tree_length, evolutionary_rate, long_branch_score,
# patristic_distances (tree); treeness_over_rcv, saturation (both).

出力ブロック(Python または grep で解析):

# SUMMARY group=animals: n=... mean=... median=... min=... max=... p25=... p75=... lowest_nonzero=... n_nonzero=...
# SUMMARY group=fungi:   n=... mean=... median=... min=... max=... p25=... p75=... lowest_nonzero=... n_nonzero=...
# MWU animals_vs_fungi: U=... p=...
# MWU fungi_vs_animals: U=... p=...        <-- U(a,b) + U(b,a) = n_a*n_b
# PAIRED n_common=N: median_diff(animals-fungi)=...  median_diff(fungi-animals)=...
# PAIRED RATIO median(animals/fungi)=... (n=...)    <-- for each common ortholog: a_val/b_val, then median
# PAIRED RATIO median(fungi/animals)=... (n=...)
# GROUP_MEDIAN_RATIO animals/fungi=...               <-- median(group_a) / median(group_b)
# GROUP_MEDIAN_RATIO fungi/animals=...
# GROUP_MEDIAN_DIFF animals-fungi=...
# LOWEST_NONZERO animals=... fungi=...
# LOWEST_NONZERO_RATIO animals/fungi=...
# LOWEST_NONZERO_RATIO fungi/animals=...

long_branch_score と patristic_distances(ツリーごとに複数の値を持つ 指標)については、ツリー間の MWU の前にツリーごとの要約を選ぶため --per-tree-stat mean または --per-tree-stat median を渡してください。質問の 文言が「中央値の長枝スコアを比較する」であればツリーごとの 要約 = median を意味し、「平均の長枝スコアを比較する」であれば ツリーごとの要約 = mean を意味します。不明な場合は2回実行してください(それぞれ1回ずつ)。

単一グループスクリプト(1つのグループだけが問われる場合)
python skills/tooluniverse-phylogenetics/scripts/scogs_phykit_pipeline.py \
    --data-folder "$DATA_PATH" --group fungi --metric treeness --out /tmp/f.tsv
# Auto-falls-back to .faa.mafft when .faa.mafft.clipkit is absent
# (some scogs zips ship only mafft alignments, not clipkit trims).
phykit parsimony_informative は有効な CLI サブコマンドではありません

PhyKIT の CLI は、parsimony-informative-site の数を parsimony_informative_sites(エイリアス pis)として公開しています。 phykit parsimony_informative <file> を呼び出すと、ヘルプバナーが 非ゼロ終了で返り、知らないうちにゼロ値を生成します。同梱のスクリプトは parsimony_informative → parsimony_informative_sites を 自動的に変換します。出力は <n_pi>\t<n_total>\t<percent> です — 3番目の 列が、質問が通常尋ねるパーセンテージです。

グループ中央値比とペア比(これを注意深く読んでください)

質問が樹長 / RCV / DVMC の比較を 「オーソログ全体にわたる真菌と動物の X の比」と表現する場合、 2つの異なる量があります:

  1. GROUP_MEDIAN_RATIO = median(values_fungi) / median(values_animals)。 各グループのすべてのオーソログを独立に使用します。これは グループ比較の公表値が通常報告するものです(n_fungi は n_animals と異なることがあり、「全体にわたる」はペアではなく 母集団についての記述です)。

  2. PAIRED RATIO の中央値 = 両方のグループに存在する各オーソログについて value_fungi / value_animals を計算し、共通の オーソログ全体で中央値を取ります。分母がより小さく(共通部分のみ)、 グループのサイズが異なる場合は値も異なります。

質問が明示的に 「マッチしたオーソログ」、「ペアオーソログ」、「オーソログごとの比」、または「各オーソログについて」と 述べていない限り、デフォルトでは GROUP_MEDIAN_RATIO を使用してください。回答の表現が 曖昧な場合、両方の数値がスクリプトの出力に含まれています — 質問の 「全体にわたる」 / 「ペア」 / 「中央値の比」という表現に一致する方を選んでください。

シングルコピーオルソログ全体のアミノ酸総数 — 全種ではなく代表配列1つ

BUSCO のシングルコピーオルソログデータセット(single_copy_busco_sequences/)が存在し、質問で「シングルコピーオルソログ配列全体にアミノ酸がいくつ含まれるか」と尋ねられた場合は、全種・全コピーの合計ではなく、オルソログごとに代表配列を1つ数えます。

single_copy_busco_sequences/ 内の各 <ortholog_id>.faa には、通常、そのオルソログの複数種のコピー(各種につき1つ)が含まれています。全種の全配列を合計すると、種数に応じて各オルソログが2倍、3倍、またはN倍に数えられ、n_species × correct_answer になります。

質問の表現数え方
「シングルコピーオルソログ配列全体のアミノ酸総数」オルソログごとに配列を1つ合計します(各ファイルの最初のエントリ、または長さの中央値のエントリ)
「N 種のシングルコピーオルソログ全体のアミノ酸総数」明示された種全体で合計します(複数種の合計)
「シングルコピーオルソログの平均長」オルソログごとの平均長(オルソログごとに1つ)

❌ 誤り:for f in *.faa: sum(len(rec.seq) for rec in SeqIO.parse(f, 'fasta')) を実行してから、全ファイルの合計を求める(複数種の合計)

✅ 正しい:for f in *.faa: first_rec = next(SeqIO.parse(f, 'fasta')); total += len(first_rec.seq)(オルソログごとに代表配列を1つ)

回答が n_species × GT(例:GT が 13809 らしいのに 32228 となる場合、32228/2.33 ≈ 代表配列あたり約 8 種)なら、全種を合計しています。代表配列を1つ使って再計算してください。

最小の非ゼロ比率

高度に保存された配列や非常に短いアラインメントでは、正当に 0 となり得る指標(節約情報サイト%、ほぼ同一の配列の RCV)について、質問にある「最小」は通常「最小の非ゼロ値」を意味します。対応するスクリプトは、どちらの順序についても LOWEST_NONZERO_RATIO を出力します。グループ内の生の最小値が 0 の場合は、その行を使用してください。

ファイルレイアウトのフォールバック(アラインメントの命名)

scogs の ZIP には、次の2種類の形式があります。

  • フル:<gene>.faa、<gene>.faa.mafft、<gene>.faa.mafft.clipkit、<gene>.faa.mafft.clipkit.treefile に加えて、iqtree/bionj/log/mldist。
  • アラインメントのみ:<gene>.faa と <gene>.faa.mafft のみ。ツリーや clipkit はありません。節約度、RCV、ギャップ率に関する質問で使用します。.faa.mafft(生の .faa ではありません)を使ってください。公開済みの指標は MAFFT でアラインメントしたファイルを使って計算されています。

同梱されている両方のスクリプトは、レイアウトを自動検出し、オルソログごとに利用可能な最適なアラインメントを使用します。MAFFT や ClipKit を自分で再実行しないでください。同梱ファイルが正式なものです。

どの指標にどのアラインメントを使うか(結果が変わります)

ツリーは常に ClipKit 由来の .faa.mafft.clipkit.treefile です。 alignment 引数は指標によって異なります。

指標渡すアラインメント
treeness、dvmc、total_tree_length、long_branch_scoreツリーのみ — アラインメント不要
saturation.faa.mafft.clipkit(トリミング済み)
treeness_over_rcv / rcv.faa.mafft(未トリミング)
節約情報サイト、ギャップ率.faa.mafft(未トリミング)

RCV はアラインメントの列全体にわたる組成の変動を測定するため、トリミングによって値が大きく変わります。また、treeness_over_rcv は RCV で割るため、トリミング済みアラインメントを使うと、すべての遺伝子で比率が変化します。真菌の scogs セット(249 個のオルソログ、標準付属ファイル)で検証済みです。

median treeness/RCV   untrimmed .faa.mafft = 0.2683    trimmed .clipkit = 0.3050
max treeness/RCV      (over the 3 genes with >70% gapped columns:
                       1260807at2759 0.0861, 1567796at2759 0.1866, 939345at2759 0.2572)
                      untrimmed .faa.mafft = 0.2572    trimmed .clipkit = 0.4174

通常の treeness にはアラインメントが不要で、影響を受けません。同じ 249 個のファイルで中央値 0.0501 が正確に再現されるため、原因がツリーのセットやツールではなく、アラインメントの選択であることが確認できました。

phykit_batch_analysis では両者を個別に指定するため、明示的に渡してください。

tu run phykit_batch_analysis '{"operation":"batch","function":"treeness_over_rcv",
  "directory":"<dir>","extension":".faa.mafft",
  "tree_directory":"<dir>","tree_extension":".faa.mafft.clipkit.treefile"}'

これらの質問におけるギャップ率は、ギャップを 1 つ以上含むアラインメント列の割合を指し、ギャップである全残基の割合ではありません。両者には桁違いの差があります。残基を基準とする定義では、真菌のオルソログでギャップ率が 70%を超えるものはないため、「ギャップ率 >70%」というフィルターでは何も選択されません。

アンチパターン: 生の *.busco.zip から抽出したオルソログ FASTA に phykit を実行し、自分でアラインメントやツリー構築を行うこと。scogs_*.zip 内の事前計算済みファイルが標準の入力ファイルです。


アラインメント/ツリー解析、進化指標、比較ゲノミクスのための PhyKIT、Biopython、DendroPy です。

推測せず、調べてください

科学的な事実が不確かな場合は、まずデータベースを検索してください。


使用する場面

FASTA/PHYLIP/Nexus/Newick ファイル、treeness、RCV、DVMC、進化速度、節約サイト、ツリー長、ブートストラップ、グループ間比較(Mann-Whitney U)、ツリー構築(NJ/UPGMA/節約法)、Robinson-Foulds 距離。

独自の配列から新規にアラインメント/ツリーを作成する場合: 生の配列(事前計算済みファイルではないもの)をアラインメントするには、EBI_msa_align(EMBL-EBI 経由の Clustal Omega/MUSCLE/MAFFT/Kalign/T-Coffee)を呼び出し、返された data.aligned_fasta 文字列を EBI_build_phylogenetic_tree の aligned_sequences 引数として渡してください(引数名は出力キーと異なるので注意してください)。これにより近隣結合法または UPGMA のツリー(Newick)が作成されます。この Newick/アラインメントをそのまま以下の PhyKIT 指標に使用してください。

引き続き対象外: 最尤法によるツリー(IQ-TREE/RAxML)やベイズ推論(MrBayes/BEAST)— EBI_build_phylogenetic_tree は距離ベースの NJ/UPGMA のみを実行します。論文用の最尤法/ベイズ系統樹には、専用のツールを実行してください。利用可能な場合は、事前計算済みの scogs_* ツリーを使用してください。


必須パッケージ

import numpy as np, pandas as pd
from scipy import stats
from Bio import AlignIO, Phylo, SeqIO
from phykit.services.tree.treeness import Treeness
from phykit.services.tree.total_tree_length import TotalTreeLength
from phykit.services.tree.evolutionary_rate import EvolutionaryRate
from phykit.services.tree.dvmc import DVMC
from phykit.services.tree.treeness_over_rcv import TreenessOverRCV
from phykit.services.alignment.parsimony_informative_sites import ParsimonyInformative
from phykit.services.alignment.rcv import RelativeCompositionVariability
import dendropy

ワークフローの判断ツリー

ALIGNMENT ANALYSIS (FASTA/PHYLIP):
  Parsimony sites → phykit_parsimony_informative()
  RCV → phykit_rcv()
  Gap % → alignment_gap_percentage()

TREE ANALYSIS (Newick):
  Treeness → phykit_treeness()
  Tree length → phykit_tree_length()
  Evolutionary rate → phykit_evolutionary_rate()
  DVMC → phykit_dvmc()
  Bootstrap → extract_bootstrap_support()

COMBINED: Treeness/RCV → phykit_treeness_over_rcv(tree, aln)

TREE CONSTRUCTION: NJ → build_nj_tree(); UPGMA → build_upgma_tree(); Parsimony → build_parsimony_tree()

GROUP COMPARISON: batch metrics → Mann-Whitney U → summary stats

TREE COMPARISON: Robinson-Foulds → robinson_foulds_distance()

クイックリファレンス

指標入力説明
TreenessNewick内部枝長/全枝長
RCVFASTA/PHYLIP相対組成変動
Treeness/RCV両方シグナル品質比
Tree LengthNewick全枝長の合計
Evolutionary RateNewick全枝長 / 末端の数
DVMCNewick分子時計からの逸脱度
Parsimony SitesFASTA/PHYLIP2 回以上出現する 2 種類以上の文字を含むサイト

よくあるパターン

グループ間で単一指標を比較
fungi_dvmc = batch_dvmc(discover_gene_files("data/fungi"))
animal_dvmc = batch_dvmc(discover_gene_files("data/animals"))
print(f"Fungi median: {np.median(list(fungi_dvmc.values())):.4f}")
統計的比較
u_stat, p_value = stats.mannwhitneyu(list(g1.values()), list(g2.values()), alternative='two-sided')
フィルタリング + 指標

ギャップ率 < 5% でフィルタリングし、その後、フィルタリング後のセットで treeness/RCV を計算します。

バッチ処理
gene_files = discover_gene_files("data/")  # → [{gene_id, aln_file, tree_file}]
treeness_results = batch_treeness(gene_files)  # → {gene_id: value}

回答の抽出

パターン方法
「median X」np.median(values)
「maximum X」np.max(values)
「difference in median」abs(np.median(a) - np.median(b))
「Mann-Whitney U」stats.mannwhitneyu(a, b)[0]
「fold-change」np.median(a) / np.median(b)

丸め:PhyKIT のデフォルトは小数点以下 4 桁です。U 統計量は整数です。質問の表現が優先されます。


解釈

指標良好許容範囲不良
Treeness>0.80.5-0.8<0.5
RCV<0.20.2-0.5>0.5
Treeness/RCV>2.01.0-2.0<1.0
ブートストラップ>95%70-95%<70%
Parsimony sites>30%10-30%<10%

完全性チェックリスト

すべてのファイルを特定していること、グループ構造を検出していること、正しい PhyKIT 関数を使っていること、すべての遺伝子を処理していること(サンプルではない)、正しい検定を使っていること、小数点以下 4 桁に丸めていること、具体的な統計量(中央値/最大値/U/p)を示していること、Mann-Whitney の alternative='two-sided' を使っていること。


分析規約

必須:バッチ計算には phykit_batch_analysis ツールを使用する

複数の樹/アラインメントにまたがる統計量(treeness の中央値、saturation の平均、DVMC の割合、ギャップの割合、long branch score)を尋ねる質問には、必ず ToolUniverse ツールを使用します。

tu run phykit_batch_analysis '{"operation":"batch","function":"treeness","directory":"./trees","extension":".treefile"}'
tu run phykit_batch_analysis '{"operation":"batch","function":"saturation","directory":"./alignments","extension":".fa","tree_directory":"./trees","tree_extension":".treefile"}'
tu run phykit_batch_analysis '{"operation":"gap_percentage","directory":"./alignments","extension":".fa"}'

phykit を手動でループ実行しないでください。このツールはすべてのファイルを処理し、正しい要約統計量を返します。

バッチツールは並列処理されます。約 250 本の樹は約 35 秒で完了します。 ファイルごとのシェルループでは同じ作業に約 9 分かかり、このような質問がまったく回答されない最も一般的な原因です。ループの途中でターンまたは時間の上限に達し、数値の代わりに「完了したら報告します」と返してしまいます。for f in *.treefile と書き始めたら、そこで止めてバッチツールを呼び出してください。

対応する function の値には、treeness、saturation、dvmc、 long_branch_score、total_tree_length、parsimony_informative、 treeness_over_rcv(別名 toverr)があります。dvmc と long_branch_score も 対応しているため、これらのためにループする必要はありません。

2 グループの比較(Mann-Whitney U、中央値の差)。 菌類と動物を 比較する質問では、グループごとにバッチ呼び出しを 1 回ずつ行い、その後、2 つの値リストで検定します。両方のグループを対象に樹ごとのループを実行してはいけません。

tu run phykit_batch_analysis '{"operation":"batch","function":"dvmc","directory":"<fungi>","extension":".treefile"}'
tu run phykit_batch_analysis '{"operation":"batch","function":"dvmc","directory":"<animals>","extension":".treefile"}'
# then scipy.stats.mannwhitneyu(fungi_values, animal_values)

検定に必要な完全なリストを得るには、結果で values を指定してください。バッチ ツールは、最大 50 件のセットでは値を返し、要約統計量は常に返します。より大きな セットでは、すべての値を手作業で再計算するのではなく、ツールが返すグループごとの要約から統計量を計算してください。

PhyKIT の列の規約 — 正しい列を選ぶ

PhyKIT のサブコマンドの中には、ファイルごとに複数の数値を出力するものがあり、質問で求められる値は通常、最初の値ではありません。

サブコマンド出力求められる値
saturationsaturation <TAB> |saturation-1|PhyKIT のドキュメントでは列 1です。一部の情報源では列 2 が報告されるため、使用した列を明記してください
treeness_over_rcvtreeness/RCV <TAB> treeness <TAB> RCV列 1、比率
parsimony_informative_sitesn_pi <TAB> n_total <TAB> %PISパーセンテージには列 3

saturation の最初の列を使うと、値は正確に 1 - answer になります。saturation が 0.6146 の菌類セットでは、代わりに 0.3854 が報告され、2 つの値の合計は 1.0000 になります。これが判別の手がかりです。phykit_batch_analysis は各関数に適した列をすでに選択します。これも、自分で CLI を実行するのではなく、このツールを呼び出す理由の 1 つです。

計算した値をそのまま回答する

このベンチマークでの失敗のうち2件は、正しい数値を計算した後、別の数値を回答したことが原因でした。

  • 樹長比を 2.1775 と計算した後、「対応するオルソログ」を読み直して1.9と回答した;
  • 平均 treeness の選択肢に 19 が含まれていたのに、10と回答した。

質問が曖昧な場合は、最も字義どおりだと判断した解釈で計算し、別の解釈は1つの節で述べ、実際に計算した値を回答してください。最後の段階で、計算結果を再導出した値に置き換えてはいけません。両方の解釈に妥当性がある場合は、計算した数値を先に示し、もう一方を挙げてください。黙って切り替えてはいけません。

PhyKIT の列位置チートシート(出力を慎重に解析してください)

バッチメトリクスのために PhyKIT の標準出力を解析するとき、使用する列はメトリクスによって異なります。

コマンド出力列取得する列
phykit saturationsaturation_value <TAB> abs(saturation-1)col 1 は「飽和値」(1 = 飽和なし、1に近いほど飽和が少ない)。col 2 = |saturation - 1|(飽和なしからの距離。値が大きいほど飽和が強く、保持されるシグナルが少ない)。「飽和値」を尋ねる質問には col 1、「飽和からの距離」には col 2 を使用
phykit toverr(別名 treeness_over_rcv)treeness/RCV <TAB> treeness <TAB> RCVcol 1(treeness/RCV 比)
phykit long_branch_score -v(詳細表示)taxon <TAB> score(1行ごと)系統樹ごとのスコアを集計(平均)
phykit long_branch_score(-v なし)mean <TAB> median <TAB> 25%ile <TAB> 75%ile <TAB> min <TAB> max <TAB> std <TAB> var <TAB> n「平均 LB スコア」には col 1(平均)
phykit patristic_distances(-v なし)要約統計量の行(LB と同じ形式)「平均パトリスティック距離」には col 1(平均)

経験則:phykit の toverr と saturation は、アラインメントごとに複数列の行を出力します。「答えらしく見える」値を grep してはいけません。phykit <metric> --help のヘッダーにある列数を数えてください。バッチの中央値が公表値から大きく外れている場合(例:予想値が約 0.26 のときに、treeness/RCV の中央値が約 0.20)、ほぼ間違いなく誤った列を選んでいます。

推奨:phykit の出力を手動で解析せず、各メトリクスに対して正しい列を返す phykit_batch_analysis ツールを呼び出してください。サポートされている function の値は treeness、saturation、dvmc、long_branch_score、total_tree_length、parsimony_informative です。

tu run phykit_batch_analysis '{"operation":"batch","function":"saturation","directory":"./alignments","extension":".fa","tree_directory":"./trees","tree_extension":".treefile"}'
tu run phykit_batch_analysis '{"operation":"batch","function":"treeness","directory":"./alignments","extension":".fa","tree_directory":"./trees","tree_extension":".treefile"}'

treeness_over_rcv(toverr / treeness/RCV 比)には対応する function がありません。直接計算するには、下記の同梱 scogs_*.py スクリプトを使用してください。

生物学的な scogs 系統樹の妥当性確認用の目安:飽和度の中央値は約0.4~0.7、treeness/RCV の中央値は約0.2~0.4、treeness の中央値は約0.05~0.15です。これらから1桁ずれた値は、列の選択が間違っていることを意味します。

同梱スクリプト:BUSCO target_orthologs の共通部分

データフォルダーに *.busco.zip ファイルと target_orthologs.txt がある場合は、同梱スクリプトを使用してください。すべての zip ファイルにわたって single_copy_busco_sequences/*.faa を手動で列挙してはいけません。

python skills/tooluniverse-phylogenetics/scripts/busco_target_orthologs.py \
  --data-folder /path/to/data

デフォルトの実行では、「総アミノ酸数」の一般的な解釈をすべて網羅した、5行の要約が出力されます。

# SUMMARY: n_targets=K, n_intersected=N (single-copy in ALL S species), intersected_total_aa=A, sum_all_aa=B
# SUMMARY group=all: intersected n=N total_aa=A, sum_all total_aa=B
# SUMMARY group=animals: sum_all total_aa=X        <-- per-group sum (animal species only)
# SUMMARY group=fungi:   sum_all total_aa=Y        <-- per-group sum (fungal species only)
適切な SUMMARY 行の選び方(よく読んでください)

質問の表現に SUMMARY 行を対応させます。

質問の表現選ぶ行理由
「すべての単一コピーオルソログ配列に含まれる合計 AA」で、データフォルダーに動物種しかない場合、または質問が単一の生物グループに言及している場合# SUMMARY group=animals: sum_all total_aa=...(または group=fungi)scogs の系統ゲノミクス解析はグループごとに実行されます。「すべて」はグループ間の和集合ではなく、そのグループ内のすべてのオルソログを指します
「すべての種で単一コピーであるオルソログ全体の合計 AA」intersected_total_aa厳密な共通部分ルールです
「種ごとの全コピーを通じた合計 AA」sum_all_aa(group=all)質問に「すべての種」と明記されている場合、またはデータフォルダーに生物グループが 1 つしかない場合に限ります

データフォルダーに動物と菌類の両方の busco zip が含まれる場合のデフォルトルール:「総アミノ酸数」に関する公表済みの回答は、ほとんどの場合、グループ間の和集合ではなく、1 つのグループ(解析グループ)を指します。group=animals: sum_all または group=fungi: sum_all を使用してください。質問に「8 種すべてで」または「菌類と動物を合わせて」と明記されていない限り、和集合の値(sum_all_aa)を選ばないでください。

このため、スクリプトは標準出力にグループごとの合計を和集合の合計より先に出力します。出力を 1 行ずつ確認し、質問が示す解析グループに合致する group=animals / group=fungi の行で止めてください。

種間の単一コピーオルソログ — 比較セットと共通部分

BUSCO の single_copy_busco_sequences/ データを数える場合は、次の 2 段階のルールに従います。

  1. まず比較セットを特定します。 データフォルダーに target_orthologs.txt(または類似の名前のサブセットリスト)がある場合、そのファイルが比較セットです。そのオルソログ ID のみを対象にします。種間の BUSCO 単一コピー配列ファイルをすべて列挙しないでください。対象リストがある場合、「すべて」が BUSCO 出力全体を意味すると決めつけないでください。

  2. 次に共通部分ルールを適用します。 種間で「単一コピーオルソログ」とは、比較セット内のすべての種で単一コピーであることを意味します。ある種の single_copy_busco_sequences/ にオルソログがない場合は、他の種には存在していても、カウントから完全に除外してください。

妥当性チェック:ある種のオルソログごとのカウントが、他の種より大幅に少ない場合(例:同じオルソログセットで、ある種は約 600 aa、他の種は 4000 aa 以上)、一部の種に存在しないオルソログがオルソログごとの平均を押し上げています。まずそれらを除外してください。

具体例。 データフォルダーに 8 種(動物 4 種、菌類 4 種)の *.busco.zip と、10 個のオルソログ ID を記載した target_orthologs.txt がある場合:

  • 誤り:8 種すべての single_copy_busco_sequences/*.faa を列挙する → 約 80 ファイル → AA を合計する → 32228 と回答する(種ごとのコピーをそれぞれ独立して扱っています)。
  • 正しい方法:対象となる 10 個の ID ごとに、8 種すべてで single_copy として出現するか確認する → 共通部分に含まれる ID のみを残す(多くの場合 5/10 — 対象 ID の一部は、ある種で複数コピーまたは欠落しています)→ 残った ID について、8 種全体の AA を合計する → 13809。
  • 「5 trees」の意味:質問に「5 trees」とあるのに treefile が 10 個見つかった場合、GT ではすべての種で単一コピーであるオルソログの共通部分が使われています。10 個すべてではありません。平均を計算する前に、共通部分を再導出してください。
サンプルではなく全データを処理する(バッチ指標)

質問がオルソログ全体の中央値/パーセンタイル/平均値を求めている場合、バッチには該当する比較セット内のオルソログをすべて含める必要があります。

  • scogs_fungi.zip には菌類のアラインメントとツリーが約 255 個、scogs_animals.zip には約 241 個含まれています。10 ファイルのサンプルから計算した中央値は、公表済みの回答ではありません。
  • phykit_batch_analysis では、必ずすべてのオルソログのファイルが含まれる展開済み scogs ディレクトリを指定してください。手動で選んだサブセットは指定しないでください。
  • 計算した RCV/treeness/DVMC の中央値が妥当性確認用の目標値から 10% を超えて乖離する場合は、まずファイル数を確認してください。サブセットだけを処理している可能性があります。
絞り込んでから計算する(計算してから絞り込まない)

「ギャップが70%を超える遺伝子のうち、Xの最大値」のような質問では、最大値を求める前に絞り込みを適用する必要があります。

# 1. Compute gap% per alignment
# 2. Keep only alignments with gap% > 70
# 3. Compute treeness/RCV ON THE FILTERED SET
# 4. Take max

すべての遺伝子について指標を計算してから最大値を取ると、全体の最大値が返されてしまい、誤りです。

動物と菌類 — 長枝スコアの集約

PhyKIT の long_branch_score -v は、分類群ごとの LB スコアを出力します(系統樹の各葉につき1行)。系統樹ごとの要約を求める場合:

  1. 系統樹ごと:phykit long_branch_score -v <tree> を実行し、分類群ごとのスコア一覧を取得します。
  2. 系統樹ごとの要約:分類群ごとのスコアの平均または中央値を使って、系統樹ごとに1つの数値にまとめます。
  3. グループごとの要約:系統樹ごとの数値を集約します(中央値/平均/MWU U 値+p 値)。

質問の表現に合わせて、系統樹ごとの要約方法を選んでください:

質問の表現…--per-tree-stat ... に指定する値
「平均長枝スコア」mean
「長枝スコアの中央値」median
「平均長枝スコア」(曖昧)両方を実行し、数値/単位に一致する方を選択します

同梱の scogs_paired_compare.py --metric long_branch_score --per-tree-stat {mean,median} は、両グループについて手順1と2をまとめて実行し、グループ間の MWU U 値+p 値を直接出力します。

よくある誤り:系統樹ごとの手順を経ずに、4種の動物と4種の菌類を直接平均すること。これでは種ごとの LB とオルソログごとの LB が混同され、差分が誤った値になります。

樹木性/RCV:適切な入力ファイルを使用する

phykit toverr(別名 treeness_over_rcv)には、アラインメントと系統樹の両方が必要です。生の .faa.mafft ではなく、トリミング済みアラインメント(*.faa.mafft.clipkit)と対応するツリーファイル(*.faa.mafft.clipkit.treefile)を使用してください。ClipKit でトリミングしたアラインメントから標準的な系統樹が作成されているため、比率を公表済みの数値に一致させるには、同じトリミング済みアラインメントで RCV を計算する必要があります。

節約法で情報を持つサイト
  • 数える前にギャップのみの列を除外してください。すべてがギャップの列は情報を持ちません。
  • 2種類以上の異なる非ギャップ文字が、それぞれ2つ以上の分類群に出現するサイトは、節約法で情報を持ちます。
  • Biopython の AlignIO または AMAS ツールを使って列を走査し、数を数えてください。
樹木性(RCV 比率)

樹木性 = 内部枝長の合計 / 系統樹全体の枝長。内部枝とは、葉(末端)に続かない枝です。

PhyKIT の使用方法

PhyKIT(pip install phykit)には、系統樹やアラインメントの統計量を計算するコマンドライン関数が用意されています。よく使われる関数:

  • phykit treeness <tree_file> — 樹木性(RCV)値を出力します
  • phykit saturation <alignment_file> -t <tree_file> — 飽和値を出力します
  • phykit dvmc <tree_file> — 分子時計からの逸脱度を出力します
  • phykit long_branch_score <tree_file> — 長枝スコア(LBS)を出力します
  • phykit alignment_length <alignment_file> — アラインメント長を出力します
  • phykit parsimony_informative <alignment_file> — 節約法で情報を持つサイトの数を出力します

複数の遺伝子系統樹/アラインメントに対して PhyKIT を実行する場合は、同梱のバッチスクリプトを使用してください:

# Treeness across all trees
python skills/tooluniverse-phylogenetics/scripts/phykit_batch.py \
  --dir scogs_fungi --function treeness --ext .treefile --stat median

# Saturation with paired alignment+tree
python skills/tooluniverse-phylogenetics/scripts/phykit_batch.py \
  --dir alignments --function saturation --tree-dir trees \
  --ext .fa --tree-ext .treefile --stat median

# Long branch score (mean per tree, then median across trees)
python skills/tooluniverse-phylogenetics/scripts/phykit_batch.py \
  --dir trees --function long_branch_score --ext .treefile \
  --per-tree-stat mean --stat median

# DVMC
python skills/tooluniverse-phylogenetics/scripts/phykit_batch.py \
  --dir trees --function dvmc --ext .treefile --stat all

# Gap percentage across all alignments
python skills/tooluniverse-phylogenetics/scripts/phykit_batch.py \
  --dir alignments --function gap_percentage --ext .fa

# Evolutionary rate (median across trees)
python skills/tooluniverse-phylogenetics/scripts/phykit_batch.py \
  --dir trees --function evolutionary_rate --ext .treefile --stat median

# Mean patristic distance per tree, then mean across trees
python skills/tooluniverse-phylogenetics/scripts/phykit_batch.py \
  --dir trees --function patristic_distances --ext .treefile --stat mean

推奨: PhyKIT を手動で実行する代わりに、phykit_batch_analysis ToolUniverse ツールを使用してください:

# Via CLI
tu run phykit_batch_analysis '{"operation":"batch","function":"treeness","directory":"/path/to/trees","extension":".treefile"}'

# Via SDK
tu.run_one_function({"name": "phykit_batch_analysis", "arguments": {"operation": "batch", "function": "saturation", "directory": "/path/to/alignments", "extension": ".fa", "tree_directory": "/path/to/trees"}})

# Gap percentage
tu run phykit_batch_analysis '{"operation":"gap_percentage","directory":"/path/to/alignments","extension":".fa"}'

主なルール:

  1. すべてのファイルを処理する — 一部だけで止めないでください。ツールが自動的にすべて処理します
  2. ギャップ率:すべてのアラインメントにおけるギャップ数の合計 / 全位置数の合計(ファイルごとの平均ではありません)
  3. 長枝スコア:各系統樹から分類群ごとのスコアを算出 → 系統樹ごとに要約(平均)→ 系統樹全体で要約(中央値)。"per_tree_stat":"mean" を使用してください
  4. 菌類と動物の比較:ファイルの順番ではなく、オルソログ ID(ファイル名のステム)で遺伝子を対応付けます。それぞれの生物のディレクトリでツールを実行してから、中央値を比較してください

参考資料

references/sequence_alignment.md、references/tree_building.md、references/parsimony_analysis.md、scripts/tree_statistics.py

コメント