Contrastive learning leverages the idea of learning embeddings such that similar data points are brought closer in the embedding space, while dissimilar points are pushed apart. Therefore, how to find these similar and dissimilar data pairs are critical in all types of contrastive learning methods.

In self-supervised contrastive learning, we define positive (similar) and negative (dissimilar) data pairs based on whether the data is from a paired sample collected simultaneously from the same subject. For a data pair {xi,yj};i, j∈{1,2, …, n} where n is the batch size.

In supervised contrastive learning for categorical response variable, positive and negative data pairs are determined based on whether the data has the same response value. Given a data pair {xi,yj} and their responses {ri,rj} where i, j∈{1,2, …, n}, the positive (or negative) pair is defined as

However, when the response variable is continuous in a supervised contrastive learning model, the traditional definition of positive and negative data pairs becomes inapplicable, as continuous response values are rarely identical. Instead, we categorize data into similar and dissimilar pairs based on the distance between their response values. Given data points {xi,yj} with continuous responses {ri,rj}, we define the distance as dij=|ri−rj|, which serves as the measure of similarity. A distance close to 0 indicates a highly similar data pair, while a larger value approaching its maximum suggests greater dissimilarity. The distance d_ij is then used in the generalized contrastive loss.

The general goal of contrastive learning is to maximize the similarity of encoded embeddings in the representation domain. To measure the similarity of two vectors, the cosine similarity is the most frequent choice. From a mini-batch point of view, suppose x_i and y_j are different omics data vectors obtained from sample i and sample j ( i,j∈{1,2,…,n}; n is the batch size), respectively, f_x and f_y are the encoders for each omics data, then the embeddings for x_i and y_j are zix=fx(xi) and zjy=fy(yj). The similarity of embeddings is defined as

where · denotes the dot product and ||·|| is the Euclidean norm of a vector.

By anchoring the embedding of one omics data z^x, we define the generalized contrastive loss to be

where wij∈[0,1] is the weight for sample pair {xi, yj}, and

is the SoftMax probability of cosine similarity of scaled embeddings and scaler τ∈ℝ+ is the temperature parameter.

By fixing the other embedding z^y as the anchor, we will get the other piece of loss Loss^y,x. For two omics data, the total contrastive loss is

Rewrite this formula in matrix form:

where ∘ denotes the Hadamard product (element-wise product), log P is applied elementwise, and (·)°(−1) denotes the Hadamard inverse. The Hadamard inverse of matrix X is defined as [X°(−1)](i,j)=1xij. Also, vector jn×1=[1, 1,…,1]T, matrix W_n×n is the weight matrix and matrix P_n×n is the SoftMax probability matrix, with P(i,j)=Pij. Here, (Wj) holds the row sums of the weights, so its Hadamard inverse normalizes each anchor by its total weight.

By this definition, we can use one generalized formula to incorporate self-supervised contrastive loss, supervised contrastive loss on categorical responses and supervised contrastive loss on continuous responses. The key lies in the definition of the weight matrix W, which captures the information of both similar and dissimilar pairs.

First, for self-supervised contrastive loss, we have Wn×n=In×n where I is an identity matrix.

Second, for supervised contrastive loss on categorical responses, W is defined as

where r is the response vector and equality function equal(·) is a true/false function that checks the element-wise equalities. Equivalently,

Third, for supervised contrastive loss on continuous responses, we propose three methods to define the weight matrix W. The underlying principle is to assign a weight to a data pair based on the distance of their responses, with the weight being negatively correlated to the response distance. Namely, the larger the distance is, the smaller the weight is and vice versa.

Suppose r_i, r_j are the responses for a data pair {xi, yj}, the norm-1 distance is dij=|ri−rj|.

Accordingly, the three weights are defined as follows.

(1) Linear weights

where the scale

and a small positive constant ϵ is to avoid zero denominator.

(2) Exponential weights

(3) Negative-log weights

where ϵ is a small positive constant to avoid zero values and we further scale wij′ to contain the weights within [0,1].

