StatQuest学习笔记14——PCA

前言——主要内容

这篇笔记是StatQuest系列视频教程的第39-43节，其中39节的内容是PCA；第40节的内容是PCA的基本思想；第41小节的内容还是PCA，此处与第39小节的内容重复，不记笔记，略过；第42节的内容是R语言与PCA；第43节的内容是Python与PCA。为了使笔记更有条理，笔记的顺序是按照视频教程40，39，42，43节的顺序进行的。

PCA的基本思想

我们先从最简单的案例入手，先看下面的一堆绿色的圆圈，如下所示：

我们可以把这些绿色圆圈当成一群正常的细胞，当然了，如果你不是学生物的，你也可以把它们理解成一群人，一群卡车，一群城市等，总之是一群事物。虽然这一群细胞看上去是一样的（至少表面上一样），但是，它们本质上也有可能不一样，例如像下面的这个样子：

这群细胞分成了三群，不过从直观上来看，你太不容易能看出这一群细胞划分成了三群。

为了能够严格地划分这些细胞，我们可以对这些细胞进行测序，确认这些细胞中有哪些基因处于活跃状态，我们就能知道这些细胞的功能，这样就知道了这些细胞的种类，当然了，如果这是一群人，我们可以测量他们的身高，血压等，如下所示：

现在我们测了这些细胞的基因，如下所示：

假设，我们此时没有那么多细胞，只有2个细胞，那么我们就可以把这两个细胞的所有基因画在一个坐标上，如下所示：

就像上图显示的这样，基因1在细胞1中的高表达，而在细胞2中低表达，而基因9正好相反，在细胞1中低表达，而在细胞2中高表达，如下所示：

从我们所列出的数据来看，细胞1和细胞2的基因表达水平很多呈相反的水平，那么我们很容易就知道，这两个细胞必定不是同一种类的细胞，如下所示：

此时，我们把问题升级一下，我们有3个细胞，现在利用二维坐标系分别绘制出这三个细胞两两之间的基因表达水平，如下所示：

但这种两两之间的比较太麻烦了，我们可以直接把它们放到一个三维坐标系中进行比较，如下所示：

在这个三维坐标系中，它的三个坐标轴分别表示细胞1，细胞2，细胞3，我们通过旋转坐标轴就可以发现每两个细胞之间的关系。

现在把问题再进一步，如果我们有4个细胞，比较直观的手段就是分别画出这4个细胞两两的平面坐标系，如下所示：

或者是画出一个四维的坐标系，像下面的这个样子：

但前面的方法就不太适用于这种情况，例如四维坐标就不太直观，此时我们可以采用主成分分析（PCA）的方法来对细胞分类。

PCA可以将不同细胞之间的关系转换到一个平面坐标系中，如下所示：

此时我们可以发现，通过PCA这种方法，可以把一群细胞中高度相关的细胞聚集起来，像下面的这个样子：

为了能够更直观地说明这种分析方法，我们可以把不同的细胞换成不同的颜色，如下所示：

此时，我们再回到最初的问题上来，通过PCA方法可以把一群细胞中高度相关的细胞给区分开来，如下所示：

PCA的结果解释

此时我们来看一下，PCA的结果具体是什么含义，在PCA图中，它的坐标轴是按重要性进行排序的。其中PC1是第一主成分轴，它的重要性要强于PC2，如下所示：

假如一个PCA长得是下面的这个样子，如下所示：

此时，如果左上角的绿点与左下角的红点的距离是d1，右下角的黄点与左下角红点的距离是d2，此时d1=d2的话，那么我们就可以认为，黄色点与红色点的差异要大于绿色点与红色点的差异（因为PC1的坐标是第一主成分），如下所示：

在我们了解了PCA之后，我们应该知道，PCA仅仅是一种划分不同数据类型的方法之一，此外，还有其他基于降维的统计学方法（dimension reduction），例如热图，t-SNE图，多维缩放（multi-dimensional scaling，MDS）等，如下所示：

PCA的计算原理

简单案例

