从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: num1 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_C2. Structure软件运行与参数设置Structure是群体结构分析的核心工具但它的参数设置对新手来说可能有点复杂。我建议从mainparams和extraparams两个配置文件开始。mainparams中这几个参数需要特别注意numreps建议设置为10000确保结果稳定burnin预热迭代次数5000是个不错的起点nummcmc正式运行的迭代次数10000次通常足够inferalpha是否推断α值对于复杂群体建议设为1popinfo是否使用先验种群信息初学者建议设为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(xV1, yV7)) geom_line(colorblue) geom_point(colorred, size3) labs(xK value, yDelta 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然后运行CLUMPPCLUMPP paramfile最后用R绘制漂亮的群体结构图library(ggplot2) library(reshape2) data - read.table(notgu.clumpp.out) data$Sample - rownames(data) melted - melt(data, id.varsSample) ggplot(melted, aes(xSample, yvalue, fillvariable)) geom_bar(statidentity, positionstack) scale_fill_brewer(paletteSet3) theme_minimal() theme(axis.text.xelement_text(angle90, hjust1))这个可视化效果比Structure自带的输出要美观得多适合直接放入论文。我曾经用这个方法发表过一篇群体遗传学文章审稿人特别称赞了图表的质量。