企业官网建设流程全解析

做群体遗传分析的人，对PopLDdecay这个工具应该都不陌生。命令行简单、算得快，LD衰减分析这块它几乎是默认选择。但我发现一个普遍现象：很多人跑完软件就停在"出图"这一步——官方脚本点一下，存张图，完事。至于.stat文件里每列是什么、bin宽度怎么回事、曲线为什么长这样、怎么把多个群体放到同一张图里，基本没深入研究过。我早几年也是这样，直到有一次要给一个物种的六个地理种群做LD衰减对比，被审稿意见逼着重做图，才把"结果可视化"这条链路从头到尾捋清楚。这篇就把我的实践整理出来，内容包括stat文件结构解析、官方绘图脚本的用法和局限、用R做自定义LD衰减图的完整代码、多群体比较的作图思路，以及我踩过的一些坑。适合刚接触群体遗传学的学生，也适合想把图件升级成发表级别的研究人员。

1. 先弄明白：LD衰减曲线和stat文件里到底有什么

1.1 LD衰减的数据是怎么算出来的

连锁不平衡（Linkage Disequilibrium，LD）刻画的是不同遗传位点等位基因的非随机组合状态。最常用的衡量指标是r²，可以通俗理解成两个SNP位点之间的"相关程度"。物理距离越近的位点，重组打断它们组合的概率越低，r²一般越高；距离拉开后，重组累积效应会让r²逐步走低。把染色体上所有SNP对按物理距离分箱、在每个bin内取平均r²，再把平均值按距离顺序连起来，就得到经典的LD衰减曲线。

这条曲线的形状不是随机长出来的，背后有明确的群体遗传学含义：有效群体大小、重组率、交配制度、历史瓶颈和扩张事件都会影响它。自交为主的物种重组机会少，LD衰减普遍慢；异交物种通常衰减快。我们常说的"LD衰减距离"，就是指r²降到某个阈值或一半时对应的物理距离，它在GWAS标记密度设计、驯化历史推断、群体结构校正等场景里都很关键。所以可视化的重点不只是"画出来"，而是"画得对、读得懂"。

PopLDdecay的处理逻辑可以概括为：读入VCF，筛选位点，两两计算r²，按距离分箱统计平均值，最后输出汇总文件。它流行是因为实现上做了不少加速优化，处理几百个样本、几十万个SNP的VCF也不会等太久。同类计算也有人用别的工具，但PopLDdecay在内存控制和输出格式的简洁性上综合体验更好，这也是我把它作为主流程工具的原因。

1.2 读懂.stat文件：三列数据与bin粒度

跑完PopLDdecay之后，最重要的产出就是.stat文件，后续所有绘图都围绕它展开。一个典型的stat文件长这样，打开后通常能看到三列：第一列是物理距离（单位kb，取bin的中点或起始位置），第二列是该bin内所有SNP对的平均r²，第三列是参与计算的SNP对数量。

Distance Mean_r^2 Number_of_SNP_pairs 10 0.42 156320 20 0.31 143877 ...

要特别提醒的是，不同版本的PopLDdecay输出的stat文件在细节上略有差别。我见过两种：一种是统一bin宽度从头到尾；另一种是同一个文件里包含两种bin粒度的数据段，近距离段用细bin（比如每10 kb一个），远距离段用粗bin（比如每100 kb一个），中间用注释行分隔。这种设计的目的是保证近端区域有足够分辨率去捕捉LD快速下降的阶段，远端区域又不会因为bin太细导致每个bin里SNP对太少、均值波动太大。

无论哪种格式，拿到文件后先head看一下，确认分隔符、表头名称和bin粒度，再决定怎么读。不要想当然用read.table默认参数直接读，少了表头或把注释行读进去，后面画图全是乱的。第三列SNP对数也建议留意：如果一个bin里只有几十对SNP，它的平均r²是非常不稳定的，这类点往往就是曲线上那些离群噪点的来源。如果你开启了-OutLD参数，还会额外得到成对LD明细文件，体积通常很大，记录所有SNP对的r²和距离。这个文件日常画LD衰减图用不上，它的主要用途是画染色体区段LD热图和单倍型区块分析，大样本全基因组跑出来的明细文件动辄几个GB，不是必要不要轻易开启。

2. 运行环节：VCF预处理和参数选择的逻辑