在这一部分，我们使用SVD（singular value decomposition，中文译名“奇异值分解”）的方法来计算PCA。我们先看一个简单的案例，在这个案例中，我们检测了6只不同小鼠的2个基因，如下所示：

其实我们可以把它再抽象化一下，把小鼠看成样本，基因看成2个变量，就像下面的这个样子：

我们还以小鼠的基因为例说明，如果我们只检测1个基因的话（Gene 1），那么我们根据基因1表达的情况，把小鼠的绘制到数轴上，如下所示：

其中小鼠1，小鼠2，小鼠3的Gene 1表达水平比较高，如下所示：

而小鼠4，小鼠5，小鼠6的Gene 1水平则较低，如下所示：

虽然这个图形比较简单，但是，我们从中还是能得到一些信息的，例如小鼠1，小鼠2和小鼠3比较接近，小鼠4，小鼠5和小鼠6比较接近，如下所示：

如果我们检测了2个基因，那么我们可以绘制一个二维坐标系，横轴是Gene 1，纵轴是Gene 2，那么小鼠1，小鼠2和小鼠3会聚在一起，小鼠4，小鼠5和小鼠6会聚在一起，如下所示：

再进一步，如果我们检测了3个基因，那么我们可以绘制三维的坐标系，如下所示：

在上图的这个3维坐系中，圆点越大，表示离你越近。

再进一步，如果我们检测了4个基因，此时我们很难绘制出四维的坐标系，那么我们就需要进行PCA分析了，PCA可以把超过4个的基因降维成二维的坐标系，如下所示：

在这个PCA的二维坐标系中，我们可以发现，小鼠4、小鼠5和小鼠6是一类，小鼠1，小鼠2和小鼠3是一类，它们的各自的基因表达模式也类似，如下所示：

PCA在对数据进行聚类（clustering）时有很大的价值，例如，经过PCA分析，在它的二维坐标轴上我们可以发现，Gene 3在x轴上对样本的区分有贡献最大（这里不太理解为什么直接就看出了Gene 3，我猜测可能是基因3在这两群小鼠中的分布比较明显，一组高，一组低），如下所示：

计算原理

为了理解PCA是如何进行计算的，我们还回最初2个基因的案例上来，如下所示：

接着是计算过程。

第一，我们根据所有的样本的这两个基因的表达情况绘制出它们的平面坐标系，如下所示：

第二，计算Gene 1表达水平的均值，就是下图中x轴上红叉的所在位置，如下所示：

第三，计算Gene 2表达水平的均值，就是下图中y轴上红叉的所在位置，如下所示：

第四，有了Gene 1和Gene 2的这两个数值，我们就得到了整体数值的中心，如下所示：

此时，我们有了这个点（蓝叉所在位置）后，我们此时就不需要再看原始数据了，如下所示：

第五，我们把所有的数据点各左下角移动，将数据的中心与原点(0,0)重合，就像下面的这个样子，从左图移成右图的样子：

此时，虽然所有的数据点都移动了，但是每个点之间的相对距离是不变的，整体数据的中心就变成了原点(0,0)，如下所示：

第六，我们找一条通过原点的直线，就像下面的这个样子，具体这条直线做什么的，先不用管：

第七，旋转这条直线（红色虚线），使这条直线尽量匹配这些数据，如下所示：

第八，最终，最匹配数据的直线就是下面的这个样子：

第1主成分的计算原理

此时，我们可能有疑问，为什么这条直线是最匹配数据的，它的计算原理是什么，那么接着看。

我们先回到最初的直线，如下所示：

为了准确地找出最佳匹配所有数据的直线，PCA会将所有数据点都役射到这条直线上来，如下所示：

此时，可以计算这些数据点到投射到这条直线上的距离，并且使这些距离最小，如下所示：

除了可以计算数据点到直线的距离最小外，还要计算所有数据点投射到这条直线上的点（图中绿叉位置所在点）到原点的距离，使这个距离最大，如下所示：

如果你不理解上述的这种描述，那么还可以这么理解，当这条直线（红色虚线）尽量匹配这些数据点时，这些数据点到这条直线的距离是在缩小的，如下所示：

