首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >转录组表达矩阵为什么需要主成分分析以及怎么做

转录组表达矩阵为什么需要主成分分析以及怎么做

作者头像
生信技能树
发布2019-08-13 15:08:06
7.3K0
发布2019-08-13 15:08:06
举报
文章被收录于专栏:生信技能树生信技能树

我们阅读量破万的综述:RNA-seq这十年(3万字长文综述)给粉丝朋友们带来了很多理解上的挑战:

所以我们开辟专栏慢慢介绍其中的一些概念性的问题,上一期:表达矩阵的归一化和标准化,去除极端值,异常值

本期我们介绍一下主成分分析,即PCA,等降维算法,因为通常RNA-seq都是测几万个基因的表达量,但实际上病人与正常人的区别并不需要这么多基因,可能仅仅是关键的几个基因的组合,就足够区分了,那么以后临床应用就简单了,仅仅是关注重要的的主成分所构成的基因即可。

其实之前我们介绍过一文看懂主成分分析,大家可以比较一下本文和之前的介绍,下面是正文:

PCA的步骤及解释

PCA大约是198x年提出的,是一种数据降维的方法。 简单将:PCA就是将高维数据通过线性变换投影到低维空间上。 PCA的指导思想是:找出最能够代表原始数据的投影方法

问题的提出

在研究生物学问题,常常希望把更多的特征(如转录组产生的表达矩阵)纳入数学模型,这些基因的表达往往存在相关性。例如人的编码基因大约有2万多,这些庞大的表达特征矩阵增加了问题的复杂性。 如果能将多个相关的特征综合为少数几个具有代表性的特征

  • 既能够代表原始特征的绝大多数信息
  • 组合后的特征互不相关,降低特征的相关性(冗余);
  • 这些特征就是主成分 寻找这些少数具有代表性特征的过程被称为主成分分析

PCA去掉了什么

在降低特征矩阵复杂性的同时,希望降维后的数据不能失真,只去除噪声冗余的数据

  • 噪音污染干扰了想听到的真正声音。假设样本中某个主要的维度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()

其实我们还留了一个悬念,就是前面的主成分,到底是由哪些基因组成呢?那些基因的重要性在该主成分的比例如何呢?

本文参与 腾讯云自媒体分享计划,分享自微信公众号。
原始发表:2019-08-09,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 生信技能树 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体分享计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 问题的提出
  • PCA实现的数学基础
  • PCA的步骤(公式推导)
  • R实现简单的PCA分析
  • 总结
  • PCA和tSNE
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档