R语言中利用prcomp函数 对基因型数据进行PCA分析

 

001、plink软件生成raw文件

root@PC1:/home/test# ls
outcome.map  outcome.ped        ## maf过滤，并将基因型数据转换为0、1、2的形式
root@PC1:/home/test# plink --file outcome --sheep --maf 0.01 --recode A 1> /dev/null; rm *.log *.nosex
root@PC1:/home/test# ls
outcome.map  outcome.ped  plink.raw
root@PC1:/home/test# sed 1d plink.raw | cut -d " " -f 7- > use.raw  ## 修改数据格式
root@PC1:/home/test# ls
outcome.map  outcome.ped  plink.raw  use.raw

 

002、R读入，并使用prcomp函数计算， 然后绘图

dat <- read.table("use.raw")                       ## 读入数据
pca <- prcomp(dat,center = TRUE,scale. = TRUE)     ## prcomp函数进行pca分析
plot(pca$x[,1], pca$x[,2], pch = 19)               ## 绘图

 

 

002、plink软件pca验证

root@PC1:/home/test# ls
outcome.map  outcome.ped                           ## plink软件计算pca
root@PC1:/home/test# plink --file outcome --sheep --pca 3 header tabs --out pca 1> /dev/null; rm *.log *.nosex
root@PC1:/home/test# ls                            ## 结果文件
outcome.map  outcome.ped  pca.eigenval  pca.eigenvec

 

R绘图：

par(mfrow = c(1,2))                                  ##分割画布
dat <- read.table("use.raw")
pca <- prcomp(dat,center = TRUE,scale. = TRUE)
plot(pca$x[,1], pca$x[,2], pch = 19, main = "R")     ## R中生成的结果

dat2 <- read.table("pca.eigenvec", header = T)
plot(dat2$PC1, dat2$PC2, pch = 19, main = "plink")   ## plink中生成的结果

 

 

为什么结果时反着呢？