有效种群大小（Ne）估计方法¶

一句话概述：有效种群大小 Ne 是群体遗传学的核心参数，反映了群体中实际参与繁殖的"等效个体数"，可通过连锁不平衡（LD）、GONE 和 PSMC 等方法从基因组数据中估计。

核心知识点速查表¶

概念	白话解释
Ne（有效种群大小）	一个理想群体要产生相同遗传漂变需要的个体数
Census size (N)	实际可数的种群个体数
Ne 通常远小于 N	因为不是所有个体都均等繁殖
LD 法	利用连锁不平衡估计近期 Ne（几代内）
GONE	利用全基因组 LD 估计历史 Ne（200代内）
PSMC	利用单个基因组的杂合度估计远古 Ne（千~百万年）
NeEstimator	基于 LD 的经典 Ne 估计软件
currentNe2	2025 年新工具，不需要遗传图谱

一、Ne 的概念（白话版）¶

比喻：一个村庄有 1000 人（N=1000），但实际只有 200 对夫妇生了后代。而且其中 10 对特别能生，贡献了大部分后代。从遗传学角度看，这个村庄"有效"只相当于一个 150 人的村庄（Ne=150）。

为什么 Ne 重要？ - Ne 决定了遗传漂变的强度（Ne 越小，漂变越强，等位基因越容易随机丢失） - Ne 决定了选择的效率（Ne 越小，弱选择效果越弱） - 保护生物学中，Ne > 50 避免近交衰退，Ne > 500 维持进化潜力 - 面试常考：Ne/N 比值通常是 0.1~0.5

二、估计方法概览¶

Ne 估计方法分类¶

方法	时间尺度	输入数据	代表工具
LD 法	近期（1~几代）	群体SNP数据	NeEstimator, currentNe2
GONE	历史（~200代）	群体SNP+遗传图谱	GONE, GONE2
PSMC	远古（千~百万年）	单个高深度基因组	PSMC, SMC++
时序法	两个时间点间	两次采样的等位基因频率	NeEstimator
杂合度法	近期	群体基因型数据	currentNe

三、方法 1：NeEstimator（LD 法估计当代 Ne）¶

安装与运行¶

# 下载 NeEstimator V2.1
# 官网：http://www.molecularfisherieslaboratory.com.au/neestimator-software/
wget http://www.molecularfisherieslaboratory.com.au/download/NeEstimator-v2.1.zip
unzip NeEstimator-v2.1.zip  # 解压

# NeEstimator 有 GUI 和命令行两种模式
# 命令行模式需要准备参数文件

# 准备输入数据（GENEPOP格式）
# GENEPOP格式示例：
cat << 'EOF' > input.gen
Title: My Population
Locus1
Locus2
Locus3
Pop
Ind1, 0101 0203 0102
Ind2, 0102 0201 0101
Ind3, 0201 0203 0202
Pop
Ind4, 0103 0201 0101
Ind5, 0102 0103 0201
EOF

# 从VCF转换为GENEPOP格式
# 用 PGDSpider 或 radiator R包转换

# 准备 NeEstimator 参数文件
cat << 'EOF' > ne_option.txt
1                        # 输入格式：1=GENEPOP
input.gen                # 输入文件名
ne_output.txt            # 输出文件名
3                        # 方法：3=LD法
0.02                     # 最低等位基因频率阈值（Pcrit）
0                        # 0=不做时序估计
EOF

# 运行 NeEstimator
./Ne2-1 i:ne_option.txt  # 运行（i:表示参数文件输入）

用 PLINK 准备 LD 数据¶

# 如果数据是VCF格式，可以先用PLINK计算LD，再手动估计Ne
# 但建议直接用 NeEstimator 或 currentNe2

# PLINK 计算 r² 用于 LD-based Ne estimation
plink --vcf input.vcf.gz \    # 输入VCF
      --r2 \                    # 计算r²
      --ld-window-r2 0 \       # 输出所有r²值（包括0）
      --ld-window 999999 \     # 最大SNP间隔（对数）
      --ld-window-kb 10000 \   # 最大物理距离10Mb
      --out ld_results         # 输出前缀

