送书 | 主成分分析PCA
主成分分析 PCA
本節作者:劉華,中國科學技術大學
版本1.0.3,更新日期:2020年6月18日
什么是PCA(Principal Component Analysis)
相關背景
在許多領域的研究與應用中,通常需要對含有多個變量的數據進行觀測,收集大量數據后進行分析尋找規律。多變量大數據集無疑會為研究和應用提供豐富的信息,但是也在一定程度上增加了數據采集的工作量。更重要的是在很多情形下,許多變量之間可能存在相關性,從而增加了問題分析的復雜性。如果分別對每個指標進行分析,分析往往是孤立的,不能完全利用數據中的信息,因此盲目減少指標會損失很多有用的信息,從而產生錯誤的結論。
因此需要找到一種合理的方法,在減少需要分析的指標同時,盡量減少原指標包含信息的損失,以達到對所收集數據進行全面分析的目的。由于各變量之間存在一定的相關關系,因此可以考慮將關系緊密的變量變成盡可能少的新變量,使這些新變量是兩兩不相關的,那么就可以用較少的綜合指標分別代表存在于各個變量中的各類信息。主成分分析與因子分析就屬于這類降維算法。
數據降維
降維就是一種對高維度特征數據預處理方法(NBT:單細胞轉錄組新降維可視化方法PHATE)。降維是將高維度的數據保留下最重要的一些特征,去除噪聲和不重要的特征,從而實現提升數據處理速度的目的。在實際的生產和應用中,降維在一定的信息損失范圍內,可以為我們節省大量的時間和成本。降維也成為應用非常廣泛的數據預處理方法。
降維具有如下一些優點:
(1) ? ?使得數據集更易使用。
(2) ? ?降低算法的計算開銷。
(3) ? ?去除噪聲。
(4) ? ?使得結果容易理解。
降維的算法有很多,比如奇異值分解(SVD)、主成分分析(PCA)、因子分析(FA)、獨立成分分析(ICA)。
PCA的概念
PCA(Principal Component Analysis),即主成分分析方法(一文看懂PCA主成分分析),是一種使用最廣泛的數據降維算法。
PCA的主要思想是將n維特征映射到k維上,這k維是全新的正交特征也被稱為主成分,是在原有n維特征的基礎上重新構造出來的k維特征。
PCA的工作就是從原始的空間中順序地找一組相互正交的坐標軸,新的坐標軸的選擇與數據本身是密切相關的。其中,第一個新坐標軸選擇是原始數據中方差最大的方向,第二個新坐標軸選取是與第一個坐標軸正交的平面中使得方差最大的,第三個軸是與第1,2個軸正交的平面中方差最大的。依次類推,可以得到n個這樣的坐標軸。通過這種方式獲得的新的坐標軸,我們發現,大部分方差都包含在前面k個坐標軸中,后面的坐標軸所含的方差幾乎為0。于是,我們可以忽略余下的坐標軸,只保留前面k個含有絕大部分方差的坐標軸。事實上,這相當于只保留包含絕大部分方差的維度特征,而忽略包含方差幾乎為0的特征維度,實現對數據特征的降維處理。
變換的步驟:
① ? ?第一步計算矩陣 X 的樣本的協方差矩陣 S(此為不標準PCA,標準PCA計算相關系數矩陣C) :
② ? ?第二步計算協方差矩陣S(或C)的特征向量 e1,e2,…,eN和特征值t = 1,2,…,N;
③ ? ?第三步投影數據到特征向量張成的空間之中。利用公式 ,其中BV值是原樣本中對應維度的值。
PCA 的目標是尋找 r ( r<n )個新變量,使它們反映事物的主要特征,壓縮原有數據矩陣的規模,將特征向量的維數降低,挑選出最少的維數來概括最重要特征。每個新變量是原有變量的線性組合,體現原有變量的綜合效果,具有一定的實際含義。這 r 個新變量稱為“主成分”,它們可以在很大程度上反映原來 n 個變量的影響,并且這些新變量是互不相關的,也是正交的。通過主成分分析,壓縮數據空間,將多元數據的特征在低維空間里直觀地表示出來。
PCA實例
PCA+特征箭頭
華大基因研究團隊發表于Microbiome的文章(Zhong et al., 2019),研究早期事件和生活方式對幼齡兒童腸道菌群和代謝的影響。在這里,我們選擇了文章中兩個PCA分析結果圖為例進行講解:
圖1c. 屬水平PCA主成分分析
Figure 1c. Genus-based principal component analysis (PCA) of children and adults.
結果:PCA分析揭示兒童和成年人腸道菌群沒有明顯區別
Principal component analysis (PCA) based on genus profiles showed no separation between Dutch children and adults.
PCA+特征箭頭+箱線圖
圖3. 多種早期事件和學齡前生活方式與腸道菌群有關。a:PCA顯示了兒童的多變量以及不同因素對PC1和PC2的主要影響。將包括早期事件和學齡前生活方式在內的18個因素進行PCA分析,其中PC1或PC2成分得分 ≥ 0.2的因素為主要影響因素。箱形圖顯示各腸型內PC1和PC2評分的總體分布(#P<0.05;Wilcoxon秩和檢驗)。
Multiple early events and pre-school lifestyle associated with the school-age gut microbiota. a PCA showing the multivariate variation of children and the major contributions of different factors to PC1 and PC2. A total of 18 factors including early events and pre-school lifestyle were subjected to PCA, and those factors with component scores for PC1 or PC2 ≥ 0.2 were shown as major contributors. Box plots showing the overall distribution of PC1 and PC2 scores within each enterotype (#P<0.05; Wilcoxon rank-sum test).
結果:我們發現母乳喂養持續時間、母親教育水平和學前飲食模式,包括攝入的蛋白質、纖維、及奶制品在PC1中貢獻最多 (15.05%),PC2中兒童的碳水化合物和脂肪總攝入量是第二個最重要的變量 (12.74%)。E3組兒童的PC1得分低于E1組,但PC2得分高于E1組(圖3a, Wilcoxon秩和檢驗,P < 0.05)。這種腸型間的差異是由PC1評分的主要因素決定的,如 E3組中較短的母乳喂養時間以及較少的膳食纖維和植物性蛋白攝入量 (Kruskal-Wallis檢驗, P < 0.05)。
We found that breastfeeding duration, educational level of mother at childbirth, and pre-school dietary patterns including intake of protein, fiber, and milk products contributed most to the variability in PC1 (15.05%, Fig. 3a), and total intake of carbohydrates and fat represented the second most important variation among children, as displayed in PC2 (12.74%, Fig. 3a). Interestingly, children in E3 exhibited lower PC1 scores but higher PC2 scores than children in E1 (Fig. 3a, Wilcoxon rank-sum test, P < 0.05). This inter-enterotype difference was governed by specific major contributors of the PC1 scores, including shorter breastfeeding duration and less intake of dietary fiber and plant-based protein in E3 as compared to the two other enterotypes (Kruskal-Wallis test, P < 0.05).
PCA+形狀分組+門著色
例3:Robert D. Finn團隊發表于Nature的文章(Almeida et al., 2019),構建了人類腸道微生物群的基因組藍圖。
Fig 5:未培養的物種具有獨特的功能。a、以已知人類參考基因組(HGR:553個基因組)和未分類物種宏基因組(UMGS:1952個基因組)的基因組特性(Genome Properties, GPs)進行PCA主成分分析,以門水平上色。
The uncultured species have a distinct functional capacity.a, Principal component analysis (PCA) based on GPs of the HGR (n = 553 genomes) and the UMGS (n = 1,952 genomes) coloured by phylum.
結果:我們使用GhostKOALA生成KEGG Orthology (KO)注釋(Pathview包:整合表達譜數據可視化KEGG通路),以跟蹤不同UMGS和HGR集合中特定功能類別的差異豐度。在全球范圍內,通過對GPs的分類組成分析,發現按門水平分類分離良好,特別是擬桿菌門和變形菌門顯示出獨特的功能特征。
In parallel, we used GhostKOALA to generate KEGG Orthology (KO) annotations to track the differential abundance of specific functional categories across the UMGS and HGR sets. Globally, by analyzing the repertoire of GPs according to the taxonomic composition, we observed a good separation by phylum (ANOSIM R = 0.42, P < 0.001), with the Bacteroidetes and Proteobacteria taxa in particular displaying very distinctive functional profiles (Fig. 5a).
總結
PCA主成分圖中坐標軸PC1/2的數值為總體差異的解釋率;圖中點代表樣品,顏色代表分組;箭頭代表原始變量,其中方向代表原始變量與主成分的相關性,長度代表原始數據對主成分的貢獻度。
做PCA,首先要構建特征/變量的協方差矩陣,然后對其特征值和特征向量進行排序,根據需要取前面最重要的部分,將后面的維數省去,可以達到降維,從而達到簡化模型或對數據進行壓縮的效果,同時最大程度的保持了原有數據的信息。
但是PCA原理主要是為了消除變量之間的相關性,并且假設這種相關性是線性的,對于非線性的依賴關系則不能得到很好的結果。同時PCA假設變量服從高斯分布,當變量不服從高斯分布(如均勻分布)時,會發生尺度縮放與旋轉。
PCA繪圖實戰
數據和代碼下載:以 https://github.com/YongxinLiu/MicrobiomeStatPlot 上的微生物組數據進行展示。
用于計算PCA 的R軟件中提供了來自不同軟件包的多個函數:
prcomp()和princomp() [內置];
PCA() [ FactoMineR包];
dudi.pca() [ ade4包];
epPCA() [ ExPosition包]。
可以通過factoextraR包和ggbiplot包來輕松提取和可視化PCA的結果。
安裝PCA分析與可視化R包
判斷每個依賴的包是否存在,沒有則安裝:
# 安裝CRAN來源R包,多個包使用循環檢測和安裝 p_list = c("FactoMineR", "dplyr", "factoextra", "ggpubr", "pca3d") for(p in p_list){if (!requireNamespace(p, quietly = TRUE))install.packages(p) }# 安裝github來源R包 suppressWarnings(suppressMessages(library(devtools))) if (!requireNamespace("ggbiplot", quietly = TRUE))install_github("vqv/ggbiplot")內置數據演示PCA繪制
library("dplyr") # 查看鳶尾花的內置數據 head(iris, n=3) # 獲得純數值表格,去除最后一行的分類型分組數據, iris <- select(iris, -Species) # prcomp函數進行PCA分析,需要解釋參數的意義???cor=T??? iris.pca <- prcomp(iris,cor = T) # 查看對象的名稱,此處返回結果中5個列表的名稱 names(iris.pca) # 對象摘要,需要解釋參數的意義???loadings = T??? summary(iris.pca,loadings = T)sdev是標準偏差; center是每列計算是減去的均值; scores即降維之后的結果,當然也可以使用predict函數,結果一樣。
我們看Proportion of Variance,即為每一個成分方差所占比例,Cumulative Proportion代表是累計比例,為Proportion of Variance的累計值,一般達到90%左右就可以代表所有數據。
library("factoextra") # 提取變量的分析結果,加載包出現在函數第一次出現前,知道函數與包的關系 get_pca_var(iris.pca)factoextra包自帶了提取變量的分析結果get_pca_var函數,其中:
coord表示用于創建散點圖的變量坐標。coord實際上就是成分載荷,指觀測變量與主成分的相關系數;
cor表示相關系數;
cos2表示因子質量,var.cos2 = var.coord * var.coord;
contrib表示包含變量對主成分的貢獻(百分比)。
圖1. PCA展示變量與主成分之間的關系,以及變量之間的關聯
菌群數據實戰
本次測試數據來自同樣來自Science:擬南芥三萜化合物特異調控根系微生物組
數據截取了3個實驗組各6個樣品的結果用于演示。數據位于Data/Science2019目錄,本次需要元數據(metadata.txt)和otu表(otutab.txt)兩個輸入文件。
貢獻率圖(scree plot)
如果我們想判斷PCA中需要多少個主成分比較好,那么可以從主成分的特征值來考慮(Kaiser-Harris準則建議保留特征值大于1的主成分);特征值表示主成分所保留的變異量(所解釋的方差);如用get_eigenvalue函數來提取特征值,結果中第一列是特征值,第二列是可解釋變異的比例,第三列是累計可解釋變異的比例。
get_eigenvalue(otu.pca)[1:3,]除了特征值大于1作為主成分個數的閾值外,還可以設置總變異的閾值(累計)作為判斷指標。
除了看表格來判斷,還可從圖形上直觀的感受下:
# 繪制崖低碎石圖(scree plot)即貢獻率圖,外層()可對保存的圖形同時預覽 (p <- fviz_eig(otu.pca, addlabels = TRUE)) ggsave(paste0("p2.pca_screen.pdf"), p, width=89, height=56, units="mm") ggsave(paste0("p2.pca_screen.png"), p, width=89, height=56, units="mm")圖2. 貢獻率圖表示各個主成分貢獻率,進而決定選擇多少主成分。可以將主成分數量限制為占總方差的比例。
如果我們想提取PCA結果中變量的信息,則可用get_pca_var():
(var <- get_pca_var(otu.pca))Quality of representation(對應var$cos2),用于展示每個變量在各個主成分中的代表性(高cos2值說明該變量在主成分中有good representation,對應在Correlation circle圖上則是接近圓周邊上;低cos2值說明該變量不能很好的代表該主成分,對應Correlation circle圖的圓心位置)。簡單的說,如果一個變量在PC1和PC2的Contributions(cos2值在各個主成分中的比例)很高的話,則說明該變量可有效解釋數據的變異,我們可以用圖形展示各個變量在PC1和PC2上的Contributions。
特征PCA圖
# 繪制變量PCA主成分分析圖 (p <- fviz_pca_var(otu.pca, col.var = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"))) ggsave(paste0("p3.pca_cos.pdf"), p, width=89, height=89, units="mm") ggsave(paste0("p3.pca_cos.png"), p, width=89, height=89, units="mm")圖3. PCA圖展示變量
可以看到,biplot圖只能用于展示變量較少的情況,當變量較多時需要進行篩選。
接下來分析觀測值,先提取出individuals信息。
(ind <- get_pca_ind(otu.pca))樣本PCA圖
然后按照上面的模式來展示下individuals的點圖,比如以cos2值來代表各個individuals點的圓圈大小
# 樣本PCA圖,點大小為cos2,形狀為21圓形,按組填充,repel避免標簽重疊 (p <- fviz_pca_ind(otu.pca, pointsize = "cos2", pointshape = 21, fill = metadata$Group, repel = TRUE)) ggsave(paste0("p4.pca_individuals.pdf"), p, width=89, height=56, units="mm") ggsave(paste0("p4.pca_individuals.png"), p, width=89, height=56, units="mm")圖4. PCA圖展示觀測值,以cos2值來代表各個individuals點的圓圈大小。
如果有分組信息,則可以將同一組的individuals圈在一起,如:
# 樣本PCA圖,只顯示點,分組著色并手動分配顏色,添加置信橢圓和圖例 (p <- fviz_pca_ind(otu.pca,geom.ind = "point", # show points only ( not "text")col.ind = metadata$Group, # color by groupspalette = c("#00AFBB", "#E7B800", "#FC4E07"),addEllipses = TRUE, # Concentration ellipseslegend.title = "Groups")) # 保存圖片,指定圖片為pdf格式方便后期修改,圖片寬89毫米,高56毫米 ggsave(paste0("p5.sample_group_ellipse.pdf"), p, width=89, height=56, units="mm") ggsave(paste0("p5.sample_group_ellipse.png"), p, width=89, height=56, units="mm")圖5. 樣本PCA圖按分組著色并添加置信橢圓
biplot樣本和特征圖
如果我們想將特征變量(vars)和樣本(individuals)同時在一張biplot圖中展示,那么就要對主要的OTUs/ASVs/分類進行排序挑選。
展示主要差異ASV與主成分的關系
# 轉換原始數據為百分比 norm <- t(t(count)/colSums(count,na=T)) * 100 # 篩選mad(median absolute deviation,中位數偏差的絕對值的中位數,衡量特異波動的方法)值大于0.5的ASV, mad.5 <- norm[apply(norm,1,mad)>0.5,] # 另一種方法:按mad值排序取前N個,如6個波動最大的ASVs # mad.5 <- head(norm[order(apply(norm,1,mad), decreasing=T),],n=6) # 計算PCA和菌與菌軸的相關性 otu.pca <- prcomp(t(mad.5)) # 繪制觀測值PCA主成分分析圖,外層()可對保存的圖形同時預覽 (p <- fviz_pca_biplot(otu.pca, col.ind = metadata$Group, palette = "jco", addEllipses = TRUE, label = "var",col.var = "black", repel = TRUE, legend.title = "Group")) # 保存圖片,指定圖片為pdf格式方便后期修改,圖片寬89毫米,高56毫米 ggsave(paste0("p6.sample_group_ellipse2.pdf"), p, width=89, height=56, units="mm") ggsave(paste0("p6.sample_group_ellipse2.png"), p, width=89, height=56, units="mm")圖6. PCA圖展示主要ASVs與主成分的關系。
我們僅用中值絕對偏差(mad)大于0.5的6個OTUs進行主成分分析,即可將三組樣品明顯分開。圖中向量長短代表差異貢獻,方向為與主成分的相關性。可以看到最長的向量ASV_2與X軸近平行,表示PC1的差異主要由此菌貢獻。其它菌與其方向相反代表OTUs間可能負相關;夾角小于90%的代表兩個OTUs有正相關。
ggbiplot包可視化PCA圖
我們也可以選擇ggbiplot包可視化PCA圖
# 加載ggbiplot并繪制觀測值PCA主成分分析圖 suppressWarnings(suppressMessages(library("ggbiplot"))) # 繪制觀測值PCA圖 ggbiplot(otu.pca, obs.scale = 1, var.scale = 1, groups = metadata$Group, ellipse = TRUE,var.axes = T) # 保存圖片,指定圖片為pdf格式方便后期修改,圖片寬89毫米,高56毫米 ggsave(paste0("p7.sample_group_ellipse3.pdf"), p, width=89, height=56, units="mm") ggsave(paste0("p7.sample_group_ellipse3.png"), p, width=89, height=56, units="mm")圖7. ggbiplot可視化PCA結果
通常情況下,需要使用PERMANOVA來檢驗不同組樣本間的微生物群落是否具有顯著差異
# 使用vegan包中的adonis函數進行PERMANOVA分析 library("vegan") otu.adonis <- adonis(t(count) ~ Group, data = metadata, permutations = 999) # 之后在繪圖代碼中將PERMANVOA結果在PCA圖中進行展示 (p1 <- p + geom_text(aes(x = - 5 , y = 4,label = paste("PERMANOVA:\n KO VS OE VS WT \n p-value = ",otu.adonis$aov.tab$`Pr(>F)`[1], sep = "")),size = 2, hjust = 0)) ggsave(paste0("p8.sample_group_ellipse4.pdf"), p1, width=120, height=100, units="mm") ggsave(paste0("p8.sample_group_ellipse4.png"), p1, width=120, height=100, units="mm")圖8. PCA結果圖添加PERMANOVA檢驗
PCA圖x/y軸添加箱線圖
在PCA圖的x和y軸添加箱線圖,可以實現進一步展示組間差異
# 需要制作PCA結果可視化的繪圖數據文件 PC1 <- otu.pca$x[,1] PC2 <- otu.pca$x[,2] otu.pca.data <- data.frame(rownames(otu.pca$x),PC1,PC2,metadata$Group) colnames(otu.pca.data) <- c("sample", "PC1", "PC2", "group") # 這里需要ggpubr包,對組間進行統計檢驗以及組合圖的拼接 library("ggpubr") # 設置比較組 my_comparisons = list(c('KO','OE'),c('OE','WT'),c('KO','WT')) # 繪制y軸為PC1值的分組箱線圖 p2 <- ggplot(otu.pca.data, aes(x = group, y= PC1, colour = group)) +geom_boxplot() +theme(panel.background = element_rect(fill = "white",colour = "white"),panel.grid = element_blank(),axis.text.y = element_blank(),legend.position="none") +xlab("") + ylab("") +stat_compare_means(comparisons = my_comparisons, label = "p.signif") # 添加顯著性檢驗 # 繪制y軸為PC2值的分組箱線圖 p3 <- ggplot(otu.pca.data, aes(x = group, y= PC2, colour = group)) +geom_boxplot(aes()) +theme(panel.background = element_rect(fill = "white",colour = "white"),panel.grid = element_blank(),axis.text.x = element_blank(),legend.position="none") +xlab("") + ylab("") +coord_flip() + # coord_flip()函數翻轉坐標軸stat_compare_means(comparisons = my_comparisons, label = "p.signif") # ggpubr::ggarrange()函數對圖進行拼接 (p4 <- ggarrange(p3, NULL, p1, p2, widths = c(5,2), heights = c(2,4), align = "hv")) ggsave(paste0("p9.sample_group_ellipse_boxplot.pdf"), p4, width=180, height=150, units="mm") ggsave(paste0("p9.sample_group_ellipse_boxplot.png"), p4, width=180, height=150, units="mm")圖9. PCA圖添加箱線圖副圖
還可以使用pca3d包對數據進行三維展示
pca3d包可視化PCA圖
我們還可以使用pca3d包對數據進行三維展示
suppressWarnings(suppressMessages(library("pca3d"))) # 繪制樣本三維PCA圖,分組著色,68%置信橢圓 (p <- pca3d(otu.pca, group = metadata$Group, show.ellipses=TRUE, ellipse.ci=0.68, show.plane=FALSE))會打開新窗口,展示三維PCA圖,而且可用鼠標托動旋轉變換觀察角度,變量p中保存了各組名、顏色 、形狀名稱和編號。
圖10. 三維PCA圖快照
注:三維PCA圖適合交互探索數據,但不能導出矢量圖,在發表文章中使用較少。
參考文獻
Huanzi Zhong, John Penders, Zhun Shi, Huahui Ren, Kaiye Cai, Chao Fang, Qiuxia Ding, Carel Thijs, Ellen E. Blaak, Coen D. A. Stehouwer, Xun Xu, Huanming Yang, Jian Wang, Jun Wang, Daisy M. A. E. Jonkers, Ad A. M. Masclee, Susanne Brix, Junhua Li, Ilja C. W. Arts & Karsten Kristiansen. (2019). Impact of early events and lifestyle on the gut microbiota and metabolic phenotypes in young school-age children. Microbiome 7, 2, doi: https://doi.org/10.1186/s40168-018-0608-z
Alexandre Almeida, Alex L. Mitchell, Miguel Boland, Samuel C. Forster, Gregory B. Gloor, Aleksandra Tarkowska, Trevor D. Lawley & Robert D. Finn. (2019). A new genomic blueprint of the human gut microbiota. Nature 568, 499-504, doi: https://doi.org/10.1038/s41586-019-0965-1
責編:劉永鑫,中科院遺傳發育所
你可能還想看
什么?你做的差異基因方法不合適?
PCA主成分分析實戰和可視化 附R代碼和測試數據
用了這么多年的PCA可視化竟然是錯的!!!
復現nature communication PCA原圖|代碼分析(一)
這篇Nature子刊文章的蛋白組學數據PCA分析竟花費了我兩天時間來重現|附全過程代碼
送書
在上周的留言送書活動中,恭喜下面這位讀者獲得書籍“Python:數據分析與可視化”,請及時與生信寶典編輯(shengxinbaodian)聯系。
本期的留言主題是:您在處理數據和學習過程中遇到了哪些值得記小本本的坑呢?歡迎轉發朋友圈并留言評論。
因有讀者反饋刷贊一事,因此更改送書規則為在留言中選擇“有緣人”獲得下面由北京大學出版社贊助的書籍(聯系小編時請附上分享截圖),內容契合留言主題且實際留言得贊越高者越可能成為“有緣人”,歡迎分享并互相監督,結果在下一期送書活動中公布:
往期精品(點擊圖片直達文字對應教程)
后臺回復“生信寶典福利第一波”或點擊閱讀原文獲取教程合集
創作挑戰賽新人創作獎勵來咯,堅持創作打卡瓜分現金大獎總結
以上是生活随笔為你收集整理的送书 | 主成分分析PCA的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 在Rstudio中点一点就出来了一个R包
- 下一篇: 微生物组-扩增子16S分析第9期(报名直