一文讀懂PCA分析 （原理、算法、解釋和可視化）

生物信息學(xué)習(xí)的正確姿勢(shì)

NGS系列文章包括NGS基礎(chǔ)、高顏值在線繪圖和分析、轉(zhuǎn)錄組分析?（Nature重磅綜述|關(guān)于RNA-seq你想知道的全在這）、ChIP-seq分析?（ChIP-seq基本分析流程）、單細(xì)胞測(cè)序分析?(重磅綜述：三萬(wàn)字長(zhǎng)文讀懂單細(xì)胞RNA測(cè)序分析的最佳實(shí)踐教程)、DNA甲基化分析、重測(cè)序分析、GEO數(shù)據(jù)挖掘（典型醫(yī)學(xué)設(shè)計(jì)實(shí)驗(yàn)GEO數(shù)據(jù)分析 (step-by-step)）、批次效應(yīng)處理等內(nèi)容。

library(knitr)
library(psych)
library(reshape2)
library(ggplot2)
library(ggbeeswarm)
library(scatterplot3d)
library(useful)
library(ggfortify)

mat_show <- function(matr) {

    printmrow <- function(x) {
        ret <- paste(paste(x,collapse = " & "),"\\\\\n")
        sprintf(ret)
    }

    align_str <- paste0('{',paste0(rep('r',ncol(matr)), collapse=""),'}')

    format_mat <- apply(matr,1,printmrow)
    add_env <- paste("\\left[\\begin{array}", align_str, 
        paste(format_mat, collapse=' '),"\\end{array}\\right]")
    return(add_env)
}

主成分分析簡(jiǎn)介

主成分分析 (PCA, principal component analysis)是一種數(shù)學(xué)降維方法, 利用正交變換 (orthogonal transformation)把一系列可能線性相關(guān)的變量轉(zhuǎn)換為一組線性不相關(guān)的新變量，也稱為主成分，從而利用新變量在更小的維度下展示數(shù)據(jù)的特征。

主成分是原有變量的線性組合，其數(shù)目不多于原始變量。組合之后，相當(dāng)于我們獲得了一批新的觀測(cè)數(shù)據(jù)，這些數(shù)據(jù)的含義不同于原有數(shù)據(jù)，但包含了之前數(shù)據(jù)的大部分特征，并且有著較低的維度，便于進(jìn)一步的分析。

在空間上，PCA可以理解為把原始數(shù)據(jù)投射到一個(gè)新的坐標(biāo)系統(tǒng)，第一主成分為第一坐標(biāo)軸，它的含義代表了原始數(shù)據(jù)中多個(gè)變量經(jīng)過(guò)某種變換得到的新變量的變化區(qū)間；第二成分為第二坐標(biāo)軸，代表了原始數(shù)據(jù)中多個(gè)變量經(jīng)過(guò)某種變換得到的第二個(gè)新變量的變化區(qū)間。這樣我們把利用原始數(shù)據(jù)解釋樣品的差異轉(zhuǎn)變?yōu)槔眯伦兞拷忉寴悠返牟町悺?/p>

這種投射方式會(huì)有很多，為了最大限度保留對(duì)原始數(shù)據(jù)的解釋，一般會(huì)用最大方差理論或最小損失理論，使得第一主成分有著最大的方差或變異數(shù) (就是說(shuō)其能盡量多的解釋原始數(shù)據(jù)的差異)；隨后的每一個(gè)主成分都與前面的主成分正交，且有著僅次于前一主成分的最大方差 (正交簡(jiǎn)單的理解就是兩個(gè)主成分空間夾角為90°，兩者之間無(wú)線性關(guān)聯(lián)，從而完成去冗余操作)。

主成分分析的意義

1. 簡(jiǎn)化運(yùn)算。

   在問(wèn)題研究中，為了全面系統(tǒng)地分析問(wèn)題，我們通常會(huì)收集眾多的影響因素也就是眾多的變量。這樣會(huì)使得研究更豐富，通常也會(huì)帶來(lái)較多的冗余數(shù)據(jù)和復(fù)雜的計(jì)算量。
   
   比如我們我們測(cè)序了100種樣品的基因表達(dá)譜借以通過(guò)分子表達(dá)水平的差異對(duì)這100種樣品進(jìn)行分類。在這個(gè)問(wèn)題中，研究的變量就是不同的基因。每個(gè)基因的表達(dá)都可以在一定程度上反應(yīng)樣品之間的差異，但某些基因之間卻有著調(diào)控、協(xié)同或拮抗的關(guān)系，表現(xiàn)為它們的表達(dá)值存在一些相關(guān)性，這就造成了統(tǒng)計(jì)數(shù)據(jù)所反映的信息存在一定程度的冗余。另外假如某些基因如持家基因在所有樣本中表達(dá)都一樣，它們對(duì)于解釋樣本的差異也沒(méi)有意義。這么多的變量在后續(xù)統(tǒng)計(jì)分析中會(huì)增大運(yùn)算量和計(jì)算復(fù)雜度，應(yīng)用PCA就可以在盡量多的保持變量所包含的信息又能維持盡量少的變量數(shù)目，幫助簡(jiǎn)化運(yùn)算和結(jié)果解釋。

2. 去除數(shù)據(jù)噪音。

   比如說(shuō)我們?cè)跇悠返闹苽溥^(guò)程中，由于不完全一致的操作，導(dǎo)致樣品的狀態(tài)有細(xì)微的改變，從而造成一些持家基因也發(fā)生了相應(yīng)的變化，但變化幅度遠(yuǎn)小于核心基因 (一般認(rèn)為噪音的方差小于信息的方差）。而PCA在降維的過(guò)程中濾去了這些變化幅度較小的噪音變化，增大了數(shù)據(jù)的信噪比。

3. 利用散點(diǎn)圖實(shí)現(xiàn)多維數(shù)據(jù)可視化。

   在上面的表達(dá)譜分析中，假如我們有1個(gè)基因，可以在線性層面對(duì)樣本進(jìn)行分類；如果我們有2個(gè)基因，可以在一個(gè)平面對(duì)樣本進(jìn)行分類；如果我們有3個(gè)基因，可以在一個(gè)立體空間對(duì)樣本進(jìn)行分類；如果有更多的基因，比如說(shuō)n個(gè)，那么每個(gè)樣品就是n維空間的一個(gè)點(diǎn)，則很難在圖形上展示樣品的分類關(guān)系。

   利用PCA分析，我們可以選取貢獻(xiàn)最大的2個(gè)或3個(gè)主成分作為數(shù)據(jù)代表用以可視化。這比直接選取三個(gè)表達(dá)變化最大的基因更能反映樣品之間的差異。（利用Pearson相關(guān)系數(shù)對(duì)樣品進(jìn)行聚類在樣品數(shù)目比較少時(shí)是一個(gè)解決辦法）

4. 發(fā)現(xiàn)隱性相關(guān)變量。

   我們?cè)诤喜⑷哂嘣甲兞康玫街鞒煞诌^(guò)程中，會(huì)發(fā)現(xiàn)某些原始變量對(duì)同一主成分有著相似的貢獻(xiàn)，也就是說(shuō)這些變量之間存在著某種相關(guān)性，為相關(guān)變量。同時(shí)也可以獲得這些變量對(duì)主成分的貢獻(xiàn)程度。對(duì)基因表達(dá)數(shù)據(jù)可以理解為發(fā)現(xiàn)了存在協(xié)同或拮抗關(guān)系的基因。

示例展示原始變量對(duì)樣品的分類

假設(shè)有一套數(shù)據(jù)集，包含100個(gè)樣品中某一基因的表達(dá)量。如下所示，每一行為一個(gè)樣品，每一列為基因的表達(dá)值。這也是做PCA分析的基本數(shù)據(jù)組織方式，每一行代表一個(gè)樣品，每一列代表一組觀察數(shù)據(jù)即一個(gè)變量。

