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

前言

主成分分析(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

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 199,175评论 5 466
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 83,674评论 2 376
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 146,151评论 0 328
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 53,597评论 1 269
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 62,505评论 5 359
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 47,969评论 1 275
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,455评论 3 390
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,118评论 0 254
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,227评论 1 294
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,213评论 2 317
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,214评论 1 328
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 32,928评论 3 316
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,512评论 3 302
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,616评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,848评论 1 255
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,228评论 2 344
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 41,772评论 2 339

推荐阅读更多精彩内容