1.数据下载以及整理

setwd("E:/001微米work_Learn/work_unzip/O2PLS/")
 
set.seed(31*12*2016)
 
if(!("test.tab" %in% list.files())){
	## If you didn't download the expression data,
	## this code will download it to the current directory (getwd())
	temp <- tempfile()
	download.file(
		"http://www.ebi.ac.uk/arrayexpress/files/E-TABM-1036/E-TABM-1036.processed.1.zip",
		temp)
	rna0 <- read.table(unzip(temp, "test.tab"), sep='\t')
	unlink(temp); rm(temp)
} else {
	## Or if you've downloaded test.tab already we simply load it
	rna0 <- read.table("test.tab", sep='\t')
	
}
rna1 <- t(rna0[-(1:2),-1])
rna2 <- matrix(as.numeric(rna1), nrow = nrow(rna1))
dimnames(rna2) <- list(as.character(unlist(rna0[1,-1])),unlist(rna0[-(1:2),1]))
rna2 <- rna2[order(row.names(rna2)), ] # Order rows according to the participant ID
 
 
 
filter_rna <- function(rna=rna, prop = 0.75){
	#calculate the maximum of gene expression per each gene and take the top
	maxGE <- apply(rna, 2, max)
	propGEmax <- quantile(maxGE, prop)
	#take the IQR of each gene and take the top genes
	IQRGE <- apply(rna, 2, IQR, na.rm=TRUE)
	propGEIQR <- quantile(IQRGE, prop)
	#selected genes/probes are the intersection of the two previous sets
	filter2 <- (intersect(which(maxGE> propGEmax), which(IQRGE> propGEIQR)))
	return(filter2)
}
rna3 <- rna2[,filter_rna(rna2)]
rm(rna0)
rm(rna1)
 
 
if(!("metabonomic_data.txt" %in% list.files())){
	temp <- tempfile()
	download.file(
		"http://msb.embopress.org/content/msb/6/1/441/DC3/embed/inline-supplementary-material-3.zip",
		temp)
	metab0 <- read.table(unzip(temp, "metabonomic_data.txt"), header = T)
	unlink(temp); rm(temp)
} else {
	## Or if you've downloaded metabonomic_data.txt already run the next line
	metab0 <- read.table("metabonomic_data.txt", header=T)
}
metab1 <- t(metab0[,-1])
colnames(metab1) <- metab0$Metabolite

缺失值处理

 
 
install.packages("VIM")
install.packages("missForest")
 
 
 
VIM::aggr(metab1, col=c('navyblue','red'), numbers=TRUE, sortVars=FALSE,
		  labels=names(data), cex.axis=.7, gap=3,
		  ylab=c("Histogram of missing data","Pattern"))
 
 
NAs_in_metab1 <- which(apply(metab1, 1, function(e) sum(is.na(e))/length(e))==1)
metab2 <- metab1[-NAs_in_metab1,]
rna4 <- rna3[-NAs_in_metab1,]
 
 
metab2.imp <- missForest::missForest(metab2, verbose = T)
metab <- scale(metab2.imp$ximp, scale=F)
rna <- scale(rna4, scale = F)
 
install.packages("gplots")
 
 
gplots::heatmap.2(cor(metab1,use = 'pair'), dendrogram='none', Rowv=F, Colv=F,trace='n',
				  breaks=seq(-1,1,length.out = 25), col=gplots::bluered)
gplots::heatmap.2(cor(metab,use = 'pair'), dendrogram='none', Rowv=F, Colv=F,trace='n',
				  breaks=seq(-1,1,length.out = 25), col=gplots::bluered)
				  
saveRDS(metab, "metab.rds")
saveRDS(rna, "rna.rds")

正式分析

注意需要数据填充和数据过滤，数据过滤IQR和Max的方法对小样本更具有稳健性
可以画 boxplot看看中位数是否相同

#install.packages("OmicsPLS")
library(OmicsPLS)
metab=readRDS("metab.rds")
rna=readRDS("rna.rds")
CV1 <- crossval_o2m_adjR2(rna, metab, 1:3, c(0,1,5,10), c(0,1,5,10),
						  nr_folds = 2, nr_cores = 4)
CV2 <- crossval_o2m(rna, metab, 1:2, 0:2, 9:11,
					nr_folds = 10, nr_cores = 4)
 
fit = o2m(rna, metab, 2, 1, 10)
fit
summary(fit)

可视化

library(illuminaHumanv3.db)
 