count <- 50
Gene1_a <- rnorm(count,5,0.5)
Gene1_b <- rnorm(count,20,0.5)
grp_a <- rep('a', count)
grp_b <- rep('b', count)
cy_data <- data.frame(Gene1 = c(Gene1_a, Gene1_b), Group=c(grp_a, grp_b))
cy_data <- as.data.frame(cy_data)
label <- c(paste0(grp_a, 1:count), paste0(grp_b, 1:count))
row.names(cy_data) <- label
library(knitr)
library(psych)
# Add additional column to data only for plotting
cy_data$Y <- rep(0,count*2)
kable(headTail(cy_data), booktabs=TRUE, caption="Expression profile for Gene1 in 100 samples")

Expression profile for Gene1 in 100 samples

	Gene1	Group	Y
a1	5.99	a	0
a2	4.66	a	0
a3	4.92	a	0
a4	4.79	a	0
…	…	NA	…
b47	20.78	b	0
b48	20.5	b	0
b49	20.46	b	0
b50	19.94	b	0

從下圖可以看出，100個(gè)樣品根據(jù)Gene1表達(dá)量的不同在橫軸上被被分為了2類，可以看做是在線性水平的分類。

library("ggplot2")
library("ggbeeswarm")

# geom_quasirandom:用于畫Jitter Plot
# theme(axis.*.y): 去除Y軸
# xlim, ylim設(shè)定坐標(biāo)軸的區(qū)間
ggplot(cy_data,aes(Gene1, Y))+geom_quasirandom(aes(color=factor(Group)), groupOnX=FALSE)+
    theme(legend.position=c(0.5,0.7)) + theme(legend.title=element_blank()) + 
    scale_fill_discrete(name="Group") + theme(axis.line.y=element_blank(), 
    axis.text.y=element_blank(), axis.ticks.y=element_blank(), 
    axis.title.y=element_blank()) + ylim(-0.5,5) + xlim(0,25)

那么如果有2個(gè)基因呢？

count <- 50
Gene2_a <- rnorm(count,5,0.2)
Gene2_b <- rnorm(count,5,0.2)

cy_data2 <- data.frame(Gene1 = c(Gene1_a, Gene1_b), Gene2 = c(Gene2_a, Gene2_b), 
        Group=c(grp_a, grp_b))
cy_data2 <- as.data.frame(cy_data2)

row.names(cy_data2) <- label

kable(headTail(cy_data2), booktabs=T, 
        caption="Expression profile for Gene1 and Gene2 in 100 samples")

Expression profile for Gene1 and Gene2 in 100 samples

	Gene1	Gene2	Group
a1	5.99	5.35	a
a2	4.66	5.31	a
a3	4.92	5.12	a
a4	4.79	5.11	a
…	…	…	NA
b47	20.78	4.95	b
b48	20.5	4.9	b
b49	20.46	4.85	b
b50	19.94	4.87	b

從下圖可以看出，100個(gè)樣品根據(jù)Gene1和Gene2的表達(dá)量的不同在坐標(biāo)軸上被被分為了2類，可以看做是在平面水平的分類。而且在這個(gè)例子中，我們可以很容易的看出Gene1對(duì)樣品分類的貢獻(xiàn)要比Gene2大，因?yàn)?code style="font-family: Consolas, Monaco, Andale Mono, monospace;line-height: 1.5;font-size: 13px;overflow-wrap: break-word;padding: 2px 4px;border-radius: 4px;margin: 0px 2px;color: rgb(233, 105, 0);background: rgb(242, 247, 251) none repeat scroll 0% 0%;border-left: 4px solid rgb(220, 230, 240);white-space: pre-wrap;font-size: 15px !important;color: rgb(217, 33, 66);border-radius: 3px;background-color: rgb(247, 247, 247);color: inherit;">Gene1在樣品間的表達(dá)差異大。

ggplot(cy_data2,aes(Gene1, Gene2))+geom_point(aes(color=factor(Group)))+
    theme(legend.position=c(0.5,0.9)) + theme(legend.title=element_blank()) + 
    ylim(0,10) + xlim(0,25)

如果有3個(gè)基因呢？

count <- 50
Gene3_a <- c(rnorm(count/2,5,0.2), rnorm(count/2,15,0.2))
Gene3_b <- c(rnorm(count/2,15,0.2), rnorm(count/2,5,0.2))

data3 <- data.frame(Gene1 = c(Gene1_a, Gene1_b), Gene2 = c(Gene2_a, Gene2_b), 
        Gene3 = c(Gene3_a, Gene3_b), Group=c(grp_a, grp_b))
data3 <- as.data.frame(data3)
data3$Group <- as.factor(data3$Group)

row.names(data3) <- label

kable(headTail(data3), booktabs=T, caption="Expression profile for 3 genes in 100 samples")

Expression profile for 3 genes in 100 samples

	Gene1	Gene2	Gene3	Group
a1	5.99	5.35	5.14	a
a2	4.66	5.31	5.05	a
a3	4.92	5.12	4.88	a
a4	4.79	5.11	4.8	a
…	…	…	…	NA
b47	20.78	4.95	4.91	b
b48	20.5	4.9	5.11	b
b49	20.46	4.85	4.95	b
b50	19.94	4.87	5.18	b

從下圖可以看出，100個(gè)樣品根據(jù)Gene1、Gene2和Gene3的表達(dá)量的不同在坐標(biāo)軸上被被分為了4類，可以看做是立體空間的分類。而且在這個(gè)例子中，我們可以很容易的看出Gene1和Gene3對(duì)樣品分類的貢獻(xiàn)要比Gene2大。

library(scatterplot3d)
colorl <- c("#E69F00", "#56B4E9")
# Extract same number of colors as the Group and same Group would have same color.
colors <- colorl[as.numeric(data3$Group)]
scatterplot3d(data3[,1:3], color=colors, xlim=c(0,25), ylim=c(0,25), zlim=c(0,25), 
    angle=55, pch=16)
legend("top", legend=levels(data3$Group), col=colorl, pch=16, xpd=T, horiz=T)

當(dāng)我們向由Gene1和Gene2構(gòu)成的X-Y平面做垂線時(shí)，可以很明顯的看出，Gene2所在的軸對(duì)樣品的分類沒(méi)有貢獻(xiàn)。因?yàn)橥渡涞?code style="font-family: Consolas, Monaco, Andale Mono, monospace;line-height: 1.5;font-size: 13px;overflow-wrap: break-word;padding: 2px 4px;border-radius: 4px;margin: 0px 2px;color: rgb(233, 105, 0);background: rgb(242, 247, 251) none repeat scroll 0% 0%;border-left: 4px solid rgb(220, 230, 240);white-space: pre-wrap;font-size: 15px !important;color: rgb(217, 33, 66);border-radius: 3px;background-color: rgb(247, 247, 247);color: inherit;">X-Y屏幕上的點(diǎn)在Y軸幾乎處于同一位置。

library(scatterplot3d)
colorl <- c("#E69F00", "#56B4E9")
colors <- colorl[as.numeric(data3$Group)]
scatterplot3d(data3[,1:3], color=colors, xlim=c(0,25), ylim=c(0,25), zlim=c(0,25),type='h')
legend("top", legend=levels(data3$Group), col=colorl, pch=16, xpd=T, horiz=T)

我們把坐標(biāo)軸做一個(gè)轉(zhuǎn)換，可以看到在由Gene1和Gene3構(gòu)成的X-Y平面上，樣品被分為了4類。Gene2對(duì)樣品的分類幾乎沒(méi)有貢獻(xiàn)，因?yàn)閹缀跛袠悠吩贕ene2維度上的值都一樣。

library(scatterplot3d)
colorl <- c("#E69F00", "#56B4E9")
colors <- colorl[as.numeric(data3$Group)]
scatterplot3d(x=data3$Gene1, y= data3$Gene3, z= data3$Gene2, color=colors, 
        xlim=c(0,25), ylim=c(0,25), zlim=c(0,25),type='h')
legend("top", legend=levels(data3$Group), col=colorl, pch=16, xpd=T, horiz=T)