同样的，这些数据点投射到直线上的点到原点的距离是加大的，如下所示：

为了更好地以数学思维来理解这个过程，此时我们先看一个点的情况，如下所示：

这个点是固定的，它到原点的距离也是固定的，如下所示：

换句话讲，当直线（红色虚线）在旋转的过程中，这个点到原点的距离是不会变的，如下所示：

当直线旋转到某一位置时，我们将这个点投射到直线上，如下所示：

此时，我们会得到一个直角，如下所示：

有了直角，我们就能得到如下的直角三角形，如下所示：

此时，我们把这个直角三角形的三条边命名为a，b,c，那么根据勾股定理，我们就会得到这个公式，即${a}^2={b}^2+{c}^2$，如下所示：

三角形的斜边a是不会变的，如下所示：

如果直角边b变大的话，如下所示：

那么直角边c就会变小，如下所示：

反之，也成立。而PCA计算的则是降低直角边b，也就是降低数据点投射到直线上的距离，如下所示：

同时，PCA还要保证投到到直线上的点到原点的距离最大，也就是说要求c边最大，如下所示：

通过上面利用直角三角形的解释，就比较容易理解PCA计算最佳匹配直线的思路了。此时，再回到起点，将所有的数据点投射到直线上，如下所示：

然后计算某个点到原点为的距离，我们把它命名为${d}_1$，如下所示：

同理，计算第2点投射到直线上的点到原点的距离，命名为${d}_2$，如下所示：

计算剩余的点投射到直线上的点到原点的距离，如下所示：

然后把这些值的平方加起来，如下所示：

此时，我们旋转直线，直到最终这个SS的距离最大（这个距离最大，也就是说数据点到直线的距离最小）为止，如下所示：

最终，得到了我们想要的直线，如下所示：

此时，这条直线就叫第一主成分（Principal Component 1，简称PC1），如下所示：

有了这个直线，它的斜率就很容易计算出来了，经计算，PC1的斜率为0.25，如下所示：

这个斜率0.25如何理解呢？它表示的意思就是，每当我们在Gene 1（x轴坐标）前进4个单位时，在这个直线上的数据点在Gene 2上就增加1个单位，如下所示：

如果按照生物学意义来理解就是，这几个样本在Gene 1坐标轴上更加分散，如下所示：

而在Gene 2上的分散程度则较小，如下所示：

我们可以根据鸡尾酒配方（cocktail recipe）来思考这个PC1（注：鸡尾酒是一种混合酒品，掺杂有不同的酒），例如为了生成PC1，我们可以加入4份的Gene 1，再加入1份的Gene 2，另外，我们可以发现，Gene 1与Gene 2的比值是4，它就说明在描述数据的分散程度方面，Gene 1这个基因更加重要，如下所示：

但在数学上，把这种鸡尾酒配方（cocktail recipe）称为是Gene 1和Gene 2的线性组合（linear combination），在统计学上，描述PC1也可以这样，即“PC1是几个变量的一个线性组合”，如下所示：

第1主成分的特征值

紧接上面的案例，继续进行，还要计算一些东西。对于PC1来说，它的斜率是0.25，那么我们在Gene 1增加4个单位，在Gene 2上就会增加1个单位，如下所示：

那么根据勾股定理，我们就可以计算出红色箭头的长度，长度为4.12，如下所示：

当我们使用SVD（全称是Singular-value decomposition，中文是奇异值分解）进行PCA计算时，我们需要把这个红色箭头定义为1个单位，如下所示：

为了将这个红色箭头定义为1个单位，那么我们就要把这个三角形的各个连长都除以4.12，如下所示：

那么，除以4.12后，各个三角形的连长就成了这个样子：

也就是说三角形的三个边长分别变成了1，0.242，0.97，不过构成PC1的线性组合，即Gene 1和Gene 2的比值还是不变的，仍旧为4，那么此时我们就可以说取0.97份的Gene 1和0.242份的Gene 2，就构成了PC1，如下所示：

此时，我们再回想一下我们计算PC1的整个过程：

