16是分析之主成分分析實戰部分

下面我們來做PCA分析，來做出很好看的圖片

之前在插幾句：非約束排序，不做統計檢驗，不能做！

#安裝包ggbiplot包

#install.packages("devtools",repo="http://cran.us.r-project.org")

#library(devtools)

#install_github("vqv/ggbiplot")

library(ggbiplot)

#修改工作目錄

setwd("E:/Shared_Folder/HG_usearch_HG")

#導入需要的作圖文件

design =read.table("map_lxdjhg.txt", header=T, row.names= 1,sep=" ")

head(design)

#讀取OTU表並進行數據整理

otu_table =read.delim("otu_table8825.txt", row.names= 1, header=T, sep=" ")

head(otu_table)

#過濾數據並排序

sub_design =design[rownames(design) %in% colnames(otu_table),]

count = otu_table[,rownames(sub_design)]

以上我們導入數據並整理數據，接下來我們開始分析PCA：

函數prcomp默認選項center=T(當然中心化是一定要做的)，但是不做標準化，我們調整參數scale. = TRUE，將方差縮放為單位方差最為合理

#進行PCA分析

otu.pca

這裡我們使用的ggbiplot，obs.scale = 1, var.scale = 1,這兩個參數表示我們使用標尺出圖，ellipse= TRUE我們添加橢圓置信區間默認是0.95置信區間，可以修改嗎？對不起，不行，寫死了，已經在函數裡面；

p=ggbiplot(otu.pca,obs.scale = 1, var.scale = 1,

groups = sub_design$SampleType,ellipse = TRUE,var.axes = T,circle = TRUE)+

scale_color_discrete(name = "") +

theme(legend.direction = "horizontal",legend.position = "top")

#

p

下面這幅圖算是將PCA的所有要素展示的非常好了，這裡有一個新的東西，之前也沒說過，就是：

平衡貢獻圓（circle of equilibrium contribution），它的半徑表示變數的向量長度對排序的平均貢獻率，可以這麼講，如果變數的箭頭線長度長於圓的半徑，則它對排序空間的貢獻大於變數的平均貢獻率；

我們使用擴增子16s，未免otu數目過多，什麼都看不清楚，自然需要我們將變數減少

#因此我們首先關注一下高丰度菌的影響

#轉換原始數據為百分比

norm =t(t(count)/colSums(count,na=T)) * 100 # normalization to total 100

head(norm)

#篩選mad值大於0.5的OTU

norm2 =norm[apply(norm,1,mean)>1,]

dim(norm2)

#重新計算PCA

otu.pca

p

groups = sub_design$SampleType,ellipse = TRUE,var.axes = T,circle = TRUE)

p

#挑選0.001丰度以上的菌

norm3 =norm[apply(norm,1,mean)>0.1,]

otu.pca1

p

groups = sub_design$SampleType,ellipse = TRUE,var.axes = T,circle = TRUE)

p

#當然我們僅僅展示樣品點坐標也是可以的只需修改var.axes = F

p

groups = sub_design$SampleType, ellipse= TRUE,var.axes = F)

#可以修改主題，我修改了一份自己覺得美觀的主題（我沒有什麼審美），大家有需要可以借鑒

p+theme_bw()+scale_colour_manual(values= mi)+labs(x=paste("PC 1 (", 61.7, "%)", sep=""),

y=paste("PC 2 (", 14.5,"%)", sep=""))+

labs(title = "PCA")+

theme(plot.title = element_text(hjust =0.5))+

guides(color=guide_legend(title =NULL),shape=guide_legend(title = NULL))+

geom_hline(yintercept=0) +geom_vline(xintercept=0)

ggbiplot是一款PCA分析結果可視化的R包工具，但是畢竟也有很多局限在裡面的，比如：

我需要修改樣品表示方式，我需要化多個置信區間，可能就不可以實現，所以將PCA分析結果提取出來，使用ggplot作圖系統來作圖會更加自由一些；

下面我們開始提取主成分分析結果

#現在我們來提取作圖坐標

#首先是樣品排序軸坐標

predict(otu.pca)

#也可以使用下面這條命令提取

yangpin

yangpin$SampleType=sub_design$SampleType

#提取荷載坐標

bianliang

#提取特徵根,這裡提供的並不是特徵值而是標準差，需要求其平方才是特徵值

eig=otu.pca$sdev

eig=eig*eig

#在這裡我設定了隨機種子，方便兩種形式圖形比較

set.seed(10)

ggplot(data=yangpin,aes(x=yangpin$PC1,y=yangpin$PC2,group=SampleType,color=SampleType))+geom_point()+

stat_ellipse(type = "t", linetype =2)+

labs(x=paste("PC 1 (", format(100 *eig[1] / sum(eig), digits=4), "%)", sep=""),

y=paste("PC 2 (", format(100 *eig[2] / sum(eig), digits=4), "%)", sep=""))+

labs(title = "PCA")

#########我們將載荷也畫上去

p+stat_ellipse(data=yangpin,aes(group=SampleType,color=SampleType),level = 0.95,type ="t", linetype = 2)+

geom_point(data=bianliang,aes(x=PC1,y=PC2),size=2)+

geom_text(data=bianliang,aes(x=PC1,y=PC2,label=spece),color="blue",hjust=1.2,vjust=-1.3,size=3)

這裡我們沒有像ggbioplot一樣畫一個箭頭，因為我沒有找到簡單的方法完成這件事，總會這樣的，如果大家有好的辦法，可以留言給我最好了，我將不勝感激。

最後還是那句話：學習永無止境，分享永不停歇！

喜歡這篇文章嗎？立刻分享出去讓更多人知道吧！