在上述例子中，我們可以很容易的區(qū)分出Gene1和Gene3可以作為分類的主成分，而Gene2則對(duì)分類沒(méi)有幫助，可以在計(jì)算中去除。

但是如果我們測(cè)序了幾萬(wàn)個(gè)基因的表達(dá)時(shí)，就很難通過(guò)肉眼去看，或者作出一個(gè)圖供我們篩選哪些基因?qū)颖痉诸愗暙I(xiàn)大。這時(shí)我們應(yīng)該怎么做呢？

其中有一個(gè)方法是，在這個(gè)基因表達(dá)矩陣中選出3個(gè)變化最大的基因做為3個(gè)主成分對(duì)樣品進(jìn)行分類。我們?cè)囼?yàn)下效果怎么樣。

# 數(shù)據(jù)集來(lái)源于 http://satijalab.org/seurat/old-get-started/
# 原始下載鏈接 http://www.broadinstitute.org/~rahuls/seurat/seurat_files_nbt.zip
# 為了保證文章的使用，文末附有數(shù)據(jù)的新下載鏈接，以防原鏈接失效
data4 <- read.table("data/HiSeq301-RSEM-linear_values.txt", header=T, row.names=1,sep="\t")
dim(data4)

## [1] 23730   301

library(useful)
kable(corner(data4,r=15,c=8), booktabs=T, caption="Gene expression matrix")

Gene expression matrix

	Hi_2338_1	Hi_2338_2	Hi_2338_3	Hi_2338_4	Hi_2338_5	Hi_2338_6	Hi_2338_7	Hi_2338_8
A1BG	9.08	0.00	0.00	1.75	0.00	0.40	0.00	0.78
A1BG-AS1	0.00	0.00	3.47	0.36	0.00	0.00	0.00	0.00
A1CF	0.00	0.05	0.00	0.00	0.00	0.00	0.00	0.00
A2LD1	0.00	0.00	0.00	0.29	0.00	9.19	0.00	0.00
A2M	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
A2M-AS1	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
A2ML1	0.10	0.00	0.14	0.00	0.00	0.00	0.00	0.00
A2MP1	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
A4GALT	0.57	0.00	0.00	0.00	0.00	0.00	0.35	0.00
A4GNT	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
AA06	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
AAAS	38.95	0.00	0.00	4.44	0.00	32.90	0.00	5.58
AACS	0.12	0.00	0.00	0.00	0.58	1.03	2.16	65.74
AACSP1	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
AADAC	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00

我們篩選變異系數(shù)最大的3個(gè)基因。在這之前我們先剔除在少于5個(gè)樣品中表達(dá)的基因和少于1000個(gè)表達(dá)的基因樣品 （這里我們把表達(dá)值不小于1的基因視為表達(dá)的基因），并把所有基因根據(jù)其在不同樣品中表達(dá)值的變異系數(shù)排序。

#去除表達(dá)值全為0的行
#data4_nonzero <- data4[rowSums(data4)!=0,]

#篩選符合要求的表達(dá)的行和列
#data4_use <- data4[apply(data4,1,function(row) sum(row>=1)>=5),]
#data4_use <- data4[,apply(data4,2,function(col) sum(col>=1)>=1000),]
data4_use <- data4[rowSums(data4>=1)>5,colSums(data4>=1)>1000]

# 對(duì)于表達(dá)譜數(shù)據(jù)，因?yàn)樯婕暗絇CR的指數(shù)擴(kuò)增，一般會(huì)取log處理
# 其它數(shù)據(jù)log處理會(huì)降低數(shù)據(jù)之間的差異，不一定適用
data4_use_log2 <- log2(data4_use+1)

dim(data4_use_log2)

## [1] 16482   301

# 計(jì)算變異系數(shù)(標(biāo)準(zhǔn)差除以平均值)度量基因表達(dá)變化幅度
#cv <- apply(data4_use_log2,1,sd)/rowMeans(data4_use_log2)
# 根據(jù)變異系數(shù)排序
#data4_use_log2 <- data4_use_log2[order(cv,decreasing = T),]

# 計(jì)算中值絕對(duì)偏差 (MAD, median absolute deviation)度量基因表達(dá)變化幅度
# 在基因表達(dá)中，盡管某些基因很小的變化會(huì)導(dǎo)致重要的生物學(xué)意義，
# 但是很小的觀察值會(huì)引入很大的背景噪音，因此也意義不大。
mads <- apply(data4_use_log2, 1, mad)
data4_use_log2 <- data4_use_log2[rev(order(mads)),]

#篩選前3行
data_var3 <- data4_use_log2[1:3,]

# 轉(zhuǎn)置矩陣使得每一行為一個(gè)樣品，每一列為一組變量
data_var3_forPCA <- t(data_var3)

dim(data_var3_forPCA)

## [1] 301   3

kable(corner(data_var3_forPCA, r=10,c=5), booktabs=TRUE, 
        caption="A table of the 3 most variable genes")

A table of the 3 most variable genes

	MT2A	ANXA1	ARHGDIB
Hi_2338_1	12.211493	11.837198	9.543283
Hi_2338_2	11.306216	8.098769	8.071623
Hi_2338_3	11.926226	10.688626	9.720535
Hi_2338_4	10.974207	9.386574	9.883376
Hi_2338_5	14.603994	10.375072	9.970379
Hi_2338_6	6.904604	11.155349	10.093510
Hi_2338_7	12.436719	10.852249	7.742882
Hi_2338_8	9.798375	9.783227	6.270716
Hi_2338_9	11.743673	9.626476	9.250251
Hi_2338_10	11.240016	11.303056	0.000000

# 獲得樣品分組信息
sample <- rownames(data_var3_forPCA)

# 把樣品名字按 <_> 分割，取出其第二部分作為樣品的組名
# lapply(X, FUC) 對(duì)列表或向量中每個(gè)元素執(zhí)行FUC操作，F(xiàn)UNC為自定義或R自帶的函數(shù)
## One better way to generate group
group <- unlist(lapply(strsplit(sample, "_"), function(x) x[2]))

##One way to generate group
#sample_split <- strsplit(sample,"_")
#group <- matrix(unlist(sample_split), ncol=3, byrow=T)[,2]
print(sample[1:4])

## [1] "Hi_2338_1" "Hi_2338_2" "Hi_2338_3" "Hi_2338_4"

print(group[1:4])

## [1] "2338" "2338" "2338" "2338"

data_var3_scatter <- as.data.frame(data_var3_forPCA)
data_var3_scatter$group <- group
kable(corner(data_var3_scatter, r=10,c=5), booktabs=TRUE, 
    caption="A table of the 3 most variable genes")

A table of the 3 most variable genes

	MT2A	ANXA1	ARHGDIB	group
Hi_2338_1	12.211493	11.837198	9.543283	2338
Hi_2338_2	11.306216	8.098769	8.071623	2338
Hi_2338_3	11.926226	10.688626	9.720535	2338
Hi_2338_4	10.974207	9.386574	9.883376	2338
Hi_2338_5	14.603994	10.375072	9.970379	2338
Hi_2338_6	6.904604	11.155349	10.093510	2338
Hi_2338_7	12.436719	10.852249	7.742882	2338
Hi_2338_8	9.798375	9.783227	6.270716	2338
Hi_2338_9	11.743673	9.626476	9.250251	2338
Hi_2338_10	11.240016	11.303056	0.000000	2338

library(reshape2)
library(ggplot2)
data_var3_melt <- melt(data_var3_scatter, id.vars=c("group"))
kable(corner(data_var3_melt, r=10,c=5), booktabs=TRUE, 
        caption="A table of the 3 most variable genes in melted format")

A table of the 3 most variable genes in melted format

group	variable	value
2338	MT2A	12.211493
2338	MT2A	11.306216
2338	MT2A	11.926226
2338	MT2A	10.974207
2338	MT2A	14.603994
2338	MT2A	6.904604
2338	MT2A	12.436719
2338	MT2A	9.798375
2338	MT2A	11.743673
2338	MT2A	11.240016

