ksrates介绍
Ksrates 是一个基于同义替换率（Ks）的系统发育分析工具，用于校正不同演化谱系之间的替换速率差异（rate heterogeneity）。它通过整合焦点物种的旁系同源基因（paralogs）和不同物种间的直系同源基因（orthologs），在系统发育框架下推断远古全基因组复制（Whole Genome Duplication，WGD）事件与物种分化事件之间的时间关系。
与传统直接比较Ks峰值的方法不同，ksrates利用多个外类群估计不同分支的演化速率，并将物种分化Ks峰校正到焦点物种的Ks尺度，从而避免由于不同谱系演化速率不同而造成WGD时间判断偏差。
GitHub：https://github.com/VIB-PSB/ksrates
使用ksrates的分析流程
Step1 数据准备
准备输入文件：
- CDS序列（多个物种）：需要统一CDS文件命名，物种简称.cds.fa
- 焦点物种GFF
- 系统发育树（Newick）：纯拓扑关系版本的进化树，移除所有的冒号、数字和支长，只保留括号和物种简称，纯粹用来表示物种间的亲缘远近顺序。
- 物种名称映射：物种简称对应的拉丁学名
可以直接使用Github上官方的配置文件脚本进行修改，也可以根据官方配置脚本的内容自己写配置文件并在主流程脚本中生成ksrates可以读取的配置文件（ksrates.conf）
Step2 焦点物种旁系Ks计算
运行
subprocess.run(
["ksrates", "paralogs-ks", abs_conf_path, "--n-threads", threads],
check=True
)得到焦点物种所有旁系同源基因Ks分布：Tgra.ks.tsv
Step3 焦点物种与所有外类群计算Ortholog Ks
cmd = [
"ksrates", "orthologs-ks",
abs_conf_path,
sp1, sp2,
"--n-threads", threads
]Step4 初始化速率校正
运行ksrates init
subprocess.run(["ksrates", "init", abs_conf_path], check=True)解析 init 生成的所需物种对列表ortholog_pairs_Tgra.tsv，补充再跑一遍ortholog Ks
什么还要再跑一遍ortholog Ks？
实际上，
init会根据系统发育树自动确定进行速率校正所需的三元组（Triplets），并列出仍需计算的物种配对。
例如
A-B
A-C
B-C我们需要补齐缺失配对
# 解析 init 生成的所需物种对列表
required_pairs = set()
with open(pair_file, 'r') as f:
for i, line in enumerate(f):
line = line.strip()
if i == 0 or not line:
continue
parts = line.split('\t')
if len(parts) >= 2:
sp1, sp2 = sorted((parts[0], parts[1]))
required_pairs.add((sp1, sp2))
required_pairs = list(required_pairs)
# 再跑一遍Ortholog Ks
print("\n" + "=" * 70)
print(f"Running MISSING pairwise ortholog Ks for trios")
print("=" * 70)
computed = set()
for sp1, sp2 in required_pairs:
sp1, sp2 = sorted((sp1, sp2))
if (sp1, sp2) in computed:
continue
computed.add((sp1, sp2))
pair_name = f"{sp1}_vs_{sp2}"
anchor1 = os.path.join(target_dir, "ortholog_distributions", f"wgd_{sp1}_{sp2}", f"{sp1}_{sp2}.ks.tsv")
anchor2 = os.path.join(target_dir, "ortholog_distributions", f"wgd_{sp2}_{sp1}", f"{sp2}_{sp1}.ks.tsv")
if (os.path.exists(anchor1) and os.path.getsize(anchor1) > 0) or \
(os.path.exists(anchor2) and os.path.getsize(anchor2) > 0):
continue
print(f"\n" + "-" * 60)
print(f"Running ortholog Ks : {sp1} vs {sp2}")
print("-" * 60)
#继续运行ksrates orthologs-ks，直到ortholog_pairs.tsv所有配对全部完成
cmd = ["ksrates", "orthologs-ks", abs_conf_path, sp1, sp2, "--n-threads", threads]
try:
subprocess.run(cmd, check=True)
except subprocess.CalledProcessError as e:
print(f"[FAILED] {pair_name}\n{e}")
failed_pairs.append(pair_name)Step5 Ks速率校正
Peak database、Ks database：注意配置数据库的路径与结果文件路径一致
peak_db = config["ksrates_databases"]["peak_database_path"]
ks_db = config["ksrates_databases"]["ks_list_database_path"]
for f in [peak_db, ks_db]:
if os.path.exists(f) and os.path.getsize(f) > 0:
print(f"[SUCCESS] {f}")
else:
print(f"[WARNING] Missing {f}")
#运行paralogs-analyses，拟合WGD峰
print("\n[INFO] Running paralog analyses")
subprocess.run(["ksrates", "paralogs-analyses", abs_conf_path], check=True)
#可视化分析，得到Ks histogram、Mixed Ks plot、Corrected Ks plot、Ks tree
print("\n[INFO] Plot mixed Ks")
subprocess.run(["ksrates", "plot-paralogs", abs_conf_path], check=True)
print("\n[INFO] Plot ortholog distributions")
subprocess.run(["ksrates", "plot-orthologs", abs_conf_path], check=True)
print("\n[INFO] Plot Ks tree")
try:
subprocess.run(["ksrates", "plot-tree", abs_conf_path], check=True)
except subprocess.CalledProcessError:
print("[WARNING] plot-tree skipped.")
if failed_pairs:
print("\nFAILED PAIRS:")
for sp in failed_pairs:
print(sp)
print("\nAll analyses finished.")