溯祖模型模拟 — msprime¶

一句话概述：msprime 是基于溯祖理论的高效群体遗传模拟器，能在秒级模拟百万样本的全基因组数据，广泛用于方法验证、统计检验和群体历史推断。

核心知识点速查表¶

概念	白话解释
溯祖（Coalescent）	从现在的样本"回溯"到共同祖先的过程
树序列（Tree Sequence）	msprime 的核心数据结构，存储谱系树
tskit	处理树序列的 Python 库
Ne	有效种群大小（决定溯祖速率）
重组率	相邻位点之间发生交换的概率
突变率	每位点每代产生突变的概率
人口学模型	描述群体大小变化、分裂、迁移的模型
SLiM	前向模拟器（正向时间），与 msprime 互补

一、溯祖理论（白话版）¶

比喻：你和你的表兄弟都来自同一个曾祖父。如果从你们"回溯"家谱，最终会在曾祖父那里"汇合"——这就是溯祖（coalescence）。

为什么用溯祖模拟？ - 前向模拟（从祖先到后代）：需要模拟整个群体所有个体所有基因组 → 超级慢 - 溯祖模拟（从后代回溯祖先）：只模拟被采样的个体的祖先谱系 → 超级快 - msprime 可以在几秒内模拟千万级样本的基因组数据

二、安装与基本用法¶

# 安装 msprime 和 tskit
pip install msprime tskit  # pip安装
# 或
conda install -c conda-forge msprime tskit  # conda安装

最简单的模拟¶

# ========== Python脚本：基本溯祖模拟 ==========
import msprime                              # 导入msprime
import tskit                                 # 导入tskit

# 模拟10个样本、100kb基因组
ts = msprime.sim_ancestry(
    samples=10,                              # 10个二倍体个体
    sequence_length=1e5,                     # 基因组长度100kb
    recombination_rate=1e-8,                 # 重组率（人类约1e-8）
    population_size=10000,                   # 有效种群大小Ne=10000
    random_seed=42                           # 随机种子（可重复）
)

# 在树序列上叠加突变
mts = msprime.sim_mutations(
    ts,                                      # 输入树序列
    rate=1e-8,                               # 突变率（人类约1.4e-8）
    random_seed=42                           # 随机种子
)

# 查看基本信息
print(f"树的数量: {mts.num_trees}")          # 区间内有多少棵不同的谱系树
print(f"变异位点数: {mts.num_sites}")        # 突变位点数量
print(f"样本数: {mts.num_samples}")          # 样本数量

# 导出为VCF（可用于下游分析）
with open("simulated.vcf", "w") as f:        # 打开输出文件
    mts.write_vcf(f)                          # 写入VCF格式

模拟复杂人口学模型¶

# ========== 模拟群体分裂+瓶颈+扩张 ==========
import msprime

# 定义人口学模型
demography = msprime.Demography()            # 创建人口学对象

# 添加群体
demography.add_population(                   # 添加群体A
    name="PopA",
    initial_size=50000                       # 当前Ne=50000
)
demography.add_population(                   # 添加群体B
    name="PopB",
    initial_size=30000                       # 当前Ne=30000
)
demography.add_population(                   # 添加祖先群体
    name="Ancestor",
    initial_size=20000                       # 祖先Ne=20000
)

# 添加群体事件（时间单位是"代"，从现在往回数）
# 1000代前：PopA和PopB从Ancestor分裂
demography.add_population_split(
    time=1000,                               # 1000代前
    derived=["PopA", "PopB"],                # 分裂出的子群体
    ancestral="Ancestor"                     # 祖先群体
)

# 500代前：PopA经历瓶颈（Ne从50000降到500，持续50代）
demography.add_population_parameters_change(
    time=500,                                # 500代前
    population="PopA",
    initial_size=500                         # 瓶颈时Ne=500
)
demography.add_population_parameters_change(
    time=550,                                # 550代前（瓶颈结束）
    population="PopA",
    initial_size=10000                       # 瓶颈前Ne=10000
)

# 模拟
ts = msprime.sim_ancestry(
    samples={"PopA": 50, "PopB": 50},        # 每个群体采50个个体
    demography=demography,                   # 人口学模型
    sequence_length=1e7,                     # 10Mb基因组
    recombination_rate=1e-8,
    random_seed=42
)

# 添加突变
mts = msprime.sim_mutations(ts, rate=1e-8, random_seed=42)

# 导出VCF（分群体）
with open("pop_split_sim.vcf", "w") as f:
    mts.write_vcf(f)

三、计算群体遗传统计量¶

# ========== 用 tskit 计算常用统计量 ==========

# 获取群体的样本列表
popA_samples = mts.samples(population=0)     # PopA的样本ID
popB_samples = mts.samples(population=1)     # PopB的样本ID

# 1. 核苷酸多样性（π）
pi_A = mts.diversity(sample_sets=[popA_samples])   # PopA的π
pi_B = mts.diversity(sample_sets=[popB_samples])   # PopB的π
print(f"PopA π = {pi_A[0]:.6f}")
print(f"PopB π = {pi_B[0]:.6f}")

