草庐IT

R数据可视化: PCA和PCoA图, 2D和3D

斗战胜佛oh 2023-03-28 原文

前言

主成分分析(Principal Components Analysis,PCA),也称主分量分析或主成分回归分析法,是一种无监督的数据降维方法。PCA通过线性变换将原始数据变换为一组各维度线性无关的表示,可用于提取数据的主要特征分量,常用于高维数据的降维。这种降维的思想首先减少数据集的维数,同时还保持数据集的对方差贡献最大的特征,最终使数据直观呈现在二维坐标系。

数据降维展示

直观上,第一主成分轴 优于 第二主成分轴,具有最大可分性。
主坐标分析(Principal Coordinates Analysis,PCoA),即经典多维标度(Classical multidimensional scaling),用于研究数据间的相似性。

【二者差异】

主成分分析(Principal components analysis,PCA)是一种统计分析、简化数据集的方法。它利用正交变换来对一系列可能相关的变量的观测值进行线性变换,从而投影为一系列线性不相关变量的值,这些不相关变量称为主成分(Principal Components)。具体地,主成分可以看做一个线性方程,其包含一系列线性系数来指示投影方向(如图)。PCA对原始数据的正则化或预处理敏感(相对缩放)。PCA是最简单的以特征量分析多元统计分布的方法。通常情况下,这种运算可以被看作是揭露数据的内部结构,从而更好的解释数据的变量的方法。


PCA示意图

主坐标分析(Principal Coordinates Analysis,PCoA),即经典多维标度(Classical multidimensional scaling),用于研究数据间的相似性。PCoA与PCA都是降低数据维度的方法,但是差异在在于PCA是基于原始矩阵,而PCoA是基于通过原始矩阵计算出的距离矩阵。因此,PCA是尽力保留数据中的变异让点的位置不改动,而PCoA是尽力保证原本的距离关系不发生改变,也就是使得原始数据间点的距离与投影中即结果中各点之间的距离尽可能相关(如图)。


PCoA示意图

如何进行PCA和PCoA分析

R中有很多包都提供了PCA和PCoA,比如常用的ade4包。本文将基于该包进行PCA和PCoA的分析,数据是自带的deug,该数据提供了104个学生9门课程的成绩(见截图)和综合评定。综合评定有以下几个等级:A+,A,B,B-,C-,D。
让我们通过PCA和PCoA来看一看这样的综合评定是否合理,是否确实依据这9门课把这104个学生合理分配到不同组(每个等级一个组)。

Deug的9门课

绘图

(1)PCA分析及作图

前文已经介绍了PCA是基于原始数据,所以直接进行PCA分析即可。相信大家都比较熟悉散点图的绘制方法,这里不再细讲,PCA分析完毕后我们直接作图展示结果。

library(ade4)
library(ggplot2)
library(RColorBrewer)
data(deug)

#PCA分析
pca<- dudi.pca(deug$tab, scal = FALSE, center = deug$cent, scan = FALSE)

#坐标轴解释量(前两轴)
pca_eig <- (pca$eig)[1:2] / sum(pca$eig)

#提取样本点坐标(前两轴)
sample_site <- data.frame({pca$li})[1:2]
sample_site$names <- rownames(sample_site)
names(sample_site)[1:2] <- c('PCA1', 'PCA2')

#以最终成绩作为分组
sample_site$level<-factor(deug$result,levels=c('A+','A','B','B-','C-','D'))

library(ggplot2)

pca_plot <- ggplot(sample_site, aes(PCA1, PCA2,color=level)) +
  theme_classic()+#去掉背景框
  geom_vline(xintercept = 0, color = 'gray', size = 0.4) + 
  geom_hline(yintercept = 0, color = 'gray', size = 0.4) +
  geom_point(size = 1.5)+  #可在这里修改点的透明度、大小
  scale_color_manual(values = brewer.pal(6,"Set2")) + #可在这里修改点的颜色
  theme(panel.grid = element_line(color = 'gray', linetype = 2, size = 0.1), 
        panel.background = element_rect(color = 'black', fill = 'transparent'), 
        legend.title=element_blank()
  )+
  labs(x = paste('PCA1: ', round(100 * pca_eig[1], 2), '%'), y = paste('PCA2: ', round(100 * pca_eig[2], 2), '%')) 