2.1 输入VCF必须处理好的三件事

PopLDdecay吃的是VCF，和大多数生信工具一样，它要求文件是排序过、bgzip压缩并建立过索引的标准格式。如果你的VCF未排序或者干脆是纯文本，建议先整理一遍再跑，否则后续按距离分箱时容易出现逻辑错误，甚至直接中断。

第二个容易忽略的点是contig信息。PopLDdecay计算距离依赖位点在染色体上的坐标，如果VCF的contig命名混乱，或者混合了不同组装的坐标，算出来的距离就是不可信的。我处理过一个外部提供的VCF，跑出来LD衰减曲线前几个点居然在上升，排查了半天才发现是contig排序的问题。现在我的习惯是：任何VCF进PopLDdecay之前，先确认排序、压缩、索引三件套齐全，再手动抽查几个位点的坐标顺序。

第三个问题是计算量。样本量特别大或者SNP特别多的时候，全基因组一次性跑内存压力会很大。最省心的做法是按染色体拆分VCF分别运行，最后合并各染色体的stat文件。拆分本身不损失信息，因为不同染色体之间本来就不应该比较LD，连锁只发生在同一条染色体内部。这一步很多教程没提，但实测对降低内存占用非常有效。

2.2 核心参数怎么选，以及背后的判断依据

PopLDdecay的核心参数不多，但每个都会直接影响stat文件和最终图形。

-MaxDist是最大两两距离阈值，默认300 kb，决定曲线横向范围。选择依据是物种的LD衰减水平：衰减快的外来异交物种，可能几十kb内r²就基本归零，300 kb够用；衰减慢的自交物种或经过强驯化的群体，r²拖到几百kb还很坚挺，这时候可以把-MaxDist调大到500甚至1000 kb，让曲线完整展示衰减过程。判断方法很简单：先看默认参数下曲线尾部是否已经降到接近0的平台，如果还在下降，说明范围不够要加大。

-MAF和-Miss是两个过滤参数，分别控制最小等位基因频率和缺失率。默认值（MAF 0.05，Miss 0.2到0.5之间）对大多数场景是合理的。需要提醒的是，过滤越严格，剩余SNP对越少，长距离bin里的SNP对数会急剧下降，尾部噪点变多。所以不是过滤越狠越好，要在数据干净和统计稳定之间找平衡。

-SubPop参数用于分组分析，输入是一个两列的分组文件，第一列是样本名，第二列是分组名。开启后软件会按组分别输出统计文件。做群体比较时我强烈建议用这个参数，而不是自己把VCF按群体拆开再分别跑，因为后者会让各群体的SNP集合不一致，r²量级不可比。

2.3 一条完整命令和运行后的检查

我常用的命令是这样的：

PopLDdecay -InVCF input.vcf.gz -OutStat species -MaxDist 300 -MAF 0.05 -Miss 0.2

跑完后先执行下面几步确认运行质量：

head -20 species.stat wc -l species.stat

看一眼行数是否正常、前三列是否完整，顺便确认bin宽度是否符合预期。如果stat文件行数少得可怜，多半是过滤参数太严格把SNP几乎滤光了，或者输入VCF本身位点数不够。这道检查很关键，很多人画出来的图只有寥寥几个点，问题就出在这步被跳过了。

如果是分组分析，命令变成：

PopLDdecay -InVCF input.vcf.gz -OutStat species -MaxDist 300 -MAF 0.05 -Miss 0.2 -SubPop group.list

输出会按组生成species.stat.1、species.stat.2等文件，顺序对应group.list里的分组顺序。画多群体对比图时，把这些文件读进去合并处理即可。

3. 从省事到精细：三条可视化路线

3.1 官方Plot_OnePop.pl适合快速验证数据

PopLDdecay自带Perl绘图脚本Plot_OnePop.pl，可以一键把stat文件画成LD衰减散点图，基本用法：

perl Plot_OnePop.pl -inStat species.stat -output species_LD -bin1 10 -bin2 100

