R語言的主成分分析（PCA）詳解和帶聚類的PCA圖繪制
最近有個老師在整理文章數(shù)據(jù)曲楚，由于分組較多棕所，想展示PCA圖瀑踢，說明不同分組的差異擂啥，公司默認給出的PCA圖比較樸素哄陶，想做的好看一些，怎么來做呢啤它。自然難不倒小編奕筐，以下幾種展示方法，讓文章圖表更炫酷（前提是樣品多变骡，分組多离赫，如果2組3重復(fù)，是不夠的）塌碌。
在RNA-seq中渊胸，主成分分析（PCA）是最常見的多元數(shù)據(jù)分析類型之一√ㄗ保基因表達定量后獲得了各樣本中所有基因的表達值信息翎猛，隨后我們通常會期望比較樣本之間在基因表達值的整體相似性或者差異程度〗邮＃基因數(shù)量成千上萬切厘，肯定不能對每個基因的表達都作個比較，這時候就要用到“降維”算法懊缺，PCA分析因此派上用場疫稿。PCA設(shè)法將N維（N=基因數(shù)量）的表達矩陣降維成只有少數(shù)幾個主成分的形式，這幾個主成分就可以代表所有基因的整體表達格局，進而據(jù)此描述樣本差異遗座。
例如這是文獻中的一些PCA分析圖舀凛。
PCA圖中，如果兩個樣本間整體基因表達值更為相似途蒋，那么它們的距離就接近猛遍；反之若整體基因表達值差異很大，則它們的距離就會更遠号坡。據(jù)此懊烤，我們就可以從中評估，不同組間樣本的基因表達是否相差更為明顯筋帖，組內(nèi)樣本的基因表達是否均勻或一致性較高奸晴，哪些處理組之間引起了相似的基因表達變化趨勢，等等日麸。
本篇教程就讓我們來學習如何使用R語言進行PCA分析。示例數(shù)據(jù)和R代碼等逮光，可添加微信公眾號”獲取代箭。
1 示例數(shù)據(jù)
“PCA.data.txt”為基因表達值矩陣。其中第一列“ensembl_id”為基因名稱涕刚，這里以ensembl id作為指代嗡综；其余各列為各樣本的名稱，記錄了RNA-seq獲得的各基因在各樣本中的表達值信息杜漠。
“group.txt”則為樣本分組文件极景，記錄了樣本所屬的不同分組。對于示例數(shù)據(jù)而言驾茴，共包含兩組盼樟，即對照組（control）和處理組（treat），每組各5個樣本锈至。
接下來晨缴，就以基因表達值矩陣文件作為輸入，展示如何進行PCA分析的一般過程峡捡。
2 R語言執(zhí)行主成分分析（PCA）
首先击碗，需要根據(jù)基因表達值進行樣本間的PCA分析，確定樣本在PCA圖中的位置们拙。
將上述基因表達值矩陣讀入到R中進行計算稍途。R語言中能夠執(zhí)行PCA分析的方法有很多，不過它們的算法都是統(tǒng)一的砚婆，隨便使用任何一個R包就可以械拍，例如這里選擇使用FactoMineR包中的PCA方法。
#讀取基因表達值矩陣
#推薦使用 log 轉(zhuǎn)化后的基因表達值，降低不同基因表達水平數(shù)量級相差過大的問題
gene <- read.delim('PCA.data.txt', row.names = 1, sep = '\t')
?
#將基因表達值矩陣作個轉(zhuǎn)置殊者，使行為樣本与境，列為基因
gene <- t(gene)
?
#我們使用 FactoMineR 包中的方法，實現(xiàn) PCA 分析和聚類添加
library(FactoMineR)
?
#樣本中基因表達值的 PCA 分析
gene.pca <- PCA(gene, ncp = 2, scale.unit = TRUE, graph = FALSE)
plot(gene.pca) #PCA 簡圖
這樣猖吴，PCA分析結(jié)果就初步得到了摔刁。從結(jié)果中我們可以看出，處理組和對照組的基因表達值整體差異還是非常明顯的海蔽，因為兩組的樣本能夠在PCA圖中區(qū)分很開共屈。
3 可視化PCA圖
上一步獲得了PCA分析結(jié)果，并觀察到明顯的組間差異党窜。但如果想把結(jié)果往文章中擺放拗引，上圖肯定是不行的，起碼要讓它好看一些幌衣。因此接下來矾削，繼續(xù)展示如何繪制一個好看的PCA圖。
例如這里選擇使用ggplot2包美化PCA圖豁护，它是一款非常出名的R語言作圖包哼凯。不過在使用ggplot2作圖之前，需要事先在上述PCA分析結(jié)果中將關(guān)鍵信息提取出楚里，例如樣本點在PCA圖中的位置信息断部，以及PCA軸的貢獻度等。
#提取樣本在 PCA 前兩軸中的坐標
pca_sample <- data.frame(gene.pca$ind$coord[ ,1:2])
head(pca_sample)
?
#提取 PCA 前兩軸的貢獻度
pca_eig1 <- round(gene.pca$eig[1,2], 2)
pca_eig2 <- round(gene.pca$eig[2,2],2 )
隨后班缎，加載ggplot2作圖包蝴光，根據(jù)提取出的樣本位置坐標以及PCA軸的貢獻度數(shù)值，繪制二維散點圖表示PCA結(jié)果达址。在同時蔑祟，我們也將樣本分組文件讀取到R中用于指定樣本的分組信息，以在圖中使用不同的顏色表示不同組的樣本苏携。
FactoMineR包也能用于繪制這類統(tǒng)計圖做瞪，以上述添加層次聚類的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ù)中包含了樣本坐標和分組信息
?
#ggplot2 繪制二維散點圖
library(ggplot2)
?
p <- ggplot(data = pca_sample, aes(x = Dim.1, y = Dim.2)) +
geom_point(aes(color = group), size = 3) + #根據(jù)樣本坐標繪制二維散點圖
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 軸貢獻度添加到坐標軸標題中
?
p
#標識樣本名稱右冻，使用 ggplot2 的拓展包 ggrepel 來完成
library(ggrepel)
?
p <- p +
geom_text_repel(aes(label = samples), size = 3, show.legend = FALSE,
box.padding = unit(0.5, 'lines'))
p
樣式是不是比最初的好看很多装蓬？
再一次，能夠很清晰地從圖中觀察到纱扭，處理組和對照組的基因表達值整體差異是非常明顯的牍帚，因為兩組的樣本能夠在PCA圖中區(qū)分很開。
4 可選添加樣本聚類
繼續(xù)乳蛾，可以選擇在PCA圖中展示“樣本聚類”暗赶。方法大致分為兩種鄙币，一種是通過樣本點的坐標擬合置信橢圓，另一種是直接將各組的邊界區(qū)樣本點以多邊形連接蹂随。不過需要注意的是十嘿，這種方法并不是真正的聚類，而是一種用于在作圖時更好地區(qū)分不同組的方法岳锁。
#添加 95% 置信橢圓绩衷，可用于表示對象分類，但只能適用于各組樣本數(shù)大于 5 個的情況
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'))
?
#多邊形連接同類別對象邊界的樣式激率，適用于各組樣本數(shù)大于 3 個的情況
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'))
這樣咳燕，一個完整的PCA分析過程就完整地展示出來了，包括輸入文件準備乒躺，如何計算PCA招盲，以及PCA圖的可視化等，還是非常簡單的對嗎嘉冒？