# Color names
LLmodule <- c("ILMN_1690209",'ILMN_1766551', 'ILMN_1749131', 'ILMN_1688423',
			  'ILMN_2102670', 'ILMN_1792323', 'ILMN_1899034', 'ILMN_1806721',
			  'ILMN_1695530', 'ILMN_1726114', 'ILMN_1751625', 'ILMN_1726114',
			  'ILMN_1753648', 'ILMN_1779043')
LLnr <- which(colnames(rna) %in% LLmodule)
rna_genenames <- select(illuminaHumanv3.db,
						keys = colnames(rna)[LLnr],
						keytype = "PROBEID", columns = "SYMBOL")[,2]
 
 
name_col <- 1 + sapply( #First sapply loops over column names
	X = colnames(metab),
	FUN = function(arg){
		crossprod(
			c(1, 1, 3, 4, 5), # Weights to be used as categories
			sapply(c("VLDL", "LDL", "IDL", "HDL","FA"), # metabolite classes
				   function(arg2){grepl(arg2, arg)} # compare class of metabolites
			)
		)
	}
)
name_col <- factor(name_col,
				   levels = c(3,2,4:6,1),
				   labels = c("VLDL", "LDL", "IDL", "HDL","FA","Other"))
 
 
library(ggplot2)
library(dplyr)
p_metab <- plot(fit, loading_name="Yj", i=1, j=2, label="c", # Plot the loadings
				alpha=0) + # set points to be 100% transparant
	##################### Add all layers ###
	theme_bw() +
	coord_fixed(ratio = 1, xlim=c(-.2,.2),ylim=c(-.2,.2)) +
	geom_point( # Set color and size
		aes(col=name_col, size = I(1+(name_col%in%c("VLDL","HDL"))),
			shape = name_col),show.legend = T) +
	theme(legend.position="right") +
	scale_color_discrete(name="Metabolite\nGroup",labels=c("VLDL", "LDL", "IDL", "HDL","FA","Other")) +
	guides(size=F) + scale_shape_discrete(name="Metabolite\nGroup",
										  labels=c("VLDL", "LDL", "IDL", "HDL","FA","Other")) +
	scale_shape_manual(name="Metabolite\nGroup", values=c(15,3,4,17,5,6)) +
	labs(title = "Metabolite joint loadings",
		 x = "First Joint Loadings", y = "Second Joint Loadings") +
	theme(plot.title = element_text(face='bold'),
		  legend.title=element_text(face='bold')) +
	geom_hline(yintercept = 0) + geom_vline(xintercept = 0)
 
alprna <- loadings(fit, "Xjoint", 1:2) %>% magrittr::raise_to_power(2) %>% rowSums
alprna[-(order(alprna,decreasing=T)[1:10])] = 0
alprna <- sign(alprna)
toprna <- which(alprna>0)
names_rna <- mapIds(illuminaHumanv3.db,
					keys = colnames(rna)[toprna],
					keytype = "PROBEID",
					column = "SYMBOL",
					multiVals = 'first')
 
 
 
p_rna <- ggplot(data.frame(x = fit$W.[, 1], y = fit$W.[, 2]),
				aes(x = x, y = y),
				#alpha = alprna,
				aes(label = NA)) +
	##################### Add all layers ###
	theme_bw() +
	coord_fixed(.8, c(-.15,.15),c(-.15,.15)) +
	geom_point(alpha = 0.5, col = 'grey') +
	geom_point(data = data.frame(x=fit$W.[LLnr,1],y=fit$W.[LLnr,2]),
			   shape = 2, col = 2, size = 2) +
	geom_text(data = data.frame(x=fit$W.[toprna,1],y=fit$W.[toprna,2]),
			  hjust = rep(c(1, 0), length.out = length(toprna)),
			  aes(label = names_rna)) +
	labs(title = "Transcript joint loadings",
		 x = "First Joint Loadings", y = "Second Joint Loadings") +
	theme(plot.title = element_text(face='bold')) +
	geom_hline(yintercept = 0) + geom_vline(xintercept = 0)
## Finally plot both plots in one figure.
(p_metab+p_rna)

注意1.
loadings函数提取出的是权重系数

edit(getAnywhere(scores.o2m))
edit(getAnywhere(loadings.o2m))

这个提取的是得分

 glimpse(s)
