从SSR到DeltaK：群体结构分析的完整流程与可视化实践

1. SSR数据准备与格式转换

第一次接触群体结构分析时，最让我头疼的就是数据格式问题。SSR（简单序列重复）作为共显性标记，其原始数据往往不符合Structure软件的要求格式。记得当时我花了整整三天时间才搞明白如何正确转换数据格式，现在我把这个踩坑经验分享给你。

Structure要求每个样本占两行：第一行包含样本名称、种群编号和等位基因1的数据，第二行则是相同的样本名称和种群编号，加上等位基因2的数据。原始SSR数据通常是一行一个样本，交替记录两个等位基因。我写了个Python脚本来自动完成这个转换：

#!/usr/bin/env python3 import sys dat = sys.argv[1] with open(dat,'r') as f: num=1 for LN in f: if num ==1: print(LN.strip()) num-=1 else: spl = LN.strip().split() cnt = len(spl) print(' '.join(spl[0:3]),' '.join([spl[i] for i in range(3,cnt,2)])) print(' '.join(spl[0:3]),' '.join([spl[i] for i in range(4,cnt,2)]))

使用这个脚本时，只需将原始数据文件作为参数传入：

python3 dat2structure.py raw_data.txt > structure_input.txt

转换后的数据格式应该是这样的：

样本1 种群1 等位基因1_A 等位基因1_B 等位基因1_C 样本1 种群1 等位基因2_A 等位基因2_B 等位基因2_C 样本2 种群2 等位基因1_A 等位基因1_B 等位基因1_C 样本2 种群2 等位基因2_A 等位基因2_B 等位基因2_C

2. Structure软件运行与参数设置

Structure是群体结构分析的核心工具，但它的参数设置对新手来说可能有点复杂。我建议从mainparams和extraparams两个配置文件开始。mainparams中这几个参数需要特别注意：

• numreps：建议设置为10000，确保结果稳定
• burnin：预热迭代次数，5000是个不错的起点
• nummcmc：正式运行的迭代次数，10000次通常足够
• inferalpha：是否推断α值，对于复杂群体建议设为1
• popinfo：是否使用先验种群信息，初学者建议设为0

运行Structure时，我习惯用shell脚本批量处理不同K值：

for K in {1..10};do (nohup structure -i structure_input.txt -m mainparams -e extraparams \ -K ${K} -o run1_k${K} &); done

这里有个实用技巧：每个K值至少运行3次，后续计算DeltaK时结果会更可靠。我曾经因为只运行一次导致结果不稳定，不得不全部重跑，浪费了一周时间。

3. DeltaK计算与最佳K值确定

Structure运行完成后，我们需要确定最佳的群体分组数K。Evanno方法通过计算DeltaK值来解决这个问题，这比直接看LnP(K)更可靠。我推荐使用structureHarvester工具来自动化这个过程。

安装structureHarvester很简单：

git clone https://github.com/dentearl/structureHarvester.git cd structureHarvester python setup.py install

使用命令行分析结果：

python structureHarvester.py --dir=./results/ --out=./summary --evanno

这个命令会生成一个包含DeltaK值的evanno.txt文件。接下来用R绘制DeltaK折线图：

library(ggplot2) data <- read.table("evanno.txt") ggplot(data, aes(x=V1, y=V7)) + geom_line(color="blue") + geom_point(color="red", size=3) + labs(x="K value", y="Delta K") + theme_minimal()

DeltaK峰值对应的K值就是最佳群体分组数。记得我第一个项目在这个环节犯过错——忽略了DeltaK而直接选择LnP(K)平台期的K值，导致后续分析全部需要重做。

4. 结果整合与可视化

得到最佳K值后，我们需要整合多次运行的重复结果。CLUMPP工具可以解决这个问题，它能对齐不同重复运行中的群体聚类。

首先准备CLUMPP的paramfile配置文件：

DATATYPE 0 INDFILE notgu.popfile OUTFILE notgu.clumpp.out MISCFILE notgu.clumpp.miscfile K 8 C 3 W 0 S 2

然后运行CLUMPP：

CLUMPP paramfile

最后，用R绘制漂亮的群体结构图：

library(ggplot2) library(reshape2) data <- read.table("notgu.clumpp.out") data$Sample <- rownames(data) melted <- melt(data, id.vars="Sample") ggplot(melted, aes(x=Sample, y=value, fill=variable)) + geom_bar(stat="identity", position="stack") + scale_fill_brewer(palette="Set3") + theme_minimal() + theme(axis.text.x=element_text(angle=90, hjust=1))

这个可视化效果比Structure自带的输出要美观得多，适合直接放入论文。我曾经用这个方法发表过一篇群体遗传学文章，审稿人特别称赞了图表的质量。