In the rest of the article, we will use all three methods to train the model and compare their performance. In general, we do not expect any single method to be universally optimal for all datasets. The best method depends on the specific distribution of the response variables and the characteristics of the data. Instead of leaving this choice open, we provide a simple validation-based rule to select the weighting method and show when each method is most robust (see Results).

The MB-SupCon-cont framework (Figure 1) is designed to integrate multi-omics data and predict both categorical and continuous responses with enhanced accuracy and interpretability. The first stage of the MB-SupCon-cont framework involves data collection. Paired omics data, such as microbiome and metabolome data, are collected along with patients’ phenotype information. In the next stage, contrastive learning is applied by training a supervised contrastive learning model to obtain encoders f_x and f_y. In this context, the generalized contrastive loss is used to accommodate different types of response data. In the final stage, prediction heads (either classifiers or regressors) are applied to the embeddings. When making new predictions, embeddings from either microbiome data or other omics data are selected based on their performance during the training phase. Additionally, unique trends related to the responses can be visualized in the lower-dimensional space.

Each omics data type has its own encoder, and the two encoders do not share weights, so each encoder can match the size and structure of its own input. Each encoder is a multilayer perceptron composed of fully connected blocks. Every block applies a linear layer, one-dimensional batch normalization, a ReLU6 activation, and dropout, followed by a final linear layer that outputs the embedding of dimension d_e. For the simulation study we varied the embedding dimension over de∈{10,20,40} to assess its effect, and we report d_e = 20 as the representative setting (Figure 2), with d_e = 10 and 40 in Supplementary Figures 1, 2. For the two real datasets we used a compact embedding of d_e = 10, at the low end of this range, which suits their smaller sample sizes. The hidden-layer widths of each encoder are listed in the Supplementary Material.

The two encoders are trained together by minimizing the generalized contrastive loss, using stochastic gradient descent with a learning rate of 0.001 and a momentum of 0.9, gradient clipping, a batch size of 32, and 1,000 epochs. Regularization is provided by dropout, batch normalization, and weight decay. For each dataset, response, and weighting scheme, the dropout rate, weight decay, and temperature were selected by the lowest validation RMSE over a grid search that used only the training and validation data. We repeated this selection across seeds and checked both how sensitive the validation RMSE was to these settings and how stable the selected values were from seed to seed. The search grid and its selected values, and the sensitivity and stability results are reported in the Supplementary Material.

We compared MB-SupCon-cont against three groups of baselines, all evaluated on the same splits and with the same regression heads. For the classical baselines, we fit (1) elastic net, (2) support vector regression (SVR), (3) random forest (RF), and (4) extreme gradient boosting (XGBoost) on the original omics data. Each baseline was tuned by 5-fold cross-validation on the combined training and validation data, with the testing data held out. The hyperparameter search grids for these four baselines are given in the Supplementary Material. As a multi-view baseline, we fit (5) canonical correlation analysis on the training data with the same number of components as the embedding dimension and passed the canonical components of each omics data type to the regression heads. To separate the contribution of the contrastive loss from that of the neural network, we trained (6) an autoencoder for each omics data type, using the same encoder architecture but a reconstruction loss in place of the contrastive loss, and used the encoder outputs as embeddings. Its dropout and weight decay were selected on its own validation reconstruction loss, so that this baseline was tuned for the objective it optimizes.

For statistical testing, we paired the 12 test RMSE values by seed for each comparison, applied a paired Wilcoxon signed-rank test, reported a 95% confidence interval for the mean difference from a paired bootstrap with 10,000 resamples, and controlled the false discovery rate (FDR) within each table using the Benjamini-Hochberg procedure (Benjamini and Hochberg, 1995).

To benchmark the efficacy of MB-SupCon-cont, we designed a simulation study. The simulation data was generated by Reverse Principal Component Analysis (Reverse PCA), a technique we developed and detailed in the Supplementary Material. It reconstructed the original omics data from its lower-dimensional representation. The choice of Reverse PCA was grounded in its simplicity and ease of implementation, making it an accessible yet robust tool for our analysis. It is assumed that for two sets of omics data, a representation space exists wherein the corresponding features from both datasets are correlated. Reverse PCA is used to reconstruct a reduced-rank approximation of the original data based on their principal components. This approach guaranteed that the vital information contained within the embeddings is preserved through the reconstruction process. In this study, we assume the generated embeddings in the representation space are corresponding to the highest-variance principal components and the response variable is linked to these embeddings.

