前言

主成分分析(Principal Components Analysis，PCA)，也稱(chēng)主分量分析或主成分回歸分析法，是一種無(wú)監(jiān)督的數(shù)據(jù)降維方法。PCA通過(guò)線(xiàn)性變換將原始數(shù)據(jù)變換為一組各維度線(xiàn)性無(wú)關(guān)的表示，可用于提取數(shù)據(jù)的主要特征分量，常用于高維數(shù)據(jù)的降維。這種降維的思想首先減少數(shù)據(jù)集的維數(shù)，同時(shí)還保持?jǐn)?shù)據(jù)集的對(duì)方差貢獻(xiàn)最大的特征，最終使數(shù)據(jù)直觀呈現(xiàn)在二維坐標(biāo)系。

數(shù)據(jù)降維展示

直觀上，第一主成分軸 優(yōu)于 第二主成分軸，具有最大可分性。
主坐標(biāo)分析（Principal Coordinates Analysis，PCoA），即經(jīng)典多維標(biāo)度（Classical multidimensional scaling），用于研究數(shù)據(jù)間的相似性。

【二者差異】

主成分分析（Principal components analysis，PCA）是一種統(tǒng)計(jì)分析、簡(jiǎn)化數(shù)據(jù)集的方法。它利用正交變換來(lái)對(duì)一系列可能相關(guān)的變量的觀測(cè)值進(jìn)行線(xiàn)性變換，從而投影為一系列線(xiàn)性不相關(guān)變量的值，這些不相關(guān)變量稱(chēng)為主成分（Principal Components）。具體地，主成分可以看做一個(gè)線(xiàn)性方程，其包含一系列線(xiàn)性系數(shù)來(lái)指示投影方向（如圖）。PCA對(duì)原始數(shù)據(jù)的正則化或預(yù)處理敏感（相對(duì)縮放）。PCA是最簡(jiǎn)單的以特征量分析多元統(tǒng)計(jì)分布的方法。通常情況下，這種運(yùn)算可以被看作是揭露數(shù)據(jù)的內(nèi)部結(jié)構(gòu)，從而更好的解釋數(shù)據(jù)的變量的方法。

PCA示意圖

主坐標(biāo)分析（Principal Coordinates Analysis，PCoA），即經(jīng)典多維標(biāo)度（Classical multidimensional scaling），用于研究數(shù)據(jù)間的相似性。PCoA與PCA都是降低數(shù)據(jù)維度的方法，但是差異在在于PCA是基于原始矩陣，而PCoA是基于通過(guò)原始矩陣計(jì)算出的距離矩陣。因此，PCA是盡力保留數(shù)據(jù)中的變異讓點(diǎn)的位置不改動(dòng)，而PCoA是盡力保證原本的距離關(guān)系不發(fā)生改變，也就是使得原始數(shù)據(jù)間點(diǎn)的距離與投影中即結(jié)果中各點(diǎn)之間的距離盡可能相關(guān)（如圖）。

PCoA示意圖

如何進(jìn)行PCA和PCoA分析

R中有很多包都提供了PCA和PCoA,比如常用的ade4包。本文將基于該包進(jìn)行PCA和PCoA的分析，數(shù)據(jù)是自帶的deug，該數(shù)據(jù)提供了104個(gè)學(xué)生9門(mén)課程的成績(jī)（見(jiàn)截圖）和綜合評(píng)定。綜合評(píng)定有以下幾個(gè)等級(jí)：A+,A,B,B-,C-,D。
讓我們通過(guò)PCA和PCoA來(lái)看一看這樣的綜合評(píng)定是否合理，是否確實(shí)依據(jù)這9門(mén)課把這104個(gè)學(xué)生合理分配到不同組（每個(gè)等級(jí)一個(gè)組）。

Deug的9門(mén)課

繪圖

（1）PCA分析及作圖

前文已經(jīng)介紹了PCA是基于原始數(shù)據(jù)，所以直接進(jìn)行PCA分析即可。相信大家都比較熟悉散點(diǎn)圖的繪制方法，這里不再細(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)

