一文看懂PCA主成分分析中介紹了PCA分析的原理和分析的意義(基本簡介如下凳怨，更多見博客)瑰艘，今天就用數(shù)據(jù)來實(shí)際操練一下。
在生信寶典公眾號(hào)后臺(tái)回復(fù)“PCA實(shí)戰(zhàn)”肤舞，獲取測(cè)試數(shù)據(jù)紫新。
一、PCA應(yīng)用
# 加載需要用到的R包library(psych)
library(reshape2)
library(ggplot2)
library(factoextra)
1. 數(shù)據(jù)初始化
# 基因表達(dá)數(shù)據(jù)
exprData <- "ehbio_salmon.DESeq2.normalized.symbol.txt"
# 非必須
sampleFile <- "sampleFile"
2. 數(shù)據(jù)讀入
# 為了保證文章的使用李剖，文末附有數(shù)據(jù)的新下載鏈接芒率，以防原鏈接失效
data <- read.table(exprData, header=T, row.names=NULL,sep="\t")
# 處理重復(fù)名字，謹(jǐn)慎處理篙顺，先找到名字重復(fù)的原因再?zèng)Q定是否需要按一下方式都保留
rownames_data <- make.names(data[,1],unique=T)
data <- data[,-1,drop=F]
rownames(data) <- rownames_data
data <- data[rowSums(data)>0,]
# 去掉方差為0 的行偶芍，這些本身沒有意義，也妨礙后續(xù)運(yùn)算
data <- data[apply(data, 1, var)!=0,]
3. 數(shù)據(jù)預(yù)處理（可選）
# 計(jì)算中值絕對(duì)偏差 (MAD, median absolute deviation)度量基因表達(dá)變化幅度
# 在基因表達(dá)中德玫，盡管某些基因很小的變化會(huì)導(dǎo)致重要的生物學(xué)意義匪蟀，
# 但是很小的觀察值會(huì)引入很大的背景噪音，因此也意義不大宰僧。
# 可以選擇根據(jù)mad值過濾材彪，取top 50， 500等做分析
mads <- apply(data, 1, mad)data <- data[rev(order(mads)),]
# 此處未選
# data <- data[1:500,]
dim(data)
4. PCA分析
# Pay attention to the format of PCA input
# Rows are samples and columns are variables
data_t <- t(data)
variableL <- ncol(data_t)
if(sampleFile != "") {
sample <- read.table(sampleFile,header = T, row.names=1,sep="\t")
data_t_m <- merge(data_t, sample, by=0)
rownames(data_t_m) <- data_t_m$Row.names
data_t <- data_t_m[,-1]
}
# By default, prcomp will centralized the data using mean.
# Normalize data for PCA by dividing each data by column standard deviation.# Often, we would normalize data.
# Only when we care about the real number changes other than the trends,
# `scale` can be set to TRUE.
# We will show the differences of scaling and un-scaling effects.
pca <- prcomp(data_t[,1:variableL], scale=T)
# sdev: standard deviation of the principle components.
# Square to get variance
# percentVar <- pca$sdev^2 / sum( pca$sdev^2)
# To check what's in pca
print(str(pca))
5. PCA結(jié)果展示
# PCA結(jié)果提取和可視化神器
# http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials/
library(factoextra)
1. 碎石圖展示每個(gè)主成分的貢獻(xiàn)
# 如果需要保持琴儿，去掉下面第1,3行的注釋
#pdf("1.pdf")
fviz_eig(pca, addlabels = TRUE)
#dev.off()
2. PCA樣品聚類信息展示
# repel=T段化，自動(dòng)調(diào)整文本位置
fviz_pca_ind(pca, repel=T)
3. 根據(jù)樣品分組上色
# 根據(jù)分組上色并繪制
fviz_pca_ind(pca, col.ind=data_t$conditions, mean.point=F, addEllipses = T, legend.title="Groups")
4. 增加不同屬性的橢圓
“convex”: plot convex hull of a set o points.
“confidence”: plot confidence ellipses around group mean points as the function coord.ellipse() [in FactoMineR].
“t”: assumes a multivariate t-distribution.
“norm”: assumes a multivariate normal distribution.
“euclid”: draws a circle with the radius equal to level, representing the euclidean distance from the center. This ellipse probably won’t appear circular unless coord_fixed() is applied.
# 根據(jù)分組上色并繪制95%置信區(qū)間
fviz_pca_ind(pca, col.ind=data_t$conditions, mean.point=F, addEllipses = T, legend.title="Groups", ellipse.type="confidence", ellipse.level=0.95)
5. 展示貢獻(xiàn)最大的變量 (基因)
1） 展示與主坐標(biāo)軸的相關(guān)性大于0.99的變量 (具體數(shù)字自己調(diào)整)
# Visualize variable with cos2 >= 0.99
fviz_pca_var(pca, select.var = list(cos2 = 0.99), repel=T, col.var = "cos2", geom.var = c("arrow", "text") )
2）展示與主坐標(biāo)軸最相關(guān)的10個(gè)變量
# Top 10 active variables with the highest cos2
fviz_pca_var(pca, select.var= list(cos2 = 10), repel=T, col.var = "contrib")
3）展示自己關(guān)心的變量（基因）與主坐標(biāo)軸的相關(guān)性分布
# Select by names
# 這里選擇的是MAD值最大的幾個(gè)基因
name <- list(name = c("FN1", "DCN", "CEMIP","CCDC80","IGFBP5","COL1A1","GREM1"))
fviz_pca_var(pca, select.var = name)
6. biPLot同時(shí)展示樣本分組和關(guān)鍵基因
# top 5 contributing individuals and variable
fviz_pca_biplot(pca,
fill.ind=data_t$conditions,
palette="joo",
addEllipses = T,
ellipse.type="confidence",
ellipse.level=0.95,
mean.point=F,
col.var="contrib",
gradient.cols = "RdYlBu",
select.var = list(contrib = 10),
ggtheme = theme_minimal())
二、PCA分析注意事項(xiàng)
一般說來凤类，在PCA之前原始數(shù)據(jù)需要中心化（centering穗泵，數(shù)值減去平均值）。中心化的方法很多谜疤，除了平均值中心化（mean-centering）外佃延，還包括其它更穩(wěn)健的方法现诀，比如中位數(shù)中心化等。
除了中心化以外履肃，定標(biāo) (Scale, 數(shù)值除以標(biāo)準(zhǔn)差) 也是數(shù)據(jù)前處理中需要考慮的一點(diǎn)仔沿。如果數(shù)據(jù)沒有定標(biāo)，則原始數(shù)據(jù)中方差大的變量對(duì)主成分的貢獻(xiàn)會(huì)很大尺棋。數(shù)據(jù)的方差與其量級(jí)成指數(shù)關(guān)系封锉，比如一組數(shù)據(jù)
(1,2,3,4)的方差是1.67，而(10,20,30,40)的方差就是167,數(shù)據(jù)變大10倍膘螟，方差放大了100倍成福。但是定標(biāo)(scale)可能會(huì)有一些負(fù)面效果，因?yàn)槎?biāo)后變量之間的權(quán)重就是變得相同荆残。如果我們的變量中有噪音的話奴艾，我們就在無形中把噪音和信息的權(quán)重變得相同，但PCA本身無法區(qū)分信號(hào)和噪音内斯。在這樣的情形下蕴潦，我們就不必做定標(biāo)。
一般而言俘闯，對(duì)于度量單位不同的指標(biāo)或是取值范圍彼此差異非常大的指標(biāo)不直接由其協(xié)方差矩陣出發(fā)進(jìn)行主成分分析潭苞，而應(yīng)該考慮對(duì)數(shù)據(jù)的標(biāo)準(zhǔn)化。比如度量單位不同真朗，有萬人此疹、萬噸、萬元蜜猾、億元秀菱，而數(shù)據(jù)間的差異性也非常大，小則幾十大則幾萬蹭睡，因此在用協(xié)方差矩陣求解主成分時(shí)存在協(xié)方差矩陣中數(shù)據(jù)的差異性很大衍菱。在后面提取主成分時(shí)發(fā)現(xiàn)，只提取了一個(gè)主成分肩豁，而此時(shí)并不能將所有的變量都解釋到脊串，這就沒有真正起到降維的作用。此時(shí)就需要對(duì)數(shù)據(jù)進(jìn)行定標(biāo)(scale)清钥，這樣提取的主成分可以覆蓋更多的變量琼锋，這就實(shí)現(xiàn)主成分分析的最終目的。但是對(duì)原始數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化后更傾向于使得各個(gè)指標(biāo)的作用在主成分分析構(gòu)成中相等祟昭。對(duì)于數(shù)據(jù)取值范圍不大或是度量單位相同的指標(biāo)進(jìn)行標(biāo)準(zhǔn)化處理后缕坎，其主成分分析的結(jié)果與仍由協(xié)方差矩陣出發(fā)求得的結(jié)果有較大區(qū)別。這是因?yàn)閷?duì)數(shù)據(jù)標(biāo)準(zhǔn)化的過程實(shí)際上就是抹殺原有變量離散程度差異的過程篡悟，標(biāo)準(zhǔn)化后方差均為1谜叹，而實(shí)際上方差是對(duì)數(shù)據(jù)信息的重要概括形式匾寝，也就是說，對(duì)原始數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化后抹殺了一部分重要信息荷腊，因此才使得標(biāo)準(zhǔn)化后各變量在主成分構(gòu)成中的作用趨于相等艳悔。因此，對(duì)同度量或是取值范圍在同量級(jí)的數(shù)據(jù)還是直接使用非定標(biāo)數(shù)據(jù)求解主成分為宜女仰。
中心化和定標(biāo)都會(huì)受數(shù)據(jù)中離群值（outliers）或者數(shù)據(jù)不均勻（比如數(shù)據(jù)被分為若干個(gè)小組）的影響猜年，應(yīng)該用更穩(wěn)健的中心化和定標(biāo)方法。
PCA也存在一些限制疾忍，例如它可以很好的解除線性相關(guān)乔外，但是對(duì)于高階相關(guān)性就沒有辦法了，對(duì)于存在高階相關(guān)性的數(shù)據(jù)一罩，可以考慮Kernel PCA袁稽，通過Kernel函數(shù)將非線性相關(guān)轉(zhuǎn)為線性相關(guān)，關(guān)于這點(diǎn)就不展開討論了擒抛。另外，PCA假設(shè)數(shù)據(jù)各主特征是分布在正交方向上补疑，如果在非正交方向上存在幾個(gè)方差較大的方向歧沪，PCA的效果就大打折扣了。
參考資料
http://stackoverflow.com/questions/22092220/plot-only-y-axis-but-nothing-else
http://www.sthda.com/english/wiki/scatterplot3d-3d-graphics-r-software-and-data-visualization
http://gastonsanchez.com/how-to/2014/01/15/Center-data-in-R/
http://www.statpower.net/Content/310/R Stuff/SampleMarkdown.html
http://www2.edu-edu.com.cn/lesson_crs78/self/02198/resource/contents/ch_05/ch_05.html
http://stackoverflow.com/questions/1249548/side-by-side-plots-with-ggplot2
R統(tǒng)計(jì)和作圖
更多閱讀
心得體會(huì)癌癥數(shù)據(jù)庫LinuxPython