这里的-bin1和-bin2要跟stat文件里两种bin粒度的实际值对应。如果stat文件从头到尾只有一个bin宽度，两个参数填同一个值就行。脚本会在输出目录生成PNG图。这条路线优点是快、无脑，适合内部检查和初步趋势判断。但局限也很明显：样式基本固定，颜色、点大小、坐标轴范围都不能灵活调整；输出是PNG位图，对期刊投稿来说分辨率不够；而且脚本依赖Perl的GD绘图模块，很多服务器默认没装，一跑就报"Can't locate GD.pm"，还得先折腾模块安装。

所以我给的建议是：官方脚本只用来做快速验证，确认数据没有大问题；真正要出图、要比较、要投稿的阶段，直接走R路线。

3.2 R加ggplot2画出可发表的定制图

我现在的主力方案是R，灵活度和可控性完全不一样。第一步读文件，注意不同版本stat文件的差异，先把列名标准化：

library(ggplot2) stat <- read.table("species.stat", header = TRUE, comment.char = "#", check.names = FALSE) # 先用 names(stat) 查看实际列名，再按下面方式统一命名 names(stat) <- c("Distance", "Mean_r2", "Npairs") stat <- stat[is.finite(stat$Mean_r2), ]

然后画散点加平滑曲线：

p <- ggplot(stat, aes(x = Distance, y = Mean_r2)) + geom_point(size = 0.8, alpha = 0.5, colour = "grey50") + geom_smooth(method = "loess", span = 0.3, se = FALSE, colour = "#B2182B", linewidth = 1) + scale_x_continuous(limits = c(0, 300), breaks = seq(0, 300, 50)) + labs(x = "Physical distance (kb)", y = expression(italic(r)^2)) + theme_classic(base_size = 14) ggsave("LD_decay.pdf", p, width = 7, height = 5)

几个值得注意的细节。geom_point的size和alpha要根据SNP对数调整：点太密就调小size、加大透明度，不然图上是一团黑。loess的span值决定平滑程度：span太小曲线跟着噪点走，span太大又抹平真实衰减趋势，我一般从0.3开始试，上下乱窜就增大，过度平滑就减小。长距离尾部噪点太严重时，直接用scale_x_continuous限幅截断，截到曲线进入平台期的地方就够，不必展示所有距离。ggsave输出PDF或SVG矢量格式，投稿时编辑字体、缩放都不会糊。

3.3 多群体LD衰减对比图的合并画法

多群体比较是文章里最常出现的场景。我的做法是把分组运行得到的多个stat文件合并进一个数据框，然后一图画多条曲线：

files <- c("A.stat", "B.stat", "C.stat") pop_names <- c("Population A", "Population B", "Population C") df_list <- lapply(seq_along(files), function(i) { d <- read.table(files[i], header = TRUE, comment.char = "#", check.names = FALSE) names(d) <- c("Distance", "Mean_r2", "Npairs") d$Population <- pop_names[i] d }) df <- do.call(rbind, df_list) ggplot(df, aes(x = Distance, y = Mean_r2, colour = Population)) + geom_point(size = 0.5, alpha = 0.3) + geom_smooth(method = "loess", span = 0.3, se = FALSE, linewidth = 1) + scale_colour_manual(values = c("#4575B4", "#F46D43", "#807DBA")) + labs(x = "Physical distance (kb)", y = expression(italic(r)^2)) + theme_classic(base_size = 14) + theme(legend.position = c(0.75, 0.8))

这里有个关键细节：如果各群体样本量差异很大，SNP对数差异也会很大，直接比较平滑曲线可能失真。画图前先检查各群体stat文件第三列的量级是否可比，如果差异悬殊，说明比较本身不公平，需要在分析层面统一样本量或过滤标准，而不是只调图。另外，分组超过五个时，我一般只画平滑曲线去掉散点，图面更干净，重点是看曲线之间的相对位置和下降速率。

4. 学会读图：曲线形态和量化指标

4.1 从曲线形状判断群体特征

LD衰减曲线不只是数学上的下降线，不同形状对应不同的群体历史，解读时要有意识。

快速衰减的群体，曲线在很短的物理距离内断崖式下降，r²从高位迅速跌到接近0。这种情况常见于有效群体大小较大、异交率高的群体，或者历史上经历过扩张的群体。重组机会多，LD被破坏得快，位点之间的关联只存在于短距离内。

