写在前面

https://github.com/vqv/ggbiplot/blob/master/README.md

前几天在《宏基因组0》微信讨论群看到了有人发了一个上面链接,点开一看居然是一条命令出帅图,真是太实用了。我立即使用本领域的OTU表上进行了测试,效果很好,现分享给大家,欢迎大家留言补充。

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) # 原始图,大家可以尝试画下,不忍直视

看,是不是一条命令就把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使用和参数

?ggbiplotggbiplot(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图表规范 2DNA提取 3 实验vs分析
  • 必备技能:1提问 2搜索 3Endnote
  • 文献阅读 1热心肠 2SemanticScholar 3geenmedical
  • 扩增子分析:1图表解读 2分析流程 3统计绘图 4群落功能 5进化树
  • 科研团队经验:1云笔记 2云协作 3公众号
  • 系列教程:1Biostar 2微生物组 3宏基因组
  • 生物科普 1肠道细菌 2生命大跃进 3细胞的暗战 4人体奥秘

写在后面

为鼓励读者交流、快速解决科研困难,我们建立了“宏基因组”专业讨论群,目前己有国内外100+ PI,800+ 一线科研人员加入。参与讨论,获得专业解答,欢迎分享此文至朋友圈,并扫码加主编好友带你入群,务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助,首先阅读《如何优雅的提问》学习解决问题思路,仍末解决群内讨论,问题不私聊,帮助同行。

学习扩增子、宏基因组科研思路和分析实战,关注“宏基因组”

点击阅读原文,跳转最新文章目录阅读
https://mp.weixin.qq.com/s/5jQspEvH5_4Xmart22gjMA

ggbiplot-最好看的PCA作图:样品PCA散点+分组椭圆+主成分丰度和相关相关推荐

  1. ggbiplot-最好看的PCA作图:样品PCA散点+分组椭圆+变量贡献与相关

    ggbiplot简介 ggbiplot是一款PCA分析结果可视化的R包工具,可以直接采用ggplot2来可视化R中基础函数prcomp() PCA分析的结果,并可以按分组着色 .分组添加不同大小椭圆. ...

  2. PCA计算原特征(指标)对主成分的贡献量/权重

    1. 用PCA反推原始特征对主成分的贡献量或权重 使用python的sklearn包中的pca函数 # -*- coding: utf-8 -*- import os import numpy as ...

  3. Python机器学习:PCA与梯度上升03求数据的主成分PCA

    梯度上升法解决最优化问题 import numpy as np import matplotlib.pyplot as plt X = np.empty((100,2)) X[:,0] = np.ra ...

  4. matlab mysvd代码解释,关于使用SVD进行PCA主成分提取的代码问题!也是必须涉及到原理的!...

    如标题所示,目的是将矩阵进行PCA分析最终得到降维后的新主成分,对于特征值特征向量的提取方法是通过svd.代码和函数名称定义如下/如图所示 function [Xm,U,L]=pca(X,K); % ...

  5. PCA主成分分析 提取主成分,过滤噪音

    前一篇提到的人脸识别中,我们在使用SVM支持向量机做人脸分类之前使用到PCA提取人脸数据中的主要成分,降低计算的维度,那么具体PCA是如何提取的呢?下文了解一下. PCA is a method to ...

  6. matlab pca svd,关于使用SVD进行PCA主成分提取的代码问题!也是必须涉及到原理的!...

    如标题所示,目的是将矩阵进行PCA分析最终得到降维后的新主成分,对于特征值特征向量的提取方法是通过svd.代码和函数名称定义如下/如图所示 function [Xm,U,L]=pca(X,K); % ...

  7. R语言主成分PCA、因子分析、聚类对地区经济研究分析重庆市经济指标

    全文下载链接:http://tecdat.cn/?p=27515 建立重庆市经济指标发展体系,以重庆市一小时经济圈作为样本,运用因子分析方法进行实证分析,在借鉴了相关评价理论和评价方法的基础上,本文提 ...

  8. 主成分分析旋转矩阵MATLAB实现,R语言高维数据的主成分pca、t-SNE算法降维与可视化分析案例报告...

    维度降低有两个主要用例:数据探索和机器学习.它对于数据探索很有用,因为维数减少到几个维度(例如2或3维)允许可视化样本.然后可以使用这种可视化来从数据获得见解(例如,检测聚类并识别异常值).对于机器学 ...

  9. R语言辅导高维数据的主成分pca、 t-SNE算法降维与可视化分析案例报告

    降低维度有两个主要用例:数据探索和机器学习.它对于数据探索很有用,因为维数减少到几个维度(例如2或3维)允许可视化样本.然后可以使用这种可视化来从数据获得见解(例如,检测聚类并识别异常值).对于机器学 ...

最新文章

  1. matrix_multiply代码解析
  2. HDU-2044-一只小蜜蜂
  3. JS代码的window.location属性详解
  4. mysql 显示前三项_详解MySQL三项实用开发知识
  5. linux基础分支,Linux基础--/etc/shadow中字段的分支和操作
  6. tensorflow2 tensorboard可视化使用
  7. java mvc引擎_SpringMvc+JavaConfig+Idea 搭建项目
  8. VB用API实现各种对话框(总结)(转载)
  9. 在Mysql中遇到关于区间范围内的索引优化
  10. 各省简称 拼音 缩写_中国各省、直辖市、自治区名称汉语拼音字母缩
  11. SnowNLP自然语言处理模块具体用法
  12. java如何动态添加数组数据_Java动态数组添加数据的方法与应用示例
  13. html表格可视化设计器,基于vue-element-ui的一款表格设计器table-making
  14. 空间命名的定义及使用:using namespace std 的用法详解
  15. java get请求 数组_GET方式请求的url参数如果是数组,该形式/base/get?foo[]=barfoo[]=baz'报错......
  16. Android 编译错误:CreateProcess error=206, 文件名或扩展名太长。
  17. 程序员的情感修养 —— 专访“非诚勿扰”牵手成功男嘉宾程序员石川
  18. Win10 批量修改文件名
  19. 无人驾驶8: 粒子滤波定位(优达学城项目)
  20. CycleGAN训练自己的数据集

热门文章

  1. 如何使用消息队列解决分布式事物?
  2. 系列文章|OKR与敏捷(二):实现全栈敏捷
  3. 网上整理的对于Rest和Restful api的理解
  4. Eclipse注释模板 注释快捷键
  5. python字典存储省份与城市_在Python中存储字典路径
  6. C++中的变量作用域介绍
  7. winform程序打包EXE三种方式
  8. oracle驱动程序包的安装失败,Maven 、oracle的jdbc的jar包下载失败
  9. 视觉传感器:3D感知算法
  10. 如何评价CVPR 2021的论文接收结果?