ggplot(data_var3_melt, aes(factor(variable),value))+ylab("Gene expression")+
    geom_violin(aes(fill=factor(group)), stat="ydensity", position="dodge",
            scale="width", trim=TRUE) +xlab(NULL)

高顏值免費(fèi)在線繪圖工具升級(jí)版來(lái)了~~~

#ggplot(data_var3_melt, aes(factor(variable),value))+ylab("Gene expression")+ 
    #geom_quasirandom(aes(color=factor(group))) +xlab(NULL)

# 根據(jù)分組數(shù)目確定顏色變量
colorA <- rainbow(length(unique(group)))

# 根據(jù)每個(gè)樣品的分組信息獲取對(duì)應(yīng)的顏色變量
colors <- colorA[as.factor(group)]

# 根據(jù)樣品分組信息獲得legend的顏色
colorl <- colorA[as.factor(unique(group))]

# 獲得PCH symbol列表
pch_l <- as.numeric(as.factor(unique(group)))
# 產(chǎn)生每個(gè)樣品的pch symbol
pch <- pch_l[as.factor(group)]

scatterplot3d(data_var3_forPCA[,1:3], color=colors, pch=pch)
legend(0,10, legend=levels(as.factor(group)), col=colorl, pch=pch_l, xpd=T, horiz=F, ncol=6)

我們看到圖中的樣品并沒(méi)有按照預(yù)先設(shè)定的標(biāo)簽完全分開(kāi)。當(dāng)然我們也可以通過(guò)其他方法篩選變異最大的三個(gè)基因，最終的分類效果不會(huì)相差很大。因?yàn)椴还茉趺春Y選，我們都只用到了3個(gè)基因的表達(dá)量。

假如我們把這個(gè)數(shù)據(jù)用PCA來(lái)分類，結(jié)果是怎樣的呢？

# Pay attention to the format of PCA input 
# Rows are samples and columns are variables
data4_use_log2_t <- t(data4_use_log2)

# Add group column for plotting
data4_use_log2_label <- as.data.frame(data4_use_log2_t)
data4_use_log2_label$group <- group

# 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(data4_use_log2_t, 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))

## List of 5
##  $ sdev    : num [1:301] 36.6 30.4 23.3 21.6 19.9 ...
##  $ rotation: num [1:16482, 1:301] -0.01133 -0.01955 -0.00199 -0.00569 -0.0204 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:16482] "MT2A" "ANXA1" "ARHGDIB" "RND3" ...
##   .. ..$ : chr [1:301] "PC1" "PC2" "PC3" "PC4" ...
##  $ center  : Named num [1:16482] 6.5 5.47 5.43 4.23 4.1 ...
##   ..- attr(*, "names")= chr [1:16482] "MT2A" "ANXA1" "ARHGDIB" "RND3" ...
##  $ scale   : Named num [1:16482] 5.62 4.96 4.78 4.13 3.98 ...
##   ..- attr(*, "names")= chr [1:16482] "MT2A" "ANXA1" "ARHGDIB" "RND3" ...
##  $ x       : num [1:301, 1:301] -37.6 -28.6 -30.6 -43.7 -11.4 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:301] "Hi_2338_1" "Hi_2338_2" "Hi_2338_3" "Hi_2338_4" ...
##   .. ..$ : chr [1:301] "PC1" "PC2" "PC3" "PC4" ...
##  - attr(*, "class")= chr "prcomp"
## NULL

從圖中可以看到，數(shù)據(jù)呈現(xiàn)了一定的分類模式 (當(dāng)然這個(gè)分類結(jié)果也不理想，我們隨后再進(jìn)一步優(yōu)化)。

library(ggfortify)
autoplot(pca, data=data4_use_log2_label, colour="group") + 
    xlab(paste0("PC1 (", round(percentVar[1]*100), "% variance)")) + 
    ylab(paste0("PC2 (", round(percentVar[2]*100), "% variance)")) + theme_bw() + 
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
    theme(legend.position="right")

采用3個(gè)主成分獲得的分類效果優(yōu)于2個(gè)主成分，因?yàn)檫@樣保留的原始信息更多。

# 根據(jù)分組數(shù)目確定顏色變量
colorA <- rainbow(length(unique(group)))

# 根據(jù)每個(gè)樣品的分組信息獲取對(duì)應(yīng)的顏色變量
colors <- colorA[as.factor(group)]

# 根據(jù)樣品分組信息獲得legend的顏色
colorl <- colorA[as.factor(unique(group))]

# 獲得PCH symbol列表
pch_l <- as.numeric(as.factor(unique(group)))
# 產(chǎn)生每個(gè)樣品的pch symbol
pch <- pch_l[as.factor(group)]

pc <- as.data.frame(pca$x)
scatterplot3d(x=pc$PC1, y=pc$PC2, z=pc$PC3, pch=pch, color=colors, 
    xlab=paste0("PC1 (", round(percentVar[1]*100), "% variance)"), 
    ylab=paste0("PC2 (", round(percentVar[2]*100), "% variance)"), 
    zlab=paste0("PC3 (", round(percentVar[3]*100), "% variance)"))

legend(-3,8, legend=levels(as.factor(group)), col=colorl, pch=pch_l, xpd=T, horiz=F, ncol=6)

PCA的實(shí)現(xiàn)原理

在上面的例子中，PCA分析不是簡(jiǎn)單地選取2個(gè)或3個(gè)變化最大的基因，而是先把原始的變量做一個(gè)評(píng)估，計(jì)算各個(gè)變量各自的變異度(方差)和兩兩變量的相關(guān)度（協(xié)方差），得到一個(gè)協(xié)方差矩陣。在這個(gè)協(xié)方差矩陣中，對(duì)角線的值為每一個(gè)變量的方差，其它值為每?jī)蓚€(gè)變量的協(xié)方差。隨后對(duì)原變量的協(xié)方差矩陣對(duì)角化處理，即求解其特征值和特征向量。原變量與特征向量的乘積（對(duì)原始變量的線性組合）即為新變量（回顧下線性代數(shù)中的矩陣乘法）；新變量的協(xié)方差矩陣為對(duì)角協(xié)方差矩陣且對(duì)角線上的方差由大到小排列；然后從新變量中選擇信息最豐富也就是方差最大的的前2個(gè)或前3個(gè)新變量也就是主成分用以可視化。下面我們一步步闡釋這是怎么做的。

我們先回憶兩個(gè)數(shù)學(xué)概念，方差和協(xié)方差。方差用來(lái)表示一組一維數(shù)據(jù)的離散程度。協(xié)方差表示2組一維數(shù)據(jù)的相關(guān)性。當(dāng)協(xié)方差為0時(shí)，表示兩組數(shù)據(jù)完全獨(dú)立。當(dāng)協(xié)方差為正時(shí)，表示一組數(shù)據(jù)增加時(shí)另外一組也會(huì)增加；當(dāng)協(xié)方差為負(fù)時(shí)表示一組數(shù)據(jù)增加時(shí)另外一組數(shù)據(jù)會(huì)降低 （與相關(guān)系數(shù)類似）。如果我們有很多組一維數(shù)據(jù)，比如很多基因的表達(dá)數(shù)據(jù)，就會(huì)得到很多協(xié)方差，這就構(gòu)成了協(xié)方差矩陣。

假如我們有一個(gè)矩陣如下，

mat <- as.data.frame(matrix(rnorm(20,0,1), nrow=4))
colnames(mat) <- paste0("Gene_", letters[1:5])
rownames(mat) <- paste0("Samp_", 1:4)
mat

##            Gene_a     Gene_b     Gene_c     Gene_d     Gene_e
## Samp_1 -1.2207554  0.7608654 -0.2921542  0.2781680 -0.5515104
## Samp_2 -0.3855339 -1.3785382  1.1008845  0.5385477  0.1083430
## Samp_3  0.5828754  0.7443152  0.2221331  1.0248715  0.2456042
## Samp_4 -1.6524972 -1.6259796  0.8235352 -0.7186743  0.3295052