原始数据；
最佳匹配直线；
单位向量计算，如下所示：

图中的红色箭头是1个长度单位，它由0.97个Gene 1和0.242个Gene 2构成，它称为PC1的奇异向量（Singular vector），或者是特征向量（Eigenvector），而每个基因的比例则被称为载荷得分（Loading Scores），如下所示：

此外，在PCA中，将SS（也就是各个点投射到最佳匹配的直线上的点到原点的距离的平方和）称为PC1的特征值（Eigenvalue），如下所示：

PC1特征向量的值（Eigenvalue for PC1）的平方根被称为PC1的奇异值（Singular Value for PC1），如下所示：

第2主成分的计算原理

前面是第1主成分（PC1）的计算，此时进行第2主成分的计算（PC2）。由于我们使用的二维坐标，PC2仅仅是一条通过了原点，并且与PC1垂直的直线，PC2并没有像PC1那样复杂的优化，如下所示：

由于PC2与PC1垂直，那么很容易地就能计算出PC2的斜率为-4，这也就是说，PC2是由-1个Gene 1和4个Gene 2构成的，如下所示：

当我们在PC2所在直线上定义一个单位长度的向量时，那么PC2就是由-0.242个Gene 1和0.97个Gene 2构成了，此时这个单位向量（蓝色箭头）就被称为PC2的奇异向量（Singular Vector）或特征向量（Eigenvector），如下所示：

而-0.242和0.97就是PC2的载荷分数，它表示的就是，基因2相当于4位的基因1，如下所示：

PC2的特征值则是数据点投到到PC2上的点到原点的距离，如下所示：

此时，PC1和PC2的计算结束，绘制最终的PCA图，如下所示：

然后旋转这个坐标，让PC1水平，PC2垂直，如下所示：

在这个新的坐标系中，图中黑色的叉就表示原始的样本6（Sample 6），如下所示：

而Sample 6位于这个点上：

同理，Sample 2在这里：

各主成分的变异

此时，再讲一些术语，先回忆一下前面的PC1的特征值和PC2的特征值，如下所示：

我们把原始的数据点投射到主成分上，如下所示：

计算投射到PC1上的点到原点的距离，如下所示：

所有的这些数值的平方和就是PC1的特征值，如下所示：

我们也可以将特征值转为PC1到原点的变异，方法就是SS除以样本数减1，如下所示：

以这个案例为例说明，假如PC1的变异（Variation）为15，PC2的变异为3，这就表明主成分的总变异为18，那么也就是说PC1在总变异中所占的比值为83%（15/18=0.83），如下所示：

而PC2占的总变异为17%（3/18=0.17），如下所示：

此时，我们把这两个主成分所占的比例画出来，这种图与条形图很像，不过在PCA中，这叫碎石图（scree plot），每个矩形表示每个PC所占总变异的比值，如下所示：

一个复杂案例

我们看一个比较复杂的案例，在这个案例中，6只小鼠，检测了3个基因，如下所示：

按照前面的思路，找到数据的中心，并将其移到原点，如下所示：

找到最佳的匹配数据的直线，也就是PC1，如下所示：

不过这个PC1有3个成分，分别是0.62个Gene 1，0.15个Gene 2和0.77个Gene 3，如下所示：

其中，Gene3 是最主要的组成部分，如下所示：

然后寻找PC2，PC2通过原点，并与PC1垂直，如下所示：

最终，PC2的成分为0.77个Gene 1，0.62个Gene 2和0.15个Gene 3，如下所示：

在PC2中，Gene 1是最主要的组成部分，如下所示：

最后，我们寻找PC3，PC3通过原点，同时与PC1和PC2垂直，如下所示：

如果我们有更多的基因，那么我们就会不断地寻找更多的主成分，也就是说不断地往这个坐标中加入与前一个主成分垂直的直线，从理论上讲，一个基因就有一个主成分（例如这个案例中，3个基因，最终找到了3个主成分），因此，主成分的数目要么跟变量相等，要么跟样本的数目相等（这个数目还是比较小的），如下所示：