pca_plot
PCA plot

整体看起来还不错,就是B-和C-的学生似乎难以区分。

(2)PCoA分析及作图

library(ade4)
library(ggplot2)
library(RColorBrewer)
library(vegan)#用于计算距离
data(deug)
tab<-deug$tab
tab.dist<-vegdist(tab,method='euclidean')#基于euclidean距离
pcoa<- dudi.pco(tab.dist, scan = FALSE,nf=3)

#坐标轴解释量(前两轴)
pcoa_eig <- (pcoa$eig)[1:2] / sum(pcoa$eig)

#提取样本点坐标(前两轴)
sample_site <- data.frame({pcoa$li})[1:2]
sample_site$names <- rownames(sample_site)
names(sample_site)[1:2] <- c('PCoA1', 'PCoA2')

#以最终成绩作为分组
sample_site$level<-factor(deug$result,levels=c('A+','A','B','B-','C-','D'))

library(ggplot2)

pcoa_plot <- ggplot(sample_site, aes(PCoA1, PCoA2,color=level)) +
  theme_classic()+#去掉背景框
  geom_vline(xintercept = 0, color = 'gray', size = 0.4) + 
  geom_hline(yintercept = 0, color = 'gray', size = 0.4) +
  geom_point(size = 1.5)+  #可在这里修改点的透明度、大小
  scale_color_manual(values = brewer.pal(6,"Set2")) + #可在这里修改点的颜色
  theme(panel.grid = element_line(color = 'gray', linetype = 2, size = 0.1), 
        panel.background = element_rect(color = 'black', fill = 'transparent'), 
        legend.title=element_blank()
  )+
  labs(x = paste('PCoA1: ', round(100 * pcoa_eig[1], 2), '%'), y = paste('PCoA2: ', round(100 * pcoa_eig[2], 2), '%')) 

pcoa_plot
PCoA plot

有时候PCA和PCoA的结果差不多,有时候某种方法能够把样本有效分开而另一种可能效果不佳,这些都要看样本数据的特性。

(3)3D图

除转录组研究以外,在16S微生物的研究中我们会根据物种丰度的文件对数据进行PCA或者PCoA分析,也是我们所说的β多样性分析。根据PCA或者PCoA的结果看感染组和对照组能否分开,以了解微生物组的总体变化情况。

β多样性分析的概念
Beta多样性指的是样本间多样性。在肠道菌群分析中,Beta多样性是衡量个体间微生物组成相似性的一个指标。通过计算样本间距离可以获得β多样性计算矩阵,后续一般会利用PCoA、进化树聚类等分析对此数值关系进行图形展示。主要基于OTU的群落比较方法,有欧式距离、bray curtis距离、Jaccard 距离,这些方法优势在于算法简单,考虑物种丰度(有无)和均度(相对丰度),但其没有考虑OTUs之间的进化关系,认为OTU之间不存在进化上的联系,每个OTU间的关系平等。另一种算法Unifrac距离法,是根据系统发生树进行比较,并根据16s的序列信息对OTU进行进化树分类, 一般有加权和非加权分析。

QIIME2中重要的Beta多样性指数:

Jaccard距离:群落差异的定性度量,即只考虑种类,不考虑丰度。

Bray-Curtis距离:群落差异的定量度量,较常用。

Unweighted UniFrac距离:包含特征之间的系统发育关系的群落差异定性度量。

Weighted UniFrac距离:包含特征之间的系统发育关系的群落差异定量度量。

R绘制三维PCoA图

解压缩通过qiime2输出的 .qza文件,获得绘图的matrix和pcoa结果文件

qiime tools export \
--input-path bray_curtis_distance_matrix.qza \
--output-path bray_curtis
解压后得到的文件

将pcoa结果整理成下表,保存为 ***_site.txt