平均值中心化 (mean centering)：中心化數(shù)據(jù)使其平均值為0

# mean-centering data for columns
# Get mean-value matrix first
mat_mean_norm <- mat - rep(colMeans(mat),rep.int(nrow(mat),ncol(mat)))
mat_mean_norm

##            Gene_a    Gene_b     Gene_c       Gene_d      Gene_e
## Samp_1 -0.5517776  1.135700 -0.7557539 -0.002560247 -0.58449593
## Samp_2  0.2834438 -1.003704  0.6372848  0.257819506  0.07535752
## Samp_3  1.2518532  1.119149 -0.2414665  0.744143242  0.21261872
## Samp_4 -0.9835194 -1.251145  0.3599356 -0.999402500  0.29651969

# mean-centering using scale for columns
scale(mat, center=T, scale=F)

##            Gene_a    Gene_b     Gene_c       Gene_d      Gene_e
## Samp_1 -0.5517776  1.135700 -0.7557539 -0.002560247 -0.58449593
## Samp_2  0.2834438 -1.003704  0.6372848  0.257819506  0.07535752
## Samp_3  1.2518532  1.119149 -0.2414665  0.744143242  0.21261872
## Samp_4 -0.9835194 -1.251145  0.3599356 -0.999402500  0.29651969
## attr(,"scaled:center")
##      Gene_a      Gene_b      Gene_c      Gene_d      Gene_e 
## -0.66897778 -0.37483430  0.46359965  0.28072824  0.03298553

中位數(shù)中心化 (median centering)：如果數(shù)據(jù)變換范圍很大或有異常值，中位數(shù)標(biāo)準(zhǔn)化效果會(huì)更好。

# median-centering data for columns
mat_median_norm <- mat - rep(apply(mat,2,median),rep.int(nrow(mat),ncol(mat)))
mat_mean_norm

##            Gene_a    Gene_b     Gene_c       Gene_d      Gene_e
## Samp_1 -0.5517776  1.135700 -0.7557539 -0.002560247 -0.58449593
## Samp_2  0.2834438 -1.003704  0.6372848  0.257819506  0.07535752
## Samp_3  1.2518532  1.119149 -0.2414665  0.744143242  0.21261872
## Samp_4 -0.9835194 -1.251145  0.3599356 -0.999402500  0.29651969

我們可以計(jì)算Gene_a的方差為 0.973082 (var(mat$Gene_a))；Gene_a和Gene_b的協(xié)方差為0.5734631。

mat中5組基因的表達(dá)值的方差計(jì)算如下：

apply(mat,2,var)

##    Gene_a    Gene_b    Gene_c    Gene_d    Gene_e 
## 0.9730820 1.7050318 0.3883852 0.5396773 0.1601483

mat中5組基因表達(dá)值的協(xié)方差計(jì)算如下：

cov(mat)

##             Gene_a     Gene_b      Gene_c      Gene_d      Gene_e
## Gene_a  0.97308195  0.5734631 -0.01954727  0.66299331  0.10613531
## Gene_b  0.57346308  1.7050318 -0.73950788  0.60717437 -0.29082853
## Gene_c -0.01954727 -0.7395079  0.38838520 -0.12438895  0.18171565
## Gene_d  0.66299331  0.6071744 -0.12438895  0.53967732 -0.03906622
## Gene_e  0.10613531 -0.2908285  0.18171565 -0.03906622  0.16014830

如果均值為0，數(shù)值矩陣的協(xié)方差矩陣為矩陣的乘積 （實(shí)際上是矩陣的轉(zhuǎn)置與其本身的乘積除以變量的維數(shù)減1）。

# Covariance matrix for Mean normalized matrix
cov(mat_mean_norm)

##             Gene_a     Gene_b      Gene_c      Gene_d      Gene_e
## Gene_a  0.97308195  0.5734631 -0.01954727  0.66299331  0.10613531
## Gene_b  0.57346308  1.7050318 -0.73950788  0.60717437 -0.29082853
## Gene_c -0.01954727 -0.7395079  0.38838520 -0.12438895  0.18171565
## Gene_d  0.66299331  0.6071744 -0.12438895  0.53967732 -0.03906622
## Gene_e  0.10613531 -0.2908285  0.18171565 -0.03906622  0.16014830

# Covariance matrix for Mean normalized matrix 
# crossprod: matrix multiplication
crossprod(as.matrix(mat_mean_norm)) / (nrow(mat_mean_norm)-1)

##             Gene_a     Gene_b      Gene_c      Gene_d      Gene_e
## Gene_a  0.97308195  0.5734631 -0.01954727  0.66299331  0.10613531
## Gene_b  0.57346308  1.7050318 -0.73950788  0.60717437 -0.29082853
## Gene_c -0.01954727 -0.7395079  0.38838520 -0.12438895  0.18171565
## Gene_d  0.66299331  0.6071744 -0.12438895  0.53967732 -0.03906622
## Gene_e  0.10613531 -0.2908285  0.18171565 -0.03906622  0.16014830

# Use %*% for matrix multiplication (slower)
t(as.matrix(mat_mean_norm)) %*% as.matrix(mat_mean_norm) / (nrow(mat_mean_norm)-1)

##             Gene_a     Gene_b      Gene_c      Gene_d      Gene_e
## Gene_a  0.97308195  0.5734631 -0.01954727  0.66299331  0.10613531
## Gene_b  0.57346308  1.7050318 -0.73950788  0.60717437 -0.29082853
## Gene_c -0.01954727 -0.7395079  0.38838520 -0.12438895  0.18171565
## Gene_d  0.66299331  0.6071744 -0.12438895  0.53967732 -0.03906622
## Gene_e  0.10613531 -0.2908285  0.18171565 -0.03906622  0.16014830

用矩陣形式書寫如下，便于理解

根據(jù)前面的描述，原始變量的協(xié)方差矩陣表示原始變量自身的方差（協(xié)方差矩陣的主對(duì)角線位置）和原始變量之間的相關(guān)程度(非主對(duì)角線位置)。如果從這些數(shù)據(jù)中篩選主成分，則要選擇方差大(主對(duì)角線值大)，且與其它已選變量之間相關(guān)性最小的變量（非主對(duì)角線值很?。Ｈ绻@些原始變量之間毫不相關(guān)，則它們的協(xié)方差矩陣在除主對(duì)角線處外其它地方的值都為0，這種矩陣成為對(duì)角矩陣。

而做PCA分析就是產(chǎn)生一組新的變量，使得新變量的協(xié)方差矩陣為對(duì)角陣，滿足上面的要求。從而達(dá)到去冗余的目的。然后再選取方差大的變量，實(shí)現(xiàn)降維和去噪。

如果正向推導(dǎo)，這種組合可能會(huì)有很多種，一一計(jì)算會(huì)比較麻煩。那反過(guò)來(lái)看呢？ 我們不去尋找這種組合，而是計(jì)算如何使原變量的協(xié)方差矩陣變?yōu)閷?duì)角陣。

數(shù)學(xué)推導(dǎo)中謹(jǐn)記的兩個(gè)概念：

1. 假設(shè): 把未求解到的變量假設(shè)出來(lái)，用符號(hào)代替；這樣有助于思考和演算

2. 逆向：如果正向推導(dǎo)求不出，不妨倒著來(lái)；盡量多的利用已有信息

前面提到，新變量(Y_m,?k)是原始變量(X_m,?n)(原始變量的協(xié)方差矩陣為(C_n,?n))的線性組合，那么假設(shè)我們找到了這么一個(gè)線性組合(命名為特征矩陣(P_n,?k))，得到一組新變量Y_m,?k?=?X_m,?nP_n,?k，并且新變量的協(xié)方差矩陣(D_k,?k)為對(duì)角陣。那么這個(gè)特征矩陣(P_n,?k)需要符合什么條件呢？

從矩陣運(yùn)算可以看出，最終的特征矩陣(P_n,?k)需要把原變量協(xié)方差矩陣(C_n,?n)轉(zhuǎn)換為對(duì)角陣(因?yàn)樾伦兞康膮f(xié)方差矩陣(D_k,?k)為對(duì)角陣)，并且對(duì)角元素從大到小排列（保證每個(gè)主成分的貢獻(xiàn)度依次降低）。

