• 利用協(xié)方差矩陣，特征值和特征向量將高緯變量投影到數個低維變量的過程；

• PCA分析的過程就是從千萬級別的SNP位點中提取關鍵信息，以便使用更少的變量就可以對樣本進行有效的刻畫和區(qū)分;

• 常用分析軟件有：R、ldak、GCTA、EIGENSOFT等；

• 其結果可以代替群體結構分析的結果，作為協(xié)方差矩陣運用于關聯(lián)分析。

  
  
  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))+ #置信區(qū)間
      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？

相互驗證，使用時選擇其一即可，當群體結構的最優(yōu)分群數較低，從進化樹和PCA結果看材料分化程度又比較高時，優(yōu)先選用PCA結果作為協(xié)方差矩陣參與關聯(lián)分析。

引用轉載請注明出處，如有錯誤敬請指出。