List of 13
 $ Comp     : int 2
 $ R2_X     : num 0.295  X的解释包括联合和正交
 $ R2_Y     : num 0.653
 $ R2_Xjoint: num 0.124
 $ R2_Yjoint: num 0.411
 $ R2_Xhat  : num 0.0144
 $ R2_Yhat  : num 0.0616
 $ R2_Xpred : num 0.116
 $ R2_Ypred : num 0.15
 $ B_T      : num [1:2, 1:2] 0.399 -0.262 -0.144 0.338
 $ B_U      : num [1:2, 1:2] 0.42 0.157 0.086 0.444
 $ flags    :List of 14
  ..$ n        : num 2
  ..$ nx       : num 1
  ..$ ny       : num 10
  ..$ stripped : logi FALSE
  ..$ highd    : logi FALSE
  ..$ ssqX     : num 332036
  ..$ ssqY     : num 68387
  ..$ varXjoint: num [1:2] 25085 16214
  ..$ varYjoint: num [1:2] 18514 9590
  ..$ varXorth : num 56639
  ..$ varYorth : num [1:10] 6687 3304 2085 1193 1089 ...
  ..$ method   : chr "O2PLS"
  ..$ call     : language o2m(X = rna, Y = metab, n = 2, nx = 1, ny = 10)
  ..$ time     : Named num 0.95
  .. ..- attr(*, "names")= chr "elapsed"
 $ digits   : num 3
 - attr(*, "class")= chr "summary.o2m"
 
> 0.295-0.124  X正交的解释
[1] 0.171
> 0.653-0.411 Y正交的解释
[1] 0.242

能告诉的信息

模型能总共解释多少X ，可以看出X Y那个解释的更多
R2_Xjoint (0.124)： X 矩阵的总方差中，由联合部分 (Joint Part) 所解释的比例。只有 12.4% 的 X 变异与 Y 共享，意味着 X 中大部分变异是 X 独有的（orthogonal）
R2_Xhat (0.0144)：由 joint 成分预测 X 的能力这里是1.44%，这是理论“重建 X”，不是 cross-validation，一般用于检查模型完整性。。
R2_Xpred (0.116)：使用Y来预测X 时，预测的决定系数R²。可以理解为Y对X的预测能力。
B _T B_U回归系数矩阵
B_T：X→Y 的联合得分系数（T 预测 U）

用 X 的 joint scores (T) 预测 Y 的 joint scores (U)

表示“X 的联合变异如何映射到 Y”

B_U：Y→X 的联合得分系数（U 预测 T）

crossval_o2m 和 crossval_o2m_adjR2 都是用于模型交叉验证的函数，但它们优化的指标不同，因此会导致模型选的成分数也不同。核心区别：

✅ 1. crossval_o2m：基于预测误差（PRESS / Q²）来选择成分

优化目标通常是 预测能力（prediction accuracy）。
评价指标常为 PRESS、Q² 或类似的预测相关的统计量。
主要用于判断 模型对新数据的预测性能。

➡ 适用于希望模型有最好预测能力的场景。

✅ 2. crossval_o2m_adjR2：基于调整后的 R²（adjusted R²）

优化目标是 解释能力（explained variance），并用 调整后的 R² 惩罚过多的成分。
adjusted R² 会随着加入过多成分而下降，因此能避免过拟合。
强调模型对原始数据变异的解释，而不是预测。

➡ 适用于更关注解释性的 O2PLS 模型选择。

⭐ 为什么会有这两种方式？

O2PLS 模型可以用于预测（如根据 X 预测 Y）
也可以用于 解释数据结构（如寻找 X 与 Y 的共享结构 vs. 各自的特有成分）

不同的任务下，所需的成分数会不同：

如果用 crossval_o2m，成分可能更多，因为更多成分提高预测能力。
如果用 crossval_o2m_adjR2，成分可能更少，因为 adjusted R² 会惩罚复杂模型。

📌 简短总结

方法	主要优化指标	强调	结果特点
crossval_o2m	Q² / PRESS（预测误差）	预测性	成分可能较多
crossval_o2m_adjR2	Adjusted R²	解释性 & 避免过拟合	成分较少、更简洁

在没有任何先验知识的情况下，选择 O2PLS 中的 联合成分数（r） 以及 X-specific 成分（X） 和 Y-specific 成分（Y） 的范围，通常会依赖于 交叉验证（cross-validation） 等方法来优化。然而，设定参数范围时的合理性需要考虑到以下几个因素：

1. 参数范围的设置：r, X 和 Y 的选择

选择 r（联合成分数）：

r = 2-5 是合理的，太少的联合成分可能导致模型不能充分捕捉到 X 和 Y 数据集之间的共变结构，而太多的联合成分可能导致过拟合，尤其是在噪声较多的情况下。
在选择联合成分时，你可以通过 交叉验证（比如 10-fold CV）来选择最佳的 r 值。通常，r 值的选择会依赖于模型解释的方差以及在测试集上的预测性能。

选择 X 和 Y 的成分数（X-specific 和 Y-specific）：

X-specific 成分（X = 1:3） 和 Y-specific 成分（Y = 1:3） 也是合理的范围，尤其是当你不确定 X 和 Y 的具体关联时，设置一个较小的范围有助于限制模型复杂性，同时通过交叉验证来确定最合适的成分数。

