R实战 | OPLS-DA（正交偏最小二乘判别分析)筛选差异变量(VIP)及其可视化

主成分分析（PCA）无监督降维方法，能有效处理高维数据。PCA对相关性小的变量不敏感PLS-DA(偏小二乘判别分析)能有效解决这个问题。OPLS-DA(正交偏最小二乘分析)结合正交信号PLS-DA筛选差异变量。

“

本分析主要用于代谢组筛选不同代谢物。

数据集

液相色谱高分辨质谱法（LTQ Orbitrap）对183名成年人的尿样进行了分析。

sacurine list 它包含三个数据矩阵：

dataMatrix样本-代谢物含量矩阵(log10转换)记录了各种代谢物在各种样本中的含量信息。共有183个样本(行)和109个代谢物(列)。

sampleMetadata年零、体重、性别等信息记录在183个样本中。

variableMetadata注释109种代谢物的细节，MSI level水平。

rm(list=ls()) #loadpackages library(ropls) #loaddata data(sacurine) #查看数据集 head(sacurine$dataMatrix[,1:2]) head(sacurine$sampleMetadata) head(sacurine$variableMetadata) #提取性别分类 genderFc=sampleMetadata[,"gender"]

>head(sacurine$dataMatrix[,1:2]) (2-methoxyethoxy)propanoicacidisomer(gamma)Glu-Leu/Ile HU_0113.0197663.888479 HU_0143.8143394.277149 HU_0153.5196914.195649 HU_0172.5621834.323760 HU_0183.7819224.629329 HU_0194.1610744.412266

>head(sacurine$sampleMetadata) agebmigender HU_0112919.75M HU_0145922.64F HU_0154222.72M HU_0174123.03M HU_0183420.96M HU_0193523.41M

OPLS-DA

#以性别为例 #通过orthoI指定正交组分数目 #orthoI=NA时，执行OPLS，合适的正交组分数通过交叉验证自动计算 oplsda=opls(dataMatrix,genderFc,predI=1,orthoI=NA)

OPLS-DA 183samplesx109variablesand1response standardscalingofpredictorsandresponse(s) R2X(cum)R2Y(cum)Q2(cum)RMSEEpreortpR2YpQ2 Total0.2750.730.6020.262120.050.05

结果中，R2X和R2Y分别表示所建模型对X和Y矩阵的解释率，Q2表示模型的预测能力，其值越接近1，模型的拟合度越好，训练集的样本就越能准确地划分为原始归属。

• Inertia(惯量)柱形图(左上)

  展示三个正交轴R2Y和Q2Y。评估正交组分是否足以通过显示累释率来评估。

• 显着性诊断(右上)

  实际和模拟模型R2Y和Q2Y值经随机排列后的散点图，模型R2Y和Q2Y(散点)大于真实值时(横线)，表示产生过拟合2。右上图，OPLS-DA模型的R2Y和Q2Y比较随机置换数据后获得的相应值。

离群点展示(左下)

展示了各样本在投影平面内以及正交投影面的距离，具有高值的样本标注出名称，表明它们与其它样本间的差异较大。颜色代表性别分组。

• x-score plot(右下)

  各样本在OPLS-DA轴中的坐标，颜色代表性别分组。

可视化

library(ggplot2)
library(ggsci)
library(tidyverse)
#提取样本在 OPLS-DA 轴上的位置
sample.score = oplsda@scoreMN %>%  #得分矩阵
  as.data.frame() %>%
  mutate(gender = sacurine[["sampleMetadata"]][["gender"]],
         o1 = oplsda@orthoScoreMN[,1]) #正交矩阵
head(sample.score)#查看

> head(sample.score)
              p1 gender         o1
HU_011 -1.582933      M -4.9806037
HU_014  1.372806      F -1.7443382
HU_015 -3.341370      M -3.4372771
HU_017 -3.590063      M -0.9794960
HU_018 -1.662716      M  0.3155845
HU_019 -2.312923      M  0.6561281