缓慢衰减的群体，曲线下降平缓，几十kb甚至几百kb之后r²依然维持在较高水平。这通常对应有效群体大小较小、自交比例高，或者经历过瓶颈与近交的群体。这类群体的LD区块很长，意味着关联分析需要的标记密度较低，但全基因组水平的统计推断要小心假阳性。驯化物种和经过强烈人工选择的群体经常表现出这种模式，因为选择会拖拽附近区域的变异，人为延长了LD。

还有一种值得警惕的形态：曲线尾部出现平台或反弹，r²没有随距离下降，反而维持在一个非零平台甚至略有上升。这通常不是真实的生物学信号，而是群体结构或样品混杂造成的伪LD——不同亚群混合在一起，会让远距离位点也表现出相关性。看到这种形态，第一反应应该是检查群体是否被正确分层，而不是急着解读。

4.2 用两个指标量化LD衰减距离

文章里光放曲线不够，审稿人通常希望看到量化的LD衰减距离。最常用的是两个指标：半衰距离（r²下降到最大值一半时对应的物理距离），以及阈值距离（r²降到0.1或0.2时对应的距离）。

用R计算很简单：

stat_fit <- stat[is.finite(stat$Mean_r2), ] r2_max <- stat_fit$Mean_r2[1] # 第一个bin的r2通常最高，比max更稳 half_dist <- stat_fit$Distance[which(stat_fit$Mean_r2 <= r2_max / 2)[1]] r2_02 <- 0.2 dist_02 <- stat_fit$Distance[which(stat_fit$Mean_r2 <= r2_02)[1]]

注意取基准值时不要用全表的max，因为曲线前段的异常噪音点可能比第一个bin的r²还高，用它算半衰距离会被带偏。我习惯直接用第一个bin的均值作为基准，更稳健。

报告这两个指标时，我一般按表格整理，配合曲线图一起呈现：

有曲线但没有数字支撑的分析，深度是不够的，这个表格能有效补上短板。

5. 常见问题速查与避坑记录

5.1 高频问题速查表

5.2 几个真金白银踩出来的坑

第一个坑是VCF没排序。有次拿到外部VCF直接跑，结果LD衰减曲线前面一段居然在上升，当时以为发现了什么特殊生物学信号，折腾了好几天想解释。后来才查出来是染色体内部SNP顺序混乱，坐标距离算得毫无意义。现在我的铁律是：任何VCF进PopLDdecay前，排序、压缩、索引三件套缺一不可，跑之前手动抽查几个位点的坐标顺序。

第二个坑是解析stat时把列搞混。有次自己写脚本，想当然以为第二列是SNP对数、第三列才是r²，结果画出来的曲线上下颠倒，数值明显不对。从那以后我养成了一个习惯：每次拿到新输出，先打印前五行核对列名，再做任何操作，哪怕这个文件已经跑过一百遍也不跳过。

第三个坑是过滤参数太贪心。有一回为了让结果更"干净"，把MAF调到0.2、缺失率卡到0.1，结果长距离bin里几乎没有SNP对了，尾部全是单点噪声，怎么平滑都救不回来。后来才意识到，过滤的每一刀都会削减SNP对数量，而SNP对数量恰恰是统计稳定性的基础。现在我会先在默认参数下跑一遍，检查SNP对数量的分布，再根据实际情况微调，绝不为了好看乱砍参数。

第四个坑和分组比较有关。做多群体对比时，我一开始图省事，各群体用各自的VCF独立跑，后来发现不同群体的变异位点集合差异很大，r²量级根本不在同一水平线上。改用-SubPop参数统一处理之后，所有群体共享同一套过滤和计算流程，可比性才真正成立。这个教训我记到现在：比较的前提是流程一致，不是画图一致。

最后分享一点个人体会。PopLDdecay出图这件事，看起来是收尾工作，实际上是最值得花时间打磨的环节之一。数据本身没有问题，图的表达质量直接决定别人能不能快速抓住群体的核心特征。我现在的标准流程是：官方脚本快速巡检一遍确认数据正常，R里做精细出图和群体比较，同时保留半衰距离和阈值距离两个数值指标配合曲线呈现。每次拿到新数据集，先跑默认参数，看stat结构、SNP对数量和曲线基本形态，再决定怎么调参数。这套流程帮我在不少项目里省下了反复返工的时间，希望对你也有用。如果你碰到stat格式差异或者曲线形态不知道怎么解释的情况，回头对照这篇的检查步骤，多半能找到问题所在。