現(xiàn)在就把求解新變量的任務(wù)轉(zhuǎn)變?yōu)榱饲蠼庠兞繀f(xié)方差矩陣的對(duì)角化問(wèn)題了。在線性代數(shù)中，矩陣對(duì)角化的問(wèn)題就是求解矩陣的特征值和特征向量的問(wèn)題。

我們舉一個(gè)例子講述怎么求解特征值和特征向量。

假設(shè)A_n,?n為n階對(duì)稱陣，如存在λ和非零向量x，使得A**x?=?λ**x，則稱λ為矩陣A_n,?n的特征值，非零向量x為為矩陣A_n,?n對(duì)應(yīng)于特征值λ的特征向量。

根據(jù)這個(gè)定義可以得出(A_n,?n???λ**E)x?=?0，由于x為非零向量，所以行列式|A???λ**E|?=?0。

由此求解出n個(gè)根λ₁,?λ₂,?…,?λ₃就是矩陣A的特征值。

回顧下行列式的計(jì)算：

• 行列式的值為行列式第一列的每一個(gè)數(shù)乘以它的余子式（余子式是行列式中除去當(dāng)前元素所在行和列之后剩下的行列式）。

• 當(dāng)行列式中存在線性相關(guān)的行或列或者有一行或一列元素全為0時(shí)，行列式的值為0。

• 上三角形行列式的值為其主對(duì)角線上元素的乘積。

• 互換行列式的兩行或兩列，行列式變號(hào)。

• 行列式的某一列（行）乘以同意書加到另一列（列）對(duì)應(yīng)元素上去，行列式不變。

簡(jiǎn)單的PCA實(shí)現(xiàn)

我們使用前面用到的數(shù)據(jù)data3來(lái)演示下如何用R函數(shù)實(shí)現(xiàn)PCA的計(jì)算，并與R中自帶的prcomp做個(gè)比較。

library(knitr)
kable(headTail(data3), booktabs=T, caption="Expression profile for 3 genes in 100 samples")

Expression profile for 3 genes in 100 samples

	Gene1	Gene2	Gene3	Group
a1	5.99	5.35	5.14	a
a2	4.66	5.31	5.05	a
a3	4.92	5.12	4.88	a
a4	4.79	5.11	4.8	a
…	…	…	…	NA
b47	20.78	4.95	4.91	b
b48	20.5	4.9	5.11	b
b49	20.46	4.85	4.95	b
b50	19.94	4.87	5.18	b

標(biāo)準(zhǔn)化數(shù)據(jù)

data3_center_scale <- scale(data3[,1:3], center=T, scale=T)
kable(headTail(data3_center_scale), booktabs=T, 
        caption="Normalized expression for 3 genes in 100 samples")

Normalized expression for 3 genes in 100 samples

	Gene1	Gene2	Gene3
a1	-0.85	1.83	-0.96
a2	-1.03	1.66	-0.98
a3	-0.99	0.72	-1.01
a4	-1.01	0.67	-1.03
…	…	…	…
b47	1.09	-0.11	-1.01
b48	1.06	-0.38	-0.97
b49	1.05	-0.61	-1
b50	0.98	-0.52	-0.95

計(jì)算協(xié)方差矩陣

data3_center_scale_cov <- cov(data3_center_scale)
kable(data3_center_scale_cov, booktabs=T, 
        caption="Covariance matrix for 3 genes in 100 samples")

Covariance matrix for 3 genes in 100 samples

	Gene1	Gene2	Gene3
Gene1	1.0000000	0.0255834	-0.0087919
Gene2	0.0255834	1.0000000	0.0553298
Gene3	-0.0087919	0.0553298	1.0000000

求解特征值和特征向量

data3_center_scale_cov_eigen <- eigen(data3_center_scale_cov)

# 特征值，從大到小排序
data3_center_scale_cov_eigen$values

## [1] 1.0580005 1.0066389 0.9353605

# 特征向量, 每一列為對(duì)應(yīng)特征值的特征向量
data3_center_scale_cov_eigen$vectors

##           [,1]        [,2]       [,3]
## [1,] 0.2192050  0.90784568  0.3574429
## [2,] 0.7223522  0.09525968 -0.6849327
## [3,] 0.6558631 -0.40834033  0.6349030

產(chǎn)生新的矩陣

pc_select = 3
label = paste0("PC",c(1:pc_select))
data3_new <- data3_center_scale %*% data3_center_scale_cov_eigen$vectors[,1:pc_select]
colnames(data3_new) <- label
kable(headTail(data3_new), booktabs=T, 
        caption="PCA generated matrix for the expression of 3 genes in 100 samples")

PCA generated matrix for the expression of 3 genes in 100 samples

	PC1	PC2	PC3
a1	0.51	-0.21	-2.17
a2	0.33	-0.37	-2.12
a3	-0.36	-0.42	-1.49
a4	-0.41	-0.43	-1.47
…	…	…	…
b47	-0.5	1.39	-0.17
b48	-0.68	1.32	0.03
b49	-0.86	1.31	0.16
b50	-0.79	1.23	0.1

比較原始數(shù)據(jù)和新產(chǎn)生的主成分對(duì)樣品的聚類

#library(scatterplot3d)
colorl <- c("#E69F00", "#56B4E9")
# Extract same number of colors as the Group and same Group would have same color.
colors <- colorl[as.numeric(data3$Group)]

# 1 row 2 columns
par(mfrow=c(1,2))

scatterplot3d(data3[,1:3], color=colors, angle=55, pch=16, main="Original data")
legend("top", legend=levels(data3$Group), col=colorl, pch=16, xpd=T, horiz=T)

scatterplot3d(data3_new, color=colors,angle=55, pch=16, main="Principle components")
legend("top", legend=levels(data3$Group), col=colorl, pch=16, xpd=T, horiz=T)

#par(mfrow=c(1,1))

利用prcomp進(jìn)行主成分分析

pca_data3 <- prcomp(data3[,1:3], center=TRUE, scale=TRUE)

#Show whats in the result returned by prcomp
str(pca_data3)

## List of 5
##  $ sdev    : num [1:3] 1.029 1.003 0.967
##  $ rotation: num [1:3, 1:3] 0.2192 0.7224 0.6559 -0.9078 -0.0953 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:3] "Gene1" "Gene2" "Gene3"
##   .. ..$ : chr [1:3] "PC1" "PC2" "PC3"
##  $ center  : Named num [1:3] 12.46 4.98 10
##   ..- attr(*, "names")= chr [1:3] "Gene1" "Gene2" "Gene3"
##  $ scale   : Named num [1:3] 7.602 0.202 5.06
##   ..- attr(*, "names")= chr [1:3] "Gene1" "Gene2" "Gene3"
##  $ x       : num [1:100, 1:3] 0.506 0.331 -0.36 -0.409 -0.27 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:100] "a1" "a2" "a3" "a4" ...
##   .. ..$ : chr [1:3] "PC1" "PC2" "PC3"
##  - attr(*, "class")= chr "prcomp"

# 新的數(shù)據(jù)，與前面計(jì)算的抑制
data3_pca_new <- pca_data3$x
kable(headTail(data3_pca_new), booktabs=T, 
        caption="PCA generated matrix usig princomp for the expression of 3 genes in 100 samples")

PCA generated matrix usig princomp for the expression of 3 genes in 100 samples

	PC1	PC2	PC3
a1	0.51	0.21	-2.17
a2	0.33	0.37	-2.12
a3	-0.36	0.42	-1.49
a4	-0.41	0.43	-1.47
…	…	…	…
b47	-0.5	-1.39	-0.17
b48	-0.68	-1.32	0.03
b49	-0.86	-1.31	0.16
b50	-0.79	-1.23	0.1

# 特征向量，與我們前面計(jì)算的一致(特征向量的符號(hào)是任意的)
pca_data3$rotation

