image.png

公眾號(hào)后臺(tái)有讀者留言問這個(gè)圖的實(shí)現(xiàn)辦法，這個(gè)圖相比于普通的PCA散點(diǎn)圖是多了一個(gè)垂直和水平的誤差線，這個(gè)如何實(shí)現(xiàn)之前還沒有嘗試過，所以查了查資料，找到了一個(gè)參考鏈接

https://cran.r-project.org/web/packages/SIBER/vignettes/Plot-SIA-ggplot2.html

按照這個(gè)參考鏈接的代碼 我們?cè)囋?。用之前提到的小麥種子數(shù)據(jù)做示例數(shù)據(jù)

加載需要用到的包

library(ggplot2)
library(ggforce)
library(tidyverse)

讀取數(shù)據(jù)

df<-read.csv("Seed_Data.csv")
df$target<-paste0("cultivar",df$target)

主成分分析

df.pca<-prcomp(df[,1:7],scale. = T)

提取主成分分析的作圖數(shù)據(jù)

df.pca$x %>% as.data.frame() %>% 
  mutate(group=df$target) -> pca.result

計(jì)算PCA結(jié)果的平均值和標(biāo)準(zhǔn)差

pca.result %>% 
  group_by(group) %>% 
  summarise(pc1m=mean(PC1),
            pc1sd=sd(PC1),
            pc2m=mean(PC2),
            pc2sd=mean(PC2)) -> pca.result.a

初步的草圖

ggplot()+
  geom_errorbar(data=pca.result.a,
                aes(x=pc1m,
                    ymin=pc2m-1.96*pc2sd,
                    ymax=pc2m+1.96*pc2sd,
                    color=group),
                width=0.2)+
  geom_errorbarh(data = pca.result.a,
                 aes(y=pc2m,
                     xmin=pc1m-1.96*pc1sd,
                     xmax=pc1m+1.96*pc1sd,
                     color=group))+
  stat_ellipse(data=pca.result,
               geom="polygon",
               aes(x=PC1,
                   y=PC2,
                   color=group,
                   fill=group),
               alpha=0.2) -> p
p

image.png

從圖中提取作圖數(shù)據(jù)

ggplot_build(p)$data[[1]] %>% 
  select(colour,x,ymin,ymax) %>% 
  pivot_longer(cols = c(ymin,ymax)) %>% 
  rename("group"="colour",
         "x" ="x",
         "y"="value") %>% 
  select(x,y,group) -> tmp.1

ggplot_build(p)$data[[2]] %>% 
  select(colour,y,xmin,xmax) %>% 
  pivot_longer(cols = c(xmin,xmax))%>% 
  rename("group"="colour",
         "y" ="y",
         "x"="value") %>% 
  select(x,y,group) -> tmp.2

tmp<-rbind(tmp.1,tmp.2)

最終作圖

ggplot()+
  geom_mark_ellipse(data=tmp,
                    aes(x=x,y=y,
                        fill=group),
                    expand = unit(0,"mm"),
                    color="white")+
  geom_errorbar(data=pca.result.a,
                aes(x=pc1m,
                    ymin=pc2m-1.96*pc2sd,
                    ymax=pc2m+1.96*pc2sd),
                width=0.2)+
  geom_errorbarh(data = pca.result.a,
                 aes(y=pc2m,
                     xmin=pc1m-1.96*pc1sd,
                     xmax=pc1m+1.96*pc1sd))+
  theme_bw()+
  labs(x="PC1 (42.8%)",y="PC2(37.9%)")+
  geom_point(data=pca.result.a,
             aes(x=pc1m,y=pc2m),
             size=3)

image.png

歡迎大家關(guān)注我的公眾號(hào)

小明的數(shù)據(jù)分析筆記本

小明的數(shù)據(jù)分析筆記本 公眾號(hào) 主要分享：1、R語言和python做數(shù)據(jù)分析和數(shù)據(jù)可視化的簡(jiǎn)單小例子；2、園藝植物相關(guān)轉(zhuǎn)錄組學(xué)、基因組學(xué)、群體遺傳學(xué)文獻(xiàn)閱讀筆記；3、生物信息學(xué)入門學(xué)習(xí)資料及自己的學(xué)習(xí)筆記！