p <- ggplot(sample.score, aes(p1, o1, color = gender)) +
  geom_hline(yintercept = 0, linetype = 'dashed', size = 0.5) + #横向虚线
  geom_vline(xintercept = 0, linetype = 'dashed', size = 0.5) +
  geom_point() +
  #geom_point(aes(-10,-10), color = 'white') +
  labs(x = 'P1(5.0%)',y = 'to1') +
  stat_ellipse(level = 0.95, linetype = 'solid', 
               size = 1, show.legend = FALSE) + #添加置信区间
  scale_color_manual(values = c('#008000','#FFA74F')) +
  theme_bw() +
  theme(legend.position = c(0.1,0.85),
        legend.title = element_blank(),
        legend.text = element_text(color = 'black',size = 12, family = 'Arial', face = 'plain'),
        panel.background = element_blank(),
        panel.grid = element_blank(),
        axis.text = element_text(color = 'black',size = 15, family = 'Arial', face = 'plain'),
        axis.title = element_text(color = 'black',size = 15, family = 'Arial', face = 'plain'),
        axis.ticks = element_line(color = 'black'))
p

差异代谢物筛选

#VIP 值帮助寻找重要的代谢物
vip <- getVipVn(oplsda)
vip_select <- vip[vip > 1]    #通常以VIP值>1作为筛选标准
head(vip_select)

vip_select <- cbind(sacurine$variableMetadata[names(vip_select), ], vip_select)
names(vip_select)[4] <- 'VIP'
vip_select <- vip_select[order(vip_select$VIP, decreasing = TRUE), ]
head(vip_select)    #带注释的代谢物，VIP>1 筛选后，并按 VIP 降序排序

> head(vip_select)   
                               msiLevel      hmdb chemicalClass
p-Anisic acid                         1 HMDB01101        AroHoM
Malic acid                            1 HMDB00156        Organi
Testosterone glucuronide              2 HMDB03193 Lipids:Steroi
Pantothenic acid                      1 HMDB00210        AliAcy
Acetylphenylalanine                   1 HMDB00512        AA-pep
alpha-N-Phenylacetyl-glutamine        1 HMDB06344        AA-pep
                                    VIP
p-Anisic acid                  2.533220
Malic acid                     2.479289
Testosterone glucuronide       2.421591
Pantothenic acid               2.165296
Acetylphenylalanine            1.988311
alpha-N-Phenylacetyl-glutamine 1.965807

#对差异代谢物进行棒棒糖图可视化
#代谢物名字太长进行转换
vip_select$cat = paste('A',1:nrow(vip_select), sep = '')
p2 <- ggplot(vip_select, aes(cat, VIP)) +
  geom_segment(aes(x = cat, xend = cat,
                   y = 0, yend = VIP)) +
  geom_point(shape = 21, size = 5, color = '#008000' ,fill = '#008000') +
  geom_point(aes(1,2.5), color = 'white') +
  geom_hline(yintercept = 1, linetype = 'dashed') +
  scale_y_continuous(expand = c(0,0)) +
  labs(x = '', y = 'VIP value') +
  theme_bw() +
  theme(legend.position = 'none',
        legend.text = element_text(color = 'black',size = 12, family = 'Arial', face = 'plain'),
        panel.background = element_blank(),
        panel.grid = element_blank(),
        axis.text = element_text(color = 'black',size = 15, family = 'Arial', face = 'plain'),
        axis.text.x = element_text(angle = 90),
        axis.title = element_text(color = 'black',size = 15, family = 'Arial', face = 'plain'),
        axis.ticks = element_line(color = 'black'),
        axis.ticks.x = element_blank())
p2

参考

1. OPLS-DA在R语言中的实现 | 小蓝哥的知识荒原 (blog4xiang.world)

2. R包ropls的偏最小二乘判别分析（PLS-DA）和正交偏最小二乘判别分析（OPLS-DA）

3. 用PLS和OPLS分析代谢组数据 - 简书 (jianshu.com)

4. ropls: PCA, PLS(-DA) and OPLS(-DA) for multivariate analysis and feature selection of omics data (bioconductor.org)

往期

1. 单组学的多变量分析|1.PCA和PLS-DA

2. 单组学的多变量分析| 2.稀疏偏最小二乘判别分析（sPLS-DA）