Algorithm 1 in the Supplementary Material gives the full implementation details of Reverse PCA. Briefly, the procedure generates two correlated omics datasets from latent embeddings and constructs a nonlinear response from linear, quadratic, and cross-omics interaction effects, and use the simulated data to evaluate MB-SupCon-cont. For each simulation setting, we generate 1,000 paired samples and set the dimensions of the two simulated omics data types to 100 and 500 features. The effect coefficients linking the latent embedding features to the response were generated from a normal distribution N(0,25), and the generated response is scaled to 0 to 100 assuming it contains “Age” information. We also considered several values for the average feature correlation, μρ∈{0.4, 0.6, 0.8}, and the embedding dimension, de∈{10,20,40}.

To test whether the advantage of MB-SupCon-cont holds beyond the Reverse PCA design, we added three more realistic regimes to the simulation. In the first regime, only a subset of the representation carries signal, and the remaining components are uninformative noise. In the second regime, only one omics layer is informative. In the third regime, the response coefficients are sparse. As a further check on the realism of the design, we validated Reverse PCA on the two real datasets by comparing the original data with their low-rank reconstruction and measuring the fraction of feature variance that the reconstruction recovered. We describe the regimes and the reverse PCA validation in full in the Supplementary Material.

To showcase the efficacy and benefits of MB-SupCon-cont, the prediction on simulated original data is compared to the prediction on learned embeddings. Our simulation demonstrated the superiority of MB-SupCon-cont in predicting continuous responses. Specifically, as the average correlation of embedding features rises, the advantage of MB-SupCon-cont grows correspondingly. At μρ=0.8, the benefit of using MB-SupCon-cont is most notable. As shown in Figure 2, in most scenarios the blue, orange and green boxes (prediction on embeddings) are significantly lower than the red box (prediction on original data). Additionally, the linear weighting method performs marginally better than the other two when predicting the simulated “Age” variable. This illustrates that MB-SupCon-cont can enhance prediction performance based on embeddings, and this enhancement increases with the average correlation coefficient of embeddings (Figure 2; Supplementary Figures 1, 2).

We repeated the simulation under the three more realistic regimes described in the Methods, which add uninformative noise components, make only one omics layer informative, or make the response coefficients sparse. In the single informative omics and sparse regimes, the embeddings still gave lower prediction error than the original data, and when only one omics layer carried signal, the informative layer produced the better predictions, as expected. In the noise regime, where most of the representation is uninformative, the advantage of the embeddings shrank to roughly a tie at the best head, as expected when there is little shared structure to exploit. This shows that the advantage of MB-SupCon-cont holds whenever the two omics layers share response-related signal and does not depend on the specific structure of the Reverse PCA design. The full results for all regimes are given in the Supplementary Material (Supplementary Figure 6; Supplementary Table 4).

We also assessed how well Reverse PCA reproduces the real data. At an embedding dimension of 20, the low-rank reconstruction recovered about 79% of the feature variance for the T2D microbiome, 49% for the T2D metabolome, 34% for the HFD microbiome, and 46% for the HFD metabolome, which shows that Reverse PCA is a reasonable generative model for real omics data. The features with the poorest reconstruction were likely to be low-abundance taxa and metabolites, which may be harder for principal component analysis (PCA) to capture accurately.

We benchmarked MB-SupCon-cont to paired microbiome and metabolome data (Zhou et al., 2019). This study enrolled 720 samples of prediabetes patients. Three continuous responses that are all associated with the disease status of diabetes patients were considered for prediction, which are Age, Body Mass Index (BMI) and Steady State Plasma Glucose (SSPG) level.

To assess the prediction performance, we employed the Root Mean Square Error (RMSE) as the evaluation metric. For a comparative perspective, predictions were also made on the original microbiome and metabolome data. The MB-SupCon-cont results used the best hyperparameter combination selected on the validation data, and the baselines were also tuned, as described in the Methods. The dataset was partitioned into a training set (70%), a validation set (15%), and a testing set (15%). The dimensions of the microbiome data are 720 × 96, and for the metabolome data, they are 720 × 724. We executed the analysis 12 times, each with a unique random seed for the training-validation-testing splits and recorded the RMSE for all prediction methods during each run.