2. 考虑数据的波动性和噪声

由于你提到的数据来源于 质谱分析（Mass Spectrometry, MS），并且是 代谢组数据和蛋白数据，质谱数据通常有较大的波动性和噪声，这对模型的选择有以下影响：

高波动性和噪声：质谱数据通常会包含很多噪声，可能有大量的 无用信号 或 技术变异。这可能导致模型捕捉到一些没有生物学意义的成分，尤其是在 X-specific 或 Y-specific 方面。为了避免噪声影响，可以考虑以下处理：
- 降噪：如果数据集中的噪声较大，可以考虑先通过 主成分分析（PCA） 或其他降维方法来减少噪声对模型的影响。
数据的波动性：质谱数据中，代谢物和蛋白质的浓度可能在不同的样本之间变化较大，这会导致 X 和 Y 之间的关系更加复杂。因此，设置 X = 1:3 和 Y = 1:3 成分数范围，可能会帮助你探索 X 和 Y 之间的潜在共变结构，而不至于过多地捕捉到由于波动性和噪声产生的误差。
r：联合成分数的范围 2-5 是合理的。通常来说，联合成分数不宜过多，因为会增加模型的复杂度和过拟合的风险。建议通过交叉验证选择最佳的 r。一般来说，在 O2PLS 中，2到3个联合成分是常见的选择，但也取决于数据本身的复杂度和噪声。
X-specific 和 Y-specific 成分数的选择（X = 1:3，Y = 1:3）：
- 对于代谢组和蛋白质数据，选择 1到3个成分 是一个合理的选择。过多的特定成分可能导致模型更难解释，而且在数据噪声较大的情况下会更容易过拟合。
- 如果你在交叉验证中发现，更多的成分（比如 X = 3，Y = 3）能有效降低误差，并且在验证集上表现更好，那么可以考虑增加这些成分数，但一般来说，1到3个成分足够捕捉数据中的主要结构。

4. 如何评估这些选择是否合理

通过以下几个步骤，你可以判断是否选择了合理的参数范围：

交叉验证：使用 交叉验证（比如 K-fold CV）来评估不同参数选择的模型性能。你可以比较不同的 r 值（联合成分数）、X-specific 和 Y-specific 成分数的模型预测误差（如 MSE 或 R²）来选择最佳的组合。
解释方差：观察 联合成分（r） 的选择是否能解释 X 和 Y 数据集的主要变异。如果你选择的 r 值过小，可能会导致模型无法充分捕捉数据的共变结构。如果选择的 r 值过大，可能会导致过拟合和对噪声的过度拟合。
生物学解释：最终，通过对结果的生物学解释，确认所选的 X-specific 和 Y-specific 成分是否与实际的生物学问题相关。如果结果符合预期或者揭示了有意义的生物学关系，那么选择的参数范围可以认为是合理的。

忽略数据的原始分布（例如，直接进行 Z-score 标准化，而不进行 log 转化 或其他特定的归一化处理）可能会对分析结果产生一定的影响。特别是在面对 偏态分布、稀疏数据 或 不同量纲的数据 时，直接使用 Z-score 标准化可能并不总是最优选择。下面详细解释忽略原始分布可能带来的影响：

1. 偏态分布（Skewed Distributions）

许多 代谢物数据、基因表达数据 和 蛋白质浓度数据 都表现出 偏态分布，即大部分数据值集中在一个范围内，而有少数极端值（例如，极低或极高的代谢物浓度或基因表达水平）。对于这种数据，直接进行 Z-score 标准化可能存在以下问题：

放大极端值的影响：Z-score 标准化会将所有数据按其均值和标准差转换，而这会使得那些本来较小的、稀有的极端值变得更加突出。对于偏态分布数据，极少数的高值或低值可能导致标准差较大，从而对 Z-score 标准化后的结果产生不成比例的影响。
偏态数据的转化需求：对于偏态分布数据，常见的做法是先进行 log 转化（或类似的平滑处理），使得数据的分布趋于对称，更适合进行后续的分析。如果忽略这种处理，直接进行 Z-score 标准化，可能会导致模型受到这些极端值的强烈影响，影响最终的分析结果。

例如，如果代谢数据中某些代谢物的浓度跨越多个数量级，而没有进行 log 转化，Z-score 标准化可能会导致这些代谢物的极高浓度占主导地位，而使得模型不能有效捕捉到整体数据的变异。

2. 稀疏数据（Sparsity）

很多基因表达数据和代谢物数据是稀疏的，尤其是基因组学和代谢组学数据中，可能有很多基因或代谢物在大部分样本中均未检测到或其表达值接近零。对于这种稀疏数据：