# Ne ≈ 1/(3 × (r²_adjusted - 1/S))
# r²_adjusted = r² - 1/(2S)，S=样本量
# 这只是粗略估计，建议用专业软件

四、方法 2：GONE（历史 Ne 变化轨迹）¶

# GONE 利用全基因组LD模式估计最近200代的Ne变化
# 优点：能看到Ne随时间的变化趋势（如瓶颈、扩张）

# 安装 GONE
git clone https://github.com/esrud/GONE.git  # 克隆仓库
cd GONE
chmod +x GONE                                  # 添加执行权限

# 准备输入文件（PLINK ped/map 格式）
plink --vcf input.vcf.gz \      # 输入VCF
      --recode \                 # 转换为ped/map格式
      --allow-extra-chr \        # 允许非标准染色体名
      --out mydata               # 输出前缀

# 准备遗传图谱文件（按染色体的遗传距离）
# 如果没有遗传图谱，可以用物理距离近似（1 Mb ≈ 1 cM）

# 运行 GONE
bash script_GONE.sh mydata       # 运行GONE

# 输出文件：
# Output_Ne_mydata — 每一代的 Ne 估计值
# 格式：代数  Ne估计值

# 可视化 Ne 变化曲线
cat << 'EOF' > plot_gone.R
ne <- read.table("Output_Ne_mydata", header=FALSE)  # 读取结果
colnames(ne) <- c("Generation", "Ne")                # 命名列

pdf("Ne_trajectory_GONE.pdf", width=10, height=6)    # 创建PDF
plot(ne$Generation, ne$Ne,                            # 画Ne变化曲线
     type="l", col="blue", lwd=2,                    # 蓝色线
     xlab="Generations ago",                          # X轴
     ylab="Effective population size (Ne)",           # Y轴
     main="Historical Ne estimated by GONE",          # 标题
     log="y")                                         # Y轴取对数
dev.off()
EOF
Rscript plot_gone.R

GONE2（2025 年更新）¶

# GONE2 是 2025 年发布的升级版
# 新增功能：
# - 可以处理群体结构（FST、迁移率、亚群数量估计）
# - 支持单倍体数据
# - 处理基因分型错误和低深度测序数据

# 安装 GONE2
git clone https://github.com/esrud/GONE2.git  # 克隆GONE2
cd GONE2
chmod +x GONE2

# 运行（用法类似 GONE）
bash script_GONE2.sh mydata

五、方法 3：PSMC（远古 Ne 变化）¶

# PSMC 用单个二倍体基因组的杂合度分布估计 Ne 随时间变化
# 适用范围：约 1万 ~ 100万年前

# 安装 PSMC
git clone https://github.com/lh3/psmc.git  # 克隆仓库
cd psmc
make                                         # 编译
cd utils
make                                         # 编译工具

# 第1步：从BAM生成一致性序列（consensus sequence）
# 需要高深度（>20X）的全基因组测序数据
samtools mpileup -C 50 -uf ref.fa sample.bam | \  # 生成pileup
    bcftools call -c | \                            # 变异检测
    vcfutils.pl vcf2fq -d 10 -D 100 | \           # 转换为fastq，深度10-100
    gzip > consensus.fq.gz                          # 压缩

# 第2步：转换格式
psmc/utils/fq2psmcfa -q20 consensus.fq.gz > input.psmcfa  # 转换为psmc输入格式

# 第3步：运行 PSMC
psmc -N25 \                      # 最大迭代次数25
     -t15 \                       # 最大时间参数
     -r5 \                        # 初始θ/ρ比值
     -p "4+25*2+4+6" \           # 时间区间模式
     -o output.psmc \             # 输出文件
     input.psmcfa                 # 输入文件

# 第4步：Bootstrap（评估置信区间）
# 先分割输入文件
psmc/utils/splitfa input.psmcfa > split.psmcfa