不过，真正在进行PCA的时候，最终的主成分数目是有限的，这个马上就要讲到。一旦计算出来了所有的主成分，那么我们就可以使用特征值来计算每个主成分的变异占总变异的比例，如下所示：

在这个案例中，PC1的变异占总变异的79%，PC2占15%，PC3占6%，把这些主成分绘制成碎石图，就是下图左图的样子，如下所示：

从碎石图中我们可以看到，PC1和PC2占了变异的绝大比例，如下所示：

这就表明了，在二维图中，我们基本上只使用PC1和PC2就能解释三维图中的数据，因为二维图中的PC1和PC2占据了整体的变异的94%，如下所示：

我们将三维图转换为二维图，画出PC1和PC2，以及原始数据点，如下所示：

将数据投射到这个新坐标的PC1上，如下所示：

再投射到PC2上，如下所示：

然后，我们把这个坐标旋转一下，方便查看，如下所示：

下面的这个点就是Sample 4对应的点，如下所示：

我们再来回顾一下这个案例，首先我们把所有的数据点绘制到三维坐标中，如下所示：

再计算主成分，如下所示：

再计算主成分的特征值与变异比例，如下所示：

最后，我们使用PC1和PC2这个二维坐标来描述原始数据，如下所示：

如果我们检测了6个样本的4个基因，那么我们并不需要绘制出4维坐标图，如下所示：

我们可以计算出它们的主成分，画出碎石图，其中我们可以发现，PC1和PC2占据了变异的90%，那么我们仅使用PC1和PC2构成的二维PCA图就可以表示原始数据，如下所示：

将原始数据分别投射到PC1和PC2上，如下所示：

下面的这两个点就对应了Sample 2，如下所示：

也就是这个点：

那么，所有的点就是这个样子：

如果我们遇到下面的这种碎石图，在这种碎石图中，PC3和PC4也占据了相当比例的总变异，那么我们仅使用2个主成分（也就是PC1和PC2）就无法精确地代表原始数据，如下所示：

不过，即使是有瑕玼的PCA二维图，也可以用于数据的分类，如下所示：

在这个图中，黑色的点仍然可以被视为比较接近的类似，如下所示：

R语言与PCA

计算PCA的工具很多，这里我们介绍一下用R来计算PCA，在R中，计算PCA的函数是prcomp()，绘图包是ggplot2，使用R，我们可以知道每个主成分的变异占总变异的比值，以及计算载荷得分（loading scores），发现哪个变量对图形的影响最大，如下所示：

生成数据集

我们先生成一个数据集，代码如下：

data.matrix <- matrix(nrow=100, ncol=10)
colnames(data.matrix) <- c(
  paste("wt",1:5,sep=""),
  paste("ko",1:5,sep=""))
rownames(data.matrix) <- paste("gene", 1:100, sep="")
for (i in 1:100){
  wt.values <- rpois(5, lambda = sample(x=10:1000, size=1))
  ko.values <- rpois(5, lambda = sample(x=10:1000, size=1))
  data.matrix[i,] <- c(wt.values, ko.values)
}
head(data.matrix)

查看一下生成的数据，如下所示：

> head(data.matrix)
      wt1 wt2 wt3  wt4 wt5 ko1 ko2 ko3 ko4 ko5
gene1 423 423 419  423 372 582 589 563 534 566
gene2 905 971 967 1031 910  95  96  96  89  96
gene3 139 177 147  176 152 408 422 419 426 456
gene4 702 712 717  691 678 534 546 594 562 515
gene5 582 572 616  569 631 315 321 314 345 330
gene6 765 853 781  792 768 316 339 317 317 294

计算PCA

使用prcomp函数来对数据集进行PCA运算，如下所示：

1	pca <- prcomp(t(data.matrix),scale=TRUE)

代码解释：由于prcomp函数在默认情况下，会把行识别为样本，列识别为基因，而我们生成的数据集data.matrix正好相反，因此用到了t函数，它表示把原来的数据集转置一下，将其行与列交换。

prcomp计算的结果含有主成分分析的各种信息，通过str()函数我们可以看一下有哪些东西，如下所示：