Z-score 标准化可能不适用：在稀疏数据中，很多变量（如基因或代谢物）的值为零或接近零，直接进行 Z-score 标准化会导致这些变量的标准差极小，可能会影响标准化后的数据表现。例如，在 代谢物丰度数据 中，很多代谢物可能大部分时间为零（如低丰度代谢物）。这种情况下，Z-score 标准化可能无法正确地反映这些变量的变异性。
中位数归一化或总和标准化：对于稀疏数据，通常会采用 中位数归一化 或 总和标准化（如相对丰度标准化）来处理数据，这样可以避免对稀疏数据的标准差过分敏感。

3. 影响结果解释和模型性能

如果忽略数据的原始分布直接进行 Z-score 标准化，可能会导致以下几方面的影响：

（a）模型稳定性与精度

如果数据存在偏态分布或高度偏斜，Z-score 标准化可能会影响模型的稳定性，因为这些数据会受到不成比例的放大或缩小。比如，如果一个代谢物在绝大部分样本中浓度极低，但某些样本的浓度非常高，Z-score 标准化后，这些高浓度样本对模型的贡献可能会过大，从而影响联合成分的提取。

（b）结果的可解释性

对于代谢组学和转录组学的数据，通常希望通过 PCA、PLS 或 O2PLS 这类方法揭示数据的潜在结构。如果没有考虑数据的原始分布，可能会导致从结果中提取的特征不符合数据的生物学特性。例如，偏态分布的数据如果不进行适当转化，可能导致模型误解一些变量的关系，从而影响最终的模型解释和生物学意义的提取。

（c）过拟合与欠拟合

过拟合：如果极端值对标准差影响过大，Z-score 标准化后模型可能会对这些异常数据点过于敏感，从而导致过拟合。
欠拟合：如果没有进行合适的数据变换（例如对于偏态分布的 log 转化），可能导致模型未能捕捉到数据中的真正结构，从而无法在验证集上表现良好，产生欠拟合。

4. 数据的自然尺度和解释

数据的原始分布反映了实际生物过程中的一些自然尺度。例如，在基因表达数据中，某些基因的表达水平可能是指数级增长的，直接进行 Z-score 标准化会忽略这些生物学上的尺度。如果你能够通过 log 转化将数据转化为更适合模型处理的形式，那么你会更好地保留这些生物学上的尺度。

5. 适当的预处理策略

为了避免这些潜在问题，在进行 O2PLS 分析之前，通常会采用以下预处理步骤：

对数转换（Log Transformation）：特别是当数据存在偏态分布时（如代谢数据），对数转换可以使数据分布更接近正态分布，从而减少偏态对模型的影响。
中位数归一化（Median Normalization）：对于基因表达数据或代谢物丰度数据，采用中位数归一化方法能够减少技术偏差，提高分析的鲁棒性。
去低变基因或代谢物：去除表达量变化非常小的基因或代谢物可以提高分析的稳定性，并减少噪声的影响。

忽略数据的原始分布，直接进行 Z-score 标准化可能会导致一些潜在问题，特别是在处理 偏态分布数据 或 稀疏数据 时。Z-score 标准化会放大或缩小某些极端值和噪声，影响最终的模型结果。
适当的预处理，如 log 转化 或 中位数归一化，通常能够更好地反映数据的原始结构，增强模型的解释性和稳定性。
O2PLS 分析前，建议仔细检查数据分布和噪声情况，选择合适的预处理方法，以获得更可靠的结果。

是的，在样本数非常少的情况下（如只有 6 个样本），设置过多的联合成分（joint components）、X-specific 成分和 Y-specific 成分会导致过拟合和计算不稳定的问题。通常来说，组件数量应该与样本数量相匹配，避免因过多的成分而增加模型的复杂度。

为什么设置为 5 个成分不合适？

过拟合的风险对噪声的影响反映到模型中造成失真

在小样本数据集上，通常选择的组件数（联合成分、X-specific 和 Y-specific）不应该超过样本数 - 1。这样可以避免模型复杂度过高，减少过拟合的可能性。

样本数和成分数之间的关系
通常建议，在 O2PLS 或 PLS 分析中， 成分数（components）应该远少于样本数。这是因为在小样本数据集上，成分数过多会使得模型在解释数据时过于复杂，容易在训练数据上找出与数据噪声匹配的模式，这些模式并不具有普适性。

理想情况下，联合成分、X-specific 成分和 Y-specific 成分的数量应该考虑以下几个方面：

