最近项目需要计算Seurat[1] (v3.1.2)进行PCA降维后选取的top20 PC对数据的解释度，但是发现Seurat自带函数计算的结果并不能直接用来计算。
RunPCA
RunPCA是seurat包自带的进行PCA降维的函数，根据高可变基因的表达矩阵进行降维，返回降维结果并存储在seurat对象中。
使用Seurat对实际数据进行降维并计算top20 PC的贡献度：
PRO <- readRDS('SP.diff_PRO.rds')
PRO <- FindVariableFeatures(PRO,selection.method = "vst", nfeatures = 2000)
PRO<- ScaleData(PRO)
PRO <- RunPCA(PRO, pc.genes = PRO@var.genes, npcs = 20, verbose = FALSE)
查看一下RunPCA后的Seurat对象中存储的PCA解释度信息：
PRO@reductions$pca@stdev
[1] 12.964063 7.052222 5.740453 4.267826 3.926931 3.588792 3.420418
[8] 3.088306 2.939566 2.801638 2.726421 2.572142 2.373437 2.228373
[15] 2.165477 2.107010 2.042534 1.992695 1.898133 1.879603
可以看到，Seurat对象只保留了PCA降维时选取的N(这里是20)维度的标准差，因此计算每个PC贡献度时无法获取分母；在此我利用两种方法，来计算top20 PC的贡献度：
- 方法1，模拟Seurat内置函数RunPCA
PRO <- readRDS('SP.diff_PRO.rds')
PRO <- ScaleData(PRO)
svd_re <- svd(PRO@assays$RNA@scale.data)
stdev_all <- svd_re$d/((length(colnames(PRO))-1)**0.5)
stdev_all[1:20]
[1] 12.964063 7.052222 5.740453 4.267826 3.926931 3.588792 3.420418
[8] 3.088306 2.939566 2.801638 2.726421 2.572142 2.373437 2.228373
[15] 2.165477 2.107010 2.042534 1.992695 1.898133 1.879603
可以看到结果与Seurat内置函数计算结果一致，且此处的stdev_all包含了所有维度的标准差信息，可以用于计算每个PC的贡献度，例如计算我所需要的top20 PC的贡献度：
stdev_all <-stdev_all**2
per_PC <- stdev_all/sum(stdev_all)
per_PC[1:20]
[1] 0.127051474 0.037596674 0.024910919 0.013769266 0.011657462 0.009736301
[7] 0.008844142 0.007210045 0.006532268 0.005933644 0.005619314 0.005001353
[13] 0.004258464 0.003753818 0.003544904 0.003356067 0.003153813 0.003001782
[19] 0.002723646 0.002670730
- 方法2，利用其它PCA降维函数进行计算，在此我利用prcomp函数进行计算，因为prcomp函数的结果会自动计算每个PC的贡献度，计算代码如下：
PRO <- readRDS('SP.diff_PRO.rds')
pro_re <- prcomp(PRO@assays$RNA@scale.data,scale=F,center=F)
summary(pro_re)
PC1 PC2 PC3 PC4 PC5 PC6
Standard deviation 35.35107 19.23036 15.65336 11.63773 10.70816 9.786102
Proportion of Variance 0.12705 0.03760 0.02491 0.01377 0.01166 0.009740
Cumulative Proportion 0.12705 0.16465 0.18956 0.20333 0.21499 0.224720
PC7 PC8 PC9 PC10 PC11 PC12 PC13
Standard deviation 9.32697 8.42135 8.01576 7.63965 7.434545 7.01385 6.47201
Proportion of Variance 0.00884 0.00721 0.00653 0.00593 0.005620 0.00500 0.00426
Cumulative Proportion 0.23357 0.24078 0.24731 0.25324 0.258860 0.26386 0.26812
PC14 PC15 PC16 PC17 PC18 PC19
Standard deviation 6.076442 5.904933 5.745503 5.569686 5.433783 5.175926
Proportion of Variance 0.003750 0.003540 0.003360 0.003150 0.003000 0.002720
Cumulative Proportion 0.271880 0.275420 0.278780 0.281930 0.284930 0.287660
PC20
Standard deviation 5.125399
Proportion of Variance 0.002670
Cumulative Proportion 0.290330
可以看到，两种方法计算的结果都是一致的，也侧面说明了方法的正确性。
top20 贡献率不高的原因
根据以上结果，得到top20 PC的贡献率为29%，这个贡献率并不高（我以往的知识中，PC1 + PC2 的贡献率>=85%），以至于我开始怀疑结果的正确性，但是最终发现结果是正确的，在此我并没有查找到相关资料来说明这一点，我只是自己推测了如下理由：
- 单细胞数据的稀疏性
- 高可变基因间的重要性差异不是很大
PCA
主成分分析方法 (Principal Component Analysis, PCA)，是一种线性的数据降维算法。PCA就是将原有的n维特征映射到新的k维上，而这k维是全新的正交特征，也被称为主成分，是在原有n维特征的基础上线性组合重新构造出来的k维特征。
PCA的原理就是从原始的n维空间中顺序地找一组相互正交的坐标轴，其中，第一个新坐标轴的方向是原始数据中方差最大的方向，第二个新坐标轴选取是与第一个坐标轴正交的平面中使得方差最大的，第三个轴是与第1,2个轴正交的平面中方差最大的，依次类推，可以得到k个这样的坐标轴。
如何在数学上实现PCA呢，实际上，通过计算数据矩阵的协方差矩阵，然后得到协方差矩阵的特征值和特征向量，选择特征值最大 (即方差最大)的k个特征值所对应的特征向量组成的矩阵即可，这样就可以将数据矩阵转换到新的空间当中，实现数据特征的降维，PCA实现的方法有两种：
- 基于特征值分解协方差矩阵实现PCA算法
- 基于SVD分解协方差矩阵实现PCA算法
