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