如何进行群体遗传多样性分析,代码和结果解读
2026/5/14 15:38:22 网站建设 项目流程

当我们谈论物种的 “生命力”,往往更关注个体的生存与繁衍,却常常忽略一个更为底层的核心——群体遗传多样性。简单来说,群体遗传多样性是指一个种群内部,所有个体所携带遗传物质的总和及其差异程度。它具体表现为基因序列的微小变异、等位基因的多样组合,如同物种的 “基因储备库”:储备越丰富,物种抵御环境变化、病虫害侵袭等风险、应对生存危机的能力也就越强,更是物种进化的基础的核心动力。无论是保护濒危物种、培育优良品种,还是解析物种适应性进化机制,群体遗传多样性分析都是不可或缺的核心手段。

近年来,随着测序成本持续降低、测序技术快速迭代,高通量测序已成为群体遗传多样性研究的核心技术,可在全基因组水平高效、精准地检测单核苷酸多态性(SNP)、插入缺失(InDel)等遗传变异,为后续多样性分析提供了海量、可靠的数据支撑。伴随数据分析方法的不断完善,目前学界公认的核心遗传多样性评价指标主要包括:核苷酸多样性(Pi)、遗传分化指数(Fst)、Tajima’s D 检验统计量。接下来,我们将系统讲解这三个核心指标的分析思路、实操步骤及结果解读,帮助快速掌握群体遗传多样性的基础分析方法。

一 近核苷酸多样性(Pi)分析

物种多样性的实质是核苷酸多样性,它衡量的是群体内不同个体在某一基因片段或全基因组上的平均核苷酸差异水平,Pi值越高,说明群体内部遗传变异越丰富,遗传多样性越高。我们通常使用VCFtools软件对该指标进行计算,主要分为两种计算方式:窗口水平计算(适用于全基因组水平的整体分析)和单SNP水平计算(适用于特定位点的精细分析)。

1.1实操代码及说明

bash # 1. 以一定的窗口步长计算Pi(全基因组常用,便于后续可视化) ./vcftools --vcf sample.vcf --out result --window-pi 20000 --window-pi-step 10000 # 代码参数说明: # ./vcftools: 软件vcftools所在的路径(若已配置环境变量,可直接写vcftools); # --vcf sample.vcf:输入的VCF格式变异文件,包含所有样本的变异信息; # --out result:输出文件的前缀,后续会生成result.windowed.pi等输出文件; # --window-pi 20000:设置计算Pi的窗口大小为20000bp(可根据物种基因组大小调整,常用10000-50000bp); # --window-pi-step 10000:设置窗口滑动步长为10000bp,即相邻两个窗口重叠10000bp,确保覆盖全基因组。 # 2. 以单SNP计算pi(针对单个变异位点,用于精细分析特定区域的变异) ./vcftools --vcf sample.vcf --out result --site-pi # 代码参数说明: # --site-pi:指定以单个SNP位点为单位计算Pi值,输出文件为result.sites.pi,包含每个SNP的Pi值及相关信息。

1.2 结果说明

窗口水平计算会生成result.windowed.pi文件,包含染色体名称、窗口起始位置、窗口终止位置、窗口内SNP数量、Pi值等信息;单SNP水平计算会生成result.sites.pi文件,包含每个SNP的染色体位置、Pi值。后续可通过R语言ggplot2等包绘制Pi值分布图,直观展示全基因组或特定区域的遗传多样性分布特征。

遗传分化指数(Fst)分析

遗传分化指数(Fst)用于衡量两个或多个群体之间的遗传差异程度,取值范围为0-1。Fst值越接近0,说明群体间遗传差异越小,基因交流越频繁;Fst值越接近1,说明群体间遗传差异越大,基因交流越少,甚至出现生殖隔离。通常以窗口步长计算Fst值,结合文本处理和可视化,可清晰展示群体间的遗传分化热点区域。

2.1实操代码及说明

bash # 1. 以一定的窗口步长计算两个群体的Fst值 ./vcftools --vcf sample.vcf --weir-fst-pop group1.txt --weir-fst-pop group2.txt --out result-group1-group2 --fst-window-size 20000 --fst-window-step 10000 # 代码参数说明: # --vcf sample.vcf:输入的两个群体合并的VCF格式变异文件(需确保两个群体的样本均包含在该文件中); # --weir-fst-pop group1.txt:指定第一个群体的样本名文件(每行一个样本名,需与VCF文件中的样本名一致); # --weir-fst-pop group2.txt:指定第二个群体的样本名文件(格式同group1.txt); # --out result-group1-group2:输出文件前缀,后续会生成result-group1-group2.windowed.weir.fst等文件; # --fst-window-size 20000:计算Fst的窗口大小为20000bp(与Pi值计算窗口一致,便于对比分析); # --fst-window-step 10000:窗口滑动步长为10000bp。 # 2. 文本处理(去除无效值,整理格式,便于后续可视化) # 去除文件第一行标题,将Fst负值替换为0(负值无生物学意义,可能由计算误差导致) sed "1d" result-group1-group2.windowed.weir.fst|awk '{if($5<0)print $1"\t"$2"\t0";else print $1"\t"$2"\t"$5}' > 0-result-group1-group2.txt # 为文件添加行号和染色体前缀(适配CMplot包可视化要求) awk '{print NR" chr"$0}' 0-result-group1-group2.txt > 1-result-group1-group2.txt

