R包ropls的偏最小二乘判別分析（PLS-DA）和正交偏最小二乘判別分析（OPLS-DA）
R包ropls的PCA、PLS-DA和OPLS-DA
在代謝組學分析中經常可以見到主成分分析（PCA）、偏最小二乘判別分析（partial least-squares discrimination analysis，PLS-DA）、正交偏最小二乘判別分析（orthogonal partial least-squares discrimination analysis，OPLS-DA）等分析方法，目的為區分樣本差異，或在海量數據中挖掘潛在標志物。PCA是最常見的基於特征分解的降維方法，關於它的相關概念可參考
PCA概述。PCA是一種無監督的模式，屬於探索性分析。與PCA不同的是，PLS-DA和OPLS-DA則是有監督的模式，屬於模型的方法。它們使用偏最小二乘回歸建立代謝物表達量與樣本類別之間的關系模型，對數據降維，這種監督模式通常可以更好地確立樣本關系，如下圖所示這樣，無監督的PCA無法很好地區分組間樣本時，而PLS-DA則實現有效分離。除了降維數據外，PLS-DA和OPLS-DA還可實現對樣品類別的預測（即用於分類），通過構建分類預測模型，可進一步用於識別更多的樣本所屬，這是探索性的PCA方法無法做到的。
PLS-DA和OPLS-DA中涉及到兩個矩陣：X矩陣為樣本-變量觀測矩陣，Y矩陣為樣本類別歸屬矩陣。通過X和Y矩陣進行建模，即通過樣本-變量關系確立樣本關系。
兩種方法相比，偏最小二乘（PLS）是一種基於預測變量和響應變量之間協方差的潛在變量回歸方法，已被證明可以有效地處理具有多共線性預測變量的數據集。正交偏最小二乘（OPLS）則分別對與響應相關且正交的預測變量的變化進行建模。將它們與判別分析結合，即分別為PLS-DA和OPLS-DA。
本篇簡介R包ropls的PLS-DA和OPLS-DA方法。
數據集
液相色譜高分辨質譜法（LTQ Orbitrap）分析了來自183位成人的尿液樣品，共計109種代謝物能夠注釋到MSI level 1或2水平。
#Bioconductor 安裝 ropls
BiocManager::install('ropls')
#加載
library(ropls)
#數據集，詳情 ?sacurine
data(sacurine)
head(sacurine$dataMatrix[ ,1:6])
head(sacurine$sampleMetadata)
head(sacurine$variableMetadata)
dataMatrix為樣本-代謝物含量矩陣，記錄了各種類型的代謝物在各樣本中的含量信息，根據質譜中代謝物的峰值強度等換算得到。共計183個樣本（行）以及109種代謝物（變量，列）。
sampleMetadata中記錄了183個樣本所來源個體的年零、體重、性別等信息。
variableMetadata為109種代謝物的注釋詳情，MSI level水平。
接下來期望根據各樣本中代謝物組成的不同，評估樣本間整體的差異，預測樣本的來源個體屬性（如性別），以及鑒定重要的代謝物等。
ropls包的PCA
首先可執行PCA探索各樣本間代謝物組成的差異。
#PCA，詳情 ?opls
sacurine.pca <- opls(x = sacurine$dataMatrix)
sacurine.pca
左上圖展示了各PCA軸的特征值，顯示前2-3軸承載了較多的方差，可通過繪制前2-3軸的排序圖觀測樣本差異。
右上方圖展示了各樣本在投影平面內以及正交投影面的距離，具有高值的樣本標注出名稱，表明它們與其它樣本間的差異較大。
左下圖為PCA前兩軸中的樣本得分，即各樣本在PC1和PC2軸中的排序坐標，可據此評估各樣本在代謝物組成上的差異。
右下圖為PCA前兩軸中的變量得分，即各代謝物變量在PC1和PC2軸中的排序坐標，邊緣處的變量表示它們在各樣本中的含量差別明顯（如在某些樣本中具有較大/較小的極端值等），即對排序空間的貢獻較大，暗示它們可能為一些重要的代謝物。
#opls() 返回結果是 S4 對象結構，可通過 @ 在其中提取所需結果
names(attributes(sacurine.pca)) #查看包含信息
head(sacurine.pca@scoreMN) #例如，提取對象得分（各樣本在 PCA 軸中的排序坐標）
head(sacurine.pca@loadingMN) #例如，提取變量得分（各代謝物在 PCA 軸中的排序坐標）
#如在排序圖中根據個體屬性（性別、年齡等）給樣本上色
par(mfrow = c(1, 2))
plot(sacurine.pca, typeVc = 'x-score', parAsColFcVn = sacurine$sampleMetadata[, 'gender'])
plot(sacurine.pca, typeVc = 'x-score', parAsColFcVn = sacurine$sampleMetadata[, 'age'])
#也可根據提取坐標自定義作圖展示，如 ggplot2
#展示前兩軸的樣本分布，並按個體性別着色
library(ggplot2)
scoreMN <- sacurine.pca@scoreMN
scoreMN <- cbind(scoreMN, sacurine$sampleMetadata)
scoreMN$samples <- rownames(scoreMN)
ggplot(scoreMN, aes(p1, p2, color = gender)) +
geom_point() +
stat_ellipse(show.legend = FALSE)
ropls包的PLS-DA
對於本文的示例，盡管通過代謝物含量的PCA能夠幫助我們識別一些重要的代謝物，以及探索樣本差異，但在評估組間區別則似乎有些欠佳：例如不同性別或年齡人群的代謝物差異在PCA圖中比較難以區分。
現在我們更換為帶監督的PLS-DA，查看代謝物對個體性別的定性響應。
#PLS-DA，詳情 ?opls
#監督分組以性別為例，orthoI = 0 時執行 PLS-DA
sacurine.plsda <- opls(x = sacurine$dataMatrix, y = sacurine$sampleMetadata[, 'gender'], orthoI = 0)
sacurine.plsda
這一步直接將全部數據集輸入構建預測模型，即將所有數據均作為訓練集使用。如果最終目的不是建模預測更多未知歸屬的數據，而是只為識別給定數據集的特征，那么這樣是完全可以的。
結果中，R2X和R2Y分別表示所建模型對X和Y矩陣的解釋率，Q2標示模型的預測能力，它們的值越接近於1表明模型的擬合度越好，訓練集的樣本越能夠被准確划分到其原始歸屬中。
性別監督的PLS-DA模型。
左上圖，展示了3個正交軸的R2Y和Q2Y。
由上圖，PLS-DA模型的R2Y和Q2Y與隨機置換數據后獲得的相應值進行比較。
左下圖，展示了各樣本在投影平面內以及正交投影面的距離，具有高值的樣本標注出名稱，表明它們與其它樣本間的差異較大。顏色代表性別分組。
右下圖，各樣本在PLS-DA軸中的坐標，顏色代表性別分組。我們可以看到，相對於上文的PCA（僅通過方差特征值分解），PLS-DA在區分組間差異時更有效（帶監督的偏最小二乘判別分析）。圖的下方還提供了R2X、R2Y等值，用於評估模型優度。
對於需要的結果提取和可視化。
#opls() 返回對象的提取方法和上文 PCA 一致
names(attributes(sacurine.plsda)) #查看包含信息
head(sacurine.plsda@scoreMN) #例如，樣本在 PLS-DA 軸上的位置
head(sacurine.plsda@loadingMN) #例如，代謝物在 PLS-DA 軸上的位置
#默認 plot() 作圖，如查看樣本差異以及幫助尋找重要的代謝物
par(mfrow = c(1, 2))
plot(sacurine.plsda, typeVc = 'x-score', parAsColFcVn = sacurine$sampleMetadata[, 'gender'])
plot(sacurine.plsda, typeVc = 'x-loading')
#其它可視化方法，提取坐標后自定義作圖（如 ggplot2 等），不再多說
類似上文的PCA，上圖右圖中指示了一些重要的代謝物變量，根據它們的擬合度貢獻推斷。
此外，還可通過變量投影重要度（Variable Importance for the Projection，VIP）衡量各代謝物組分含量對樣本分類判別的影響強度和解釋能力，輔助標志代謝物的篩選。通常以VIP值>1作為篩選標准。
#VIP 值幫助尋找重要的代謝物
vipVn <- getVipVn(sacurine.plsda)
vipVn_select <- vipVn[vipVn > 1] #這里也通過 VIP>1 篩選下
head(vipVn_select)
vipVn_select <- cbind(sacurine$variableMetadata[names(vipVn_select), ], vipVn_select)
names(vipVn_select)[4] <- 'VIP'
vipVn_select <- vipVn_select[order(vipVn_select$VIP, decreasing = TRUE), ]
head(vipVn_select) #帶注釋的代謝物，VIP>1 篩選后，並按 VIP 降序排序
如上所述，建模的另一目的為預測更多未知歸屬的數據。
如上過程直接使用所有數據構建了預測模型，目的僅為識別已知數據集的分類特征以及鑒定重要的代謝物，模型自身的擬合優度通過R2X、R2Y、Q等值評估。當然作為擴展，如果存在更多的其它樣本，則也可以使用構建好的模型進一步預測它們的分類（這里為根據代謝物組成判斷樣本來源個體的性別），測試模型性能以及確定更多樣本的歸屬。
由於沒有更多樣本了，因此接下來為了演示這種預測功能，將原數據集分為兩部分，一部分數據運行PLS-DA構建模型，並通過另一部分數據測試擬合模型。
#帶模型預測性能評估的執行命令
#如下示例，奇數行的數據用於構建 PLS-DA，偶數行的數據測試 PLS-DA
#此外還可以自定義截取數據子集構建訓練集或測試集等
sacurine.plsda <- opls(x = sacurine$dataMatrix, y = sacurine$sampleMetadata[, 'gender'], orthoI = 0, subset = 'odd')
sacurine.plsda
#檢查關於訓練子集的預測
trainVi <- getSubsetVi(sacurine.plsda)
table(sacurine$sampleMetadata[, 'gender'][trainVi], fitted(sacurine.plsda))
#計算測試子集的性能
table(sacurine$sampleMetadata[, 'gender'][-trainVi], predict(sacurine.plsda, sacurine$dataMatrix[-trainVi, ]))
在這個示例中，使用了一半的數據用於構建PLS-DA，作為預測模型區分樣本分組；另一半數據使用構建好的PLS-DA作測試。結果顯示，構建好的預測模型能夠很好地通過訓練集樣本的代謝物含量信息，對樣本歸類；使用該預測模型擬合測試集，也能實現較好的分類。
ropls包的OPLS-DA
但是PLS-DA容易出現過擬合（overfitting）問題。所謂過擬合，即通過訓練集建立了一個預測模型，它在訓練集上表現出色，但通過測試集測試時卻表現不佳。過擬合是機器學習中的一個常見問題，主要出現在具有比樣本數量更多的變量數量的數據集的分析中。
相較於PLS-DA，OPLS-DA可以更好地避免過擬合現象，但與PLS-DA相比通常沒有預測性能優勢的提升。如果PLS-DA模型尚可，仍推薦PLS-DA。
盡管本文的示例數據沒有出現過擬合現象，但不影響繼續展示OPLS-DA方法。我們再通過OPLS-DA查看代謝物對個體性別的定性響應，仍然直接使用所有數據作為輸入。
#OPLS-DA，詳情 ?opls
#監督分組以性別為例，orthoI = NA 時執行 OPLS-DA
sacurine.oplsda <- opls(x = sacurine$dataMatrix, y = sacurine$sampleMetadata[, 'gender'], predI = 1, orthoI = NA)
sacurine.oplsda
性別監督的OPLS-DA模型。
左上、右上、左下圖的含義參考上文PLS-DA。
右下圖，橫坐標為預測組分，縱坐標為正交組分。並在下方展示了用於評估模型優度的統計量。
其它內容，如數據提取、可視化、重要變量評估、模型性能測試等，參考上文PLS-DA。
參考資料
https://www.bioconductor.org/packages/devel/bioc/vignettes/ropls/inst/doc/ropls-vignette.html#partial-least-squares-pls-and-pls-da
https://cn.bing.com/search?q=PLSDA%20OPLSDA&qs=n&form=QBRE&sp=-1&pq=plsda%20oplsda&sc=2-12&sk=&cvid=7BEA53B98EE647C7A3E8CEBE30A6CAB8