R语言之ggbiplot-PCA作图:样品PCA散点+分组椭圆+主成分丰度

20次阅读
没有评论

共计 3279 个字符,预计需要花费 9 分钟才能阅读完成。

ggbiplot 简介

ggbiplot 是一款 PCA 分析结果可视化的 R 包工具,可以直接采用 ggplot2 来可视化 R 中基础函数 prcomp() PCA 分析的结果,并可以按分组着色、分组添加不同大小椭圆、主成分与原始变量相关与贡献度向量等。

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 安装和官方示例

R 包,建议在 Rstudio 中使用更方便

# 安装包,安装过请跳过

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

library(devtools) 

install_github(“vqv/ggbiplot”)
# 最简单帅气的例子

 data(wine) 

wine.pca <- prcomp(wine, scale. = TRUE) # 演示样式  

ggbiplot(wine.pca, obs.scale = 1, var.scale = 1, groups = wine.class, ellipse = TRUE, circle = TRUE) 

+ scale_color_discrete(name = ”) 

+ theme(legend.direction = ‘horizontal’, legend.position = ‘top’) # 基本样式  

plot(wine.pca$x) # 原始图,大家可以尝试画下,不忍直视

R 语言之 ggbiplot-PCA 作图:样品 PCA 散点 + 分组椭圆 + 主成分丰度

看,是不是一条命令就把 prcomp()得到的主成分分析 PCA 的结果展示的淋漓尽致。是不是瞬间有了高分文章的 B 格。

主要结果说明:

  • 坐标轴 PC1/ 2 的数值为总体差异的解释率;
  • 图中点代表样品,颜色代表分组,图例在顶部有三组;
  • 椭圆代表分组按默认 68% 的置信区间加的核心区域,便于观察组间是否分开;
  • 箭头代表原始变量,其中方向代表原始变量与主成分的相关性,长度代表原始数据对主成分的贡献度。

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)

R 语言之 ggbiplot-PCA 作图:样品 PCA 散点 + 分组椭圆 + 主成分丰度

可以看到三个组在第一主轴上明显分开。

展示主要差异菌与主成分的关系

# 显著高丰度菌的影响  

# 转换原始数据为百分比  

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)

R 语言之 ggbiplot-PCA 作图:样品 PCA 散点 + 分组椭圆 + 主成分丰度

我们仅用中值绝对偏差 (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 # 标签是否缩写

正文完
 0