各个主成分的可解释方差,然后选取 累计可解释方差在 80-90% 的前 k 个主成分。但是图里都是零点零零几?
基于可解释方差(推荐)
除了基于 Tracy–Widom statistics 检验主成分的显著性外,还可以根据每个主成分的可解释方差计算。一般,选取累计解释 80-90% 的前 k 个主成分就足够。Plink、GCTA 等工具不能输出各个主成分的可解释方差,要这个信息的话可以用 vcfR、SNPRelate、bigsnpr、pcadapt 等 R 包。
SNPRelate 的并行计算速度比较快,以它为例,计算 PCA 并且得到可解释方差:
# from shiyanhe and zhaozhuji.net
# 从 Bioconductor 安装 SNPRelate 包和它依赖的 gdsfmt 包
if (!requireNamespace("BiocManager", quietly=TRUE))
install.packages("BiocManager")
BiocManager::install("gdsfmt")
BiocManager::install("SNPRelate")
# 加载 gdsfmt 和 SNPRelate 包
library(gdsfmt)
library(SNPRelate)
# 输入 PLINK 文件路径
bed.fn <- "/your_folder/your_plink_file.bed"
fam.fn <- "/your_folder/your_plink_file.fam"
bim.fn <- "/your_folder/your_plink_file.bim"
# 将 PLINK 文件转为 GDS 文件
snpgdsBED2GDS(bed.fn, fam.fn, bim.fn, "test.gds")
# 读取 GDS 文件
genofile <- snpgdsOpen("test.gds")
# 根据 LD 过滤 SNPs,阈值根据需要设定
set.seed(1000)
snpset <- snpgdsLDpruning(genofile, ld.threshold=0.2)
# 选择 SNP pruning 后要保留的 SNP
snpset.id <- unlist(unname(snpset))
# 计算 PCA,num.thread 是并行的线程数
pca <- snpgdsPCA(genofile, snp.id=snpset.id, num.thread=10)
# 以百分比形式输出 variance proportion
print(pca$varprop*100)
# 绘制前 30 个主成分的碎石图
# from shiyanhe and zhaozhuji.net
library(ggplot2)
K= 30
qplot(x = 1:K, y = (pca$varprop[1:K]), col = "red", xlab = "PC", ylab = "Proportion of explained variance") +
geom_line() + guides(colour = FALSE) +
ggtitle(paste("Scree Plot - K =", K))
各个主成分的可解释方差,然后选取 P 值显著且累计可解释方差在 80-90% 的前 k 个主成分。但是图里都是零点零零几?这到底怎么看?有没有人解答一下