ggbiplot-最好看的PCA作圖：樣品PCA散點+分組橢圓+主成分丰度和相關
寫在前面
前幾天在《巨集基因組0》微信討論群看到了有人發了一個上面連結，點開一看居然是一條命令出帥圖，真是太實用了。我立即使用本領域的OTU表上進行了測試，效果很好，現分享給大家，歡迎大家留言補充。
ggbiplot簡介
ggbiplot是一款PCA分析結果視覺化的R包工具，可以直接採用ggplot2來視覺化R中基礎函式prcomp() PCA分析的結果，並可以按分組著色 、分組新增不同大小橢圓、主成分與原始變數相關與貢獻度向量等。
An implementation of the biplot using ggplot2. The package provides two functions: ggscreeplot()
ggbiplot(). ggbiplot aims to be a drop-in replacement for the built-in R function biplot.princomp() with extended functionality for labeling groups, drawing a correlation circle, and adding Normal probability ellipsoids.
ggbiplot安裝和官方示例
R包，建議在Rstudio中使用更方便
# 安裝包，安裝過請跳過
install.packages("devtools" 
看，是不是一條命令就把prcomp()得到的主成分分析PCA的結果展示的淋漓盡致。是不是瞬間有了高分文章的B格。
主要結果說明:
- 座標軸PC1/2的數值為總體差異的解釋率；
- 圖中點代表樣品，顏色代表分組，圖例在頂部有三組；
- 橢圓代表分組按預設68%的置信區間加的核心區域，便於觀察組間是否分開；
- 箭頭代表原始變數，其中方向代表原始變數與主成分的相關性，長度代表原始資料對主成分的貢獻度。
更詳細的PCA原理、推導、圖解，請跳轉《一文看懂PCA主成分分析》，再點選閱讀原文。重點關注——PCA結果解釋部分。
OTU表實戰
本實戰，基於本公眾號之前釋出文章 《擴增子分析教程-3統計繪圖-衝擊高分文章》中提供測試資料的OTU表、實驗設計和物種註釋資訊，需要者可前往下載。
PCA分析OTU表和視覺化
# 菌群資料實戰
# 讀入實驗設計
design = read.table("design.txt", header=T, row.names= 1, sep="\t")
# 讀取OTU表
otu_table = read.delim("otu_table.txt", row.names= 1, header=T, sep="\t")
# 過濾資料並排序
idx = rownames(design) %in% colnames(otu_table)
sub_design = design[idx,]
count = otu_table[, rownames(sub_design)]
# 基於OTU表PCA分析
otu.pca <- prcomp(t(count), scale. = TRUE)
# 繪製PCA圖，並按組新增橢圓
ggbiplot(otu.pca, obs.scale = 1, var.scale = 1,
groups = sub_design$genotype, ellipse = TRUE,var.axes = F)
可以看到三個組在第一主軸上明顯分開。
展示主要差異菌與主成分的關係
# 顯著高丰度菌的影響
# 轉換原始資料為百分比
norm = t(t(count)/colSums(count,na=T)) * 100 # normalization to total 100
# 篩選mad值大於0.5的OTU
mad.5 = norm[apply(norm,1,mad)>0.5,]
# 另一種方法：按mad值排序取前6波動最大的OTUs
mad.5 = head(norm[order(apply(norm,1,mad), decreasing=T),],n=6)
# 計算PCA和菌與菌軸的相關性
otu.pca <- prcomp(t(mad.5))
ggbiplot(otu.pca, obs.scale = 1, var.scale = 1,
groups = sub_design$genotype, ellipse = TRUE,var.axes = T)
我們僅用中值絕對偏差(mad)最大的6個OTUs進行主成分分析，即可將三組樣品明顯分開。圖中向量長短代表差異貢獻，方向為與主成分的相關性。可以看到最長的向量Streptophyta與X軸近平行，表示PC1的差異主要由此菌貢獻。其它菌與其方向相反代表OTUs間可能負相關；夾角小於90%的代表兩個OTUs有正相關。
ggbiplot官網簡介
ggbiplot
An implementation of the biplot using ggplot2. The package provides two functions: ggscreeplot() and ggbiplot(). ggbiplot aims to be a drop-in replacement for the built-in R function biplot.princomp() with extended functionality for labeling groups, drawing a correlation circle, and adding Normal probability ellipsoids.
ggbiplot使用和引數
?ggbiplot
ggbiplot(pcobj, choices = 1:2, scale = 1, pc.biplot =
TRUE, obs.scale = 1 - scale, var.scale = scale, groups =
NULL, ellipse = FALSE, ellipse.prob = 0.68, labels =
NULL, labels.size = 3, alpha = 1, var.axes = TRUE, circle
= FALSE, circle.prob = 0.69, varname.size = 3,
varname.adjust = 1.5, varname.abbrev = FALSE, ...)
pcobj # prcomp()或princomp()返回結果
choices # 選擇軸，預設1：2
scale # covariance biplot (scale = 1), form biplot (scale = 0). When scale = 1, the inner product between the variables approximates the covariance and the distance between the points approximates the Mahalanobis distance.
obs.scale # 標準化觀測值
var.scale # 標準化變異
pc.biplot # 相容 biplot.princomp()
groups # 組資訊，並按組上色
ellipse # 新增組橢圓
ellipse.prob # 置信區間
labels # 向量名稱
labels.size # 名稱大小
alpha # 點透明度 (0 = TRUEransparent, 1 = opaque)
circle # 繪製相關環(only applies when prcomp was called with scale = TRUE and when var.scale = 1)
var.axes # 繪製變數線-菌相關
varname.size # 變數名大小
varname.adjust # 標籤與箭頭距離 >= 1 means farther from the arrow
varname.abbrev # 標籤是否縮寫猜你喜歡
- 一文讀懂：1微生物組 2進化樹 3預測群落功能
- 擴增子分析：1圖表解讀 2分析流程 3統計繪圖 4群落功能 5進化樹
- 科研團隊經驗：1雲筆記 2雲協作 3公眾號
- 系列教程：1Biostar 2微生物組 3巨集基因組
寫在後面
為鼓勵讀者交流、快速解決科研困難，我們建立了“巨集基因組”專業討論群，目前己有國內外100+ PI，800+ 一線科研人員加入。參與討論，獲得專業解答，歡迎分享此文至朋友圈，並掃碼加主編好友帶你入群，務必備註“姓名-單位-研究方向-職稱/年級”。技術問題尋求幫助，首先閱讀《如何優雅的提問》學習解決問題思路，仍末解決群內討論，問題不私聊，幫助同行。
學習擴增子、巨集基因組科研思路和分析實戰，關注“巨集基因組”