> str(pca)
List of 5
 $ sdev    : num [1:10] 9.34 1.85 1.6 1.45 1.24 ...
 $ rotation: num [1:100, 1:10] 0.104 -0.107 0.106 -0.103 -0.106 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:100] "gene1" "gene2" "gene3" "gene4" ...
  .. ..$ : chr [1:10] "PC1" "PC2" "PC3" "PC4" ...
 $ center  : Named num [1:100] 489 526 292 625 460 ...
  ..- attr(*, "names")= chr [1:100] "gene1" "gene2" "gene3" "gene4" ...
 $ scale   : Named num [1:100] 84.2 455.8 142.2 82.1 143.2 ...
  ..- attr(*, "names")= chr [1:100] "gene1" "gene2" "gene3" "gene4" ...
 $ x       : num [1:10, 1:10] -8.92 -9.15 -8.71 -8.75 -8.77 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:10] "wt1" "wt2" "wt3" "wt4" ...
  .. ..$ : chr [1:10] "PC1" "PC2" "PC3" "PC4" ...
 - attr(*, "class")= chr "prcomp"

从这个结果中我们可以发现，有sdev，rotation，center，scale，x。其中x含有主成分，我们查看一下，如下所示：

我们可以发现，x中有10个主成分（这是因为有10个样本），第1主成分占据了原始数据总变异的大部分，第2主成分占据了总变异的次大部分等。为了绘制二维图，我们通常使用PC1和PC2（个别情况下也会使用PC2和PC3），如下所示：

我们此时，使用x的前2列（即PC1和PC2）绘制一个二维图，如下所示：

plot(pca)

其中，pca$x[,1]是PC1，pca$x[,2]是PC2，如下所示：

10个样本，分别位于图形的左右两侧，此时我们再看一下PC1占的总变异的比例，如下所示：

1 2	pca.var <- pca$sdev^2 pca.var.per <- round(pca.var/sum(pca.var)*100,1)

结果如下所示：

> pca.var
 [1] 8.727308e+01 3.432450e+00 2.546956e+00 2.094293e+00 1.534746e+00 1.230408e+00
 [7] 8.259253e-01 6.435822e-01 4.185583e-01 1.134741e-29
> pca.var.per
 [1] 87.3  3.4  2.5  2.1  1.5  1.2  0.8  0.6  0.4  0.0

绘制出各个主成分对应的碎石图，如下所示：

barplot(pca.var.per, 
        main="Scree Plot", 
        xlab="Principal Component", 
        ylab="Percent Variation")

上图中最长的那个就是PC1所占的总变异的比例，也就是说PC1就能在很大程度上解释两组数据的区分，我们在PC1和PC2的二维图上也能看出来，这两组数据之间主要是由PC1分开的，如下所示：

ggplot2

ggplot2也能绘制类似的图形，先把数据转换为ggplot2能识别的格式，如下所示：

pca.data <- data.frame(Sample=rownames(pca$x),
                       X=pca$x[,1],
                       Y=pca$x[,2])
# one column with the sample ids
# Tow columns for the X and Y coordinates for each sample
pca.data

生成的数据如下所示：

> pca.data
    Sample         X          Y
wt1    wt1 -8.920074  0.9049750
wt2    wt2 -9.154298 -0.8467998
wt3    wt3 -8.712942 -1.0756417
wt4    wt4 -8.746978  2.6207842
wt5    wt5 -8.771551 -1.5756187
ko1    ko1  8.854899  2.7689355
ko2    ko2  8.976103 -1.7855946
ko3    ko3  9.077341  0.2534441
ko4    ko4  8.767547  1.2557291
ko5    ko5  8.629953 -2.5202131

接着使用ggplot绘图，如下所示：

ggplot(data=pca.data, aes(x=X, y=Y, label=Sample))+
  geom_text() +
  xlab(paste("PC1 - ", pca.var.per[1], "%", sep=""))+
  ylab(paste("PC2 - ", pca.var.per[2], "%", sep=""))+
  theme_bw()+
  ggtitle("My PCA Graph")