# 2. Fst
fst = mts.Fst(sample_sets=[popA_samples, popB_samples])  # 计算Fst
print(f"Fst = {fst[0]:.4f}")

# 3. Tajima's D
tajD_A = mts.Tajimas_D(sample_sets=[popA_samples])   # PopA的Tajima's D
print(f"PopA Tajima's D = {tajD_A[0]:.4f}")
# D < 0 → 群体扩张或正向选择
# D > 0 → 群体收缩或平衡选择
# D ≈ 0 → 中性进化

# 4. 沿基因组的窗口统计量
windows = list(range(0, int(1e7), 100000)) + [int(1e7)]  # 每100kb一个窗口
pi_windows = mts.diversity(
    sample_sets=[popA_samples],
    windows=windows                          # 指定窗口
)

# 5. 位点频率谱（SFS）
afs = mts.allele_frequency_spectrum(
    sample_sets=[popA_samples],              # 指定样本
    polarised=True                           # 极化（有祖先状态）
)
print("SFS:", afs)

四、模拟自然选择（结合 SLiM）¶

# msprime 本身只模拟中性进化
# 需要选择的场景可以用 msprime 初始化 + SLiM 前向模拟

# 简单方法：用 msprime 模拟大部分基因组（中性），
# 然后在感兴趣的区域用 SLiM 模拟选择

# 或者用 stdpopsim（标准化的群体遗传模拟）
# pip install stdpopsim
import stdpopsim                             # 导入stdpopsim

# 使用人类的标准模型
species = stdpopsim.get_species("HomSap")    # 人类
model = species.get_demographic_model("OutOfAfrica_3G09")  # OOA模型
contig = species.get_contig("chr22")         # 22号染色体
samples = {"YRI": 50, "CEU": 50, "CHB": 50}  # 三个群体各50人

engine = stdpopsim.get_engine("msprime")     # 使用msprime引擎
ts = engine.simulate(model, contig, samples)  # 模拟

print(f"模拟了 {ts.num_sites} 个变异位点")

五、常见报错与解决¶

报错信息	原因	解决方法
`ValueError: population not found`	群体名拼写错误	检查 `demography.add_population` 中的名称
`MemoryError`	基因组太大或样本太多	分染色体模拟，或减少样本量
`populations must be inactive`	人口学事件时间顺序错误	确保事件按时间从近到远排列
`recombination_rate must be >= 0`	重组率参数有误	检查参数值是否合理
VCF 中样本名都是 tsk_0、tsk_1	msprime 默认命名	用 `individual_names` 参数指定名称

六、面试高频问题¶

Q1：溯祖模拟和前向模拟的区别？

溯祖模拟从现在回溯到祖先（只跟踪被采样个体的谱系），非常快，适合中性进化。前向模拟从祖先到后代（模拟整个群体），很慢但可以模拟选择、频率依赖性适合度等复杂场景。msprime 属于溯祖，SLiM 属于前向。

Q2：msprime 在实际研究中怎么用？

主要用于：(1) 验证统计方法在已知模型下的表现；(2) 做零假设模拟（中性模型下的 Fst/π 分布）；(3) 估计统计功效（如检测选择的能力）；(4) ABC 推断（近似贝叶斯计算）。

Q3：树序列（tree sequence）有什么优势？

树序列是一种高效的数据结构，存储的是谱系树而不是基因型矩阵。它比 VCF 节省几个数量级的存储空间，并且可以直接计算各种统计量（π、Fst、SFS 等），速度比从 VCF 计算快得多。

七、速查表¶

# === msprime 速查 ===

# 基本模拟
ts = msprime.sim_ancestry(samples=100, sequence_length=1e6, 
                           recombination_rate=1e-8, population_size=10000)
mts = msprime.sim_mutations(ts, rate=1e-8)

# 导出VCF
with open("out.vcf", "w") as f: mts.write_vcf(f)

# 人口学模型
dem = msprime.Demography()
dem.add_population(name="A", initial_size=10000)
dem.add_population(name="B", initial_size=10000)
dem.add_population(name="AB", initial_size=20000)
dem.add_population_split(time=1000, derived=["A","B"], ancestral="AB")

# 常用统计量（tskit）
ts.diversity(sample_sets=[samples])        # π
ts.Fst(sample_sets=[pop1, pop2])           # Fst
ts.Tajimas_D(sample_sets=[samples])        # Tajima's D
ts.allele_frequency_spectrum()             # SFS

# 人类参数参考：
# Ne ≈ 10,000
# 突变率 μ ≈ 1.4e-8 /bp/代
# 重组率 r ≈ 1e-8 /bp/代
# 世代时间 ≈ 25-30 年

参考资料：msprime 1.0 (Baumdicker et al. 2022, Genetics)、tskit文档、stdpopsim (Adrion et al. 2020)、Johnson et al. 2024 (Mol Ecol Res)