#提取樣本點(diǎn)坐標(biāo)（前兩軸）
sample_site <- data.frame({pca$li})[1:2]
sample_site$names <- rownames(sample_site)
names(sample_site)[1:2] <- c('PCA1', 'PCA2')

#以最終成績(jī)作為分組
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)+  #可在這里修改點(diǎn)的透明度、大小
  scale_color_manual(values = brewer.pal(6,"Set2")) + #可在這里修改點(diǎn)的顏色
  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

整體看起來(lái)還不錯(cuò)，就是B-和C-的學(xué)生似乎難以區(qū)分。

（2）PCoA分析及作圖

library(ade4)
library(ggplot2)
library(RColorBrewer)
library(vegan)#用于計(jì)算距離
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)

#提取樣本點(diǎn)坐標(biāo)（前兩軸）
sample_site <- data.frame({pcoa$li})[1:2]
sample_site$names <- rownames(sample_site)
names(sample_site)[1:2] <- c('PCoA1', 'PCoA2')

#以最終成績(jī)作為分組
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)+  #可在這里修改點(diǎn)的透明度、大小
  scale_color_manual(values = brewer.pal(6,"Set2")) + #可在這里修改點(diǎn)的顏色
  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

有時(shí)候PCA和PCoA的結(jié)果差不多，有時(shí)候某種方法能夠把樣本有效分開(kāi)而另一種可能效果不佳，這些都要看樣本數(shù)據(jù)的特性。

（3）3D圖

除轉(zhuǎn)錄組研究以外，在16S微生物的研究中我們會(huì)根據(jù)物種豐度的文件對(duì)數(shù)據(jù)進(jìn)行PCA或者PCoA分析，也是我們所說(shuō)的β多樣性分析。根據(jù)PCA或者PCoA的結(jié)果看感染組和對(duì)照組能否分開(kāi)，以了解微生物組的總體變化情況。

β多樣性分析的概念
Beta多樣性指的是樣本間多樣性。在腸道菌群分析中，Beta多樣性是衡量個(gè)體間微生物組成相似性的一個(gè)指標(biāo)。通過(guò)計(jì)算樣本間距離可以獲得β多樣性計(jì)算矩陣，后續(xù)一般會(huì)利用PCoA、進(jìn)化樹(shù)聚類(lèi)等分析對(duì)此數(shù)值關(guān)系進(jìn)行圖形展示。主要基于OTU的群落比較方法，有歐式距離、bray curtis距離、Jaccard 距離，這些方法優(yōu)勢(shì)在于算法簡(jiǎn)單，考慮物種豐度（有無(wú)）和均度（相對(duì)豐度），但其沒(méi)有考慮OTUs之間的進(jìn)化關(guān)系，認(rèn)為OTU之間不存在進(jìn)化上的聯(lián)系，每個(gè)OTU間的關(guān)系平等。另一種算法Unifrac距離法，是根據(jù)系統(tǒng)發(fā)生樹(shù)進(jìn)行比較，并根據(jù)16s的序列信息對(duì)OTU進(jìn)行進(jìn)化樹(shù)分類(lèi)， 一般有加權(quán)和非加權(quán)分析。

QIIME2中重要的Beta多樣性指數(shù)：

Jaccard距離：群落差異的定性度量，即只考慮種類(lèi)，不考慮豐度。

Bray-Curtis距離：群落差異的定量度量，較常用。

Unweighted UniFrac距離：包含特征之間的系統(tǒng)發(fā)育關(guān)系的群落差異定性度量。

Weighted UniFrac距離：包含特征之間的系統(tǒng)發(fā)育關(guān)系的群落差異定量度量。

R繪制三維PCoA圖

解壓縮通過(guò)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%'))

3D PCoA

注意沒(méi)有l(wèi)egend，需要AI加入。
后期需要繼續(xù)摸索，其實(shí)可以加legend的，只是目前自己的技術(shù)做不到。。。

PCA思想解析：
http://m.itdecent.cn/p/09bae5cbdc53