从这个图形中，我们可以发现，X轴上标记了PC1所占总变异的比例为87.3%。Y轴上标记了PC2所占总变异的比例为3.4%。

loading scores计算

prcomp()函数计算的结果中含有loading score，即rotation这一部分的内容，如下所示：

1	loading_score <- pca$rotation[,1]

由于PC1占所了总变异的大部分，因此我们只需要看PC1中的loading score，如下所示：

由于loading_score有负数，因此我们需要用abs()函数取它们的绝对值，如下所示：

取完绝对值，对所有的loading score进行排序，从高到低，代码如下：

1 2	gene_scores <- abs(loading_score) gene_score_ranked <- sort(gene_scores, decreasing = TRUE)

此时，还要获排名前10的loading score基因的名称，如下所示：

1 2	top_10_genes <- names(gene_score_ranked[1:10]) top_10_genes

结果如下：

> top_10_genes <- names(gene_score_ranked[1:10])
> top_10_genes
 [1] "gene10" "gene98" "gene34" "gene94" "gene95" "gene99" "gene76" "gene47" "gene49"
[10] "gene72"

最后，我们看一下这些基因的loading score的数值，如下所示：

1	pca$rotation[top_10_genes,1]

结果如下所示：

> pca$rotation[top_10_genes,1]
    gene10     gene98     gene34     gene94     gene95     gene99     gene76     gene47 
 0.1069736  0.1069637  0.1069535  0.1069478 -0.1069281  0.1069267  0.1069243 -0.1069209 
    gene49     gene72 
 0.1069161  0.1069037

Python与PCA

利用Python也可以进行PCA分析，在这一部分中，我们会介绍如何这些内容：

利用Python生成一些我们要分析的数据；
利用sklearn模块中的PCA()来进行分析，其中sklearn是scikit-learn的简称，它是一个机器学习的模块，支持包括分类、回归、降维和聚类四大机器学习算法。还包含了特征提取、数据处理和模型评估三大模块；
计算每个主成分占总变异的比例；
利用matplotlib来绘制PCA图；
计算变量的loading scores。

如下所示：

第一步，导入各种函数与包

在进行PCA分析之前，Python需要导入各种包与函数，如果没有这些包，可以查一下如何下载，我使用的是Anaconda发行版，里面已经内置了各种有关数据分析的包，不用再下载各种包，直接导入即可，如下所示：

import pandas as pd 
# import pandas package,
# pandas is short for "panel data"
import numpy as np
# numpy package will generate random numbers 
import random as rd
# random package will generate an example dataset
from sklearn.decomposition import PCA
from sklearn import preprocessing
# usr a function PCA() from sklearn.decomposition 
# PCA() use to perform PCA
# preprocessing uses to scaling the data before perform PCA
import matplotlib.pyplot as plt
# draw PCA graphd

第二步：生成模拟的数据

这一步与R中生成模拟数据集的目的一样，就是生成一个数据集，用于PCA分析，如下所示：

genes = ['gene' + str(i) for i in range(1, 101)]
wt = ['wt' + str(i) for i in range(1,6)]
ko = ['ko' + str(i) for i in range(1,6)]
data = pd.DataFrame(columns=[*wt, *ko], index=genes)
for gene in data.index:
    data.loc[gene, 'wt1':'wt5'] = np.random.poisson(lam=rd.randrange(10,1000),size=5)
    data.loc[gene, 'ko1':'ko5'] = np.random.poisson(lam=rd.randrange(10,1000),size=5)
print(data.head())
print(data.shape)
# generate simulation dataset

第三，进行PCA分析

在进行PCA分析之前，先对数据进行中心化（我的理解是Z转换，转换后的数据均值是0，标准差是1），如下所示：

1	scaled_data = preprocessing.scale(data.T)

需要注意的是，上述代码中使用了data.T，这是对表达的数据集进行转换，将其行与列交换，跟R中类似。此处我们使用的sklearn包中的preprocessing函数，我们还可以使用另外的函数进行转换，就是StandardScaler().fit_transform(data.T)。第二个函数在机器学习使用更加广泛，如下所示：

