03高通量测序-PCA中的主要概念

2021-01-14 本文已影响0人不到7不改名

为什么要使用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轴上对样本的区分有贡献最大

å�¾ç��

计算PCA

我们还回最初2个基因的案例上来：

第一，我们根据所有的样本的这两个基因的表达情况绘制出它们的平面坐标系
第二，计算Gene 1表达水平的均值，就是下图中x轴上红叉的所在位置
第三，计算Gene 2表达水平的均值，就是下图中y轴上红叉的所在位置
第四，有了Gene 1和Gene 2的这两个数值，我们就得到了整体数值的中心
第五，我们把所有的数据点往左下角移动，将数据的中心与原点(0,0)重合，就像下面的这个样子，从左图移成右图的样子,此时，虽然所有的数据点都移动了，但是每个点之间的相对距离是不变的，整体数据的中心就变成了原点(0,0)。

å�¾ç��

第六，我们找一条通过原点的直线，就像下面的这个样子，具体这条直线做什么的，先不用管
第七，旋转这条直线（红色虚线），使这条直线尽量匹配这些数据
第八，最终，最匹配数据的直线就是下面的这个样子：
å�¾ç��

第1主成分的计算

此时，我们可能有疑问，为什么这条直线是最匹配数据的，它的计算原理是什么，那么接着看。我们先回到最初的直线，为了准确地找出最佳匹配所有数据的直线，PCA会将所有数据点都映射到这条直线上来，此时，可以计算这些数据点到投射到这条直线上的距离，并且使这些距离最小，除了可以计算数据点到直线的距离最小外，还要计算所有数据点投射到这条直线上的点（图中绿叉位置所在点）到原点的距离，使这个距离最大。

å�¾ç��

通过勾股定理，我们可以知道，因为a不变，当数据点到直线的距离最短时(b)，投影点到原点的距离最大(c)

å�¾ç��

我们计算投影点到原点的距离，我们把它我们把它命名为d1，d2......计算剩余的投影点到原点的距离。然后把这些值的平方加起来称为SS(distances)。我们旋转直线，直到SS的值最大,此时数据点到直线的距离最短。这条直线就叫第一主成分（Principal Component 1，简称PC1）。

å�¾ç��

第1主成分的特征值

对于PC1，斜率为0.25。也就是说基因1增加4个单位，基因2增加1个单位。计算出红色箭头的长度为4.12。

å�¾ç��

当你用SVD(singular value decomposition,奇异值分解)进行PCA时，红色箭头的长度=1，我们所要做的就是把这个三角形缩小到红箭头是1个单位时，只需每边除以4.12。

image-20210101132637041.png

三个边长分别变成了1，0.242，0.97，但，Gene1/Gene2仍然等于4。此时我们可以说PC1由0.97的Gene1和0.242的Gene2构成。

image-20210101132838793.png

我们回顾一下计算过程：

原始数据
最佳拟合的直线
单位向量计算

这个单位向量由基因1的0.97和基因2的0.242部分组成，称为PC1的“奇异向量”(Singular Vector)或“特征向量”(Eigenvector)，每个基因的比例则被称为载荷得分（Loading Scores）。原始数据的投影点到原点的距离的平方SS被称为PC1的特征值(Eigenvalue)。PC1的特征值的平方根叫做PC1的奇异值(Singular value)。

image-20210101134323147.png

第2成分的计算

因为是一个二维图，PC2只是一条垂直于PC1的穿过原点的直线，没有任何进一步的优化要做。由于PC2与PC1垂直，所以斜率为-4，也就是说PC2由-1份Gene1和4份Gene2组成。如果我们对所有东西进行缩放，得到一个单位向量，PC2由-0.242个Gene 1和0.97个Gene 2构成，称为PC2的“奇异向量”(Singular Vector)或“特征向量”(Eigenvector)。每个基因的比例则被称为载荷得分（Loading Scores），告诉我们，就基因值如何投射到PC2上而言，Gene2的重要性是Gene1的4倍

image-20210101135414086.png

最后，PC2的特征值是投影点到原点的距离的平方和。

image-20210101135622196.png

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

图片

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

图片

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

图片

而Sample 6位于这个点上：

图片

同理，Sample 2在这里：

图片

计算每个PC变化百分比和碎石图

我们可以将特征值转化为PC1到原点的变异，通过除以样本大小减1：n-1。这个例子，假设PC1的变异为15,PC2的变异为3。这意味着PCs的变化是15 +3= 18。PC1占了PCs变异的15 / 18= 0.83 = 83%。PC2占了3/18= 0.17= 17%的PCs变异。

image-20210101140855916.png

碎石图(scree plot)是用图形表示每个PC所占的变异百分比。

image-20210101141023866.png

稍微复杂的例子

PCA有3个变量(在这种情况下，3个基因)几乎等同于2个变量

image-20210101141148982.png

然后找到经过原点的最佳拟合直线，和之前一样，最佳拟合线是PC1。PC1现在有三种成分：0.62的Gene1、0.15的Gene2和 0.77的Gene3，Gene3 是最主要的组成部分。然后求出PC2，它经过原点并垂直于PC1。PC2现在有三种成分：0.77的Gene1、0.62的Gene2和 0.15的Gene3，Gene1 是最主要的组成部分。然后，我们找到了PC3，这条最合适的直线，它通过原点并垂直于PC1和PC2。如果我们有更多的基因，我们就会通过添加垂线和旋转它们来不断寻找越来越多的主成分，理论上，每个基因(或变量)都有一个。但在实际操作中，PC的数量不是变量的数量或者样本的数量，取其中较小的一个。一旦你找出了所有的主成分，你可以使用特征值(即SS(距离))来确定每个PC的变化比例。在这个例子中，PC1=79%，PC2=15%,PC3=6% ,PC1和PC2占了变异的绝大比例。这就表明了，在二维图中，我们基本上只使用PC1和PC2就能解释三维图中的数据，因为二维图中的PC1和PC2占据了整体的变异的94%，