联合成分（joint components）：通常建议联合成分数少于样本数，并且可以通过交叉验证来选择最优值。常见的选择范围可能是 1 至 3 个联合成分。
X-specific 和 Y-specific 成分：也不宜设置过多，尤其是在样本数小的情况下。通常来说，X-specific 和 Y-specific 成分的数量应与联合成分相似，或者稍少一些。可以设置为 1 到 2 个，而不是过多。
3. 计算不稳定性
在小样本和高维数据的情况下，计算不稳定性是另一个潜在问题。成分数过多可能导致协方差矩阵的估计不稳定，因为数据的维度（变量数）远高于样本数，导致计算时出现数值问题或无法收敛。因此，在这种情况下，选择合适的成分数非常重要，以确保模型在计算上能够稳定运行。
4. 样本数量对 O2PLS 的影响
O2PLS 是一种多变量统计方法，通过提取联合成分、X-specific 成分和 Y-specific 成分来分析 X 和 Y 之间的关系。在样本数非常少的情况下，通常不建议选择过多的成分，因为这会使得模型面临高维低样本数（HDLSS）问题，从而导致结果不稳定，难以提取出可靠的结构信息。
推荐的做法
减少成分数：对于只有 6 个样本的情况，建议将联合成分（joint components）和 X-specific、Y-specific 成分数限制在较小的范围内。一般来说，选择 1 到 3 个成分会更加合理。如果选择了 5 个成分，可能导致模型过于复杂，不易泛化。

通过交叉验证选择成分数：使用交叉验证（CV）来优化联合成分和特异性成分数的选择。交叉验证能够帮助找到最优的成分数，以避免过拟合。可以通过设置不同的成分范围（例如 1:3）并选择交叉验证得分最佳的组合。

增加样本数：如果可能的话，增加样本数是解决样本数小于变量数问题的根本方法。如果无法增加样本，建议谨慎选择模型的复杂度（即减少成分数）并优化数据预处理。

简化模型结构：在小样本数据集的情况下，简化模型结构（例如减少成分数、使用正则化等）能够帮助提高模型的稳定性，避免过拟合。

**B_T 的对角元素为零** 在 O2PLS 输出中的含义是，**在特定的联合成分上，代谢物（X）和基因表达（Y）数据之间没有显著的共享变异**。
	这个解释既可以表示 **没有共享变异**，也可以表示 **共享变异不足**，具体情况取决于数据的实际结构和模型设定。
 
### 1. **没有共享变异**
   - 如果在某些成分上，**B_T** 的对角元素为零，通常意味着在这些成分上，代谢物数据和基因表达数据之间 **几乎没有共同的变异**。
   - 换句话说，**X** 和 **Y** 的变异是完全独立的，代谢物和基因表达之间没有线性关系，或者说，选择的成分没有捕捉到这两者之间的任何共享信息。这时，模型不能有效地在这些成分上捕捉到它们之间的相关性。
 
### 2. **共享变异不足**
   - 另一方面，**B_T 的对角元素为零** 也可能表示，虽然 **X** 和 **Y** 之间存在某种共享的变异，但这些变异的 **解释能力不足**，即它们对整体变异的贡献非常小。 
   - 这种情况下，代谢物和基因表达数据之间确实有共享变异，但它们对目标变量的影响不足以显著区分或解释联合成分。
   - **共享变异不足** 通常会表现在联合成分上，虽然 **B_T** 有非零的元素（即有共享变异），但这些共享变异的幅度较小，可能无法有效捕捉 **X** 和 **Y** 之间的显著关系。
 
### 3. **模型的影响**
   - **B_T 对角元素为零** 可能会影响你对模型的解释和理解：
     - 如果没有共享变异（完全为零），则说明代谢物和基因表达之间几乎没有直接的关系，或者选择的成分未能有效捕捉两者之间的相关性。
     - 如果共享变异不足（接近零），可能意味着模型的表现不够好，需要调整模型的成分数目，或者改进数据的预处理（如去除低变基因、标准化等）来更好地捕捉两者之间的关系。
 
### 4. **如何判断是“没有共享变异”还是“共享变异不足”？**
   - 需要结合模型的 **其他输出**（如 **T**、**H**、**R²** 和 **MSE** 等）来进一步评估：
     - 如果模型的预测效果很好（即 **R²** 较高），但某些联合成分的 **B_T 对角元素为零**，这可能表示这些成分对模型贡献较小，或只影响 **X** 或 **Y** 中某一部分的变异。
     - 如果模型的预测效果较差（即 **MSE** 很高），那么可能是因为 **X** 和 **Y** 之间确实没有显著的共享变异，或者选取的成分数目不足。
 