此外，我们需要注意的是，sklearn中利用preprocessing.scale计算变异与R中利用scale()和prcomp()计算变异的公式略有不同，区别就在于前者的分子是测量值的数目，后者是这个数目减1，如下所示：

后者的分子是测量值的数目减1，这会导致最终的结果变大，这是一种无偏含估计，如下所示：

不过这两种方法都不影响PCA的分析，如下所示：

它的局限就在于，它们会对最终的绘图造成一定的影响，因为最终绘图中的坐标来源于scaled后的数值乘以loading scores，如下所示：

继续Python计算PCA，此时我们要创建一个PCA对象，如下所示：

1	pca = PCA()

这里要提一下，PCA不仅是一个函数，能够返回计算结果，在sklearn中，它还用于创建一个对象，使用一个数据集进行训练，进行对另外一个数据集进行计算。不过在这里，我们只用PCA计算这个数据集（我们不涉及机器学习的环境设置），其他的步骤非常繁琐，但也由此，我们也能了解一下机器学习的内容，如下所示：

计算PCA的所有结果，例如loading scores和每个主成分的变量占总变异的比例，如下所示：

1	pca.fit(scaled_data)

结果如下所示：

1
2
3

>>> pca.fit(scaled_data)
PCA(copy=True, iterated_power='auto', n_components=None, random_state=None,
  svd_solver='auto', tol=0.0, whiten=False)

接着，基于loading scores和scaled data来生成PCA图形中的坐标系，如下所示：

1	pca_data = pca.transform(scaled_data)

第四，绘图

首先，绘制一个碎石图（Scree plot），看一下应该使用多少个主成分比较合适，如下所示：

1	per_var = np.round(pca.explained_variance_ratio_*100,decimals=1)

这段代码主要是用于计算每个主成分所占总变量的比例，接着添加坐标轴的标签，如下所示：

1	labels = ['PC' + str(x) for x in range(1, len(per_var)+1)]

其次，使用matplotlib来绘制相应的图形，如下所示：

plt.bar(x=range(1,len(per_var)+1), height=per_var, tick_label=labels)
plt.ylabel('Percentage of Explained Variance')
plt.xlabel('Principal Component')
plt.title('Scree Plot')
plt.show()

图形如下所示：

最后，我们绘制PCA图，我们使用pd.DataFrame函数把新的坐标系放到一个矩阵中来，它的行表示样本名称，列表示PC的标签，如下所示：

1	pca_df = pd.DataFrame(pca_data, index=[wt, ko], columns=labels)

继续绘制含有标题名与坐标标签的散点图，如下所示：

plt.scatter(pca_df.PC1,pca_df.PC2)
plt.title('My PCA Graph')
plt.xlabel('PC1 - {0}%'.format(per_var[0]))
plt.ylabel('PC2 - {0}%'.format(per_var[1]))

最后，把样本名也添加上去，如下所示：

for sample in pca_df.index:
    plt.annotate(sample, (pca_df.PC1.loc[sample],pca_df.PC2.loc[sample]))
plt.show()

最终的结果图如下所示：

从这张图中我们可以看到，wt组的数据集中在右侧，ko组都集中在了左侧。

计算loading score

我们还可以计算一下PC1的loading scores，从而找出哪些基因在不同的分类方面发挥了主要作用，如下所示：

首先，我们要创建一个pandas的series对象，这个对象用于计算loading scores，如下所示：

loading_scores = pd.Series(pca.components_[0], index=genes)
# in python, first  number in index is 0, not 1
sorted_loading_scores = loading_scores.abs().sort_values(ascending=False)
# achieve absolute value using abs() function and sort them
top_10_genes = sorted_loading_scores[0:10].index.values
print(loading_scores[top_10_genes])

最终结果如下所示：

>>> print(loading_scores[top_10_genes])
gene31   -0.107133
gene32    0.107131
gene34    0.107070
gene10   -0.107065
gene45    0.107065
gene28   -0.107064
gene30   -0.107062
gene39    0.107053
gene93   -0.107041
gene79    0.107036
dtype: float64s

计算结束。