In Table 1, the mean RMSE values averaged over the 12 runs are contrasted across the various prediction methods. For predictions based on MB-SupCon-cont embeddings, each weighting method is shown with its best performing prediction head, and the full per-head results are given in the Supplementary Material. Predictions based on metabolome embeddings achieved the lowest RMSE for all three responses. The best average test RMSE for Age was achieved with linear weights (4.2379), for BMI with negative-log weights (1.5221), and for SSPG with negative-log weights (20.0849). Each of these was significantly more accurate than the best tuned baseline (paired Wilcoxon signed-rank test, p = 0.0005 for Age, p = 0.0093 for BMI, and p = 0.0010 for SSPG; FDR-adjusted p-value of at most 0.014). The canonical correlation analysis and autoencoder baselines were clearly worse (FDR-adjusted p-value of at most 0.004). The autoencoder uses the same encoder architecture, was tuned on its own reconstruction objective, and differs from MB-SupCon-cont only in the loss function, so this shows that the improvement comes from the contrastive learning formulation rather than from the neural network architecture alone. For the microbiome data, the tuned classical baselines had lower average RMSE than MB-SupCon-cont for all three responses (with only BMI significant; FDR-adjusted p = 0.023), which indicates that in this study the metabolome carries the stronger integrative signal. A boxplot comparison of all prediction methods for predicting three responses is also presented in Figure 3.

For comparison, PCA was performed on both (1) the original microbiome and metabolome data, and (2) the MB-SupCon-cont microbiome and metabolome embeddings, each shown at the weighting method that performed best for that omics and response. From Figure 4, both Panel A (metabolome) and Panel B (microbiome) show that the MB-SupCon-cont embeddings (second rows) exhibit a much more distinctive varying trend in the lower-dimensional space compared to the original metabolome and microbiome data (first rows). This indicates that MB-SupCon-cont embeddings provide more accurate response-related representations of the omics data and enhance the visualization of underlying data patterns. The first two principal components of the metabolome embeddings explain a large proportion of the variance (about 0.83 for Age, 0.59 for BMI, and 0.64 for SSPG), compared with about 0.14 to 0.16 for the original metabolome data.

To further validate the superiority of MB-SupCon-cont in predicting continuous responses, we applied it to another independent multi-omics murine model studying the effect of a high-fat diet on the microbial and metabolic changes in mice (Morton et al., 2019). This model utilized 434 paired microbiome and metabolome data for integration. In contrast to the previous T2D study, the HFD study presented substantially higher dimensionality in the paired omics data, with the microbiome data dimension being 434 × 913 and the metabolome data dimension being 434 × 11,978. We still split the dataset into training (70%), validation (15%), and testing data (15%). The analysis was conducted 12 times using different random seeds. As in the T2D study, the MB-SupCon-cont results used the best hyperparameter combination selected on the validation data, and the baselines were also tuned. Responses Age and Weight are selected as the prediction tasks.

Supplementary Table 1 and Supplementary Figure 3 have shown the comparison of prediction RMSEs on testing data (15% of all data) across different methods. Similar to the T2D study, predictions based on metabolome embeddings achieved the lowest average RMSE for both Age and Weight. The best average test RMSE for Age was achieved with negative-log weights (7.0759) and for Weight with linear weights (6.3413). For Age, this was significantly lower than the best tuned baseline (paired Wilcoxon signed-rank test, p = 0.0005; FDR-adjusted p-value is 0.0010), while for Weight the improvement was small and not statistically significant (p = 0.3804). The canonical correlation analysis and autoencoder baselines were again clearly worse for both responses (FDR-adjusted p-value of at most 0.0010). Meanwhile, the color variations corresponding to different response values are clearly distinguishable in the lower-dimensional PCA scatter plots (PC2 vs. PC1) of metabolome and microbiome embeddings (Supplementary Figures 4, 5). The first two principal components of these embeddings explain a large share of their variance (about 0.98 for Weight and 0.52 for Age in the metabolome embeddings, and about 0.97 for Weight and 1.00 for Age in the microbiome embeddings), compared with about 0.16 for the original metabolome data and 0.09 for the original microbiome data. This study reaffirms the stable and consistent performance of MB-SupCon-cont.

