基因家族鉴定与进化¶

一句话概述：基因家族是由一个祖先基因经复制和分化产生的一组序列相似、功能相关的基因，鉴定基因家族并分析其进化是比较基因组学的核心任务。

核心知识点速查表¶

概念	说明
基因家族	来源于共同祖先的一组同源基因（白话：同一个老祖宗分化出来的基因群）
直系同源(Ortholog)	物种分化产生的同源基因（功能通常保守）
旁系同源(Paralog)	基因复制产生的同源基因（功能可能分化）
OrthoFinder	基因家族鉴定的主流工具
CAFE	基因家族扩张/收缩分析工具
MCL聚类	Markov Cluster Algorithm，蛋白质序列聚类算法
Ka/Ks	非同义/同义替换比，判断选择压力
系统发育	用进化树展示基因/物种间的亲缘关系

一、基因家族鉴定流程¶

1.1 OrthoFinder鉴定直系同源基因组¶

# === 安装OrthoFinder ===
conda install -c bioconda orthofinder    # conda安装

# === 准备输入文件 ===
# 每个物种一个蛋白质FASTA文件，放在同一目录
mkdir proteomes/                          # 创建输入目录
cp species1.pep.fa proteomes/             # 复制物种1蛋白质
cp species2.pep.fa proteomes/             # 复制物种2蛋白质
cp species3.pep.fa proteomes/             # 复制物种3蛋白质

# === 运行OrthoFinder ===
orthofinder \
    -f proteomes/ \        # 输入目录（包含所有物种的蛋白质FASTA）
    -t 16 \                # 比对用的线程数
    -a 8 \                 # 分析用的线程数
    -M msa \               # 用多序列比对建树（更准确）
    -S diamond              # 用DIAMOND加速序列搜索（比BLAST快500倍）

# === OrthoFinder 结果文件 ===
# Results_*/
# ├── Orthogroups/
# │   ├── Orthogroups.tsv            # ★每个基因家族包含的基因
# │   ├── Orthogroups_UnassignedGenes.tsv  # 未分组的孤儿基因
# │   └── Orthogroups.GeneCount.tsv  # ★每个物种每个家族的基因数
# ├── Orthologues/                    # 物种对之间的直系同源关系
# ├── Gene_Trees/                     # 每个基因家族的进化树
# ├── Species_Tree/                   # 物种树
# ├── Comparative_Genomics_Statistics/
# │   └── Statistics_Overall.tsv     # ★总体统计
# └── Phylogenetic_Hierarchical_Orthogroups/
#     └── N0.tsv                     # 层次化直系同源组

1.2 特定基因家族的HMM搜索¶

# === 使用HMMER搜索特定基因家族 ===
# 场景：已知某基因家族的HMM模型，想在新基因组中找该家族成员

# 第1步：下载Pfam HMM模型
wget https://www.ebi.ac.uk/interpro/wwwapi/entry/pfam/PF00069?annotation=hmm  # 激酶域
# 或从Pfam数据库下载对应家族的HMM文件

# 第2步：用hmmsearch搜索
hmmsearch \
    --tblout kinase_hits.tbl \     # 表格输出
    --domtblout kinase_domains.tbl \  # 结构域表格输出
    -E 1e-5 \                      # E-value阈值
    --cpu 8 \                       # CPU线程数
    PF00069.hmm \                   # HMM模型文件
    species.pep.fa                   # 蛋白质序列文件

# 第3步：提取候选基因
awk '$1!~/^#/ {print $1}' kinase_hits.tbl | \   # 提取基因ID
    sort -u > kinase_genes.list                   # 去重

二、基因家族扩张/收缩分析（CAFE）¶

# === 安装CAFE5 ===
conda install -c bioconda cafe    # conda安装

# === 准备输入文件 ===
# 1. 基因家族计数文件（从OrthoFinder结果转换）
# 格式：Family_ID\tspecies1\tspecies2\tspecies3...
# 注意：过滤掉只在1个物种中出现的家族和超大家族（>100基因）

# 从OrthoFinder结果生成CAFE输入
python3 << 'PYEOF'
import pandas as pd

# 读取OrthoFinder基因计数
df = pd.read_csv("Orthogroups.GeneCount.tsv", sep="\t")  # 读取计数文件
df = df.drop(columns=["Total"])                            # 删除Total列

# 过滤：至少2个物种有基因，且最大基因数<100
mask = (df.iloc[:, 1:] > 0).sum(axis=1) >= 2              # 至少2个物种有
mask &= df.iloc[:, 1:].max(axis=1) < 100                  # 最大数<100
df_filtered = df[mask]

# 添加CAFE格式的描述列
df_filtered.insert(0, "Desc", "(null)")                    # 添加描述列

