利用协方差矩阵，特征值和特征向量将高纬变量投影到数个低维变量的过程；
PCA分析的过程就是从千万级别的SNP位点中提取关键信息，以便使用更少的变量就可以对样本进行有效的刻画和区分;
常用分析软件有：R、ldak、GCTA、EIGENSOFT等；
-
其结果可以代替群体结构分析的结果，作为协方差矩阵运用于关联分析。
Wang et al., 2013, Nature Communications
1. 下载及安装
1.1 下载地址
https://cnsgenomics.com/software/gcta/#Download
1.2 安装
$ unzip gcta_1.92.0beta3.zip
# 调用
$ ./gcta64
2. 主成分计算
2.1 数据格式整理
GCTA识别的.bin, .N.bin和.id 格式文件，需要在plink中进行整理。
$ cd your/path/plink1.9
# 将ped和map格式转成bed
$ ./plink --file genotype.id --allow-extra-chr --make-bed --out genotype.id --autosome-num 27
# 将bed转成GCTA需要的 .grm.bin, .grm.N.bin和.grm.id 格式文件
$ ./plink --file genotype.id --allow-extra-chr --make-grm-bin --out genotype.id --autosome-num 27
2.2 计算
$ cd your/path/gcta_1.92.0beta3
$ ./gcta64 --grm genotype.id --pca 20 --out pca
--pca 20：保留前20个主成分
生成.eigenval和.eigenvec两个文件
.eigenval 特征值
.eigenvec 特征向量
将两个文件下载到本地
3. 在R中进行展示
3.1 整理结果文件
> setwd("D:/数据/GWAS/主成分")
> eigenvec <- read.table("pca.eigenvec", quote="\"", comment.char="")
> colnames(eigenvec)<-c("FID","sample",paste0("PC",1:20))
> write.table(eigenvec[2:ncol(eigenvec)],file = "pca.eigenvector.xls",sep = "\t",row.names = F,col.names = T,quote = F)
整理后的eigenvec
> eigenval <- read.table("pca.eigenval", quote="\"", comment.char="")
> pcs<-paste0("pc",1:nrow(eigenval))
> eigenval[nrow(eigenval),1]<-0
# 计算解释率
> percentage<-eigenval$V1/sum(eigenval$V1)*100
> eigenval_df<-as.data.frame(cbind(pcs,eigenval[,1],percentage),stringsAsFactors = F)
> names(eigenval_df)<-c("pcs","variance","proportation")
> eigenval_df$variance<-as.numeric(eigenval_df$variance)
> eigenval_df$proportation<-as.numeric(eigenval_df$proportation)
> write.table(eigenval_df,file = "pca.eigenvalue.xls",sep = "\t",quote = F,row.names = F,col.names = T)
> head(eigenval_df)
pcs variance proportation
1 pc1 8.16174 4.842549
2 pc2 6.17792 3.665503
3 pc3 3.44002 2.041043
4 pc4 3.31091 1.964439
5 pc5 2.97959 1.767860
6 pc6 2.68970 1.595861
3.2 可视化
> eigvec <-read.table("pca.eigenvector.xls",header = T)
> eigval <- read.table("pca.eigenvalue.xls",header = T)
# 整理pop文件，order为序号，group为群体结构分群结果(group分组；color每个组的颜色，pch每个组点的形状）
> popinfo <- read.csv( "pca_pop.csv")
> head(popinfo)
order exp_id vcf_id group color pch
1 1 1 1 G2 #9ACD32 16
2 2 2 2 G2 #9ACD32 16
3 3 3 3 G2 #9ACD32 16
4 4 4 4 G2 #9ACD32 16
5 5 5 5 G2 #9ACD32 16
6 6 6 6 G2 #9ACD32 16
> pop <- unique(popinfo[,4:6])
> print(pop)
group color pch
1 G2 #9ACD32 16
55 G1 #FF4500 15
62 G3 #6495ED 17
2D图 (PC1&PC2为例)
> library("ggplot2")
> group=popinfo$group
> pch=popinfo$pch
> color=popinfo$color
> pdf("pca_pc1 vs. pc2.pdf", width = 8,height = 6)
> ggplot(data = eigvec, aes(x = PC1, y = PC2, group = group)) +
geom_point(alpha = 1,col=color,pch=pch)+
stat_ellipse(geom = "polygon",alpha=0.5,level = 0.95,aes(fill=group))+ #置信区间
geom_hline(yintercept = 0, linetype = "dashed", color = "grey",size=1)+
geom_vline(xintercept = 0, linetype = "dashed", color = "grey",size=1)+
theme_bw()+
theme(panel.grid = element_blank())+
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.5))+ #边框
theme(legend.text = element_text(size = 16, face = "bold.italic"),legend.title = element_blank())+ #图例
xlab(paste0("PC1 (", round(eigval[eigval$pcs == "pc1", 3], 2), "%)"))+ ylab(paste0("PC2 (", round(eigval[eigval$pcs == "pc2", 3], 2), "%)")) +
theme(axis.title.x = element_text(face = "bold", size = 18, colour = "black"),axis.title.y = element_text(face = "bold", size = 18, colour = "black"),axis.text.x = element_text(size = 14,face = "bold", colour = "black"),axis.text.y = element_text(face = "bold", size = 14, colour = "black"))
> dev.off()
PCA_2D
3D图
> library(scatterplot3d)
> library(grDevices）
> x_lab <- paste0("PC1 (", round(eigval[eigval$pcs == "pc1", 3], 2), "%)")
> y_lab <- paste0("PC2 (", round(eigval[eigval$pcs == "pc2", 3], 2), "%)")
> z_lab <- paste0("PC3 (", round(eigval[eigval$pcs == "pc3", 3], 2), "%)")
> pdf("pca_3d.pdf",width = 8,height = 6)
> scatterplot3d(x =eigvec$PC1, y = eigvec$PC2, z=eigvec$PC3,
color = color,pch = pch,
xlab=x_lab, ylab=y_lab, zlab=z_lab,
angle=45,type = "p",
mar = c(3.5,3.5,3.5,6),
cex.symbols = 1.5, cex.axis = 1, cex.lab = 1.5,
font.axis = 2, font.lab = 2, scale.y = 1)
> legend("topright", legend = pop$group,
col = pop$color,
pch = pop$pch,
inset = -0.1, xpd = TRUE, horiz = TRUE)
> dev.off()
pca_3D
群体结构 or PCA？
相互验证，使用时选择其一即可，当群体结构的最优分群数较低，从进化树和PCA结果看材料分化程度又比较高时，优先选用PCA结果作为协方差矩阵参与关联分析。
引用转载请注明出处，如有错误敬请指出。