2.2 数据可视化(R语言CMplot包)

CMplot包是全基因组关联分析及遗传多样性分析中常用的可视化工具,可快速绘制曼哈顿图、QQ图等,此处用于展示两个群体的Fst值分布,直观呈现遗传分化热点区域。

r # 加载CMplot包(若未安装,先执行install.packages("CMplot")) library(CMplot) # 读取处理后的Fst数据 g<- read.table(file ="1-result-group1-group2.txt",sep = "",header = F) # 绘制Fst值曼哈顿图 CMplot(g, plot.type="m", # 绘制曼哈顿图 LOG10=F, # 不进行LOG10转换(Fst值本身范围0-1,无需转换) cex=c(0.5,0.5,0.5), # 调整点、坐标轴标签、标题的字体大小 col=c("#073763","#999999"), # 设置点的颜色 ylab="(FST)_g1-g2", # Y轴标签,标注两个群体的名称 ylim=c(-1.5,15), # Y轴范围(根据实际Fst值调整,避免超出范围) threshold=NULL, # 不设置显著阈值线(可根据需求添加) chr.den.col=NULL, # 不显示染色体密度图 file="pdf", # 输出文件格式 file.name="FST", # 输出文件名称 dpi=300, # 图片分辨率(300dpi为学术常用) file.output=TRUE, # 输出图片文件 verbose=TRUE) # 显示运行日志

2.3结果解读

绘制的曼哈顿图中,每个点代表一个窗口的Fst值,横坐标为染色体位置,纵坐标为Fst值。Fst值较高的区域(通常Fst>0.2)可视为群体间的遗传分化热点区域,可能是受自然选择、地理隔离等因素影响,导致群体间基因交流受阻,形成遗传差异。

三 Tajima’s D 检验分析

Tajima’s D 检验是基于等位基因频率分布,判断群体是否经历选择压力、瓶颈效应、种群扩张等进化事件的重要指标。Tajima’s D值的取值范围无明确界限,通常认为:D>0时,可能存在平衡选择或种群收缩,等位基因频率趋于均匀;D<0时,可能存在正选择或种群扩张,低频等位基因频率升高;D≈0时,群体符合中性进化模型,未受到明显的选择压力或种群波动。

3.1 实操代码及说明

bash # 以一定步长计算Tajima’s D值 ./vcftools --vcf sample.vcf --TajimaD 10000 --out sample # 代码参数说明: # --TajimaD 10000:设置计算Tajima’s D的窗口大小为10000bp(可根据物种基因组大小调整); # --out sample:输出文件前缀,生成sample.TajimaD.txt文件,包含染色体名称、窗口起始位置、窗口终止位置、Tajima’s D值等信息。

3.2数据可视化(R语言ggplot2包)

ggplot2包具有强大的绘图功能,可绘制Tajima’s D值的折线图,按染色体分面展示,清晰呈现不同染色体的进化特征。

r # 加载ggplot2包(若未安装,先执行install.packages("ggplot2")) library(ggplot2) # 读取Tajima’s D数据(需替换为实际文件路径) mydata<-read.table(file ="F:/研究生数据/TajimaD/FL.TajimaD.txt",header = T) # 去除缺失值(避免绘图报错) mydata1<-na.omit(mydata) # 绘制Tajima’s D值折线图(按染色体分面) ggplot(mydata1, aes(x=BIN_START/1000000,y=TajimaD,group=factor(CHROM),colour=CHROM))+ geom_line()+ # 绘制折线图 facet_wrap(~mydata1$CHROM)+ # 按染色体分面,每个染色体单独展示 xlab("Physical distance(Mb)")+ # X轴标签,单位转换为Mb(便于阅读) ylab("Tajima's D")+ # Y轴标签 theme(legend.position = "none") # 隐藏图例(染色体颜色无实际意义,可根据需求保留) # 可选:随机抽取1000个数据点(用于快速预览,避免数据量过大导致绘图卡顿) a<-mydata1[sample(nrow(mydata1), 1000), ]

3.3 结果解读

分面折线图中,每个子图代表一条染色体,横坐标为物理距离(Mb),纵坐标为Tajima’s D值。若某一区域的Tajima’s D值显著偏离0,需结合研究背景分析可能的进化原因:例如,连续的D<0区域可能提示该区域受到正选择,与物种的适应性进化相关;而D>0区域可能提示该区域存在平衡选择,维持了群体的遗传多样性。

四 注意事项与总结

群体遗传多样性分析是解析物种进化、指导物种保护和品种改良的核心手段,本文围绕核苷酸多样性(Pi)、遗传分化指数(Fst)、Tajima’s D 检验三个核心指标,详细讲解了从数据计算、文本处理到可视化的完整实操流程,基于VCFtools和R语言实现了分析的可复现性。通过上述分析,可清晰掌握群体的遗传多样性水平、群体间的遗传分化程度及群体的进化历史,为后续的研究工作提供坚实的数据分析支撑。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询