# 保存
df_filtered.to_csv("cafe_input.tsv", sep="\t", index=False)  # 输出
print(f"保留 {len(df_filtered)} 个基因家族")
PYEOF

# === 2. 超度量物种树（带分歧时间） ===
# 使用r8s或MCMCtree从物种树推断分歧时间
# 树的格式：((sp1:10,sp2:10):5,sp3:15);  （数字是百万年）

# === 运行CAFE5 ===
cafe5 \
    -i cafe_input.tsv \           # 基因家族计数文件
    -t species_tree.nwk \         # 超度量物种树（带分歧时间）
    -p \                          # 计算p-value
    -k 2 \                        # 假设2个lambda值（进化速率）
    -o cafe_output                 # 输出目录

# === 结果解读 ===
# cafe_output/
# ├── Base_asr.tre          # 祖先基因家族大小重建
# ├── Base_family_results.txt  # ★每个家族的扩张/收缩和p-value
# ├── Base_branch_probabilities.txt  # 每个分支的变化概率
# └── Base_clade_results.txt  # 每个支系的显著变化家族

# === 提取显著扩张/收缩的基因家族 ===
# 找出p-value < 0.05的家族
awk '$2 < 0.05' cafe_output/Base_family_results.txt | \
    head -20   # 查看前20个显著变化的家族

三、系统发育分析¶

# === 单基因家族建树流程 ===

# 第1步：多序列比对
mafft \
    --auto \                    # 自动选择最佳策略
    --thread 8 \                # 线程数
    gene_family.fa > aligned.fa  # 输入→输出

# 第2步：修剪比对（去除低质量区域）
trimal \
    -in aligned.fa \            # 输入比对结果
    -out trimmed.fa \           # 输出修剪结果
    -automated1                  # 自动选择最佳修剪策略

# 第3步：建树（最大似然法）
iqtree2 \
    -s trimmed.fa \             # 输入修剪后的比对
    -m MFP \                    # 自动选择最佳模型(ModelFinder Plus)
    -bb 1000 \                  # 1000次ultrafast bootstrap
    -alrt 1000 \                # 1000次SH-aLRT检验
    -nt AUTO                     # 自动选择线程数

# 产出：trimmed.fa.treefile（Newick格式进化树）

# === Ka/Ks分析（选择压力检测） ===
# Ka/Ks < 1: 纯化选择（功能保守）
# Ka/Ks = 1: 中性进化
# Ka/Ks > 1: 正选择（功能分化）

# 使用KaKs_Calculator
KaKs_Calculator \
    -i gene_pairs.axt \         # 成对比对文件（AXT格式）
    -o kaks_results.txt \       # 输出结果
    -m YN                        # Yang-Nielsen方法

# 使用PAML的codeml进行位点模型检验（更严格）
# 需要准备ctl配置文件，设置model=0(M0)和model=2(M2a)

四、面试高频考点¶

Q1: Ortholog和Paralog的区别？¶

Ortholog：物种分化事件产生，A物种的基因X和B物种的基因X'。白话：不同物种里的"对应基因"
Paralog：基因复制事件产生，同一物种内的基因X和基因X'。白话：同一物种里的"兄弟基因"
面试关键：Ortholog通常功能保守，Paralog功能可能分化（新功能化/亚功能化）

Q2: 如何判断基因家族扩张有生物学意义？¶

CAFE统计检验p-value < 0.05
扩张基因的表达模式（是否在特定组织/条件下高表达）
扩张基因的功能富集（GO/KEGG分析）
正选择信号（Ka/Ks > 1的位点）
与表型/适应性的关联

常见报错与解决¶

报错	原因	解决方案
OrthoFinder内存不足	物种太多或蛋白太多	用`-S diamond`加速，减少物种数
CAFE报错`unreliable`	输入树不是超度量树	用r8s/MCMCtree校正分歧时间
IQ-TREE `alignment too short`	序列太短	降低bootstrap次数或换NJ方法
HMMER无结果	E-value阈值太严格	放宽到1e-3再手动过滤

速查表¶

# === 基因家族分析一键流程 ===
# 1. OrthoFinder鉴定基因家族
orthofinder -f proteomes/ -t 16 -M msa -S diamond

# 2. CAFE分析扩张/收缩
cafe5 -i cafe_input.tsv -t tree.nwk -p -o cafe_out

# 3. 特定家族建树
mafft --auto family.fa > aln.fa
trimal -in aln.fa -out trim.fa -automated1
iqtree2 -s trim.fa -m MFP -bb 1000

# 4. 选择压力分析
KaKs_Calculator -i pairs.axt -o kaks.txt -m YN