##             PC1         PC2        PC3
## Gene1 0.2192050 -0.90784568  0.3574429
## Gene2 0.7223522 -0.09525968 -0.6849327
## Gene3 0.6558631  0.40834033  0.6349030

比較手動(dòng)實(shí)現(xiàn)的PCA與prcomp實(shí)現(xiàn)的PCA的聚類結(jié)果

#library(scatterplot3d)
colorl <- c("#E69F00", "#56B4E9")
# Extract same number of colors as the Group and same Group would have same color.
colors <- colorl[as.numeric(data3$Group)]

# 1 row 2 columns
par(mfrow=c(1,2))

scatterplot3d(data3_new, color=colors,angle=55, pch=16, main="PCA by steps")
legend("top", legend=levels(data3$Group), col=colorl, pch=16, xpd=T, horiz=T)

scatterplot3d(data3_pca_new, color=colors,angle=55, pch=16, main="PCA by prcomp")
legend("top", legend=levels(data3$Group), col=colorl, pch=16, xpd=T, horiz=T)

#par(mfrow=c(1,1))

自定義PCA計(jì)算函數(shù)

ct_PCA <- function(data, center=TRUE, scale=TRUE){
  data_norm <- scale(data, center=center, scale=scale)
  data_norm_cov <- crossprod(as.matrix(data_norm)) / (nrow(data_norm)-1)
  data_eigen <- eigen(data_norm_cov)

  rotation <- data_eigen$vectors
  label <- paste0('PC', c(1:ncol(rotation)))
  colnames(rotation) <- label
  sdev <- sqrt(data_eigen$values)
  data_new <- data_norm %*% rotation
  colnames(data_new) <- label
  ct_pca <- list('rotation'=rotation, 'x'=data_new, 'sdev'=sdev)
  return(ct_pca)
}

比較有無(wú)scale對(duì)聚類的影響，從圖中可以看到，如果不對(duì)數(shù)據(jù)進(jìn)行scale處理，樣品的聚類結(jié)果更像原始數(shù)據(jù)，本身數(shù)值大的基因?qū)χ鞒煞值呢暙I(xiàn)會(huì)大。如果關(guān)注的是每個(gè)變量自身的實(shí)際方差對(duì)樣品分類的貢獻(xiàn)，則不應(yīng)該SCALE；如果關(guān)注的是變量的相對(duì)大小對(duì)樣品分類的貢獻(xiàn)，則應(yīng)該SCALE，以防數(shù)值高的變量導(dǎo)入的大方差引入的偏見(jiàn)。

data3_pca_noscale_step = ct_PCA(data3[,1:3], center=TRUE, scale=FALSE)

# 特征向量
data3_pca_noscale_step$rotation

##                PC1         PC2           PC3
## [1,]  0.9999446062 0.010502677 -0.0006915646
## [2,]  0.0006682513 0.002223103  0.9999973056
## [3,] -0.0105041862 0.999942374 -0.0022159613

# 新變量
data3_pca_noscale_pc <- data3_pca_noscale_step$x

#library(scatterplot3d)
colorl <- c("#E69F00", "#56B4E9")
# Extract same number of colors as the Group and same Group would have same color.
colors <- colorl[as.numeric(data3$Group)]

# 1 row 2 columns
par(mfrow=c(2,2))

scatterplot3d(data3[,c(1,3,2)], color=colors, angle=55, pch=16, main="Original data")
legend("top", legend=levels(data3$Group), col=colorl, pch=16, xpd=T, horiz=T)

scatterplot3d(data3_pca_noscale_pc, color=colors,angle=55, pch=16, main="PCA (no scale)")
legend("top", legend=levels(data3$Group), col=colorl, pch=16, xpd=T, horiz=T)

scatterplot3d(data3_center_scale[,c(1,3,2)], color=colors, angle=55, pch=16, 
        main="Original data (scale)")
legend("top", legend=levels(data3$Group), col=colorl, pch=16, xpd=T, horiz=T)

scatterplot3d(data3_new, color=colors,angle=55, pch=16, main="PCA (scale)")
legend("top", legend=levels(data3$Group), col=colorl, pch=16, xpd=T, horiz=T)

#par(mfrow=c(1,1))

PCA結(jié)果解釋

prcomp函數(shù)會(huì)返回主成分的標(biāo)準(zhǔn)差、特征向量和主成分構(gòu)成的新矩陣。接下來(lái)，探索下不同主成分對(duì)數(shù)據(jù)差異的貢獻(xiàn)和主成分與原始變量的關(guān)系。

• 主成分的平方為為特征值,其含義為每個(gè)主成分可以解釋的數(shù)據(jù)差異，計(jì)算方式為eigenvalues = (pca$sdev)^2

• 每個(gè)主成分可以解釋的數(shù)據(jù)差異的比例為percent_var = eigenvalues*100/sum(eigenvalues)

• 可以使用summary(pca)獲取以上兩條信息。

這兩個(gè)信息可以判斷主成分分析的質(zhì)量：

• 成功的降維需要保證在前幾個(gè)為數(shù)不多的主成分對(duì)數(shù)據(jù)差異的解釋可以達(dá)到80-90%。

指導(dǎo)選擇主成分的數(shù)目：

• 選擇的主成分足以解釋的總方差大于80% (方差比例碎石圖)

• 從前面的協(xié)方差矩陣可以看到，自動(dòng)定標(biāo)(scale)的變量的方差為1 (協(xié)方差矩陣對(duì)角線的值)。待選擇的主成分應(yīng)該是那些方差大于1的主成分，即其解釋的方差大于原始變量（特征值碎石圖，方差大于1，特征值也會(huì)大于1，反之亦然）。

鑒定核心變量和變量間的隱性關(guān)系:

• 原始變量與主成分的相關(guān)性Variable correlation with PCs (var.cor) = loadings * sdev

• 原始數(shù)據(jù)對(duì)主成分的貢獻(xiàn)度 var.cor^2 / (total var.cor^2)

在測(cè)試數(shù)據(jù)中，scale后，三個(gè)主成分對(duì)數(shù)據(jù)差異的貢獻(xiàn)度大都在30%左右，而未scale的數(shù)據(jù)，三個(gè)主成分對(duì)數(shù)據(jù)差異的貢獻(xiàn)度相差很大。這是因?yàn)槿齻€(gè)基因由于自身表達(dá)量級(jí)所引起的方差的差異導(dǎo)致它們各自對(duì)數(shù)據(jù)的權(quán)重差異，從而使主成分偏向于數(shù)值大的變量。

PCA應(yīng)用于測(cè)試數(shù)據(jù)

前面用到一組比較大的測(cè)試數(shù)據(jù)集，并做了PCA分析，現(xiàn)在測(cè)試不同的處理對(duì)結(jié)果的影響。

首先回顧下我們用到的數(shù)據(jù)。

library("gridExtra")
# Pay attention to the format of PCA input 
# Rows are samples and columns are variables
# data4_use_log2_t <- t(data4_use_log2)

# Add group column for plotting
# data4_use_log2_label <- as.data.frame(data4_use_log2_t)
# data4_use_log2_label$group <- group

kable(corner(data4_use_log2_label), digits=3, caption="Single cell gene expression data")

Single cell gene expression data

	MT2A	ANXA1	ARHGDIB	RND3	BLVRB
Hi_2338_1	12.211	11.837	9.543	9.762	9.162
Hi_2338_2	11.306	8.099	8.072	10.107	9.423
Hi_2338_3	11.926	10.689	9.721	9.388	9.694
Hi_2338_4	10.974	9.387	9.883	8.792	9.767
Hi_2338_5	14.604	10.375	9.970	8.815	9.350

比較對(duì)數(shù)運(yùn)算和scale對(duì)樣品分類的影響。