整理后文件
#载入数据
setwd("D:/your_work_path/")
pca_site <- read.delim('***_site.txt', sep = '\t', stringsAsFactors = FALSE)

#scatterplot3d包
library(scatterplot3d) 
#调整角度,保存
scatterplot3d(pca_site$PC1, pca_site$PC2, pca_site$PC3, color = rep(c('#f94144', '#f9c74f', '#5390d9'), c(8, 6, 6)), angle=40,cex.symbols=1.5,cex.axis=0.8, pch = rep(rep(c(15,16,15,16,15,16), c(4, 4, 3, 3, 3, 3))),xlab = paste('PCoA1: 15%'), ylab = paste('PCoA2: 12%'), zlab = paste('PCoA3: 8%'))
3D PCoA

注意没有legend,需要AI加入。
后期需要继续摸索,其实可以加legend的,只是目前自己的技术做不到。。。

PCA思想解析:
https://www.jianshu.com/p/09bae5cbdc53

有关R数据可视化: PCA和PCoA图, 2D和3D的更多相关文章

  1. ruby - 解析 RDFa、微数据等的最佳方式是什么,使用统一的模式/词汇(例如 schema.org)存储和显示信息 - 2

    我主要使用Ruby来执行此操作,但到目前为止我的攻击计划如下:使用gemsrdf、rdf-rdfa和rdf-microdata或mida来解析给定任何URI的数据。我认为最好映射到像schema.org这样的统一模式,例如使用这个yaml文件,它试图描述数据词汇表和opengraph到schema.org之间的转换:#SchemaXtoschema.orgconversion#data-vocabularyDV:name:namestreet-address:streetAddressregion:addressRegionlocality:addressLocalityphoto:i

  2. ruby - Ruby 有 `Pair` 数据类型吗? - 2

    有时我需要处理键/值数据。我不喜欢使用数组,因为它们在大小上没有限制(很容易不小心添加超过2个项目,而且您最终需要稍后验证大小)。此外,0和1的索引变成了魔数(MagicNumber),并且在传达含义方面做得很差(“当我说0时,我的意思是head...”)。散列也不合适,因为可能会不小心添加额外的条目。我写了下面的类来解决这个问题:classPairattr_accessor:head,:taildefinitialize(h,t)@head,@tail=h,tendend它工作得很好并且解决了问题,但我很想知道:Ruby标准库是否已经带有这样一个类? 最佳

  3. ruby - Ruby 中的波形可视化 - 2

    我即将开始一个将录制和编辑音频文件的项目,我正在寻找一个好的库(最好是Ruby,但会考虑Java或.NET以外的任何库)以进行实时可视化波形。有人知道我应该从哪里开始搜索吗? 最佳答案 要流入浏览器的数据量很大。Flash或Flex图表可能是唯一能提高内存效率的解决方案。Javascript图表往往会因大型数据集而崩溃。 关于ruby-Ruby中的波形可视化,我们在StackOverflow上找到一个类似的问题: https://stackoverflow.c

  4. ruby - 我如何添加二进制数据来遏制 POST - 2

    我正在尝试使用Curbgem执行以下POST以解析云curl-XPOST\-H"X-Parse-Application-Id:PARSE_APP_ID"\-H"X-Parse-REST-API-Key:PARSE_API_KEY"\-H"Content-Type:image/jpeg"\--data-binary'@myPicture.jpg'\https://api.parse.com/1/files/pic.jpg用这个:curl=Curl::Easy.new("https://api.parse.com/1/files/lion.jpg")curl.multipart_form_

  5. 世界前沿3D开发引擎HOOPS全面讲解——集3D数据读取、3D图形渲染、3D数据发布于一体的全新3D应用开发工具 - 2

    无论您是想搭建桌面端、WEB端或者移动端APP应用,HOOPSPlatform组件都可以为您提供弹性的3D集成架构,同时,由工业领域3D技术专家组成的HOOPS技术团队也能为您提供技术支持服务。如果您的客户期望有一种在多个平台(桌面/WEB/APP,而且某些客户端是“瘦”客户端)快速、方便地将数据接入到3D应用系统的解决方案,并且当访问数据时,在各个平台上的性能和用户体验保持一致,HOOPSPlatform将帮助您完成。利用HOOPSPlatform,您可以开发在任何环境下的3D基础应用架构。HOOPSPlatform可以帮您打造3D创新型产品,HOOPSSDK包含的技术有:快速且准确的CAD

  6. FOHEART H1数据手套驱动Optitrack光学动捕双手运动(Unity3D) - 2

    本教程将在Unity3D中混合Optitrack与数据手套的数据流,在人体运动的基础上,添加双手手指部分的运动。双手手背的角度仍由Optitrack提供,数据手套提供双手手指的角度。 01  客户端软件分别安装MotiveBody与MotionVenus并校准人体与数据手套。MotiveBodyMotionVenus数据手套使用、校准流程参照:https://gitee.com/foheart_1/foheart-h1-data-summary.git02  数据转发打开MotiveBody软件的Streaming,开始向Unity3D广播数据;MotionVenus中设置->选项选择Unit

  7. 使用canal同步MySQL数据到ES - 2

    文章目录一、概述简介原理模块二、配置Mysql使用版本环境要求1.操作系统2.mysql要求三、配置canal-server离线下载在线下载上传解压修改配置单机配置集群配置分库分表配置1.修改全局配置2.实例配置垂直分库水平分库3.修改group-instance.xml4.启动监听四、配置canal-adapter1修改启动配置2配置映射文件3启动ES数据同步查询所有订阅同步数据同步开关启动4.验证五、配置canal-admin一、概述简介canal是Alibaba旗下的一款开源项目,Java开发。基于数据库增量日志解析,提供增量数据订阅&消费。Git地址:https://github.co

  8. Unity 3D 制作开关门动画,旋转门制作,推拉门制作,门把手动画制作 - 2

    Unity自动旋转动画1.开门需要门把手先动,门再动2.关门需要门先动,门把手再动3.中途播放过程中不可以再次进行操作觉得太复杂?查看我的文章开关门简易进阶版效果:如果这个门可以直接打开的话,就不需要放置"门把手"如果门把手还有钥匙需要旋转,那就可以把钥匙放在门把手的"门把手",理论上是可以无限套娃的可调整参数有:角度,反向,轴向,速度运行时点击Test进行测试自己写的代码比较垃圾,命名与结构比较拉,高手轻点喷,新手有类似的需求可以拿去做参考上代码usingSystem.Collections;usingSystem.Collections.Generic;usingUnityEngine;u

  9. [Vuforia]二.3D物体识别 - 2

    之前说过10之后的版本没有3dScan了,所以还是9.8的版本或者之前更早的版本。 3d物体扫描需要先下载扫描的APK进行扫面。首先要在手机上装一个扫描程序,扫描现实中的三维物体,然后上传高通官网,在下载成UnityPackage类型让Unity能够使用这个扫描程序可以从高通官网上进行下载,是一个安卓程序。点到Tools往下滑,找到VuforiaObjectScanner下载后解压数据线连接手机,将apk文件拷入手机安装然后刚才解压文件中的Media文件夹打开,两个PDF图打印第一张A4-ObjectScanningTarget.pdf,主要是用来辅助扫描的。好了,接下来就是扫描三维物体。将瓶

  10. ruby-on-rails - 创建 ruby​​ 数据库时惰性符号绑定(bind)失败 - 2

    我正在尝试在Rails上安装ruby​​,到目前为止一切都已安装,但是当我尝试使用rakedb:create创建数据库时,我收到一个奇怪的错误:dyld:lazysymbolbindingfailed:Symbolnotfound:_mysql_get_client_infoReferencedfrom:/Library/Ruby/Gems/1.8/gems/mysql2-0.3.11/lib/mysql2/mysql2.bundleExpectedin:flatnamespacedyld:Symbolnotfound:_mysql_get_client_infoReferencedf

随机推荐