### 5. **如何优化模型？**
   - **增加联合成分数量**：如果模型的 **B_T 对角元素为零**，可以尝试增加更多的联合成分。通过调整成分数目，模型可能能够捕捉到更多的共享变异。
   - **重新选择合适的成分数**：通过 **交叉验证** 或 **模型选择指标**（如 **AIC**、**BIC** 或 **R²**）来优化选择的联合成分数，确保模型能够最大程度地提取代谢物和基因表达之间的共享变异。
   - **数据预处理**：如果数据本身存在很大的噪音或不均衡，可以尝试进行更好的预处理，比如去除低变基因、数据标准化等，可能会帮助模型更好地捕捉两者之间的关系。
 
### 6. **实际场景中的应用**
   - 在实际的代谢物与基因表达分析中，代谢物和基因之间的共享变异可能并非是 **一对一** 或简单线性关系，因此一些 **B_T 对角元素为零** 是可以接受的，尤其是在复杂的生物学系统中，代谢物和基因的关系可能是间接的、非线性的或者通过其他因素调节的。
   - 例如，某些代谢物可能通过转录因子或其他中介因子间接影响基因表达，因此它们的直接共享变异可能不显著，但间接的影响仍然存在。
 
### 总结：
**B_T 的对角元素为零** 的意义通常可以解释为在这些成分上 **代谢物和基因表达之间没有显著的共享变异**，但它也可以表明 **共享变异不足**，即变异存在，但贡献较小。在 O2PLS 模型中，选择合适的联合成分数量是关键，过少的成分数可能无法捕捉所有的共享变异，而数据预处理和模型评估方法有助于提高模型的解释力和预测能力。

 哪个结果更可信？
标准化 (scale = T) 的结果通常更可信，因为它消除了不同变量尺度（例如代谢物和基因表达的量纲差异）对模型的影响。标准化后的数据让所有变量的方差都变为 1，避免了某些变量因变异过大或过小而在模型中占主导地位，这有助于捕捉到数据中的潜在共享变异。
不标准化 (scale = F) 的结果可能不够可靠，特别是当 X 和 Y 数据中包含具有不同尺度或方差的变量时，这会导致模型未能有效捕捉到共享变异，从而可能出现 B_T 的对角元素为零 的情况。
4. 为什么 scale = T 更有意义？
在大多数的 多元数据分析方法（例如 PCA、PLS、O2PLS）中，标准化是一个常见的预处理步骤，因为数据集中的变量通常具有不同的量纲和方差。标准化能确保每个变量对模型的贡献是平等的，并且更好地捕捉到变量间的相对关系。
特别是在代谢物数据和基因表达数据中，数据的量纲（如基因表达值、代谢物浓度等）和方差可能差异较大，未经标准化的数据可能导致某些变量的影响过大，甚至掩盖了代谢物和基因表达之间的潜在关系。
5. 总结

指标	全称	含义
R²_X	Explained variance in X	模型（包括联合 + 正交部分）对 X 的总解释能力
R²_Xjoint	Explained variance in X-joint	模型中联合成分对 X 的解释能力（X 中和 Y 有关的那部分）
R²_Xhat	Reconstructed variance in X	用所有模型成分（联合 + 正交）重构 X 的能力
R²_Xpred	Predicted variance in X	用 Y 中的联合成分预测 X 中联合部分的能力（可预测的部分）

        R²_X（X 的总解释度）
        ├── R²_Xjoint（X 中的联合部分：与 Y 有关）
        │     └── R²_Xpred（可由 Y 联合成分预测的 X 联合成分）
        └── R²_Xorthogonal（X 中的特异部分，仅模型解释，Y 无关）

也可以理解为：

R²_X = R²_Xjoint + R²_Xorthogonal（正交部分）

R²_Xhat ≈ R²_X（模型对 X 的重构能力，理论上两者接近）

R²_Xpred ≤ R²_Xjoint（Y 可预测的那部分）

R²_X 是最泛的，表示模型对整个 X 的解释能力，包括：
- 联合成分解释的部分（与 Y 有共享结构）
- 正交成分解释的部分（X 特有的变异）
R²_Xjoint 表示模型在 X 中提取出的、与 Y 共享的部分。
- 这部分说明了 X 和 Y 有多少信息是互相关联的。
R²_Xpred 是 Y 的联合成分对 X 联合部分的预测能力。
- 代表了 X 能够被 Y 预测的部分，是交叉预测能力的重要指标。
- 越高，说明 X 和 Y 的耦合越强。
R²_Xhat 通常是模型重构后的 X 与原始 X 的相关性。
- 如果模型拟合良好，R²_Xhat 应该接近 R²_X。

✅ 举个简单例子

假设你得到以下结果：

