如何用 GCTA软件进行PCA分析并进行可视化

大家好，我是邓飞，今天介绍一下如何用GCTA计算PCA，并使用R语言进行可视化。

一、GCTA 做 PCA 的核心步骤

1. 数据准备与质控：需输入 PLINK 格式（.bed/.bim/.fam）的基因型数据，先过滤低质量 SNP（缺失率 > 0.05）和个体（缺失率 > 0.1），排除亲缘关系过近的个体。
2. 构建遗传相关矩阵（GRM）：用--make-grm命令，基于全基因组 SNP 计算所有个体间的遗传相似度矩阵，仅保留常染色体 SNP 避免性别偏差。
3. 执行 PCA 分析：用--pca [k]命令（k 指定需提取的主成分个数，常用 10-20），基于 GRM 进行特征值分解。
4. 结果输出与解读：生成.eigenvec（个体 - 主成分坐标）和.eigenval（各主成分方差解释占比）文件，可用于绘制 PCA 散点图或作为协变量纳入后续分析。

二、GCTA 做 PCA 的核心优势

• 遗传针对性强：直接基于遗传标记计算相似度，主成分能精准反映个体间遗传差异，可有效校正 GWAS 中的种群分层偏倚。
• 计算效率高：优化了 GRM 计算算法，支持大样本（数万个体）和高密度 SNP（百万级）数据，运行速度优于传统 PCA 工具（如 EIGENSOFT）。
• 兼容性好：与 PLINK 数据格式无缝衔接，输出结果可直接导入 R、Python 等软件进行可视化，且能与 GCTA 的其他功能（如遗传力估算）联动。

三、GCTA 做 PCA 的的步骤

使用--make-grm-alg进行设置，0位Yang的方法，1位Van的方法。

「Yang的方法：」

GRM = sum{[(xij- 2pi)*(xik- 2pi)] / [2pi(1-pi)]}/N

「Van的方法：」

GRM = sum[(xij- 2pi)(xik- 2pi)] / sum[2pi(1-pi)]

1 Yang的方法

gcta64 --bfile ../test --make-grm --make-grm-alg 0 --out kinship_yang
gcta64 --grm kinship_yang --pca 20 --out pca_re

结果生成：

pca_re.eigenval  pca_re.eigenvec  pca_re.log

2 Van的方法

这里，将0变为1.

gcta64 --bfile ../test --make-grm --make-grm-alg 1 --out kinship_yang
gcta64 --grm kinship_yang --pca 20 --out pca_re

结果生成：

pca_re.eigenval  pca_re.eigenvec  pca_re.log

PCA可视化

这里，先对数据进行处理，计算每个主成分解释百分比，以及前几个PCA的累计百分比。

pcaal = fread("pca_re.eigenval")
head(pcaal)
pcaal$por = pcaal$V1/sum(pcaal$V1)
pcaal$cumula = cumsum(pcaal$por)
pcaal$index = as.factor(1:dim(pcaal)[1])
head(pcaal)

碎石图（折线图）

这里，选择前10个主成分。

# 选择最佳的PCA个数：碎石折线图
pcaal[1:10,] %>%
  ggplot(aes(x=index,y=por, group=1))+
  geom_point(size=4)+
  geom_line()+
  labs(title="Scree plot: PCA")

碎石图（条形图）

这里，选择前10个主成分。

# 选择最佳的PCA个数：碎石条形图
pcaal[1:10,] %>%
  ggplot(aes(x=index,y=por))+
  geom_col()+
  labs(title="Scree plot: PCA on scaled data")

PCA可视化

「代码：」

ggplot(pcaec,aes(x = V3,y = V4)) + geom_point() + 
  xlab(paste0("PC1 ",round(pcaal$por[1],4),"%")) + 
  ylab(paste0("PC2 ",round(pcaal$por[2],4),"%"))

PCA分析拓展

1，PCA分析，可以根据分组，绘制置信区间

分组PCA

2，PCA分析中，可以将PCA的百分比和累计百分比绘制到一张图上面。