# 运行100次bootstrap
seq 100 | xargs -i -P 4 \        # 并行4个任务
    psmc -N25 -t15 -r5 -b \       # -b 表示bootstrap模式
         -p "4+25*2+4+6" \
         -o round-{}.psmc \
         split.psmcfa

# 第5步：合并结果并画图
cat output.psmc round-*.psmc > combined.psmc  # 合并所有结果

psmc/utils/psmc_plot.pl \
    -g 5 \                        # 世代时间（年/代）
    -u 2.5e-8 \                   # 突变率（每位点每代）
    -p plot_output \              # 输出前缀
    combined.psmc                 # 输入合并后的文件

# 输出 plot_output.pdf — Ne 随时间变化的曲线
# X轴：年（对数刻度），Y轴：Ne

六、方法比较与选择指南¶

特征	NeEstimator	GONE/GONE2	PSMC
时间尺度	当代（几代）	历史（~200代）	远古（千~百万年）
输入数据	群体 SNP	群体 SNP + 图谱	单个高深度基因组
最少样本量	≥25 个体	≥50 个体	1 个体
需要遗传图谱	不需要	需要（GONE2可选）	不需要
群体结构影响	敏感	较敏感	不敏感
保护遗传学适用	最佳	适用	不适用于当代 Ne

七、常见报错与解决¶

报错信息	原因	解决方法
`Ne = Infinity`	样本量太小或SNP太少	增加样本量（≥25）或 SNP 数量
`Negative r² values`	数据质量差或样本太少	过滤低质量SNP，增加样本
GONE 结果 Ne 值异常大	存在群体结构或基因流	用 GONE2 处理结构
PSMC 曲线不平滑	基因组覆盖度不够	需要 >20X 深度
`Too few chromosomes`	GONE 需要多条染色体	确保有多条染色体的数据
NeEstimator 和 GONE 结果差很大	估计的时间尺度不同	正常，两者回答不同的问题

八、面试高频问题¶

Q1：Ne 和 N（普查种群大小）的关系？

Ne 通常远小于 N。比值 Ne/N 约 0.1~0.5。原因：不均等的性比（一个雄性配多个雌性）、世代重叠、种群大小波动（瓶颈效应使 Ne 接近最小值的调和平均数）。

Q2：Ne=50 和 Ne=500 的意义？

50/500 法则（Frankham et al.）：Ne>50 可避免短期近交衰退；Ne>500 可维持长期进化适应能力。这是保护生物学中的重要阈值。

Q3：PSMC 和 GONE 能互相替代吗？

不能。PSMC 看远古（万~百万年），GONE 看近期（几百代）。两者时间尺度不重叠或只有少量重叠，应该联合使用以获得完整的历史 Ne 变化图景。

九、速查表¶

# === Ne 估计速查 ===

# 1. NeEstimator（当代Ne）
./Ne2-1 i:ne_option.txt

# 2. GONE（历史Ne，~200代）
plink --vcf input.vcf.gz --recode --out mydata
bash script_GONE.sh mydata

# 3. GONE2（2025新版，处理结构）
bash script_GONE2.sh mydata

# 4. PSMC（远古Ne，单基因组）
samtools mpileup -C 50 -uf ref.fa sample.bam | bcftools call -c | vcfutils.pl vcf2fq -d 10 -D 100 | gzip > consensus.fq.gz
fq2psmcfa -q20 consensus.fq.gz > input.psmcfa
psmc -N25 -t15 -r5 -p "4+25*2+4+6" -o output.psmc input.psmcfa
psmc_plot.pl -g 5 -u 2.5e-8 -p plot output.psmc

# Ne 重要阈值：
# Ne > 50  → 避免近交衰退
# Ne > 500 → 维持进化潜力
# Ne/N ≈ 0.1~0.5（经验值）

参考资料：GONE2/currentNe2 (Santiago et al. 2025, Nature Communications)、NeEstimator V2 (Do et al. 2014)、PSMC (Li & Durbin 2011)、Waples 2024 "The Idiot's Guide to Ne"