前言
主成分分析(Principal Components Analysis,PCA)端壳,也稱主分量分析或主成分回歸分析法告丢,是一種無監(jiān)督的數(shù)據(jù)降維方法。PCA通過線性變換將原始數(shù)據(jù)變換為一組各維度線性無關(guān)的表示损谦,可用于提取數(shù)據(jù)的主要特征分量岖免,常用于高維數(shù)據(jù)的降維。這種降維的思想首先減少數(shù)據(jù)集的維數(shù)照捡,同時還保持?jǐn)?shù)據(jù)集的對方差貢獻最大的特征颅湘,最終使數(shù)據(jù)直觀呈現(xiàn)在二維坐標(biāo)系。
直觀上栗精,第一主成分軸 優(yōu)于 第二主成分軸闯参,具有最大可分性瞻鹏。
主坐標(biāo)分析(Principal Coordinates Analysis,PCoA)鹿寨,即經(jīng)典多維標(biāo)度(Classical multidimensional scaling)新博,用于研究數(shù)據(jù)間的相似性。
【二者差異】
主成分分析(Principal components analysis脚草,PCA)是一種統(tǒng)計分析赫悄、簡化數(shù)據(jù)集的方法。它利用正交變換來對一系列可能相關(guān)的變量的觀測值進行線性變換馏慨,從而投影為一系列線性不相關(guān)變量的值埂淮,這些不相關(guān)變量稱為主成分(Principal Components)。具體地写隶,主成分可以看做一個線性方程倔撞,其包含一系列線性系數(shù)來指示投影方向(如圖)。PCA對原始數(shù)據(jù)的正則化或預(yù)處理敏感(相對縮放)樟澜。PCA是最簡單的以特征量分析多元統(tǒng)計分布的方法误窖。通常情況下,這種運算可以被看作是揭露數(shù)據(jù)的內(nèi)部結(jié)構(gòu)秩贰,從而更好的解釋數(shù)據(jù)的變量的方法。
主坐標(biāo)分析(Principal Coordinates Analysis柔吼,PCoA)毒费,即經(jīng)典多維標(biāo)度(Classical multidimensional scaling),用于研究數(shù)據(jù)間的相似性愈魏。PCoA與PCA都是降低數(shù)據(jù)維度的方法觅玻,但是差異在在于PCA是基于原始矩陣,而PCoA是基于通過原始矩陣計算出的距離矩陣培漏。因此溪厘,PCA是盡力保留數(shù)據(jù)中的變異讓點的位置不改動,而PCoA是盡力保證原本的距離關(guān)系不發(fā)生改變牌柄,也就是使得原始數(shù)據(jù)間點的距離與投影中即結(jié)果中各點之間的距離盡可能相關(guān)(如圖)畸悬。
如何進行PCA和PCoA分析
R中有很多包都提供了PCA和PCoA,比如常用的ade4包。本文將基于該包進行PCA和PCoA的分析珊佣,數(shù)據(jù)是自帶的deug蹋宦,該數(shù)據(jù)提供了104個學(xué)生9門課程的成績(見截圖)和綜合評定。綜合評定有以下幾個等級:A+,A,B,B-,C-,D咒锻。
讓我們通過PCA和PCoA來看一看這樣的綜合評定是否合理冷冗,是否確實依據(jù)這9門課把這104個學(xué)生合理分配到不同組(每個等級一個組)。
繪圖
(1)PCA分析及作圖
前文已經(jīng)介紹了PCA是基于原始數(shù)據(jù)惑艇,所以直接進行PCA分析即可蒿辙。相信大家都比較熟悉散點圖的繪制方法,這里不再細(xì)講,PCA分析完畢后我們直接作圖展示結(jié)果思灌。
library(ade4)
library(ggplot2)
library(RColorBrewer)
data(deug)
#PCA分析
pca<- dudi.pca(deug$tab, scal = FALSE, center = deug$cent, scan = FALSE)
#坐標(biāo)軸解釋量(前兩軸)
pca_eig <- (pca$eig)[1:2] / sum(pca$eig)
#提取樣本點坐標(biāo)(前兩軸)
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
整體看起來還不錯,就是B-和C-的學(xué)生似乎難以區(qū)分习瑰。
(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)
#坐標(biāo)軸解釋量(前兩軸)
pcoa_eig <- (pcoa$eig)[1:2] / sum(pcoa$eig)
#提取樣本點坐標(biāo)(前兩軸)
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
有時候PCA和PCoA的結(jié)果差不多,有時候某種方法能夠把樣本有效分開而另一種可能效果不佳甜奄,這些都要看樣本數(shù)據(jù)的特性柠横。
(3)3D圖
除轉(zhuǎn)錄組研究以外,在16S微生物的研究中我們會根據(jù)物種豐度的文件對數(shù)據(jù)進行PCA或者PCoA分析课兄,也是我們所說的β多樣性分析牍氛。根據(jù)PCA或者PCoA的結(jié)果看感染組和對照組能否分開,以了解微生物組的總體變化情況烟阐。
β多樣性分析的概念
Beta多樣性指的是樣本間多樣性搬俊。在腸道菌群分析中,Beta多樣性是衡量個體間微生物組成相似性的一個指標(biāo)蜒茄。通過計算樣本間距離可以獲得β多樣性計算矩陣唉擂,后續(xù)一般會利用PCoA、進化樹聚類等分析對此數(shù)值關(guān)系進行圖形展示檀葛。主要基于OTU的群落比較方法玩祟,有歐式距離、bray curtis距離屿聋、Jaccard 距離空扎,這些方法優(yōu)勢在于算法簡單,考慮物種豐度(有無)和均度(相對豐度)润讥,但其沒有考慮OTUs之間的進化關(guān)系转锈,認(rèn)為OTU之間不存在進化上的聯(lián)系,每個OTU間的關(guān)系平等楚殿。另一種算法Unifrac距離法撮慨,是根據(jù)系統(tǒng)發(fā)生樹進行比較,并根據(jù)16s的序列信息對OTU進行進化樹分類勒魔, 一般有加權(quán)和非加權(quán)分析甫煞。
QIIME2中重要的Beta多樣性指數(shù):
Jaccard距離:群落差異的定性度量,即只考慮種類冠绢,不考慮豐度抚吠。
Bray-Curtis距離:群落差異的定量度量,較常用弟胀。
Unweighted UniFrac距離:包含特征之間的系統(tǒng)發(fā)育關(guān)系的群落差異定性度量楷力。
Weighted UniFrac距離:包含特征之間的系統(tǒng)發(fā)育關(guān)系的群落差異定量度量喊式。
R繪制三維PCoA圖
解壓縮通過qiime2輸出的 .qza文件,獲得繪圖的matrix和pcoa結(jié)果文件
qiime tools export \
--input-path bray_curtis_distance_matrix.qza \
--output-path bray_curtis
將pcoa結(jié)果整理成下表萧朝,保存為 ***_site.txt
#載入數(shù)據(jù)
setwd("D:/your_work_path/")
pca_site <- read.delim('***_site.txt', sep = '\t', stringsAsFactors = FALSE)
#scatterplot3d包
library(scatterplot3d)
#調(diào)整角度岔留,保存
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%'))
注意沒有l(wèi)egend,需要AI加入检柬。
后期需要繼續(xù)摸索献联,其實可以加legend的,只是目前自己的技術(shù)做不到何址。里逆。。