本节内容属于思农数据工程产物,欢迎大家批评指正。
opls函数介绍
紧接着我们对于其在r中的实现进行学习,即基于opls函数,其基本形式如下,需要重点关注的参数为predI、orthoI、fig.pdfC;其中predI表示用来展示与Y相关的正交变量数量,最多为10个,设置为NA表示其根据情况自动调整,默认为NA;orthoI表示正交垂直变量数量即用来展示与Y不相关的差异的主成分个数,其数值决定此函数执行何种分析,当它设置为0时,表示进行PLS-DA分析,设置为非0整数时表示进行OPLS-DA分析,其最多为9个,NA代表自动调整,默认为0;fig.pdfC和info.txtc表示可以设置图片和文件显示和存储途径,none代表不进行显示,interactive代表存储交互环境中,即r环境中,myfile.pdf/txt代表指定输出路径,默认为interactive。
特别注意:设置导出为pdf格式且为特定导出路径后,若出现rstudio中plot里后续运行其他可视化代码不出图的情况,可先运行dev.off()退出导出pdf,后续可视化即可恢复正常;
opls(
x,
y = NULL,
predI = NA,
orthoI = 0,
algoC = c("default", "nipals", "svd")[1],
crossvalI = 7,
log10L = FALSE,
permI = 20,
scaleC = c("none", "center", "pareto", "standard")[4],
subset = NULL,
plotSubC = NA,
fig.pdfC = c("none", "interactive", "myfile.pdf")[2],
info.txtC = c("none", "interactive", "myfile.txt")[2]
)
不同分析实现代码:
#PCA,仅需导入代谢物相对丰度组成矩阵
opls(data)
#PLS-DA,需导入代谢物相对丰度组成矩阵、分组文件
opls(data,group)
#OPLS-DA,导入代谢物相对丰度组成矩阵、分组文件的同时需调整orthoI非0,正整数或NA,这里以NA举例
opls(data,group,orthoI=NA)
OPLS-DA结果展示与解释
对于opls-da几种可视化方式进行学习,同时了解其图片背后的生态学意义,主要参考
https://mp.weixin.qq.com/s/XASA-pOlTccMK1YS9le2dQ、https://mp.weixin.qq.com/s/R7cjdpQmqImfSbnymezLjg
热图设置分组标签颜色
pheatmap:color <- list(Group=c(Group0="#ED7576",Group1="#86B1D4"),xx=c().....)
patchwork拼图报错
主要原因ggplot2更新至3.5.0导致patchwork不能正常使用,需将其退回至版本3.4.4,方可正常使用。
ggClusterNet::Robustness.Random.removal函数调整
在进行鲁棒性运算时报以上错误,初始时不知所云,故对函数进行逐行检查,发现了潜在的一些问题,并对其进行修改。原函数部分代码:
我们可以注意到14行物种的平均相对丰度是基于原始的otutab来进行计算的,而网络矩阵是基于相关矩阵的来进行筛选的,同时更新的物种的相对丰度向量也是基于相关矩阵来进行筛选的,这可能会导致一个问题,就是如果原始otutab和相关性矩阵物种顺序是否一致,如若不一致,将会导致sp.ra2中的物种无法和网络矩阵相匹配,我们对这一可能的问题进行检查:
首先检查otutab物种排列顺序:
其次检查相关性矩阵排列顺序:
显然,随着对前500物种的筛选,导致相关性物种矩阵顺序发生改变,因此sp.ra2=sp.ra[colSums(abs(cor))>0]这一代码将会导致错误的保留结果,无法与网络矩阵对应,因此,我们调整otutab的物种顺序与相关性矩阵一致后再计算sp.ra,同时因为相关性矩阵中保留物种并非全部物种,因此需保留原始otutab进行相对丰度计算,故对以上函数部分作出以下修改:
#去除第14行
#25行后加入以下代码
otutab1 <- otutab[,rownames(cor)]
sp.ra <- colMeans(otutab1)/mean(rowSums(otutab))
#继续运行剩下代码
#修改后 sum(row.names(network.raw) == names(sp.ra2))输出结果为网络矩阵行数
#加权鲁棒性结果更正
作者:思农生信团队