6. GWAS：主成分分析——GCTA

• 利用协方差矩阵，特征值和特征向量将高纬变量投影到数个低维变量的过程；

• 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 /gss1/home/fzhang/sunwei/plink1.9
# 将ped和map格式转成bed
./plink --file root12 --allow-extra-chr --make-bed --out root12 --autosome-num 27
# 将bed转成GCTA需要的 .grm.bin, .grm.N.bin和.grm.id 格式文件
./plink --file root12 --allow-extra-chr --make-grm-bin --out root.id --autosome-num 27

2.2 计算

cd /gss1/home/fzhang/sunwei/gcta_1.92.0beta3
./gcta64 --grm root.id --pca 20 --out pcaroot

--pca 20：保留前20个主成分
生成.eigenval和.eigenvec两个文件
.eigenval 特征值
.eigenvec 特征向量
改名称为eigenval.txt和eigenvec.txt，下载到本地

3. 在R中进行展示

3.1 整理结果文件

> setwd("D:/研究生/数据/GWAS/主成分")
> eigenvec <- read.table("eigenvec.txt", 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("eigenval.txt", 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结果作为协方差矩阵参与关联分析。

引用转载请注明出处，如有错误敬请指出。