ct_pca_2d_plot <- function(pca, data_with_label, labelName='group', title='PCA') {
  # sdev: standard deviation of the principle components.
  # Square to get variance
  percentVar <- pca$sdev^2 / sum( pca$sdev^2)
  #data <- data_with_label
  #data[labelName] <- factor(unlist(data[labelName]))
  level <- length(unique(unlist(data_with_label[labelName])))
  shapes = (1:level)%%30  # maximum allow 30 types of symbols
  p = autoplot(pca, data=data_with_label, colour=labelName, shape=labelName) + 
      scale_shape_manual(values=shapes) +
      xlab(paste0("PC1 (", round(percentVar[1]*100), "% variance)")) + 
      ylab(paste0("PC2 (", round(percentVar[2]*100), "% variance)")) + 
      theme_bw() + theme(legend.position="right") + labs(title=title) +
      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
  return(p)
}

# 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.
data4_use_t <- t(data4_use)
ori_scale_pca_test <- prcomp(data4_use_t, scale=T)
ori_no_scale_pca_test <- prcomp(data4_use_t, scale=F)
log_scale_pca_test <- prcomp(data4_use_log2_t, scale=T)
log_no_scale_pca_test <- prcomp(data4_use_log2_t, scale=F)

ori_scale_pca_plot = ct_pca_2d_plot(ori_scale_pca_test, data4_use_log2_label, 
        title="Scaled original data")
ori_no_scale_pca_plot = ct_pca_2d_plot(ori_no_scale_pca_test, data4_use_log2_label, 
        title="Un-scaled original data")
log_scale_pca_plot = ct_pca_2d_plot(log_scale_pca_test, data4_use_log2_label, 
        title="Scaled log transformed data")
log_no_scale_pca_plot = ct_pca_2d_plot(log_no_scale_pca_test, data4_use_log2_label, 
        title="Un-scaled log transformed data")

a__ <- grid.arrange(ori_scale_pca_plot, ori_no_scale_pca_plot, log_scale_pca_plot, 
        log_no_scale_pca_plot, ncol=2)

如果首先提取500個(gè)變化最大的基因，再執(zhí)行PCA分析會(huì)怎樣呢？

data4_use_mad <- apply(data4_use, 1, mad)
data4_use_mad_top500 <- t(data4_use[rev(order(data4_use_mad))[1:500],])

data4_use_log2_mad <- apply(data4_use_log2, 1, mad)
data4_use_log2_mad_top500 <- t(data4_use_log2[rev(order(data4_use_log2_mad))[1:500],])

ori_scale_pca_top500 <- prcomp(data4_use_mad_top500, scale=T)
ori_no_scale_pca_top500 <- prcomp(data4_use_mad_top500, scale=F)
log_scale_pca_top500 <- prcomp(data4_use_log2_mad_top500, scale=T)
log_no_scale_pca_top500 <- prcomp(data4_use_log2_mad_top500, scale=F)

ori_scale_pca_plot_t5 = ct_pca_2d_plot(ori_scale_pca_top500, data4_use_log2_label, 
        title="Scaled original data")
ori_no_scale_pca_plot_t5 = ct_pca_2d_plot(ori_no_scale_pca_top500, 
        data4_use_log2_label, title="Un-scaled original data")
log_scale_pca_plot_t5 = ct_pca_2d_plot(log_scale_pca_top500, 
        data4_use_log2_label, title="Scaled log transformed data")
log_no_scale_pca_plot_t5 = ct_pca_2d_plot(log_no_scale_pca_top500, 
        data4_use_log2_label, title="Un-scaled log transformed data")

a__ <- grid.arrange(ori_scale_pca_plot_t5, ori_no_scale_pca_plot_t5, 
                    log_scale_pca_plot_t5, log_no_scale_pca_plot_t5, ncol=2)

### PCA注意事項(xiàng)

1. 一般說(shuō)來(lái)，在PCA之前原始數(shù)據(jù)需要中心化（centering，數(shù)值減去平均值）。中心化的方法很多，除了平均值中心化（mean-centering）外，還包括其它更穩(wěn)健的方法，比如中位數(shù)中心化等。

2. 除了中心化以外，定標(biāo) (Scale, 數(shù)值除以標(biāo)準(zhǔn)差) 也是數(shù)據(jù)前處理中需要考慮的一點(diǎn)。如果數(shù)據(jù)沒(méi)有定標(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倍。

3. 但是定標(biāo)(scale)可能會(huì)有一些負(fù)面效果，因?yàn)槎?biāo)后變量之間的權(quán)重就是變得相同。如果我們的變量中有噪音的話，我們就在無(wú)形中把噪音和信息的權(quán)重變得相同，但PCA本身無(wú)法區(qū)分信號(hào)和噪音。在這樣的情形下，我們就不必做定標(biāo)。

4. 一般而言，對(duì)于度量單位不同的指標(biāo)或是取值范圍彼此差異非常大的指標(biāo)不直接由其協(xié)方差矩陣出發(fā)進(jìn)行主成分分析，而應(yīng)該考慮對(duì)數(shù)據(jù)的標(biāo)準(zhǔn)化。比如度量單位不同，有萬(wàn)人、萬(wàn)噸、萬(wàn)元、億元，而數(shù)據(jù)間的差異性也非常大，小則幾十大則幾萬(wàn)，因此在用協(xié)方差矩陣求解主成分時(shí)存在協(xié)方差矩陣中數(shù)據(jù)的差異性很大。在后面提取主成分時(shí)發(fā)現(xiàn)，只提取了一個(gè)主成分，而此時(shí)并不能將所有的變量都解釋到，這就沒(méi)有真正起到降維的作用。此時(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)化的過(guò)程實(shí)際上就是抹殺原有變量離散程度差異的過(guò)程，標(biāo)準(zhǔn)化后方差均為1，而實(shí)際上方差是對(duì)數(shù)據(jù)信息的重要概括形式，也就是說(shuō)，對(duì)原始數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化后抹殺了一部分重要信息，因此才使得標(biāo)準(zhǔn)化后各變量在主成分構(gòu)成中的作用趨于相等。因此，對(duì)同度量或是取值范圍在同量級(jí)的數(shù)據(jù)還是直接使用非定標(biāo)數(shù)據(jù)求解主成分為宜。

5. 中心化和定標(biāo)都會(huì)受數(shù)據(jù)中離群值（outliers）或者數(shù)據(jù)不均勻（比如數(shù)據(jù)被分為若干個(gè)小組）的影響，應(yīng)該用更穩(wěn)健的中心化和定標(biāo)方法。

6. PCA也存在一些限制，例如它可以很好的解除線性相關(guān)，但是對(duì)于高階相關(guān)性就沒(méi)有辦法了，對(duì)于存在高階相關(guān)性的數(shù)據(jù)，可以考慮Kernel PCA，通過(guò)Kernel函數(shù)將非線性相關(guān)轉(zhuǎn)為線性相關(guān)，關(guān)于這點(diǎn)就不展開(kāi)討論了。另外，PCA假設(shè)數(shù)據(jù)各主特征是分布在正交方向上，如果在非正交方向上存在幾個(gè)方差較大的方向，PCA的效果就大打折扣了。

參考資料

• https://www.zhihu.com/question/20998460

• PCA 教程1

• PCA 文字化描述

• pca1

• ggplot2 axis

• scatterplot3D

• 穩(wěn)健PCA

• http://www.nlpca.org/pca_principal_component_analysis.html

• Data centering

• Sample R markdown

• 矩陣特征值，對(duì)稱矩陣的對(duì)角化

• Detail usage and visualization of prcomp result

• ggplot2 side by side plot

PCA主成分分析實(shí)戰(zhàn)和可視化 | 附R代碼和測(cè)試數(shù)據(jù)

用了這么多年的PCA可視化竟然是錯(cuò)的?。?！

這篇Nature子刊文章的蛋白組學(xué)數(shù)據(jù)PCA分析竟花費(fèi)了我兩天時(shí)間來(lái)重現(xiàn)|附全過(guò)程代碼

復(fù)現(xiàn)nature communication PCA原圖|代碼分析（一）

本文節(jié)選自：這個(gè)為生信學(xué)習(xí)和生信作圖打造的開(kāi)源R教程真香?。?！