We defined three weighting schemes, which raises two questions: how to choose among them, and whether all schemes are worth keeping. Our results show that the best scheme depends on the response and can even differ between the two omics. In the T2D study, the lowest error on the metabolome came from linear weights for Age and negative-log weights for BMI and SSPG, while on the microbiome it came from linear weights for Age and SSPG and exponential weights for BMI. In the HFD study, the lowest error on the metabolome came from negative-log weights for Age and linear weights for Weight, while on the microbiome linear weights were best for both Age and Weight (Table 1 and Supplementary Table 1). Thus, the best scheme depends on the response and should be selected from the data rather than fixed in advance.

To make this choice automatic, we select the scheme by validation performance. For each dataset and response, we take the three schemes at their best prediction heads and keep the one with the lowest validation RMSE. We then evaluate the selected scheme on the held-out test data. This rule selected the best performing scheme on the test set in nine of the ten (90%) scenarios across the two studies. In the one exception, HFD Weight, the scheme it chose was only about 0.2 RMSE behind the best. The validation-based rule is therefore a simple and dependable way to make the choice.

The schemes also differ in how they tolerate measurement noise in the response. To assess their robustness, we added Gaussian noise to the training and validation responses at levels ranging from 0 up to 1 standard deviation (SD) of the response (a noise fraction of 0 to 1), while keeping the test response clean. As shown in Figure 5, the three schemes degraded at different rates as noise increased, and thus their relative ranking could change depending on response ranges. For responses with a narrow range, such as T2D Age and BMI, the ranking of the three schemes remained largely unchanged as the noise increased. For responses with a wide range, such as T2D SSPG and HFD Weight, the ranking changed under heavy noise. Also, the exponential and negative-log weighting schemes became the most robust under heavier noise because they placed more weight on pairs with similar responses. For example, on the metabolome embeddings with the prediction head fixed at random forest, we varied the weighting scheme and the noise level (Figure 5). Measured by the increase in mean test RMSE as the response noise rose from 0 to 1 SD, the exponential scheme was the most robust for T2D SSPG, with its error rising by only 16.4 compared with 33.5 for linear weights. For HFD Weight, the negative-log scheme was the most robust, with its error increasing 1.7 compared with 3.9 for linear weights. In both cases, at the highest noise level, the more robust scheme also had the lowest validation RMSE and was therefore selected by the validation-based rule.

These results lead to a simple guideline. The weighting scheme is not a fixed part of the method but a data-driven choice that the validation-based rule can make reliably. The exponential and negative-log schemes earn their place because they are more robust than linear weights when a wide-range response may be measured with error. The three schemes therefore form a useful family rather than redundant options.

The principal component scatter plots suggest that the MB-SupCon-cont embeddings organize samples by the response. To test this directly, we measured how closely the geometry of the embedding space matches the geometry of the response. For each embedding we computed the pairwise Euclidean distances between samples, and did the same for the response, then measured the correlation between the two distance matrices with a Mantel test (Mantel, 1967) using 999 permutations. A high Mantel correlation means that samples close in the embedding space also have similar responses.

The Mantel test quantitatively confirms that distances in the best-scheme embeddings track the response distances much more closely than distances in the original data do. In the T2D study the Mantel correlations were about 0.27 to 0.82 across the two omics and three responses, and in the HFD study about 0.54 to 0.83 across the two omics and both responses (all significant with p = 0.001), compared with about 0.00 to 0.21 for PCA of the original data. For the autoencoder, although tuned on its own reconstruction objective, it reached only about 0.01 to 0.20 and was significant in only four of the ten cases. This confirms that the MB-SupCon-cont embeddings encode response structure that neither linear nor deep unsupervised representation recovers (Figure 6; Supplementary Table 11).

