R語(yǔ)言的主成分分析（PCA）詳解和帶聚類(lèi)的PCA圖繪制

最近有個(gè)老師在整理文章數(shù)據(jù)，由于分組較多，想展示PCA圖，說(shuō)明不同分組的差異，公司默認(rèn)給出的PCA圖比較樸素，想做的好看一些，怎么來(lái)做呢。自然難不倒小編，以下幾種展示方法，讓文章圖表更炫酷（前提是樣品多，分組多，如果2組3重復(fù)，是不夠的）。
在RNA-seq中，主成分分析（PCA）是最常見(jiàn)的多元數(shù)據(jù)分析類(lèi)型之一。基因表達(dá)定量后獲得了各樣本中所有基因的表達(dá)值信息，隨后我們通常會(huì)期望比較樣本之間在基因表達(dá)值的整體相似性或者差異程度?；驍?shù)量成千上萬(wàn)，肯定不能對(duì)每個(gè)基因的表達(dá)都作個(gè)比較，這時(shí)候就要用到“降維”算法，PCA分析因此派上用場(chǎng)。PCA設(shè)法將N維（N=基因數(shù)量）的表達(dá)矩陣降維成只有少數(shù)幾個(gè)主成分的形式，這幾個(gè)主成分就可以代表所有基因的整體表達(dá)格局，進(jìn)而據(jù)此描述樣本差異。

例如這是文獻(xiàn)中的一些PCA分析圖。

PCA圖中，如果兩個(gè)樣本間整體基因表達(dá)值更為相似，那么它們的距離就接近；反之若整體基因表達(dá)值差異很大，則它們的距離就會(huì)更遠(yuǎn)。據(jù)此，我們就可以從中評(píng)估，不同組間樣本的基因表達(dá)是否相差更為明顯，組內(nèi)樣本的基因表達(dá)是否均勻或一致性較高，哪些處理組之間引起了相似的基因表達(dá)變化趨勢(shì)，等等。

本篇教程就讓我們來(lái)學(xué)習(xí)如何使用R語(yǔ)言進(jìn)行PCA分析。示例數(shù)據(jù)和R代碼等，可添加微信公眾號(hào)”獲取。

1 示例數(shù)據(jù)

“PCA.data.txt”為基因表達(dá)值矩陣。其中第一列“ensembl_id”為基因名稱，這里以ensembl id作為指代；其余各列為各樣本的名稱，記錄了RNA-seq獲得的各基因在各樣本中的表達(dá)值信息。

“group.txt”則為樣本分組文件，記錄了樣本所屬的不同分組。對(duì)于示例數(shù)據(jù)而言，共包含兩組，即對(duì)照組（control）和處理組（treat），每組各5個(gè)樣本。

接下來(lái)，就以基因表達(dá)值矩陣文件作為輸入，展示如何進(jìn)行PCA分析的一般過(guò)程。

2 R語(yǔ)言執(zhí)行主成分分析（PCA）

首先，需要根據(jù)基因表達(dá)值進(jìn)行樣本間的PCA分析，確定樣本在PCA圖中的位置。

將上述基因表達(dá)值矩陣讀入到R中進(jìn)行計(jì)算。R語(yǔ)言中能夠執(zhí)行PCA分析的方法有很多，不過(guò)它們的算法都是統(tǒng)一的，隨便使用任何一個(gè)R包就可以，例如這里選擇使用FactoMineR包中的PCA方法。

#讀取基因表達(dá)值矩陣
#推薦使用 log 轉(zhuǎn)化后的基因表達(dá)值，降低不同基因表達(dá)水平數(shù)量級(jí)相差過(guò)大的問(wèn)題
gene <- read.delim('PCA.data.txt', row.names = 1, sep = '\t')
?
#將基因表達(dá)值矩陣作個(gè)轉(zhuǎn)置，使行為樣本，列為基因
gene <- t(gene)
?
#我們使用 FactoMineR 包中的方法，實(shí)現(xiàn) PCA 分析和聚類(lèi)添加
library(FactoMineR)
?
#樣本中基因表達(dá)值的 PCA 分析
gene.pca <- PCA(gene, ncp = 2, scale.unit = TRUE, graph = FALSE)
plot(gene.pca)  #PCA 簡(jiǎn)圖

這樣，PCA分析結(jié)果就初步得到了。從結(jié)果中我們可以看出，處理組和對(duì)照組的基因表達(dá)值整體差異還是非常明顯的，因?yàn)閮山M的樣本能夠在PCA圖中區(qū)分很開(kāi)。

3 可視化PCA圖

上一步獲得了PCA分析結(jié)果，并觀察到明顯的組間差異。但如果想把結(jié)果往文章中擺放，上圖肯定是不行的，起碼要讓它好看一些。因此接下來(lái)，繼續(xù)展示如何繪制一個(gè)好看的PCA圖。

