数据集
液相色谱高分辨质谱法（LTQ Orbitrap）分析了来自183位成人的尿液样品。
sacurine list 包含了三个数据矩阵：
dataMatrix为样本-代谢物含量矩阵(log10转换过)，记录了各种类型的代谢物在各样本中的含量信息。共计183个样本（行）以及109种代谢物（列）。
sampleMetadata中记录了183个样本所来源个体的年零、体重、性别等信息。
variableMetadata为109种代谢物的注释详情，MSI level水平。
rm(list = ls()) # load packages library(ropls) # load data data(sacurine) #查看数据集 head(sacurine$dataMatrix[ ,1:2]) head(sacurine$sampleMetadata) head(sacurine$variableMetadata) #提取性别分类 genderFc = sampleMetadata[, "gender"]
> head(sacurine$dataMatrix[ ,1:2]) (2-methoxyethoxy)propanoic acid isomer (gamma)Glu-Leu/Ile HU_011 3.019766 3.888479 HU_014 3.814339 4.277149 HU_015 3.519691 4.195649 HU_017 2.562183 4.323760 HU_018 3.781922 4.629329 HU_019 4.161074 4.412266
> head(sacurine$sampleMetadata) age bmi gender HU_011 29 19.75 M HU_014 59 22.64 F HU_015 42 22.72 M HU_017 41 23.03 M HU_018 34 20.96 M HU_019 35 23.41 M
OPLS-DA
# 分组以性别为例 # 通过orthoI指定正交组分数目 # orthoI = NA时，执行OPLS，并通过交叉验证自动计算适合的正交组分数 oplsda = opls(dataMatrix, genderFc, predI = 1, orthoI = NA)
OPLS-DA 183 samples x 109 variables and 1 response standard scaling of predictors and response(s) R2X(cum) R2Y(cum) Q2(cum) RMSEE pre ort pR2Y pQ2 Total 0.275 0.73 0.602 0.262 1 2 0.05 0.05
Snipaste_2021-10-28_21-32-57
结果中，R2X和R2Y分别表示所建模型对X和Y矩阵的解释率，Q2表示模型的预测能力，它们的值越接近于1表明模型的拟合度越好，训练集的样本越能够被准确划分到其原始归属中。
Inertia(惯量)柱形图(左上)
展示了3个正交轴的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
Snipaste_2021-10-28_22-49-44
差异代谢物筛选
#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
Snipaste_2021-10-28_23-35-09
参考
OPLS-DA在R语言中的实现 | 小蓝哥的知识荒原 (blog4xiang.world)
R包ropls的偏最小二乘判别分析（PLS-DA）和正交偏最小二乘判别分析（OPLS-DA）
用PLS和OPLS分析代谢组数据 - 简书 (jianshu.com)
ropls: PCA, PLS(-DA) and OPLS(-DA) for multivariate analysis and feature selection of omics data (bioconductor.org)
