转录组表达矩阵为什么需要主成分分析以及怎么做
PCA大约是198x年提出的,是一种数据降维的方法。
简单将:PCA就是将高维数据通过线性变换投影到低维空间上。
PCA的指导思想是:找出最能够代表原始数据的投影方法。
问题的提出
在研究生物学问题,常常希望把更多的特征(如转录组产生的表达矩阵)纳入数学模型,这些基因的表达往往存在相关性。例如人的编码基因大约有2万多,这些庞大的表达特征矩阵增加了问题的复杂性。
如果能将多个相关的特征综合为少数几个具有代表性的特征:
既能够代表原始特征的绝大多数信息;
组合后的特征互不相关,降低特征的相关性(冗余);
这些特征就是主成分
寻找这些少数具有代表性特征的过程被称为主成分分析。
在降低特征矩阵复杂性的同时,希望降维后的数据不能失真,只去除噪声或冗余的数据
噪音污染干扰了想听到的真正声音。假设样本中某个主要的维度1(基因的表达),是想要听到的“真正声音”,本来含有的“信息强度”(方差)是很大的,但由于维度1和其他维度(相关的基因表达)有着千丝万缕的相关性,受到这些相关维度的干扰,维度1所包含的信息强度被削弱了。希望能通过PCA处理,是维度1与其他维度的相关性尽可能减弱,从而恢复维度1的应有的信息,让模型“听的更清晰”。
冗余就是多余,有没有不影响结果。在样本中有些维度,在所有的样本中的变化都不明显(有些基因的表达在不同样本中没有差异),极端时在所有样本中该维度的值都相等,该维度的方差接近于零。这些维度对区分不同的样本起不到丝毫作用,通过PCA去掉这些维度。
PCA实现的数学基础
想要降噪和去冗余,首先的首先,需要将这两种标准用数值表示,即计算各维度间的相关性和方差。相关性越高,噪声越大,方差越小,冗余性越高。
有没有一种数据指标能够同时描述这两种信息呢?
协方差矩阵可以度量维度与维度之间的关系,矩阵对角线上的值是各个维度上的方差(信息),其他值是两两维度间的协方差(相关性)。
让不同维度之间相关性最小,即让协方差矩阵中的非对角线元素基本为零--矩阵对角化(线性代数)。
对角化后的矩阵,对角线上是协方差矩阵的特征值,既可表示各维度的新方差,又可代表各维度去躁后本身该具有的信息强度,完成了去躁声的工作。
对角化后的矩阵,对角线上较小的新方差就是要去除的冗余维度,只保留较大的新方差(信息)的维度,完成去冗余的工作。
PCA的步骤(公式推导)
1.形成样本矩阵,样本中心化
假设一个样本集X,里面有N个样本,每个样本的维度为d
将这些样本组织成样本矩阵的形式,即每行为一个样本,每一列为一个维度(如基因表达量),得到样本矩阵S:
将样本矩阵进行中心化,即保证每个维度的均值为零,让矩阵的每一列除以减去对应的均值即可。
2.计算样本矩阵的协方差矩阵
3.对协方差矩阵进行特征值分解,选取最大的p个特征值对应的特征向量组成投影矩阵
对角化协方差矩阵C,矩阵C是对称矩阵,对称矩阵对角化就是找到一个正交矩阵P。
我相信你肯定会看不下去了这些公式,所以建议你读之前我们介绍过:一文看懂主成分分析,跟着一步步代码的形式理解PCA的过程即可,空间想象力不是每个人都有的!
R实现简单的PCA分析
R包含有很多实现PCA分析的函数,区别主要在于特征值的分解方法不同。
# 使用对角化的协方差矩阵进行降维
> set.seed(0.1234)
> S <- matrix(sample(1:50, 30,replace = TRUE), nrow = 10, ncol = 3)
> S
[,1] [,2] [,3]
[1,] 14 21 41
[2,] 4 21 25
[3,] 39 42 46
[4,] 1 46 37
[5,] 34 10 37
[6,] 23 7 34
[7,] 43 9 42
[8,] 14 15 25
[9,] 18 21 44
[10,] 33 37 15
> C <- cov(S)
> C
[,1] [,2] [,3]
[1,] 211.56667 -38.633333 30.244444
[2,] -38.63333 198.100000 -7.044444
[3,] 30.24444 -7.044444 99.377778
> E <- eigen(C)
> E
eigen() decomposition
$values
[1] 249.24458 168.12906 91.67081
$vectors
[,1] [,2] [,3]
[1,] 0.7718771 0.5833028 -0.2529104
[2,] -0.6084512 0.7931044 -0.0277946
[3,] 0.1843717 0.1753377 0.9670904
> S1 <- S %*% E$vectors
> S1
[,1] [,2] [,3]
[1,] 5.588043 32.01027 35.526273
[2,] -5.080675 23.37184 22.581931
[3,] 13.029353 64.12472 33.455278
[4,] -20.395126 43.55360 34.250882
[5,] 26.981061 34.25083 26.905444
[6,] 19.762652 24.92917 26.869571
[7,] 35.458264 39.58414 29.492497
[8,] 6.288803 24.44625 20.219595
[9,] 9.228666 34.86950 37.415903
[10,] 5.724824 51.22392 5.131912
## 函数princomp()的频谱分解方法
> P1 <- princomp(S)
> P1$scores
Comp.1 Comp.2 Comp.3
[1,] -4.0705437 -5.226150 -8.3413444
[2,] -14.7392614 -13.864580 4.6029976
[3,] 3.3707665 26.888299 -6.2703497
[4,] -30.0537124 6.317173 -7.0659530
[5,] 17.3224743 -2.985594 0.2794844
[6,] 10.1040650 -12.307251 0.3153573
[7,] 25.7996776 2.347714 -2.3075686
[8,] -3.3697834 -12.790179 6.9653339
[9,] -0.4299204 -2.366926 -10.2309740
[10,] -3.9337621 13.987494 22.0530165
## 函数prcomp()使用奇异值分解(SVD)进行PCA
> P1 <- princomp(S)
> P1$scores
Comp.1 Comp.2 Comp.3
[1,] -4.0705437 -5.226150 -8.3413444
[2,] -14.7392614 -13.864580 4.6029976
[3,] 3.3707665 26.888299 -6.2703497
[4,] -30.0537124 6.317173 -7.0659530
[5,] 17.3224743 -2.985594 0.2794844
[6,] 10.1040650 -12.307251 0.3153573
[7,] 25.7996776 2.347714 -2.3075686
[8,] -3.3697834 -12.790179 6.9653339
[9,] -0.4299204 -2.366926 -10.2309740
[10,] -3.9337621 13.987494 22.0530165
## FactoMineR的PCA函数的奇异值分解(SVD)进行PCA
> P3 <- PCA(S)
> P3$svd
$vs
[1] 1.1434171 0.9747107 0.8617055
$U
[,1] [,2] [,3]
[1,] 0.02711376 0.3436771 -0.98936747
[2,] -1.19988150 -0.8303852 -0.60529318
[3,] 0.65511153 1.9414306 1.05205820
[4,] -1.55185618 1.4638329 -0.58679291
[5,] 1.04510421 -0.5465876 0.08082268
[6,] 0.52340656 -0.9510399 -0.52810832
[7,] 1.71281093 -0.2285315 0.28689861
[8,] -0.57453078 -1.1590963 -0.21489903
[9,] 0.34808947 0.5668166 -0.92957243
[10,] -0.98536801 -0.6001167 2.43425385
$V
[,1] [,2] [,3]
[1,] 0.6750282 0.01900612 0.7375471
[2,] -0.5026055 0.74367327 0.4408376
[3,] 0.5401155 0.66827302 -0.5115530
总结
PCA最终目的就是降噪和消灭冗余维度,以使降低维度的同时保存数据原有的特征不失真。
PCA常用数学方法是协方差矩阵对角化和奇异值分解。
PCA只是一种常用的降维方法,针对不同的数据集,应当选取适合的降维方法来得到最优的结果。
常用的降维方法还有LDA,tSNE等,下面我们简单看看代码的区别
PCA和tSNE
首先加载必备的包:
set.seed(123456789)
#set.seed()产生随机数
#用于设定随机数种子,一个特定的种子可以产生一个特定的伪随机序列,这个函数的主要目的,
#是让你的模拟能够可重复出现,因为很多时候我们需要取随机数,但这段代码再跑一次的时候,
#结果就不一样了,如果需要重复出现同样的模拟结果的话,就可以用set.seed()。
library(pheatmap)
library(Rtsne)
library(ggfortify)
library(mvtnorm)
下面代码演示一个随机矩阵生成
ng=500
nc=20
a1=rnorm(ng*nc);dim(a1)=c(ng,nc)
a2=rnorm(ng*nc)+3;dim(a2)=c(ng,nc)
a3=cbind(a1,a2)
colnames(a3)=c(paste0('cell_01_',1:nc),paste0('cell_02_',1:nc))
rownames(a3)=paste('gene_',1:ng,sep = '')
pheatmap(a3)
可以看到,很明显,我们的细胞可以区分成为2类,而且并不需要那么多基因来区分他们,所以我们需要PCA分析,找到决定性的因素,通常是前几个主成分即可。
a3=t(a3);dim(a3) ## PCA分析
pca_dat <- prcomp(a3, scale. = TRUE)
p=autoplot(pca_dat) + theme_classic() + ggtitle('PCA plot')
print(p)
# 这个时候细胞被区分开,而且是很明显的一个主成分。
人造矩阵,使用tSNE当然也可以区分,这里看不出tSNE的优势,不过在单细胞领域,就不一样了,很多时候都是PCA区分效果不行,但是tSNE就很棒。
set.seed(42)
tsne_out <- Rtsne(a3,pca=FALSE,perplexity=10,theta=0.0) # Run TSNE
tsnes=tsne_out$Y
colnames(tsnes) <- c("tSNE1", "tSNE2")
ggplot(tsnes, aes(x = tSNE1, y = tSNE2))+ geom_point()
其实我们还留了一个悬念,就是前面的主成分,到底是由哪些基因组成呢?那些基因的重要性在该主成分的比例如何呢?