例如這里選擇使用ggplot2包美化PCA圖，它是一款非常出名的R語(yǔ)言作圖包。不過(guò)在使用ggplot2作圖之前，需要事先在上述PCA分析結(jié)果中將關(guān)鍵信息提取出，例如樣本點(diǎn)在PCA圖中的位置信息，以及PCA軸的貢獻(xiàn)度等。

#提取樣本在 PCA 前兩軸中的坐標(biāo)
pca_sample <- data.frame(gene.pca$ind$coord[ ,1:2])
head(pca_sample)
?
#提取 PCA 前兩軸的貢獻(xiàn)度
pca_eig1 <- round(gene.pca$eig[1,2], 2)
pca_eig2 <- round(gene.pca$eig[2,2],2 )

隨后，加載ggplot2作圖包，根據(jù)提取出的樣本位置坐標(biāo)以及PCA軸的貢獻(xiàn)度數(shù)值，繪制二維散點(diǎn)圖表示PCA結(jié)果。在同時(shí)，我們也將樣本分組文件讀取到R中用于指定樣本的分組信息，以在圖中使用不同的顏色表示不同組的樣本。

FactoMineR包也能用于繪制這類(lèi)統(tǒng)計(jì)圖，以上述添加層次聚類(lèi)的PCA繼續(xù)。

#讀取并合并樣本分組信息
group <- read.delim('group.txt', row.names = 1, sep = '\t', check.names = FALSE)
group <- group[rownames(pca_sample), ]
pca_sample <- cbind(pca_sample, group)
?
pca_sample$samples <- rownames(pca_sample)
head(pca_sample)  #作圖數(shù)據(jù)中包含了樣本坐標(biāo)和分組信息
?
#ggplot2 繪制二維散點(diǎn)圖
library(ggplot2)
?
p <- ggplot(data = pca_sample, aes(x = Dim.1, y = Dim.2)) +
geom_point(aes(color = group), size = 3) +  #根據(jù)樣本坐標(biāo)繪制二維散點(diǎn)圖
scale_color_manual(values = c('orange', 'purple')) +  #自定義顏色
theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), 
    legend.key = element_rect(fill = 'transparent')) +  #去除背景和網(wǎng)格線
labs(x =  paste('PCA1:', pca_eig1, '%'), y = paste('PCA2:', pca_eig2, '%'), color = '')  #將 PCA 軸貢獻(xiàn)度添加到坐標(biāo)軸標(biāo)題中
?
p

#標(biāo)識(shí)樣本名稱，使用 ggplot2 的拓展包 ggrepel 來(lái)完成
library(ggrepel)
?
p <- p + 
geom_text_repel(aes(label = samples), size = 3, show.legend = FALSE, 
    box.padding = unit(0.5, 'lines'))

p

樣式是不是比最初的好看很多？

再一次，能夠很清晰地從圖中觀察到，處理組和對(duì)照組的基因表達(dá)值整體差異是非常明顯的，因?yàn)閮山M的樣本能夠在PCA圖中區(qū)分很開(kāi)。

4 可選添加樣本聚類(lèi)

繼續(xù)，可以選擇在PCA圖中展示“樣本聚類(lèi)”。方法大致分為兩種，一種是通過(guò)樣本點(diǎn)的坐標(biāo)擬合置信橢圓，另一種是直接將各組的邊界區(qū)樣本點(diǎn)以多邊形連接。不過(guò)需要注意的是，這種方法并不是真正的聚類(lèi)，而是一種用于在作圖時(shí)更好地區(qū)分不同組的方法。

#添加 95% 置信橢圓，可用于表示對(duì)象分類(lèi)，但只能適用于各組樣本數(shù)大于 5 個(gè)的情況
p + stat_ellipse(aes(color = group), level = 0.95, show.legend = FALSE)
?
p + stat_ellipse(aes(fill = group), geom = 'polygon', level = 0.95, alpha = 0.1, show.legend = FALSE) +
scale_fill_manual(values = c('orange', 'purple'))
?
#多邊形連接同類(lèi)別對(duì)象邊界的樣式，適用于各組樣本數(shù)大于 3 個(gè)的情況
library(plyr)
cluster_border <- ddply(pca_sample, 'group', function(df) df[chull(df[[1]], df[[2]]), ])
?
p + geom_polygon(data = cluster_border, aes(color = group), fill = NA, show.legend = FALSE)
?
p + geom_polygon(data = cluster_border, aes(fill = group), alpha = 0.2, show.legend = FALSE) +
scale_fill_manual(values = c('orange', 'purple'))

這樣，一個(gè)完整的PCA分析過(guò)程就完整地展示出來(lái)了，包括輸入文件準(zhǔn)備，如何計(jì)算PCA，以及PCA圖的可視化等，還是非常簡(jiǎn)單的對(duì)嗎？