指标	数值
R²_X	0.85
R²_Xjoint	0.30
R²_Xpred	0.12
R²_Xhat	0.83

解释：

模型可以解释 X 中 85% 的变异（R²_X）
其中只有 30% 是和 Y 有关的（R²_Xjoint）
Y 只能预测 X 的 12%（R²_Xpred），说明虽然 X 和 Y 有一定共享信息，但预测能力并不强
R²_Xhat ≈ R²_X，模型的拟合还不错

✅ 总结

指标	指代	作用
R²_X	总解释度	模型对 X 的整体拟合能力
R²_Xjoint	联合解释度	X 中与 Y 相关联的部分
R²_Xpred	可预测度	Y 能预测 X 的那一部分
R²_Xhat	重构度	用模型重构 X 的拟合程度，通常 ≈ R²_X

缩写	全称	含义
SSQ_X	Sum of Squares of X	X 数据集的总变异（平方和）
SSQ_Y	Sum of Squares of Y	Y 数据集的总变异（平方和）
VAR_Xjoint	Variance in X - joint	X 中由联合成分（与 Y 有共享信息）解释的变异量
VAR_Xorth	Variance in X - orthogonal	X 中由正交成分（X 特有信息）解释的变异量

1. SSQ_X / SSQ_Y

是原始数据矩阵 X、Y 中各列平方和的总和，即 数据中所有变量的总变异量。
通常作为归一化或计算 R² 的分母。

📌 SSQ_X = VAR_Xjoint + VAR_Xorth + 残差（未能解释的部分）

2. VAR_Xjoint（联合变异）

这部分表示 O2PLS 模型提取出的 X 中与 Y 有共同变化趋势的那部分信息。
对应于联合成分 T 的贡献。

3. VAR_Xorth（正交变异）

表示 X 数据中只在 X 自身内部变化、但与 Y 无关的那部分。
对应于正交成分 T_X。

✅ 三、一个类比（形象点）

假设 X 是转录组数据，Y 是代谢组数据：

SSQ_X 就是你所有基因表达量的总波动强度。
VAR_Xjoint 是那些和代谢水平明显“联动”的基因表达变化（比如某个代谢物升高时某些基因总是被激活）。
VAR_Xorth 是那些只在转录组自己内部变化，但和代谢没什么关系的表达波动（比如一些 housekeeping gene）。

✅ 四、与 R² 的关系

你可以把它们理解为 R² 的“源头”：

$R_{X}^{2} = \frac{V A R _{X j o in t} + V A R _{X or t h}}{S S Q _{X}}$

如果只考虑联合部分的解释能力：
$R_{X j o in t}^{2} = \frac{V A R _{X j o in t}}{S S Q _{X}}$

同理适用于 Y。

✅ 五、实际应用中怎么看这些数值？

假设你看到以下模型输出：

指标	数值
SSQ_X	100
VAR_Xjoint	20
VAR_Xorth	65

你可以得出：

联合成分解释了 X 总变异的 20%
正交成分解释了 X 的 65%
剩下 15 是误差项（模型没解释到的）

这说明大部分变异是 X 自身的结构，不是与 Y 共享的 → X 和 Y 的关联不强。

这段 R 代码的目的是在一个转录组（或其他组学）数据集中 筛选出变异性和表达量都比较高的基因，以便后续分析（比如 O2PLS、聚类、可视化等）。下面是逐步解释：

✅ 数据背景假设

dat.omics2 是一个 样本 × 基因 的矩阵或数据框，行为样本，列为基因表达值。
你希望筛选出在所有样本中：
- 表达量高（最大表达量高）
- 表达变异大（IQR 大）

2. X first joint loading 代表什么

O2PLS 的 first joint loading（例如 X 端的 $p_{1}$ ）表示 X 变量在第一联合成分方向上的贡献，而这个联合成分是通过最大化与 Y 联合得分的协方差来求的。
所以：
- 它不是单纯的“X 自身方差贡献”
- 它是“X 变量对 X 中那部分能解释 Y 的成分”的贡献
如果一个转录本在 $p_{1}$ 中系数很大，这意味着它在与 Y 高度相关的方向上贡献大，而不是单纯对 X 自身重要。
同理，Y 端的 $q_{1}$ 表示 Y 变量在与 X 最相关方向上的贡献。

某个转录的 loading 大，说明对自生的影响大还是对 Y 的代谢物的影响大？

在 O2PLS 的 first joint loading 中，这个“大”是指它在解释 X 那部分与 Y 相关的信号时的重要性，所以更贴近“对 Y 的代谢物影响大”的含义，而不是 X 自身总方差的贡献。

lvzhaopeng的数字花园

探索

生物 O2PLS官方文档