引言拿到测序数据第一件事为什么是PCA做群体遗传学、医学遗传学或者生态基因组学大家手里最终都会有几套基因型数据可能是全基因组重测序可能是芯片分型也可能是简化基因组。数据一到手很多人第一反应就是上STRUCTURE或者ADMIXTURE去猜群体构成但我个人的习惯恰恰相反我永远先跑PCA也就是主成分分析。先跑PCA有四个现实理由它速度极快、几乎不需要先验假设、能一眼看出样本有没有分群趋势更重要的是能帮你快速识别离群样本和可能的样本标签错误。这些信息在后边的群体分层校正、亲缘关系分析、系统发育推断里全都用得上。这篇博文就围绕“群体PCA分析及可视化”这条主线从数据准备、质控、工具选型、参数设置一直聊到R语言出图把整套流程拆开揉碎讲清楚顺便把我踩过的坑也一并交代了。这篇内容适合谁只要你的工作涉及基因组测序数据的群体结构分析无论是做作物群体进化、家养动物驯化历史还是人类群体遗传学都可以照着流程走一遍。小白也能跟得上我会从文件格式和底层概念讲起。1. 群体PCA的基本思路与选型逻辑1.1 从数学直觉理解PCA在群体中看到的东西PCA本身是个通用降维算法但在遗传数据上有它特殊的“翻译”方式。每个个体在几十上百万个SNP位点上的基因型可以理解成一个超级高维的向量比如一个个体有50万个SNP那这个个体就是50万维空间里的一个点所有个体共同构成了一个高维点云。这50万个维度不是独立的很多位点之间由于连锁不平衡存在强相关而这些相关性的根本来源很大程度来自群体历史。假设一个群体经历了隔离分化隔离时间内两个子群体各自积累不同的突变组合等位基因频率朝着不同方向漂移体现在数据上就是高维空间里样本点沿少数几个主轴明显分开。PCA做的正是从高维协方差结构里抽取特征向量最大的几个方向把个体的相对位置映射到这些方向上形成一个可以画在二维平面上的散点图。有趣的是做PCA时不需要预先指定“有几个群体”也不需要在输入文件里写群体标签它纯粹是一个无监督探索过程。这一点在实际应用里非常重要因为我们通常对群体的先验划分并不确定特别是野生材料、地方品种、混杂个体混在一起的时候强加标签反而容易让分析变得先入为主。所以我常跟同事说PCA是群体遗传学里最诚实的一项可视化技术它只基于数据本身说话。1.2 不同分析场景下PCA扮演的角色不同研究场景PCA解决的核心问题不太一样。在群体遗传结构研究里PCA主要用来确认分群模式是否与地理分布、物种信息、表型类型吻合常见的结论形式是“几个主成分把样本分成几大簇簇间差异对应物种/群体分化”。在医学或动植物全基因组关联分析GWAS里PCA的主要用途则是“主成分作为协变量纳入统计模型校正群体分层导致的假阳性”这个场景下PC的轴本身不需要解释只要它在统计上能代表潜在结构就行。在生态学里PCA常用于分析环境梯度对应的微地理结构比如高度隔离的小种群在连续环境梯度上是否呈现出渐变模式。选型逻辑也受这些场景影响。如果只是快速检查样本聚类和离群点PLINK1.9甚至不需要太复杂的参数一句话命令就能跑。如果要做严谨的群体遗传学分析特别是参考了1000 Genomes这样的国际项目工作流程很多人会选用EIGENSOFT的smartpca因为它自带一些针对遗传数据设计的细节处理比如自动去除极端离群样本后再迭代计算只是输入输出都比较老派需要用特定的格式。后面我会给出一套结合PLINK和EIGENSOFT的完整操作。1.3 PCA与ADMIXTURE、系统发育树是互补关系新手容易陷入一个误区觉得做一种分析就可以知一切了。PCA、ADMIXTURE、NJ树三者在信息层次上完全是互补的。PCA擅长在连续二维空间里展示个体间的相似关系但对于某些历史事件比如一次非常古老的群体混合PCA往往只能看到“个体散成一条线”此时ADMIXTURE可以量化祖先组分比例提供更直观的混合解释。系统进化树则呈现的是离散的分支关系适合展示谱系亲缘层级但对连续渐变的结构不如PCA平滑。在实际项目里我的习惯是先用PCA做全局概览确认样本没有标注错误、离群样本都被识别出来接着用ADMIXTURE跑1到10个K值看区间稳定情况最后在系统发育树上验证主要分支。三者结论互相印证共识才敢写进论文。这里先明确PCA的定位后面每一步才能踩在点上。2. 数据准备与质量控制出图前最该花时间的环节2.1 输入格式与标准处理链路PCA分析并不直接吃FASTQ或者VCF里的多元位点信息绝大多数工具输入都是二进制PLINK格式也就是bed/bim/fam三个文件或者从VCF转换而来。如果是从头开始的测序数据整个流程可以整理成一条固定链路FASTQ比对到参考基因组得到BAM文件再用GATK或bcftools进行变异检测得到VCF然后用PLINK把VCF转成bed格式。我给一个常用的转换命令适用于PLINK1.9和PLINK2.xplink --vcf input.vcf.gz --make-bed --out project --double-id --allow-extra-chr 0--double-id的原因是VCF里有些样本名重复出现在各个位点加上这个参数可以在转换时消除编号重复匹配问题。--allow-extra-chr 0则是因为很多非人类物种的染色体命名可能带有数字前缀或者scaffold默认参数没包含这些染色体名会报警。这一步转换后的bed文件中存的是0/1/2计数的等位基因拷贝数0代表参考纯合1代表杂合2代表变异纯合。很多人在拿到公共数据库数据的时候下载到的已经是PLINK格式那就可以直接跳过转换进入质控阶段。不管是哪种来源都要在PCA之前做质控这个环节决定了结果的可信度上限。2.2 按位点和样本双重过滤基因型数据质控分为位点层面和样本层面缺一不可。位点层面核心参数是缺失率--geno、次等位基因频率--maf和哈迪温伯格平衡检验--hwe样本层面核心参数是缺失率--mind。基因组变异的质控阈值在不同场景下差别不小我给出一套群体结构分析常用的配置并用注释说明每个参数的意图plink --bfile project --geno 0.1 --mind 0.1 --maf 0.05 --hwe 0.000001 --make-bed --out project_qc这里的--geno 0.1表示删除在超过10%样本中缺失的SNP位点--mind 0.1表示删除基因型缺失超过10%的样本这两个参数是为了避免缺失率过高带来的系统偏差。--maf 0.05表示只保留次等位基因频率大于5%的位点极小频率的位点在PCA中容易产生单点驱动的假信号。--hwe 0.000001是严格的哈迪温伯格过滤目的其实是尽量剔除明显由测序或基因分型错误导致的位点。这里有个特别重要的实操心得在做群体PCA时哈迪温伯格过滤的阈值必须结合样本分组来看如果把多个结构差异很大的群体放在一起跑全样本的HWE那些导致分群的位点反而可能带有极强的HWE偏离而被错误过滤。更稳妥的做法是先不做HWE过滤跑一次PCA看清分群结构后再按已知群体分组分别做HWE过滤或者干脆只过滤极端的HWE偏差把阈值放到1e-10甚至更低。我在实际项目中通常在有明确分组信息时报--hwe 0.000001没有分组信息时先跳过HWE。2.3 LD剪枝防止连锁位点主导主成分连锁不平衡是指两个位点由于物理距离近、重组率低而在群体中具有相关性的等位基因组合在PCA中如果包含大量连锁位点这些区域联合起来会形成巨大的冗余信号掩盖真正的群体演化关系。比如一个染色体倒位区段里可能有几千个高度连锁的位点它们的信号在PCA中会被过度放大以至成簇的样本不是依据全基因组相似性聚合而是单纯因为某几个重组冷点区域的一致性被聚在一起。因此所有严谨的PCA分析跑之前都会做LD剪枝即选取相互之间连锁程度较低的位点集合plink --bfile project_qc --indep-pairwise 50 5 0.2 --out project_prune plink --bfile project_qc --extract project_prune.prune.in --make-bed --out project_pruned--indep-pairwise 50 5 0.2的含义是以50个SNP为窗口每次滑动5个SNP计算窗口内两两SNP的LD剔除r²大于0.2的位点。窗口大小和r²阈值直接影响保留位点数目。r²阈值设到0.2是群体分层分析中比较常规的平衡点既保留足够的独立信号又不引入太多连锁冗余。如果是做亲缘关系估计或者近交系数阈值可能要更严比如0.05或0.1。这一步做完后还要检查一下保留的位点数如果剪枝之后只剩下几万位点在几十人的小群体里能用的信号就比较勉强了如果是全基因组重测序数据剪枝后保留十几万到几十万位点是很正常的此时PCA的目视结果会干净很多。2.4 样本查重与亲缘关系初步筛查群体PCA一个容易被忽视的坑是样本重复或者存在极近亲缘关系个体。同卵双胞胎或者高度近交的样本在PCA图里会重叠成同一个点如果样本标签恰好把两个批次搞混还会把原本无关联的群体错误连接。处理办法是用PLINK2的KING或GCTA的GRM做亲缘关系估计然后剔除掉PI_HAT大于0.25的一级亲属样本。plink2 --bfile project_pruned --make-king-table --out project_king awk $3 0.25 {print $1, $2} project_king.kin0 related_pair.txt实际操作中大批量公共数据经常有重复样本残留这步骤能帮你省掉很多后续解释上的麻烦。注意这一步要在LD剪枝后的数据上做因为连锁位点会严重扭曲亲缘关系估计。3. 核心工具选型与参数实测对比3.1 用PLINK1.9快速跑PCA整套流程PLINK是群体遗传学里最常用的工具之一1.9版自带的--pca能够同时输出特征向量投影坐标和特征值。在LD剪枝后的数据上直接运行plink --bfile project_pruned --pca 10 --out project_pca--pca 10表示计算前10个主成分对于绝大多数群体结构研究保守估计前10个PC已经覆盖了绝大多数可解释的群体结构。输出包括两个关键文件project_pca.eigenvec是样本在前10个PC的坐标值格式是FID、IID、PC1、PC2……PC10project_pca.eigenval是每个PC对应的特征值特征值除以总特征值可以得到每个PC解释的方差比例。PLINK2也在做同样的事命令基本一致plink2 --bfile project_pruned --pca 10 approx --out project_pca2approx表示使用近似算法在样本和位点较多时速度提升很明显实测结果与精确算法差别通常很小。PLINK2的eigenvec文件格式略微不同默认多了一列列头后续处理要注意。PLINK的优势在于依赖少、跨平台、速度快缺点在于它对离群样本不敏感全部样本无论有多极端都会被纳入计算所以跑完一次PCA后必须自己看图剔除离群值。3.2 用EIGENSOFT smartpca做更接近群体的分析很多人在论文中看到的PCA图标注的是EIGENSOFT因为它就是1000 Genomes项目里群体PCA的主力工具。smartpca对数据格式有要求需要EIGENSTRAT格式的三种文件后缀为.geno的基因型矩阵.snp的位点信息.ind的个体信息。从PLINK格式转成EIGENSTRAT格式可以用convertf命令配置文件写法如下genotypename: project_pruned.bed snpname: project_pruned.bim indivname: project_pruned.fam outputformat: EIGENSTRAT genotypeoutname: project.pruned.geno snpoutname: project.pruned.snp indivoutname: project.pruned.ind转换完之后创建smartpca参数文件这里贴一份我在实践中反复调整过的版本genotypename: project.pruned.geno snpname: project.pruned.snp indivname: project.pruned.ind evecoutname: project.pca.evec evaloutname: project.pca.eval outlinemode: 2 numoutevec: 10 numoutlieriter: 0 nsnpldregress: 0outlinemode: 2的作用是把离群样本信息单独写一个文件而不是直接把它从主输出里剔除。numoutlieriter我习惯设置为0这一步是很多教程会讲但我不推荐默认开启的迭代去离群步骤。原因在后面第5.3节细说简单讲就是好几次我用了它的默认去离群迭代结果把一些真正有生物学意义的边缘群体给当成噪声剔掉了这个默认的“智能”机制在野生群体上经常自作主张。跑完后evec文件第一行有个约定俗成的无用行正式分析时删掉即可。EIGENSOFT最大的特点是它在计算的时候会自动把每个PC轴做标准化并且输出每一轮迭代中离群样本的信息。如果你的数据是不同测序深度、不同实验室来源合并出来的smartpca的结果里能够暴露更多样本层面的批次效应。3.3 用GCTA跑PCA处理大样本场景GCTA本身是作为复杂性状遗传力估计工具开发的但它提供的--pca功能在处理特别大的样本量时非常稳定因为它基于基因关系矩阵GRM直接做特征分解。命令如下gcta64 --bfile project_pruned --autosome --make-grm --out project_grm --thread-num 10 gcta64 --grm project_grm --pca 10 --out project_gcta_pca这两步的思路先从基因型计算全样本的GRM这是一个N乘N的矩阵代表个体间全基因组相似度然后对这个矩阵做特征值分解得到的特征向量就是每个样本在新空间里的坐标。当样本数达到数万的时候PLINK的PCA可能会因为内存占用过高而变慢GCTA的GRM流程配合作业调度器分块计算更从容。当然PLINK2也有--pca的近似算法在万级样本上实测速度不差具体选哪个更多看个人习惯和工作流环境。3.4 工具对比与选型建议工具输入格式擅长场景注意点PLINK1.9/2.0bed/bim/famVCF自动转常规群体结构快速探索、离群样本初筛不会主动识别离群样本需结合图判断EIGENSOFT smartpcaEIGENSTRAT正式论文级群体PCA含离群值记录格式转换繁琐参数文件有坑GCTAbed/bim/fam数万级样本PCAGRM需复用需要两步操作生成GRM耗内存我的建议是日常探索阶段用PLINK正式分析阶段输出也用PLINK的eigenvec配合R作图如果审稿人要求用EIGENSOFT复算一遍再补跑也不迟。但有一点要留意不同工具算出的PC1、PC2方向可能呈镜像翻转这是特征向量不定号的正常现象不是分析出错。4. 绘图的完整流程与实用进阶技巧4.1 用R语言ggplot2绘制干净的PCA散点图可视化是整个流程里最有成就感的部分也是踩坑比例最高的部分。经过前面几轮质控和PCA跑完现在手上有project_pca.eigenvec和project_pca.eigenval两个文件。推荐直接用R语言读取因为后续要做颜色映射、添加椭圆、调整图例标签都比其他工具灵活很多。library(tidyverse) eigenvec - read_table(project_pca.eigenvec, col_names TRUE) eigenval - read_table(project_pca.eigenval, col_names eigenvalue) # 向量顺序在PLINK2里可能涉及是否含列头读取时留意 pve - eigenval$eigenvalue / sum(eigenval$eigenvalue) * 100 data_plot - eigenvec %% mutate(population case_when( str_detect(#IID, ^W) ~ West, str_detect(#IID, ^E) ~ East, TRUE ~ South )) ggplot(data_plot, aes(x PC1, y PC2, color population)) geom_point(size 2, alpha 0.8) stat_ellipse(aes(fill population), geom polygon, alpha 0.1, level 0.95) labs(x paste0(PC1 (, round(pve[1], 2), %)), y paste0(PC2 (, round(pve[2], 2), %))) theme_bw(base_size 14) theme(legend.position right)这里pve计算的是主成分解释的方差比例放在坐标轴标签里才符合论文规范。很多初学者容易忽略这个只写PC1和PC2读者完全不知道这张图到底解释了多少变异。群体标签的映射可以根据自己的样本命名规则修改我的习惯是写一个干净的映射表明确样本ID与群体名的对应关系不推荐在脚本里硬编码一堆正则。stat_ellipse加的95%置信椭圆在群体数量多、重叠严重的时候反而会让画面挤在一起这时候可以选择不加椭圆或只给某些重点群体加。如果PCA图上某几个群体分不开而你现在没有先验标签就不要强行着色先用离散点看真实分布。4.2 PC坐标轴的正负向问题与批量出图根据自己的实际经验刚跑出来的PCA有时会出现明明该分在右边的群体跑到左边或者第一主成分和第二主成分对应关系整个反过来的情况。这是PCA本身固有的一个特征特征向量的方向符号没有确定意义PC1指向“左”还是“右”并不重要重要的是不同样本在PC1上的相对距离。作图展示时为了让视觉更直观可以自定义乘以-1把坐标翻转本质上并没有改变样本间的相对结构。在这类图里我会习惯性把PC1-PC10全部画出来生成一个大图矩阵方便检查后边几个PC是否存在微结构或者批次效应library(patchwork) plot_list - list() for (i in 1:9) { pc_x - paste0(PC, i) pc_y - paste0(PC, i 1) p - ggplot(data_plot, aes(x .data[[pc_x]], y .data[[pc_y]], color population)) geom_point(size 1.8, alpha 0.7) labs(x paste0(pc_x, (, round(pve[i], 2), %)), y paste0(pc_y, (, round(pve[i 1], 2), %))) theme_bw(base_size 11) plot_list[[i]] - p } wrap_plots(plot_list, ncol 3)第3个或第4个PC如果暴露出清晰的亚群分离但PC1、PC2看不出来这种情况在存在地理隔离和环境适应梯度时很常见。如果PC10以后还能看到明显分群往往意味着数据质量有问题比如样本污染、DNA混池、测序深度变异带来的批次效应或者近缘物种混入。4.3 碎石图、三维图与交互式输出除了散点图特征值碎石图在结果展示里也用得到。特征值下降坡度能直观提示你需要保留多少个主成分一般来说特征值趋于平台的点就是有意义信号与噪声的分界。基础绘制方式非常简单pve_df - tibble(PC 1:length(pve), PVE pve) ggplot(pve_df, aes(x PC, y PVE)) geom_bar(stat identity, fill steelblue, alpha 0.8) geom_line() labs(x Principal Component, y Variance Explained (%)) theme_bw()三维图在群体样本量少、结构单纯时其实不太常用但在探索多条地理支系同时分化、样本沿着三方向分离的场景里一个三维旋转视图比三张二维图更容易看清立体结构。R语言的plotly包可以基于已经算好的PC坐标直接生成可交互的三维图鼠标拖拽旋转非常直观library(plotly) plot_ly( data_plot, x ~PC1, y ~PC2, z ~PC3, color ~population, type scatter3d, mode markers, marker list(size 4) )这种交互式HTML图适合在组会汇报时现场演示但论文投稿阶段还是老老实实输出高清静态图。保存出图时建议用ggsave搭配dpi 300、宽高按实际需要设置PDF格式作为矢量图也不错。4.4 分析脚本的整理与会话管理我自己的绘图脚本习惯是用sessionInfo()把R版本和所有依赖包版本打印到日志文件里保证日后复现时不会因为包版本更新导致结果差异。多个人协作的时候这个做法能让图表质量分析可追溯。此外如果要做多个不同数据集、不同剪枝阈值的PCA建议给每个版本单独建一个干净的目录而不是全部堆在同一个文件夹里否则后边的图和中间文件极易混淆。我见过有同事把十几个版本的PCA输出放在同一目录文件名只差一个下划线最后提交论文时选错图的事故并不少见。命名规则从第一天就定好能省很多事。5. 常见问题与排查技巧实录5.1 PCA结果出现镜像、正负号翻转怎么办镜像、正负号翻转是PCA的固有现象因为特征分解得到的特征向量在正方向上没有数学约束比如PC1的载荷向量取反之后仍然是同一个主成分。因此同一份数据用PLINK、EIGENSOFT、R的prcomp跑出来的结果在坐标正负方向上可能不一致这就是为什么有些论文里的图左侧群体在你的图里变成右侧而不是数据错误。要判断两个结果是否一致正确做法是比较距离矩阵而不是直接比较坐标正负。也就是说把样本间的PC距离矩阵存下来跑一个两两比较如果距离矩阵高度吻合说明算法方向上虽然翻转但结构相同可以直接在可视化时统一方向。R里用dist()就能做到。dist_mat1 - dist(eigenvec_plink[, 3:12]) dist_mat2 - dist(eigenvec_eigensoft[, 2:11]) cor(as.vector(dist_mat1), as.vector(dist_mat2))如果相关系数接近1那么放心用任意一版出图。5.2 样本量不均造成的假“分群”PCA对样本量非常敏感。一个群体如果包含200个个体另一个群体只有10个个体哪怕它们在真实遗传距离上差异不太大PCA也倾向于把样本量大的群体向远离原点方向“拉”因为协方差矩阵中的方差主要由大群体决定。这种情况下的PC1解释的很大一部分变异其实只是样本量配比造成的并不完全反映深层的分群历史。针对这个问题的办法有三个第一在条件允许时尽量平衡各群体的样本量每个群体至少保持20到30个个体以上第二在做正式分析时对比有和没有大样本群体时PCA的差异假如去掉大群体后原本的“分群”马上消失说明区分度很大程度来自样本量而非遗传差异第三结合ADMIXTURE的交叉验证结果检查如果K值对应的组分比例没有稳定的停滞点那就要小心PCA里的某些分群是伪信号。实际工作中很多野生群体的样本量天然不均衡我能给的建议是不要靠肉眼崇拜一个漂亮的分离图多问一句这分离到底稳不稳。5.3 要不要自动剔除离群样本再跑一版针对“第一次PCA发现有离群样本该不该剔除”这个问题我的回答是离群样本要区分情况。如果是明显的样本污染、标签错误或者与已知群体完全不相干的近缘物种混入应该剔除。但如果离群样本来自一个新群体或者来自一个遗传分化极深的地方品种此时将其剔除可能会抹掉真实生物学信息。EIGENSOFT默认的离群迭代剔除会把偏离超过6个标准差的样本去掉然后重新计算PCA。这个机制在大规模公共数据清洗中有用比如1000 Genomes那种标准化流程。但在我处理一些地方品种和野生群体时这个自动剔除逻辑几次把真正稀有的边缘群体删掉了害得我后来重跑分析才发现重要群体不在图里。稳妥的流程是先保留全部样本跑一次PCA手动查看离群样本对应的原始样本记录和测序质量指标人工决定去留再跑第二次。手动标记比盲目套用算法可靠得多。5.4 PCA结果与实际分组不一致时怎么处理有时PCA的聚类结构和地理分布、个体表型或者已知谱系关系不一致别急着怀疑数据。优先级最高的检查是数据个体编号是否与表型记录对齐只要样本编号错位一个整张图可能都是垃圾。其次检查是否出现了批次效应比如两个批次的DNA提取时间不同、测序批次不同这些技术偏差在PCA里也会造成人为聚类典型特征是PC上形成与实验批次严格吻合的簇而这种簇的生物信息无法解释。处理这种批次效应需要先做样本重测或使用统计方法校正。如果技术因素排除了聚类结果与地理不一致也可能是基因流、驯化后再杂交、人为迁移等因素的真实体现。做群体遗传学就该明白PCA图只是一面镜子它反映的是数据中真实的遗传相似性真实历史往往比想象更复杂。5.5 方向轴不可解释的噪声与PCA结果稳定性验证PCA结果稳定性是很容易被忽视的问题。评价PCA稳定性比较常见的做法是多次抽样或自助法比较结果。比如从数据中随机抽取80%的位点重复做一次PCA再用Procrustes分析看看两次坐标之间的一致性。如果翻来覆去结果差别很大说明PCA结构并不稳健不能用来支撑任何群体历史结论。R里可以用vegan::procrustes()完成这个比较。res1 - prcomp(t(geno_matrix[, sample_1]), scale. TRUE) res2 - prcomp(t(geno_matrix[, sample_2]), scale. TRUE) proc_test - vegan::procrustes(res1$x[, 1:5], res2$x[, 1:5]) summary(proc_test)如果Procrustes模拟的平方残差过高需要对上游数据做更严格的质量过滤或者考虑是否真存在足以支撑分群的遗传信号。经验总结与一点个人体会我在这里最后再分享一个实际工作里的经验做群体PCA时多花一点时间在数据清洗和质控上永远比在作图阶段花时间调整调色板、加椭圆线要值得多。质控做得扎实PCA结果的可解释性会好很多出现虚假阳性分群的概率也会低很多。另外一个容易被忽略的小技巧是跑完PCA后建议把前几版的质控参数和样本数量完整记录下来尤其是在论文审稿阶段被问到“为什么样本数从1万变成9000”的时候一份清晰的记录能省下大量解释时间。我习惯把每次PLINK命令、参数阈值、过滤前后的样本数和位点数都写进一个简单的流程记录文件里并把对应的图形都归档好。这个习惯帮我在后边的投稿回复里避开了很多麻烦。希望这篇从原理到实操再到排查技巧的完整记录能帮你在自己的群体数据里少走一些弯路。如果你在实践过程中遇到别的奇怪现象欢迎回来一起讨论。