Principal Component Analysis (PCA)
1 PCA原理
1.1 背景
假设我们测序了9个基因在2个细胞中的表达量并绘制成图。细胞1和细胞2呈负相关,基因在两个细胞中的表达量呈现相反关系(即在细胞1中高表达的基因在细胞2中低表达,反之亦然)。
假如我们测序了基因在3个细胞中的表达量。我们需要画3张二维图来理解细胞之间的关系。
或者画一张三维图。
当我们测序了基因在更多细胞中的表达量时,就很难再用更多的2维图来理解细胞之间的关系了。
画一张n维图也不能解决这个问题。
相反,我们会用一张主成分分析图来帮助理解细胞之间的关系。主成分分析图根据细胞之间的相关性绘图。
在图上,相关性高的点(细胞)会聚在一起。
1.2 用single value decomposition(SVD)法理解主成分分析
以2个基因在6个老鼠样本中的表达量为例,将数据绘制在二维图上。
第一步,分别计算基因1和基因2的均值(红叉),然后找到中心点,即所有数据的均值点(蓝叉)。
移动所有点直到中心点与原点重合(这样的移动并不改变数据的相对大小)。
第二步,用一条经过原点的直线来拟合所有的点。
所有点到原点的距离是固定的,这些点到拟合线的距离(拟合距离)越大,则其投影点到原点的距离(投影距离)越小。
到拟合线的距离越小,则其投影点到原点的距离越大。
这个关系可以用勾股定理来解释:a2=b2+c2,a不变;b越大,c越小;b越小,c越大。
分别将这6个点的投影距离标记为d1,d2,...,d6,然后计算这些投影距离的平方和,即sum ofsquared distance(ssd)=d12+d22+d32+d42+d52+d62。当SS获得最大值时(投影距离最大,拟合距离最小),这条线即为最佳拟合线。也即,最佳拟合直线符合最小二乘法原则。
这条直线称为第一主成分(PC1),它的斜率是0.25。
这意味着每当PC1沿基因1向右前进4个单位时,PC1沿基因2向上就增加1个单位。
同时也意味着,数据沿基因1分散地更开,沿基因2分散的程度较小;在描述数据沿PC1变化方面,基因1比基因2更重要;基因1对PC1的作用是基因2的4倍。
第三步,根据勾股定理计算出斜边的长度。
三条边都除以斜边长度,使得斜边长度为1,以达到标准化。
斜边是一个单位长的向量,称为PC1的eigenvector。它由0.97个基因1和0.242个基因2组成。这两个比例(0.97和0.242)称为PC1的loading scores。
投影距离的平方和(ssd)称为PC1的eigenvalue。ssd的平方根称为PC1的singular value。
第四步,第二主成分是经过原点,且与PC1垂直的直线。
每当PC2沿基因1向左前进一个1,则沿基因2向上前进4。
计算PC2的eigenvector和loading scores。
基因1和基因2的比例(-0.242:0.97)表明,在描述数据沿PC2变化时,基因2比基因1重要4倍。
计算PC2的ssd,即eigenvalue。
第五步,转动PC1和PC2,使之变成水平和垂直,得到数据沿PC1和PC2变化的平面图。
用每个主成分的ssd除以n-1(samplesize-1,即自由度),得到每个主成分所能解释的数据的变异(variation)。
然后计算每个主成分所能解释的变异占总的变异的比例,并绘制出条形图。
当有三个基因时,计算方法类似,要求三个主成分两两垂直。
当计算多个基因时亦如是,要求所有主成分两两垂直。
第六步,当前两个主成分所能解释的变异超过90%时,可以用PC1v.s. PC2的二维图来表示高维数据。
当PC1和PC2所能解释的变异占数据总变异的比例较低时,用PC图描述的数据准确性会下降。但是仍能大致看出哪些数据相关性比较强(变异比较相似),聚在了一起。
1.3 注意事项
1.3.1 PC1比PC2解释的变异多
从上面的图中可以看出,PC1是解释变异最多的主成分,也是最重要的主成分,PC2次之。换句话说,在PC1 vs PC2图中,PC1单位长度代表的变异比PC2的大。所以在下图中,尽管红点到黄点和到绿点的距离相同,但是红点和黄点之间的差异要比红点和绿点之间的差异大。
1.3.2 确保变量的数量级一致
如下图,math比reading多一个数量级,math的取值是0-100,reading的取值是0-10。
计算PC1后发现math的loading score是reading的10倍,即math对PC1的重要性是reading的10倍。
而这一切仅仅是因为math比reading多一个数量级,即math的取值是reading的10倍。
如果将math的值都除以10,重新绘图。
重新计算主成分后发现,math和reading有相同的loading score,在解释数据沿PC1变化时有相同的重要性。在实际操作中也可以让math和reading分别除以自己的标准差,以达到标准化的效果。
1.3.3 计算主成分之前确保数据的中心与原点重合
如果使用的R包或者其它算法,没有自动中心化的功能,就要自己事先手动中心化。否则,计算出的PC1会穿过原点,但不一定会穿过中心点,不符合预期。
1.3.4 主成分的个数如何确定
case1:还是以之前的数据为例。
经过分析只能找到两个主成分。如果有第三个主成分,它必须同时与PC1和PC2垂直。在二维平面里,这是不可能的。
case2:假设math和reading 100%相关(即呈现函数关系)。
所有点都在PC1上,它们在PC2的投影距离是0。
也就是说PC1解释了100%的数据变异,PC1占数据总变异的比例为100%。在这种情况下只有一个主成分,即PC1。
case3:假设现在只有2个样本/学生/列,2个变量/行。2个点的中心与2个点共线。
由于三点共线,所以也只有一个主成分。
case4:假设现在有3个变量/行,2个样本/列,并用三维图表示两个样本点。
由于只有2个点,且2个点的中心(原点)与2个点共线,所以只有1个主成分。
case5:假设现在有3个变量/行,3个样本/列。
由于3点确定一个平面,而一个平面只能有两个主成分。所以这种情况下,数据最多有两个主成分。
总结:PCA的数量不能超过行和列中的最小值。比如有5 sample,10个基因,PCA最多有5个。
1.4 以基因/行/variable对第一、二主成分的影响(loadingscores)作为横纵坐标将样本/列/sample绘制在二维图上;以样本/列/sample对第一、二主成分的影响(loadingscores)作为横纵坐标将基因/行/variable绘制在二维图上
第一步,先回顾一下主成分分析的过程。以13个基因在2个细胞中的表达为例。
用一条经过所有点的中心(均值点)的直线拟合数据,使得所有点的投影距离平方和最大。该直线为第一主成分(PC1)。
经过均值点并与PC1垂直的直线为第二主成分。
旋转坐标轴,直到PC1水平,PC2垂直。
数据沿水平方向(PC1)变化较大,沿垂直方向(PC2)次之。
也就是说PC1能够解释的变异最多,占数据总变异的比例最大,PC2次之。
上图是根据基因在2个细胞中的表达差异,将基因绘制在PC图上。下面讨论如何根据基因的表达差异将细胞绘制在PC图上。
第二步,主成分是用来解释变异的,其方向代表了数据变异分布的方向。但是,不同点的变异大小是不同的,越分散的点相较于中心点变异越大,越集中的点变异越小。以PC1为例,a和f分散在PC1的两端,变异比在中间的点大。所以它们在PC1中的比重较大,对PC1的影响也较大。这种影响可以用loading scores量化,影响越大,loading scores绝对值越大,符号代表方向。
同样的,b、c和d分布在PC2的两端,对PC2的影响较大;a和f分布在PC2的中间,对PC2 的影响较小。对PC2的影响同样可以量化。
第三步,利用每个基因在每个细胞中的表达量(read
count)和每个基因对不同主成分的影响(loading scores)来计算每个细胞的每个主成分值。以下图为例,细胞1的PC1值=基因a在细胞1中的read数*基因a对PC1的影响+基因b在细胞1中的read数*基因b对PC1的影响+…。
经过计算得到细胞1的所有主成分值。
细胞1的PC1和PC2值即是该细胞在PC1 vs
PC2图上的横纵坐标。
同理,计算出细胞2的PC1和PC2值并绘制在图上。
如果有更多的细胞,依次计算出每个细胞的PC1和PC2值并绘制上图上。如果两个细胞的基因表达情况相似,它们的PC1和PC2值也会相近,在图上的距离也较小,聚类时会被聚在一起。
基于这样的原理,大量的细胞根据基因的表达相似性,在PC图上被分成了不同的群体。
基因在不同细胞中的表达差异越大,就越有可能分散在主成分的两端,对主成分的影响也就越大(loading
scores越大),在计算细胞的PC值时所起的作用就越大。由于其在不同细胞的表达量不同,计算出的细胞的PC值差异就越大,细胞就更容易分散在PC图上不同的区域。如下图,对PC1影响较大的基因(在PC1上的loading
scores较大)是分开Dermal cells和Neural
cells的主要因素。
类似地,对PC2影响较大的基因是分开Blood
cells和Dermal cells的主要因素。
2 如何用R分析主成分
2.1 用基础R包分析主成分
data.matrix <- matrix(nrow=100, ncol=10)
data.matrix
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8][,9] [,10]
[1,] NA NA NA NA NA NA NA NA NA NA
[2,] NA NA NA NA NA NA NA NA NA NA
[3,] NA NA NA NA NA NA NA NA NA NA
……………
[98,] NA NA NA NA NA NA NA NA NA NA
[99,] NA NA NA NA NA NA NA NA NA NA
[100,] NA NA NA NA NA NA NA NA NA NA
#we make an empty matrix where we will measure 10 samples with 100 genesin each sample
colnames(data.matrix) <-c(paste0("wt",1:5),paste0("ko",1:5))
#add column names
data.matrix
wt1 wt2 wt3 wt4 wt5 ko1 ko2 ko3 ko4 ko5
[1,] NA NA NA NA NA NA NA NA NA NA
[2,] NA NA NA NA NA NA NA NA NA NA
[3,] NA NA NA NA NA NA NA NA NA NA
……………
[98,] NA NA NA NA NA NA NA NA NA NA
[99,] NA NA NA NA NA NA NA NA NA NA
[100,] NA NA NA NA NA NA NA NA NA NA
#the first 5 samples are "wt" (wildtype) while the last 5 are"ko" (knock-out)
rownames(data.matrix) <- c(paste0("gene",1:100))
#add row names
data.matrix
wt1 wt2 wt3 wt4 wt5 ko1 ko2 ko3 ko4 ko5
gene1 NA NA NA NA NA NA NA NA NA NA
gene2 NA NA NA NA NA NA NA NA NA NA
gene3 NA NA NA NA NA NA NA NA NA NA
…………….
gene98 NA NA NA NA NA NA NA NA NA NA
gene99 NA NA NA NA NA NA NA NA NA NA
gene100 NA NA NA NA NA NA NA NA NA NA
#each row represent each gene with name“gene1”…”gene100”
set.seed(1)
#set
up pseudo random number so that the sampling results can be repeated
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)
}
#here we give the genes mimic read count
data.matrix[1:3,]
wt1 wt2 wt3 wt4 wt5 ko1 ko2 ko3 ko4 ko5
gene1 850 820 857 800859 285 326 298 272 275
gene2 519 493 493 514512 330 349 336 298 344
gene3 128 113 130 127
128 664 645 706 669 709
#look at the mimic read counts of the first three lines
data.cpm <- log2(edgeR::cpm(data.matrix)+0.001)
#normalize the read count matrix to eliminate the difference of librarysize
#scale by log2(); adding 0.001 is to avoid negative infiniteresulting from log2(0)
data.cpm[1:3,]
wt1 wt2 wt3 wt4 wt5 ko1 ko2 ko3
gene1 14.0831014.02084 14.08855 13.99467 14.09358 12.48940 12.68098 12.55031
gene2 13.3714213.28686 13.29090 13.35648 13.34713 12.70087 12.77932 12.72343
gene3 11.3522411.16208 11.36823 11.33996 11.34754 13.70948 13.66530 13.79453
ko4 ko5
gene1 12.4261112.43326
gene2 12.5577912.75618
gene3 13.72435
13.79945
#look the first 3 lines to check the cpm values
pca <- prcomp(t(data.cpm),scale=T)
#do PCA on the cpm matrix
#by default, prcomp() expects each row is a sample and eachcolumn is a gene, so here t() is used to transpose the matrix
str(pca)
List of 5
$ sdev : num [1:10] 9.39 1.73 1.59 1.38 1.3 ...
$ rotation: num [1:100, 1:10] -0.1062 -0.105 0.1063 0.0492 0.0236 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:100] "gene1" "gene2" "gene3" "gene4" ...
.. ..$ : chr [1:10] "PC1" "PC2" "PC3" "PC4" ...
$ center : Named num [1:100] 564 419 402 795 572 ...
..- attr(*, "names")= chr [1:100] "gene1" "gene2" "gene3" "gene4" ...
$ scale : Named num [1:100] 288.7 93.5 292.3 22.4 29.2 ...
..- attr(*, "names")= chr [1:100] "gene1" "gene2" "gene3" "gene4" ...
$ x : num [1:10, 1:10] -9 -8.77 -8.96 -8.74 -9.08 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:10] "wt1" "wt2" "wt3" "wt4" ...
.. ..$ : chr [1:10] "PC1" "PC2" "PC3" "PC4" ...
- attr(*, "class")= chr "prcomp"
#pca is a list containing five elements named "sdev","rotation","center","scale" and "x"
pca$x
PC1 PC2 PC3 PC4 PC5 PC6 PC7
wt1 -9.000989 3.72325653 -1.8147624 0.1641587 -0.8558589 0.1191858 -0.3392712
wt2 -8.773236-1.62190849 -0.4101916 -2.2986849 1.2299475 0.2951986 -1.1995011
wt3 -8.955221-1.21169551 0.2153501 -0.3003204-1.1449331 1.6320785 1.3622772
wt4 -8.743347-0.85969323 -0.6059461 0.4116646 1.0018345 -1.1691670 0.6797577
wt5 -9.076914-0.02029153 2.6147249 1.9486497 -0.1832338 -0.8993623 -0.4845155
ko1 8.687179 -1.31832099 -1.2880958 2.3550991 1.1849067 0.9787371 -0.3018722
ko2 8.855259 -1.69954777 -1.4201004 -0.4829262-2.2148896 -1.4662056 -0.0192424
ko3 9.170681 1.58466276 0.4076585-1.0480894 1.8907847 -0.6902354 1.0601891
ko4 8.938537 0.86692510 -0.4963450 0.1681750-0.1323754 0.6884019 -0.4703884
ko5 8.898051 0.55661312 2.7977079 -0.9177261-0.7761825 0.5113684 -0.2874330
PC8 PC9 PC10
wt1 0.539269863 -0.1969555 3.471182e-15
wt2 -0.068458996-0.3043779 3.134645e-15
wt3 -0.291868629-0.2587362 2.910866e-15
wt4 0.424196728 1.3031354 2.140649e-15
wt5 -0.612200296-0.5304322 3.058317e-15
ko1 0.721201453 -0.3555572 3.060052e-15
ko2 0.004468956 -0.3708382 3.020154e-15
ko3 -0.268497039-0.7181874 2.770353e-15
ko4 -1.569391808 0.8949513 3.441691e-15
ko5 1.121279769 0.5369978 2.593412e-15
#x contains the PCA values
pca_value <- as.data.frame(pca$x)
#convert x into a data.frame and save it in pca_value for theconvenience of adding group information later
group <- as.factor(c(rep("wt",5),rep("ko",5)))
#create a factor vector to storegroup information
pca_value$group <- group
#add a new column to store the group information for the convenience ofadding colors on the graph later
pca_value[1:3,]
PC1 PC2 PC3 PC4 PC5 PC6 PC7
wt1 -9.000989 3.723257 -1.8147624 0.1641587 -0.8558589 0.1191858 -0.3392712
wt2 -8.773236-1.621908 -0.4101916 -2.2986849 1.2299475 0.2951986 -1.1995011
wt3 -8.955221-1.211696 0.2153501 -0.3003204-1.1449331 1.6320785 1.3622772
PC8 PC9 PC10 group
wt1 0.5392699 -0.1969555 3.471182e-15 wt
wt2 -0.0684590-0.3043779 3.134645e-15 wt
wt3 -0.2918686-0.2587362 2.910866e-15 wt
#check the first 3 lines to ensure the group information is added
plot(pca_value[,1:2],col=pca_value$group)
#use the first two PCs to draw the graph because they accountfor the most variation in the original data
#five “wt” samples cluster on the left side of PC1 and five “ko”samples cluster on the right side
pca$sdev
[1] 9.393123e+001.728837e+00 1.585614e+00 1.379070e+00 1.295212e+00 1.014360e+00
[7] 7.906336e-017.636112e-01 6.705670e-01 2.995201e-15
#sdev store the standard deviation each PC accounts for
pca.var <- pca$sdev^2
#squre “sdev” to calculate the variation each PC accounts for
pca.var.per <- round(pca.var/sum(pca.var)*100,1)
pca.var.per
[1] 88.2 3.0 2.5 1.9 1.7 1.0 0.6 0.6 0.4 0.0
#calcualte the percentage of variation each PC accounts for in the total
variation of the original data
barplot(pca.var.per,main="Scree Plot",xlab="Principal
Component",ylab="Precent Variation",ylim=c(0,100),names=colnames(pca$x))
#PC1 account for almost all variation, so there is a big differencebetween the two clusters
#the first two PCs account for over 90% of the variance, so it isrecommend to use the first two PCs to represent the whole variance
#screeplot(pca,ylim=c(0,100)) can draw the same figure
library(ggplot2)
#load ggplot2 package to re-draw the graph in ggplo2
pca.data <-
data.frame(Sample=rownames(pca$x),X=pca$x[,1],Y=pca$x[,2])
#format the data the way ggplot2 likes it
head(pca.data)
Sample X Y
wt1 wt1 -9.000989 3.72325653
wt2 wt2 -8.773236 -1.62190849
wt3 wt3 -8.955221 -1.21169551
wt4 wt4 -8.743347 -0.85969323
wt5 wt5 -9.076914 -0.02029153
ko1 ko1 8.687179 -1.31832099
#look at the first 6 lines of the new format
ggplot(pca.data,aes(x=X,y=Y))+
geom_text(aes(label=Sample))+
xlab(paste0("PC1-",pca.var.per[1],"%"))+
ylab(paste0("PC2-",pca.var.per[2],"%"))
#the result is the same as in the basic graph where five “wt” samplescluster one the left side while five “ko” samples cluster one the right side
pca$rotation[1:3,]
PC1 PC2 PC3 PC4 PC5 PC6
gene1 -0.10617570.002900561 -0.001113391 -0.0009577261 -0.029006816 -0.01596844
gene2 -0.10497220.014801401 0.009529656 0.0025598752 -0.024295101 -0.10285015
gene3 0.1062584 0.018544811 0.026482083 -0.0003765003 0.008739291 0.01372661
PC7 PC8 PC9 PC10
gene1 0.02447527-0.01712789 -0.083049168 -0.10197924
gene2 0.03105548 0.15132017 -0.048253454 0.02736120
gene3 0.02407890 0.02355022 0.001263962 -0.04512896
#rotation contains the loading scores of each gene i.e. theirinfluence on each PC.
loading_scores <- pca$rotation[,1]
#get the loading scores for PC1
gene_scores <- abs(loading_scores)
#we want to see which genes have the largest influence on PC1i.e. the genes that have largest absolute loading scores
gene_score_ranked <- sort(gene_scores,decreasing=T)
#rank the absolute loading scores from high to low
top_10_genes <- names(gene_score_ranked[1:10])
#get the top 10 genes that influence PC1 most
pca$rotation[top_10_genes,1]
gene98 gene65 gene71 gene14 gene90 gene28 gene25
0.1064237 0.1064129 -0.1063941 0.1063901 0.1063806 0.1063521 0.1063429
gene49 gene12 gene54
0.1063260 -0.1063250 -0.1063234
#the genes that havepositive values push “ko” samples to the right side of PC1 while those havingnegative values push “wt” samples to the right side of PC1
2.2 用FactoMineR和factoextra包分析主成分
library(FactoMineR)
library(factoextra)
#use these two package to do PCA
pca <- PCA(t(data.cpm), graph = FALSE)
#do PCA on the data matrix
#by default, PCA() expects each row is a sample and eachcolumn is a gene, so here t() is used to transpose the matrix
#by default, the parameter “graph” is true, here we setit false because the graph is hard to understand and very ugly (as shown in thefollowing), so we only want the pca results
class(pca)
[1] "PCA" "list"
names(dat.pca)
[1] "eig" "var" "ind" "svd" "call"
#the result of PCA is a listcontaining 5 elements
pca$eig
eigenvalue percentage of variancecumulative percentage of variance
comp 188.3913754 88.3913754 88.39138
comp 2 2.9364305 2.9364305 91.32781
comp 3 2.6361069 2.6361069 93.96391
comp 4 1.7642427 1.7642427 95.72816
comp 5 1.6438299 1.6438299 97.37199
comp 6 1.0381624 1.0381624 98.41015
comp 7 0.6423611 0.6423611 99.05251
comp 8 0.5536955 0.5536955 99.60620
comp 9 0.3937957 0.3937957 100.00000
#”eig” contains the percentageand accumulative percentage of variance each PC accounts for
barplot(pca$eig[,1],main="Scree Plot",xlab="Principal
Component",ylab="Precent Variation",ylim=c(0,100),names=rownames(pca$eig))
#PC1 account for almost all variation, so there is a big differencebetween the two clusters
#the first two PCs account for over 90% of the variance, so it isrecommend to use the first two PCs to represent the whole variance
#the result is the same as we did using prcomp()
fviz_pca_ind(pca,
geom.ind ="point",#represent each sample with a point
col.ind = group,#colored according to
group
addEllipses = TRUE,
legend.title ="Groups"
)
#plot the pca results in a more clear graph, the same as whatwe saw in ggplot2 figure
2.3 主成分分析的局限
PCA本质上是将方差最大的方向作为主要特征,并且在各个正交方向上将数据“离相关”,也就是让它们在不同正交方向上没有相关性。
2.3.1 不适用于非线性的相关性
PCA是一种线性降维算法,它可以很好地拆解线性相关,但是对于特征与特征之间的复杂多项式关系(非线性的高阶相关性)就没有办法了。对于存在高阶相关性的数据,可以考虑Kernel PCA,通过Kernel函数将非线性相关转为线性相关。但是,不是所有的非线性相关都能转化成线性相关,在这种情况下使用PCA可能会导致欠拟合的情形发生。
2.3.2 不能解释非正交方向上的变异
另外,PCA假设数据各主要变化是分布在正交方向上,如果在非正交方向上存在几个方差较大的方向,PCA的效果就大打折扣了。
2.3.3 不适用于PC1和PC2解释变异较少的情况
在实际的数据分析中,PC1和PC2所能解释的变异占数据总变异的比例往往很低,只用PC1和PC2来展示数据,更多维度的信息会被丢失。而且,代表细胞的点在空间上的分布比较散,不容易得到对这些细胞簇的直观的感受(如下图)。
这是因为,线性降维算法的一个主要问题是当它们集中将不相似的数据点放置在较低维度区域时,数据点相距较远。但是为了在低维、非线性曲面的流形上表示高维数据,我们也需要把相似的数据点放在一起展示,这不是线性降维算法所能做的。
参考资料:
PCA降维原理
//www.greatytc.com/p/ab89a675c00a
降维算法之PCA(主成分分析)--无监督
//www.greatytc.com/p/4d14546d020a
StatQuest_ Principal Component Analysis (PCA), Step-by-Step
https://www.bilibili.com/video/BV1sb411s7Ld/?spm_id_from=333.337.search-card.all.click&vd_source=ede39469cc5ecafb5afeaa8bbd3a2985
Principal Component Analysis (PCA) clearly explained (2015) StatQuest
https://www.bilibili.com/video/BV1WB4y1U7io/?spm_id_from=333.337.search-card.all.click&vd_source=ede39469cc5ecafb5afeaa8bbd3a2985
StatQuest PCA - Practical Tips
https://www.bilibili.com/video/BV1NF411g79V/?spm_id_from=333.337.search-card.all.click&vd_source=ede39469cc5ecafb5afeaa8bbd3a2985
StatQuest_ PCA main ideas in only 5 minutes!!!
https://www.bilibili.com/video/BV1st411q7iY/?spm_id_from=333.337.search-card.all.click&vd_source=ede39469cc5ecafb5afeaa8bbd3a2985
StatQuest_ PCA in R
https://www.bilibili.com/video/BV1st411q7r4/?spm_id_from=333.337.search-card.all.click&vd_source=ede39469cc5ecafb5afeaa8bbd3a2985