We also measured the effective dimensionality of each embedding using the participation ratio of its singular value spectrum, where a low value means the embedding is concentrated in a few directions and a high value means it is more spread out. The three weighting schemes produced clearly different geometries. The linear scheme gave the most concentrated embeddings, with a participation ratio around 1 to 2, while the exponential and negative-log schemes gave more spread-out embeddings, with a participation ratio up to about 7. The linear scheme was strongest on most metabolome cells (Mantel correlation about 0.82 to 0.92 on the T2D metabolome), while the negative-log scheme was strongest on some microbiome cells, such as the T2D microbiome for Age and BMI (about 0.62 and 0.61, compared with about 0.27 and 0.34 for linear weights). These results show that the weighting schemes produce genuinely different embedding geometries, not just minor variations, which gives a geometric explanation for why different schemes can be preferred on different data.

As MB-SupCon-cont embeddings are supervised by the response, we next examined which microbes and metabolites they encode. Features were ranked by their association with the response using both the embedding and the raw data, and the resulting two rankings were compared. For the first ranking, we define the response axis as the direction in the embedding space along which the response changes fastest, obtained by regressing the response on the embedding and keeping the fitted values. The regression is fitted on the training samples and applied to the testing samples, so the axis is never evaluated on the samples that defined it. Every feature is then ranked by its Spearman correlation with this axis. The second ranking comes from the raw data, where every feature is ranked by its Spearman correlation with the response itself. Both rankings use the testing samples only and are averaged over the 12 splits. Table 2 reports the top five microbes along each response axis, with both correlations and the number of splits in which each feature recurs. The full ranked lists for both microbes and metabolites are given in Supplementary Table 12.

The two rankings largely agree. The response axis follows the held-out response with a correlation of 0.52 to 0.93 across the studies, omics data types, and responses, and the two rankings share 9 to 17 of their top 20 features, against 0.03 to 4.2 expected by chance. Three controls were passed through the same pipeline on the same samples: the same encoder retrained on permuted responses (shuffled-response encoder), an autoencoder, and principal component analysis of the raw data. The shuffled-response encoder reached a correlation of at most 0.34 and shared at most 12 of the top 20 features, and the two unsupervised representations reached at most 0.62 and 14 of 20. MB-SupCon-cont was ahead of all three in every one of the ten cases (Figure 7; Supplementary Table 13). The shuffled-response encoder is the informative comparison, because it differs from MB-SupCon-cont only in that its supervision has been destroyed. The biological signal captured by the embedding is driven by the response and not from the omics data alone. The leading features also recur across splits, and the strongest of them entered the top 20 in all 12. For these features, the correlation with the axis exceeded the direct correlation in every case, which indicates that the embedding retains the response signal they carry rather than diluting it.

The leading features are also consistent with known biology. In the T2D study, the metabolite leading the Age axis was dehydroepiandrosterone sulfate (DHEA-S), matching its well-established decline over the human lifespan (Orentreich et al., 1984). The BMI axis was led by Odoribacter, a butyrate-producing genus, which decreased as BMI increased (Gomez-Arango et al., 2016). Along the SSPG axis, a measure of insulin resistance, unclassified Clostridiales taxa decreased, which agrees with a report from the same prediabetes cohort that linked a decline in the order Clostridiales to lower insulin sensitivity (Zhou et al., 2020). In the HFD study, members of the family Muribaculaceae (S24-7) decreased along the Weight axis while Erysipelotrichaceae increased, in agreement with previous reports on diet-induced changes in the mouse gut microbiome (Daniel et al., 2014; Fleissner et al., 2010). On the metabolome side of the same axis most of the leading features were never identified, and the highest ranked compounds that do carry an identification were the bile acids glycodeoxycholic acid and chenodeoxycholic acid, both of which increased with Weight, which is consistent with the evidence that high-fat diets can alter host bile acid composition (Devkota et al., 2012).

These results show that MB-SupCon-cont embeddings are built from biologically meaningful features and can still be read in terms of individual microbes and metabolites. As the encoder is supervised by the response, they recover the known response biology. The embeddings are therefore interpretable, and they can be read as a summary of the microbes and metabolites that carry the response signal.