R語言給PCA加個小圈圈

Hello，小夥伴們，今天這篇推文上周發過一次，但因為部分文字內容有誤沒有審查準確，現在修改後再推一次。

小夥伴們，在遇到組學實驗數據分析得時候，是少不了繪製PCA圖的，但是除了常規的PCA圖以外，往往也需要在我們的流程結果的PCA上展現組內樣品的分布範圍：

像這樣的

這樣的

其實，本質上這些圓圈是在PCA圖的基礎上，加了一個置信區間。但是我之前都是在ppt上手動加的。這樣總覺得哪裡不太對，萬一這樣和老師說，到時候發文章被審稿人問起怎麼辦呢?

所以這種分析要求，最好的方法是提交給分析。另外我們公司R語言培訓中也有專門的PCA的圖形繪製。其中也涉及了用ggplot2來繪製PCA加圈的圖。例如我們這篇文章的封面圖片就是用ggplot2畫的。但是ggplot功能強大，也就意味著參數多啊，所以還是我們的ggord簡單易上手，另外還可以加入變量的」小線條「。

所以老師如果想自己動手，豐衣足食的話。別急，我們有簡化版的PCA繪製教程。用R語言繪製完成後，再用PS和AI適當修改一下就可以大功告成了。

PCA降維分析，不僅可以展現我們的樣品重複性，還可以用來發現內在的生物學規律。所以這次我們就圍繞這兩個主題開展今天文章的內容。

首先今天我們用的R包是ggord（https://github.com/fawda123/ggord）

因為這個包並不放在CRAN中，所以直接用install.packages(「ggord」)是不可行的。

所以我們需要在github上下載。

#首先安裝devtools

install.packages('devtools')

library(devtools)

#我們再安裝ggord

install_github('fawda123/ggord')

library(ggord)

好的，接下來就要進入正題了，以我們的流程中的表達量總表為例

列名為樣品名稱，分為A,B,C,D四組，每組4個生物學重複。

行名為基因ID

數值是RPKM值

library(ggord)

library(ggplot2)#我們需要一些ggplot2功能，如果需要ggsave保存的話

getwd()#檢查一下路徑

setwd("C:/Users/woney/Desktop")#將工作目錄改成桌面

pca_data=read.table("all.txt",header=T,sep="\t",row=1)#讀取表格

pca_data=t(as.matrix(pca_data))#將表達量總表反轉，PCA函數計算的是展現的第一列的ID

pca_group=factor(c(rep('A',4),rep('B',4), rep('C', '4'), rep('D', '4')))#分組，後續可以自己設置對應的組名，rep('組名','樣品數')

pca1=prcomp(pca_data,center=TRUE,retx=T)#prcomp是R自帶的pca分析函數

p1 <- ggord(pca1, grp_in = pca_group, arrow=0, vec_ext =0,txt=NULL)

噠噠噠，出圖如下

     
然後，我們輸入?ggord就可以查看具體的參數了。根據幫助文檔的參數就可以自行修改。

我簡單列舉一下我做的參數修改。

變成多邊形

p2 <-ggord(pca1, grp_in = pca_group, arrow=0,vec_ext = 0,txt=NULL,ellipse=F,ellipse_pro=0.68,hull=T)

修改置信區間，以及將橢圓改成空心

p3<-ggord(pca1, grp_in = pca_group, arrow=0,vec_ext =0,txt=NULL,ellipse=T,ellipse_pro=0.68,poly = F)

不展示點

p4 <-ggord(pca1, grp_in = pca_group, arrow=0,vec_ext =0,txt=NULL,ellipse=F,ellipse_pro=0.68,hull=T,obslab =TRUE,repel=TRUE)

此外，ggord是基於ggplot2基礎上繪圖，那麼理論上ggplot2的圖層語法是可以用的。

比如加上標籤

library(ggplot2)

p5=p1+geom_text(aes(label=rownames(pca_data)),vjust=1.5,colour="red")

或者theme（）主題都是可以的。總之ggord設置了大量的參數，都可以修改圖形，其中作者也寫了一些例子（https://github.com/fawda123/ggord）或者參考幫助文檔中的參數（?ggord），可以幫助我們達成我們的目的。

那麼第二點，我們如果想做成有線條和箭頭來展現我們數據中的一些變量之間的關係呢。其實呢，在PCA中，最原始的就是會把其中的一些變量給展現出來，但是平時的分析中，我們都不會畫出來，大家可以想像一下，當幾萬的基因代表的線條展現出來的時候。圖形長什麼樣子。

通過某個變量所代表的線條在PC1和PC2上的投影，我們可以看出這個變量對樣本分離的貢獻度，線條越長，代表投影越大，影響越顯著，同時兩個線條之間的夾角，可以理解為兩個變量之間的相關性，夾角小於90度，可以認為兩個變量正相關，大於90度，可以認為兩個變量負相關。

 

那麼具體又該如何去畫呢，其實用ggord是很簡單去實現的。（其實ggord可以畫很多圖，有CCA，RDA，MAC，LDA等等）

 

我們以16s測序為對象，我們的項目中會有一個表格對各個樣本的α-多樣性指數值數匯總，存於 04.alpha_diversity\summary.alpha_diversity.xls

如果像展現環境因子與多樣性指數之間的關係，可以把環境的理化因子放入其中。如下表

數據整理完，就可以開始畫圖了。

library(ggord)

library(ggplot2)

getwd()#檢查一下路徑

setwd("C:/Users/woney/Desktop")#將工作目錄改成桌面

pca_data=read.table("summary.alpha_diversity.xls",header=T,sep="\t",row= 1)#讀取表格不需將表達量總表反轉，PCA函數計算的是展現的第一列的ID

pca_group=factor(c(rep('SN','3'),rep('SNP','3'), rep('SNPK','3')))#分組，後續可以自己設置對應的組名

pca1=prcomp(pca_data,center=TRUE,retx=T)#prcomp是R自帶的pca分析函數

p1 <-ggord(pca1, grp_in = pca_group, arrow=0.3, vec_ext =700,txt=5,repel=T)

 

參數也是自己可以修改，結果如下：

我們簡單的可以推測，其中PH和TN影響最大，鹽度與α-多樣性指數呈負相關，溫度呈現正相關。

今天的內容就到這裡啦~

參考文獻

https://github.com/fawda123/ggord

https://doi.org/10.1016/j.soilbio.2016.04.005