Lehtiö et al., 2021: Multimodal data integration

In modern -omics experiments, it is common to measure multiple modalities at once. It is a central challenge to integrate these complementary measurements to a unified view of the biological state of the specimen. alphapepttools seamlessly integrates MS-proteomics data with other modalities via the shared scverse infrastructure.

Study background

In this tutorial, we explore the integration of multiple modalities in a clinical setting by reanalysing the dataset by Lehtiö et al., 2021. In the original publication, the authors recorded a deep proteome of non-small cell lung cancer (NSCLC) samples from 141 patients with fractionated DDA-TMT proteomics. They identified 6 proteomic cancer subtypes and associated them with histological subtypes, immune infiltration, and clinical outcomes. Notably, they compared the proteomic subtypes post hoc with matched measurements of RNA microarray + methylation experiments from a previous publication. Here, we demonstrate how these matched -omics measurements can instead be jointly modelled in a MOFA+ analysis.

References

1. Proteomics Data: Lehtiö, J. et al. Proteogenomics of non-small cell lung cancer reveals molecular subtypes associated with specific therapeutic targets and immune-evasion mechanisms. Nat Cancer 2, 1224–1242 (2021). PXD020191

2. RNA microarray data Karlsson, A. et al. Gene Expression Profiling of Large Cell Lung Cancer Links Transcriptional Phenotypes to the New Histological WHO 2015 Classification. J Thorac Oncol 12, 1257–1267 (2017).

3. Methylation data Karlsson, A. et al. Genome-wide DNA Methylation Analysis of Lung Carcinoma Reveals One Neuroendocrine and Four Adenocarcinoma Epitypes Associated with Patient Outcome. Clin Cancer Res 20, 6127–6140 (2014).

1. Proteomics Analysis

To run this tutorial, you need to install the scverse packages decoupler, muon, mofapy2, and mofax in your local environment

# ! pip install muon mofapy2 mofax decoupler

import anndata as ad
import alphapepttools as apt
import pandas as pd
from alphabase.pg_reader.keys import PGCols
import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import silhouette_score
import seaborn as sns
import scanpy as sc

import matplotlib.patches as mpatches

# Multimodal integration
import mudata as md
import muon

# Interpretation of data
import decoupler as dc

1.1. Helper functions

We define a few helper functions that we will reuse in the rest of the tutorial:

def mvalue(adata: ad.AnnData, *, alpha: float = 0.01, layer: str | None = None, copy: bool = False) -> ad.AnnData:
    r"""Conversion of methylation beta values to M values

    .. math::

        M = \log_{2}{\frac{\beta + \alpha}{(1-\beta) + \alpha}}


    Parameters
    ----------
    adata
        AnnData object
    alpha
        Pseudocount added to the data to prevent infinite values
    layer
        Layer to act on. If `None`, uses adata.X
    copy
        If `True`, returns a copy of the anndata object, else modifies the object inplace

    References
    ----------
    - Du, P., Zhang, X., Huang, CC. et al.
    Comparison of Beta-value and M-value methods for quantifying methylation levels by microarray analysis.
    BMC Bioinformatics 11, 587 (2010). https://doi.org/10.1186/1471-2105-11-587
    - Ricard Argelaguet (GitHub): https://github.com/bioFAM/MOFA/issues/44

    """
    adata = adata.copy() if copy else adata
    data = adata.X if layer is None else adata.layers[layer]

    data = np.log2((data + alpha) / ((1 - data) + alpha))

    if layer is None:
        adata.X = data
    else:
        adata.layers[layer] = data

    return adata if copy else None

# Shared color scheme for the omics modalities, reused across loading and enrichment plots
palette = {"proteins": "#2b7bba", "rna": "#705eaa", "methylation": "#2c944c"}


def plot_modality_loadings(
    data: pd.DataFrame,
    *,
    score: str,
    feature: str,
    palette: dict,
    modality_column: str = "modality",
    xlabel: str = "Score",
    ax: plt.Axes | None = None,
) -> plt.Axes:
    """Horizontal bar plot of per-feature scores, colored by omics modality.

    Used to display either MOFA+ factor loadings or gene set enrichment scores. Each
    bar corresponds to a single feature (a gene or a pathway) and is colored by the
    modality it originates from.

    Parameters
    ----------
    data
        Long-format dataframe with one row per feature.
    score
        Column holding the numeric value drawn as bar length (x-axis).
    feature
        Column holding the feature labels shown on the y-axis.
    palette
        Mapping of modality name to color.
    modality_column
        Column holding the modality of each feature.
    xlabel
        Label of the x-axis.
    ax
        Axes to draw on. If `None`, a new figure and axes are created.

    Notes
    -----
    Feature labels that occur more than once in `data` (e.g. a pathway or gene shared
    across modalities) are rendered in bold to highlight cross-modal agreement.
    """
    data = data.sort_values(by=score)
    colors = data[modality_column].map(palette)

    if ax is None:
        _, ax = plt.subplots(figsize=(4, 8))

    ax.barh(range(len(data)), data[score], color=colors)
    ax.set_yticks(range(len(data)))
    labels = ax.set_yticklabels(data[feature], fontsize=8)

    # Highlight features that occur in more than one row (e.g. shared across modalities)
    feature_counts = data[feature].value_counts()
    for label, name in zip(labels, data[feature], strict=True):
        if feature_counts[name] > 1:
            label.set_fontweight("bold")
    ax.axvline(0, color="black", linewidth=0.8)  # zero reference line
    ax.set_xlabel(xlabel)

    # Legend mapping each color back to its modality
    handles = [mpatches.Patch(color=palette[modality], label=modality) for modality in data[modality_column].unique()]
    ax.legend(handles=handles, title="Modality", bbox_to_anchor=(1.01, 1), loc="upper left")

    ax.spines[["top", "right"]].set_visible(False)
    return ax

1.2. Read MS-proteomics data

We will first focus on the MS-proteomics data. We downloaded the gene-level aggregated data from the corresponding pride repository PXD020191.

In a first step, we will ingest and preprocess the MS-proteomics data separately from the other modalities.

report_path = apt.data.get_data("lehtioe2021_pg")

column_mapping = {
    PGCols.GENES: "Protein accession",
    PGCols.DESCRIPTION: "Description",
    "protein_coverage": "Coverage",
    "n_proteins": "Amount Proteins",
    PGCols.PROTEINS: "Proteins in group",
}
measurement_regex = r"^set\d{2}_tmt10plex_\d{3}[NC]?$"

adata = apt.io.read_pg_table(
    report_path,
    search_engine="diann",
    column_mapping=column_mapping,
    measurement_regex=measurement_regex,
)
adata

/Users/lucas-diedrich/Documents/Projects/scverse/alphatools/programming/alphatools/docs/notebooks/studies/gene-symbols_table.txt already exists (46.51789855957031 MB)

AnnData object with n_obs × n_vars = 170 × 14058
    var: 'description', 'protein_coverage', 'n_proteins', 'proteins'

We obtained 170 samples with a total of 14058 unique protein groups.

Now, we read the available clinical metadata, which was extracted from the supplementary information of the original study.

metadata_path = apt.data.get_data("lehtioe2021_metadata")
metadata = pd.read_csv(metadata_path, sep="\t", index_col="tmt_set_and_label")
metadata

/Users/lucas-diedrich/Documents/Projects/scverse/alphatools/programming/alphatools/docs/notebooks/studies/lethioe2021_metadata.tsv already exists (0.04194831848144531 MB)

	sample_id	proteome_subtype	age	sex	previous_histology	who2015_histology	ac_mrna_subtype_tcga	sqcc_mrna_subtype_tcga	nsclc_mrna_subtype_karlsson_et_al._2017	smokingstatus	...	zfhx3_mut	tmb	cancer_testis_antigens	ncps	purity	methylation_scores_overall	methylation_scores_promoter	methylationpresent	panelseqpresent	transcriptomicspresent
tmt_set_and_label
set13_tmt10plex_129N	LU85	1	67.0	F	AC	AC	Terminal respiratory unit	NaN	5.0	Never_smoker	...	0.0	5.86458	0	8	0.30	0.491897	0.158030	True	True	True
set13_tmt10plex_128C	LU277	1	80.0	F	AC	AC	Terminal respiratory unit	NaN	5.0	Never_smoker	...	0.0	2.93229	2	5	0.25	0.525326	0.176840	True	True	True
set07_tmt10plex_126	LU376	1	72.0	M	AC	AC	Terminal respiratory unit	NaN	9.0	Never_smoker	...	0.0	2.93229	0	2	0.18	0.510334	0.170185	True	True	True
set11_tmt10plex_129C	LU1115	1	80.0	M	AC	AC	Terminal respiratory unit	NaN	5.0	Never_smoker	...	0.0	5.86458	4	5	0.31	0.505183	0.163569	True	True	True
set12_tmt10plex_128N	LU804	1	73.0	F	AC	AC	Terminal respiratory unit	NaN	5.0	Never_smoker	...	0.0	5.27812	1	7	0.38	0.533692	0.181161	True	True	True
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
set04_tmt10plex_128C	LU745	6	70.0	M	SqCC	SqCC	NaN	Classical	2.0	Former	...	0.0	8.79687	2	22	0.62	NaN	NaN	False	True	True
set14_tmt10plex_128C	LU1020	6	80.0	M	SqCC	SqCC	NaN	Classical	2.0	Former	...	0.0	24.04480	5	13	0.57	0.486081	0.154912	True	True	True
set06_tmt10plex_129C	LU372	6	69.0	M	SqCC	SqCC	NaN	Classical	2.0	Current	...	0.0	15.24790	9	27	0.68	0.437679	0.136954	True	True	True
set05_tmt10plex_128C	LU1040	6	71.0	F	SqCC	SqCC	NaN	Classical	2.0	Former	...	0.0	19.93960	4	13	0.37	0.510035	0.171375	True	True	True
set08_tmt10plex_130N	LU398	6	77.0	F	SqCC	SqCC	NaN	Classical	2.0	Former	...	0.0	12.90210	3	24	0.36	0.458176	0.153726	True	True	True

141 rows × 57 columns

We add the metadata to the anndata object with alphapepttools.pp.add_metadata and extract batch information (tmt_plex) from the TMT sample id.

adata = apt.pp.add_metadata(adata, metadata, axis=0, keep_data_shape=False)
adata.obs = adata.obs.reset_index(names="tmt_sample_id").set_index("sample_id")
adata.obs["tmt_plex"] = adata.obs["tmt_sample_id"].str.extract("set([0-9]+)_")[0].tolist()
adata.obs["tmt_channel"] = adata.obs["tmt_sample_id"].str.extract("_([A-Za-z0-9]+)$")[0].tolist()
adata

AnnData object with n_obs × n_vars = 141 × 14058
    obs: 'tmt_sample_id', 'proteome_subtype', 'age', 'sex', 'previous_histology', 'who2015_histology', 'ac_mrna_subtype_tcga', 'sqcc_mrna_subtype_tcga', 'nsclc_mrna_subtype_karlsson_et_al._2017', 'smokingstatus', 'smokingstatus_binary', 'stage_general', 'stage_detailed', 'os_years', 'osbin', 'relapse_years', 'relapse_bin', 'adjuvant_therapy', 'pdl1_ihc', 'cd3_stromal_ihc', 'cd8_stromal_ihc', 'stroma_estimate', 'immune_estimate', 'arid1a_mut', 'arid2_mut', 'asxl1_mut', 'atm_mut', 'blm_mut', 'braf_mut', 'cdkn2a_mut', 'egfr_mut', 'fat1_mut', 'keap1_mut', 'kmt2c_mut', 'kmt2d_mut', 'kras_mut', 'nf1_mut', 'nfe2l2_mut', 'notch2_mut', 'pik3ca_mut', 'pold1_mut', 'pten_mut', 'rb1_mut', 'rbm10_mut', 'smarca4_mut', 'stk11_mut', 'tp53_mut', 'zfhx3_mut', 'tmb', 'cancer_testis_antigens', 'ncps', 'purity', 'methylation_scores_overall', 'methylation_scores_promoter', 'methylationpresent', 'panelseqpresent', 'transcriptomicspresent', 'tmt_plex', 'tmt_channel'
    var: 'description', 'protein_coverage', 'n_proteins', 'proteins'

The study design included 16 rotating in-set references, and one full reference set (Set17). These samples were used for normalization and QC, but do not contain biological information (i.e. they are absent from the metadata table). They are removed by the apt.pp.add_metadata call, yielding 141 biological samples.

1.3. Quality control + Preprocessing

To evaluate between-run variation, we check the missingness patterns of each batch. Note that we do not look at sample-level aggregates of intensity values as the data has already been transformed by the search engine (see methods of the original paper) and aggregated intensity values are difficult to interpret at this stage.

apt.metrics.fraction_complete(adata, axis="obs", column="fraction_complete_obs")
apt.metrics.fraction_complete(adata, axis="var", column="fraction_complete_var")

fig, axm = apt.pl.create_figure(1, 2, figsize=(8, 3), width_ratios=(1, 4))

apt.pl.violinplot(ax=axm[0], data=adata, direct_columns=["fraction_complete_obs"])
apt.pl.label_axes(ax=axm[0], xlabel=None, ylabel="Feature completeness", title="Feature completeness")
axm[0].set_ylim(0.0, 1.0)


apt.pl.boxplot(ax=axm[1], data=adata, grouping_column="tmt_plex", value_column="fraction_complete_obs")

apt.pl.label_axes(
    ax=axm[1], xlabel="TMT Plex (Batch)", ylabel="Feature completeness", title="Feature completeness per batch"
)
axm[1].set_ylim(0.0, 1.0)

plt.tight_layout()
plt.show()

We observe a low variability of the feature missingness between samples and between batches, indicating a good comparability of the samples between batches.

We also look at the feature-wise missingness: Most features are 100% complete, suggesting deep proteome profiling

# Histogram of detection frequency combined with a rank plot
fig, axm = apt.pl.create_figure(1, 1, figsize=(8, 4))

# Histogram
ax = axm.next()
apt.pl.histogram(
    data=adata,
    value_column="fraction_complete_var",
    bins=50,
    ax=ax,
    color=apt.pl.BaseColors.get("green"),
    hist_kwargs={"histtype": "stepfilled", "edgecolor": "black"},
)

apt.pl.label_axes(
    ax,
    xlabel="Fraction of samples with detection",
    ylabel="Number of proteins",
    title="Protein Detection Frequency Distribution",
    enumeration="A",
)

Following the original study, we drop all incomplete features.

adata = apt.pp.filter_data_completeness(adata, max_missing=0, action="drop")

Transformation and centering

We follow the original preprocessing steps and perform log2-transformation with alphapepttools.pp.nanlog and median centering of the features. We store the intermediate results in separate layers of the anndata object.

adata.layers["raw"] = adata.X.copy()

apt.pp.nanlog(adata)
adata.layers["log2"] = adata.X.copy()

adata.X = (adata.X - np.nanmedian(adata.X, axis=0)).copy()
adata.layers["median_centered"] = adata.X.copy()

Visual inspection

We check the data qualitatively and investigate whether histological subtypes separate, representing expected biological variation.

apt.tl.pca(adata)

fig, axm = apt.pl.create_figure(1, 2, figsize=(8, 4))


apt.pl.plot_pca(adata, y_column=2, ax=axm[0], color_map_column="tmt_plex", legend="auto")
apt.pl.label_axes(ax=axm[0], title="PCA (TMT plex)")


apt.pl.plot_pca(adata, y_column=2, ax=axm[1], color_map_column="who2015_histology", legend="auto")
apt.pl.label_axes(ax=axm[1], title="PCA (Histological subtype)")

plt.tight_layout()

Indeed, as expected, the data separates based on histological subtypes and not TMT runs.

2. Multimodal Analysis | Integrate with complementary modalities

In the original paper, the proteomics data was clustered and the proteomic phenotypes were associated post-hoc with complementary RNA and methylation microarray data.

Here, we try to model these different modalities all at once with a MOFA+ analysis. MOFA+ is a probabilistic latent factor model that tries to identify shared axes of variation in a multi-modal dataset. Similar to other matrix factorization techniques like principal component analysis, MOFA+ assumes that the joint observed -omics measurements \(X\) of \(P\) features in \(N\) samples can be explained by a small number (\(K\)) of unobserved, latent variables, that are associated with specific biological processes in the sample:

\[ X^{N\times P} = A^{N\times K} \times L^{K\times P} + \varepsilon \]

MOFA+ returns an estimate of the activity \(A\) of these latent factors in the individual samples as a \(N\) sample \(\times\) \(K\) latent factors matrix and the factor loadings matrix \(L\) of the shape \(K\) factors \(\times\) \(P\) features. As all modalities are jointly modelled, the factor loadings contain estimates of the model about which features are coregulated across all -omics layers.

References

• MOFA+ Argelaguet, R. et al. MOFA+: a statistical framework for comprehensive integration of multi-modal single-cell data. Genome Biol 21, 111 (2020).

• muon: Bredikhin, D., Kats, I. & Stegle, O. MUON: multimodal omics analysis framework. Genome Biology 23, 42 (2022).

2.1. Load complementary modalities

First, we load the other available -omics measurements. Note that we obtained the preprocessed data from the supplementary information of the original paper and used anndata.read_excel to convert the datasets into anndata objects.

RNA

We download the data and read the anndata object

rna_path = apt.data.get_data("lehtioe2021_rna")
rna = ad.read_h5ad(rna_path)
rna

/Users/lucas-diedrich/Documents/Projects/scverse/alphatools/programming/alphatools/docs/notebooks/studies/karlsson2014_nslc_rna.h5ad already exists (13.766532897949219 MB)

AnnData object with n_obs × n_vars = 118 × 14548
    obs: 'proteome_subtype', 'age', 'sex', 'previous_histology', 'who2015_histology', 'ac_mrna_subtype_tcga', 'sqcc_mrna_subtype_tcga', 'nsclc_mrna_subtype_karlsson_et_al._2017', 'smokingstatus', 'smokingstatus_binary', 'stage_general', 'stage_detailed', 'os_years', 'osbin', 'relapse_years', 'relapse_bin', 'adjuvant_therapy', 'pdl1_ihc', 'cd3_stromal_ihc', 'cd8_stromal_ihc', 'stroma_estimate', 'immune_estimate', 'arid1a_mut', 'arid2_mut', 'asxl1_mut', 'atm_mut', 'blm_mut', 'braf_mut', 'cdkn2a_mut', 'egfr_mut', 'fat1_mut', 'keap1_mut', 'kmt2c_mut', 'kmt2d_mut', 'kras_mut', 'nf1_mut', 'nfe2l2_mut', 'notch2_mut', 'pik3ca_mut', 'pold1_mut', 'pten_mut', 'rb1_mut', 'rbm10_mut', 'smarca4_mut', 'stk11_mut', 'tp53_mut', 'zfhx3_mut', 'tmb', 'cancer_testis_antigens', 'ncps', 'purity', 'methylation_scores_overall', 'methylation_scores_promoter', 'methylationpresent', 'panelseqpresent', 'transcriptomicspresent'

Methylation

Next, we load the methylation data. We need to convert the available \(\beta\) values (representing the fraction of methylated CpG sites with methylation among all measurements for a certain site) to \(M\)-values, to match MOFA+’s expectation that the data is approximately normally distributed. See this recommendation by the authors

methylation_path = apt.data.get_data("lehtioe2021_methylation")
methylation = ad.read_h5ad(methylation_path)

methylation.layers["beta_values"] = methylation.X.copy()
# Convert beta values to M-values
# See recommendation by original authors: https://github.com/bioFAM/MOFA/issues/44#issuecomment-526280692
mvalue(methylation, alpha=0.01, layer=None, copy=False)
methylation

/Users/lucas-diedrich/Documents/Projects/scverse/alphatools/programming/alphatools/docs/notebooks/studies/karlsson2014_nslc_methylation.h5ad already exists (12.006929397583008 MB)

AnnData object with n_obs × n_vars = 113 × 11937
    obs: 'proteome_subtype', 'age', 'sex', 'previous_histology', 'who2015_histology', 'ac_mrna_subtype_tcga', 'sqcc_mrna_subtype_tcga', 'nsclc_mrna_subtype_karlsson_et_al._2017', 'smokingstatus', 'smokingstatus_binary', 'stage_general', 'stage_detailed', 'os_years', 'osbin', 'relapse_years', 'relapse_bin', 'adjuvant_therapy', 'pdl1_ihc', 'cd3_stromal_ihc', 'cd8_stromal_ihc', 'stroma_estimate', 'immune_estimate', 'arid1a_mut', 'arid2_mut', 'asxl1_mut', 'atm_mut', 'blm_mut', 'braf_mut', 'cdkn2a_mut', 'egfr_mut', 'fat1_mut', 'keap1_mut', 'kmt2c_mut', 'kmt2d_mut', 'kras_mut', 'nf1_mut', 'nfe2l2_mut', 'notch2_mut', 'pik3ca_mut', 'pold1_mut', 'pten_mut', 'rb1_mut', 'rbm10_mut', 'smarca4_mut', 'stk11_mut', 'tp53_mut', 'zfhx3_mut', 'tmb', 'cancer_testis_antigens', 'ncps', 'purity', 'methylation_scores_overall', 'methylation_scores_promoter', 'methylationpresent', 'panelseqpresent', 'transcriptomicspresent'
    var: 'probeID', 'ID'
    layers: 'beta_values'

2.2. Create mudata object

To jointly represent the cohort across all -omics layers, we create a mudata object. mudata binds multiple modalities, represented as anndata objects, into a single container. Each modality fits into a .mod slot of the object, and the resulting data structure behaves very similarly to the familiar anndata containers.

mdata = md.MuData({"proteins": adata, "rna": rna, "methylation": methylation})

# We remove columns with missing values
mdata.obs = mdata.obs.merge(metadata.set_index("sample_id"), left_index=True, right_index=True).dropna(axis=1)
mdata.var["modality"] = (
    pd.DataFrame.from_dict({key: mdata.varmap[key].squeeze() for key in mdata.varmap}).idxmax(axis=1).tolist()
)

mdata

MuData object with n_obs × n_vars = 141 × 36278
  obs:	'proteome_subtype', 'age', 'sex', 'previous_histology', 'who2015_histology', 'smokingstatus', 'smokingstatus_binary', 'os_years', 'osbin', 'stroma_estimate', 'immune_estimate', 'cancer_testis_antigens', 'ncps', 'methylationpresent', 'panelseqpresent', 'transcriptomicspresent'
  var:	'modality'
  3 modalities
    proteins:	141 x 9793
      obs:	'tmt_sample_id', 'proteome_subtype', 'age', 'sex', 'previous_histology', 'who2015_histology', 'ac_mrna_subtype_tcga', 'sqcc_mrna_subtype_tcga', 'nsclc_mrna_subtype_karlsson_et_al._2017', 'smokingstatus', 'smokingstatus_binary', 'stage_general', 'stage_detailed', 'os_years', 'osbin', 'relapse_years', 'relapse_bin', 'adjuvant_therapy', 'pdl1_ihc', 'cd3_stromal_ihc', 'cd8_stromal_ihc', 'stroma_estimate', 'immune_estimate', 'arid1a_mut', 'arid2_mut', 'asxl1_mut', 'atm_mut', 'blm_mut', 'braf_mut', 'cdkn2a_mut', 'egfr_mut', 'fat1_mut', 'keap1_mut', 'kmt2c_mut', 'kmt2d_mut', 'kras_mut', 'nf1_mut', 'nfe2l2_mut', 'notch2_mut', 'pik3ca_mut', 'pold1_mut', 'pten_mut', 'rb1_mut', 'rbm10_mut', 'smarca4_mut', 'stk11_mut', 'tp53_mut', 'zfhx3_mut', 'tmb', 'cancer_testis_antigens', 'ncps', 'purity', 'methylation_scores_overall', 'methylation_scores_promoter', 'methylationpresent', 'panelseqpresent', 'transcriptomicspresent', 'tmt_plex', 'tmt_channel', 'fraction_complete_obs'
      var:	'description', 'protein_coverage', 'n_proteins', 'proteins', 'fraction_complete_var'
      uns:	'variance_pca_obs'
      obsm:	'X_pca_obs'
      varm:	'PCs_pca_obs'
      layers:	'raw', 'log2', 'median_centered'
    rna:	118 x 14548
      obs:	'proteome_subtype', 'age', 'sex', 'previous_histology', 'who2015_histology', 'ac_mrna_subtype_tcga', 'sqcc_mrna_subtype_tcga', 'nsclc_mrna_subtype_karlsson_et_al._2017', 'smokingstatus', 'smokingstatus_binary', 'stage_general', 'stage_detailed', 'os_years', 'osbin', 'relapse_years', 'relapse_bin', 'adjuvant_therapy', 'pdl1_ihc', 'cd3_stromal_ihc', 'cd8_stromal_ihc', 'stroma_estimate', 'immune_estimate', 'arid1a_mut', 'arid2_mut', 'asxl1_mut', 'atm_mut', 'blm_mut', 'braf_mut', 'cdkn2a_mut', 'egfr_mut', 'fat1_mut', 'keap1_mut', 'kmt2c_mut', 'kmt2d_mut', 'kras_mut', 'nf1_mut', 'nfe2l2_mut', 'notch2_mut', 'pik3ca_mut', 'pold1_mut', 'pten_mut', 'rb1_mut', 'rbm10_mut', 'smarca4_mut', 'stk11_mut', 'tp53_mut', 'zfhx3_mut', 'tmb', 'cancer_testis_antigens', 'ncps', 'purity', 'methylation_scores_overall', 'methylation_scores_promoter', 'methylationpresent', 'panelseqpresent', 'transcriptomicspresent'
    methylation:	113 x 11937
      obs:	'proteome_subtype', 'age', 'sex', 'previous_histology', 'who2015_histology', 'ac_mrna_subtype_tcga', 'sqcc_mrna_subtype_tcga', 'nsclc_mrna_subtype_karlsson_et_al._2017', 'smokingstatus', 'smokingstatus_binary', 'stage_general', 'stage_detailed', 'os_years', 'osbin', 'relapse_years', 'relapse_bin', 'adjuvant_therapy', 'pdl1_ihc', 'cd3_stromal_ihc', 'cd8_stromal_ihc', 'stroma_estimate', 'immune_estimate', 'arid1a_mut', 'arid2_mut', 'asxl1_mut', 'atm_mut', 'blm_mut', 'braf_mut', 'cdkn2a_mut', 'egfr_mut', 'fat1_mut', 'keap1_mut', 'kmt2c_mut', 'kmt2d_mut', 'kras_mut', 'nf1_mut', 'nfe2l2_mut', 'notch2_mut', 'pik3ca_mut', 'pold1_mut', 'pten_mut', 'rb1_mut', 'rbm10_mut', 'smarca4_mut', 'stk11_mut', 'tp53_mut', 'zfhx3_mut', 'tmb', 'cancer_testis_antigens', 'ncps', 'purity', 'methylation_scores_overall', 'methylation_scores_promoter', 'methylationpresent', 'panelseqpresent', 'transcriptomicspresent'
      var:	'probeID', 'ID'
      layers:	'beta_values'

2.3. Run MOFA+ factor analysis

MOFA+ can be run with the muon framework. Here, we define some parameters of the MOFA+ object, largely following the recommendations of the authors and default settings.

In addition, we model all samples, even if a modality is missing "use_obs": "union" to leverage the full dataset, and assume that all modalities are normally distributed ("likelihoods": "gaussian"). The assumption of a gaussian distribution is acceptable as we log-transformed the data beforehand.

In total, we model 25 factors ("n_factors": 25), which is in line with the assumption that at most 25 different latent processes explain the observed -omics measurements. By setting ard_factors to True, MOFA+ automatically regularizes the usage of the factors, i.e. the activities of unnecessary factors are shrunk toward zero.

mofa_kwargs = {
    "use_obs": "union",  # Use all available samples
    "n_factors": 25,  # 25 factors
    # Use defaults as recommended by the authors
    "scale_groups": True,
    "scale_views": True,
    "ard_factors": True,
    "ard_weights": True,
    "spikeslab_weights": True,
    "spikeslab_factors": True,
    "likelihoods": "gaussian",
    # To reduce runtime, set convergence mode to "fast"
    "convergence_mode": "slow",
    # Fitting parameters - set seed to 1 for reproducibility
    "seed": 1,
}

# Runs 5min on MacOS + M4
# ard_weights: True, as only a multiple views are provided
muon.tl.mofa(mdata, **mofa_kwargs, outfile="mofa-multimodal.model.hdf5")

        #########################################################
        ###           __  __  ____  ______                    ### 
        ###          |  \/  |/ __ \|  ____/\    _             ### 
        ###          | \  / | |  | | |__ /  \ _| |_           ### 
        ###          | |\/| | |  | |  __/ /\ \_   _|          ###
        ###          | |  | | |__| | | / ____ \|_|            ###
        ###          |_|  |_|\____/|_|/_/    \_\              ###
        ###                                                   ### 
        ######################################################### 
         


Scaling views to unit variance...

Scaling groups to unit variance...

Loaded view='proteins' group='group1' with N=141 samples and D=9793 features...
Loaded view='rna' group='group1' with N=141 samples and D=14548 features...
Loaded view='methylation' group='group1' with N=141 samples and D=11937 features...


Model options:
- Automatic Relevance Determination prior on the factors: True
- Automatic Relevance Determination prior on the weights: True
- Spike-and-slab prior on the factors: True
- Spike-and-slab prior on the weights: True
Likelihoods:
- View 0 (proteins): gaussian
- View 1 (rna): gaussian
- View 2 (methylation): gaussian




######################################
## Training the model with seed 1 ##
######################################



Converged!



#######################
## Training finished ##
#######################


Warning: Output file mofa-multimodal.model.hdf5 already exists, it will be replaced
Saving model in mofa-multimodal.model.hdf5...
Saved MOFA embeddings in .obsm['X_mofa'] slot and their loadings in .varm['LFs'].

2.4. Exploration

We can explore the fitted MOFA+ model

Plot explained variance as proxy for relevance

MOFA+ provides an estimate of how much variance is explained by each modality and each factor which we can extract and plot. This gives us an idea of which factors explain a relevant amount of the data and which factors we should prioritize in the downstream analysis.

explained_variance = (
    pd.DataFrame.from_dict(mdata.uns["mofa"]["variance"], orient="index") / 100
)  # (MOFA+ reports in percent, we scale to floats)

fig, axs = plt.subplots(
    2,
    2,
    figsize=(15, 4),
    gridspec_kw={"height_ratios": [0.4, 0.6], "width_ratios": [0.97, 0.03], "hspace": 0},
    squeeze=True,
)

explained_variance.T.plot.bar(ax=axs[0, 0], color=["#572c92", "#9b76ce", "#decff4"], stacked=True)
axs[0, 0].set_xticks([])
axs[0, 0].set_yticks(np.arange(0, 0.31, 0.05), minor=True)
axs[0, 0].grid(axis="y", alpha=0.5, which="minor")
axs[0, 0].grid(axis="y", alpha=0.5, which="major")

axs[0, 0].spines[["top", "right"]].set_visible(False)
axs[0, 0].set_ylabel("Total variance")
axs[0, 1].set_axis_off()


sns.heatmap(
    explained_variance,
    square=True,
    annot=True,
    fmt=".2f",
    ax=axs[1, 0],
    cbar_ax=axs[1, 1],
    cbar_kws={"shrink": 0.1, "label": "Variance explained"},
    vmin=0,
    vmax=0.2,
    cmap="Purples",
)

plt.show()

We see that the first ~7 factors each explain more than 5% of variance. We can also see that some factors share variance across all modalities (most prominently factor 1), while in other cases, the explained variance is dominated by certain modalities. For example, factors 0, 2, and 3 are dominated by the RNA and Protein modalities, while factors 4 and 5 are only relevant for the methylation measurements.

Factor activity analysis

We can plot the MOFA+ embeddings similar to a PCA. Note that we plot the first and third components of the MOFA+ embedding, as the 2nd component does not capture variation aligned with histological differences:

fig, axm = apt.pl.create_figure(1, 1, figsize=(5, 4))

muon.pl.embedding(
    mdata,
    basis="X_mofa",
    color="who2015_histology",
    dimensions=(0, 2),
    palette=apt.pl.BasePalettes().get("qualitative"),
    size=200,
    ax=axm[0],
    title="",
    show=False,
)

apt.pl.label_axes(title="MOFA+ (Histological subtype)", xlabel="MOFA+ Factor 0", ylabel="MOFA+ Factor 2", ax=axm[0])

We see that factor 2 strongly separates the subtypes squamous cell carcinoma (SqCC), adenocarcinoma (AC), and Large Cell Neuroendocrine Carcinoma (LCNEC).

Investigate factor loadings

In the following, we will investigate factor 2 of the MOFA+ model in more detail, given that it separates the 3 major histological subtypes in the cohort. For that, we first extract the MOFA+ - inferred gene weights per modality from the mudata object:

mofa_loadings = pd.DataFrame(
    mdata.varm["LFs"], index=pd.MultiIndex.from_frame(mdata.var.reset_index(names="genes")[["genes", "modality"]])
)

mofa_loadings

		0	1	2	3	4	5	6	7	8	9	...	15	16	17	18	19	20	21	22	23	24
genes	modality
A1BG	proteins	0.172913	-0.026706	-0.002021	-0.098686	0.001982	0.000210	-0.012240	-0.144932	-0.014722	0.000722	...	0.681463	-0.175271	0.038054	-0.040847	-1.087811	0.005504	0.025088	0.007520	-0.050735	0.003166
A2M	proteins	0.340784	-0.023044	0.035473	-0.272198	0.010998	0.001424	0.002820	-0.002830	-0.016485	-0.001102	...	0.891250	0.042361	0.016262	-0.066233	-0.575053	0.016547	-0.010618	0.005945	-0.078280	0.004734
A2ML1	proteins	-0.132505	-0.132383	0.575635	-0.289802	-0.002140	-0.000554	-0.217328	-0.050090	0.000430	0.002430	...	-0.024155	-0.000864	0.120288	0.015771	-0.046177	-0.016202	-0.122171	0.006784	0.011952	0.009444
A4GALT	proteins	-0.066090	0.127836	0.352663	-0.090918	0.016027	0.000483	-0.208224	-0.338341	-0.007403	0.005031	...	-0.000371	0.233258	0.494228	0.062048	-1.614278	-0.048914	0.305061	-0.020099	0.178623	-0.005935
AAAS	proteins	-0.232066	0.069589	-0.001847	0.032239	0.018148	0.002200	-0.004632	0.000892	0.000540	-0.002735	...	-0.019154	-0.132013	-0.010862	0.053626	-0.354351	0.002360	-0.036245	0.035460	0.020794	-0.000196
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
ZXDC	methylation	0.000017	-0.000212	-0.009404	0.008360	-0.290754	0.002572	-0.000100	-0.000274	-0.000365	0.011620	...	0.000092	-0.002573	0.004897	-0.000717	-0.000716	-0.009946	-0.001911	-0.000028	0.016411	-0.153422
ZYG11B	methylation	-0.003426	0.001652	0.002152	0.081731	-0.000691	-0.028402	0.000542	-0.000204	-0.000411	-0.060362	...	0.004852	-0.001320	-0.002273	-0.001021	-0.000299	0.073336	-0.000647	0.000009	-0.032413	-0.009731
ZYX	methylation	0.000226	0.037661	-0.008152	0.008459	-1.956825	0.023051	0.011551	0.000038	-0.000253	-0.085083	...	0.000873	-0.003955	0.029490	-0.001608	-0.000294	0.013852	-0.013685	0.000009	0.038148	-0.102275
ZZEF1	methylation	0.000446	-0.009386	-0.005098	0.000230	-0.000111	-0.000078	-0.000112	0.000010	-0.001361	0.427192	...	0.000463	0.000119	0.004909	-0.000308	0.003108	-0.005345	-0.006716	0.000003	0.093664	-0.482233
ZZZ3	methylation	0.003696	-0.000245	0.002400	0.004755	-0.003714	0.003177	0.000238	-0.001424	-0.000038	0.298336	...	-0.003184	0.007405	-0.001292	0.001310	-0.000815	-0.016864	0.003076	-0.000004	-0.037704	0.209354

36278 rows × 25 columns

We can first investigate, which genes are particularly important for factor 2 by plotting the top gene weights in this factor:

fig, ax = plt.subplots(figsize=(4, 8))

plot_modality_loadings(
    mofa_loadings.nlargest(30, 2).reset_index(),
    score=2,
    feature="genes",
    palette=palette,
    xlabel="Factor 2 loading",
    ax=ax,
)

plt.tight_layout()
plt.show()

Both RNA and protein features dominate the top loadings (features that are associated with high factor activity, i.e. Squamous Cell Carcinoma). Interestingly, the top loading is associated with TP63 (RNA; Rank 26 for TP63-protein), a commonly used immunohistochemistry marker for squamous cell carcinoma (e.g. Rossi et al,2009).

Functional annotation with decoupler.mt.gsea

Next, we will run gene set enrichment analysis via the decoupler framework to functionally annotate the inferred factors.

For that, we will first download the hallmarks geneset resource with decoupler and clean it up:

# recommendation by Reimand et al, 2019 https://doi.org/10.1038/s41596-018-0103-9 (there: 350)
MIN_PATHWAY_SIZE = 5
MAX_PATHWAY_SIZE = 350
PATHWAY_DATABASE = "hallmark"

resource = dc.op.resource("MSigDB", organism="human")

network = (
    resource.loc[resource["collection"] == PATHWAY_DATABASE]
    # .groupby("geneset")
    # .filter(lambda pathway: (len(pathway) >= MIN_PATHWAY_SIZE) & (len(pathway) <= MAX_PATHWAY_SIZE))
    .drop_duplicates()
    # Decoupler expects the network as bipartite graph in which genesets (sources) point to genes (targets)
    .rename(columns={"genesymbol": "target", "geneset": "source"})
)

Then, we run the gene set enrichment analysis on every modality separately:

results = []
for mod, loadings in mofa_loadings.groupby("modality", observed=True):
    # Decoupler expects a (factor x genes) dataframe without duplicates
    decoupler_input = (
        loadings.reset_index(level="genes").drop_duplicates(subset="genes", keep=False).set_index("genes").T
    )
    scores, pvals = dc.mt.gsea(data=decoupler_input, net=network)

    # Generate a long dataframe that contains factor number, pathway, enrichment score, FDR, and indicator of modality
    results.append(
        pd.merge(  # noqa: PD015
            scores.melt(ignore_index=False, var_name="pathway", value_name="score").reset_index(names="factor"),
            pvals.melt(ignore_index=False, var_name="pathway", value_name="pval_adj").reset_index(names="factor"),
            on=["factor", "pathway"],
        ).assign(modality=mod)
    )

# Concatenate results
gsea_results = pd.concat(results, axis=0)
gsea_results

	factor	pathway	score	pval_adj	modality
0	0	HALLMARK_ADIPOGENESIS	-1.084259	0.784876	methylation
1	1	HALLMARK_ADIPOGENESIS	-1.162857	0.212362	methylation
2	2	HALLMARK_ADIPOGENESIS	-1.031550	0.923214	methylation
3	3	HALLMARK_ADIPOGENESIS	-0.948506	0.856833	methylation
4	4	HALLMARK_ADIPOGENESIS	-1.014372	1.000000	methylation
...	...	...	...	...	...
1245	20	HALLMARK_XENOBIOTIC_METABOLISM	1.583415	0.000000	rna
1246	21	HALLMARK_XENOBIOTIC_METABOLISM	1.429975	0.060680	rna
1247	22	HALLMARK_XENOBIOTIC_METABOLISM	-1.101423	0.596774	rna
1248	23	HALLMARK_XENOBIOTIC_METABOLISM	0.815107	1.000000	rna
1249	24	HALLMARK_XENOBIOTIC_METABOLISM	1.157962	0.965313	rna

3750 rows × 5 columns

We now select which pathways to display. Two parameters control this:

• FACTOR — the MOFA+ factor whose loadings were tested for enrichment. We use factor "2", as it separates the histological subtypes.

• SIGNIFICANCE_THRESHOLD — the largest adjusted p-value (FDR) a pathway may have to count as significant. We keep pathways at FDR ≤ 0.01.

We subset the enrichment results to the significant pathways of this factor and plot them:

# Parameters
FACTOR = "2"
SIGNIFICANCE_THRESHOLD = 1e-2

significant_pathways = gsea_results.loc[
    (gsea_results["factor"] == FACTOR) & (gsea_results["pval_adj"] <= SIGNIFICANCE_THRESHOLD)
].sort_values(by="score")

fig, ax = plt.subplots(figsize=(4, 8))

plot_modality_loadings(
    significant_pathways,
    score="score",
    feature="pathway",
    palette=palette,
    xlabel="Enrichment score",
    ax=ax,
)

plt.tight_layout()
plt.show()

Most pathways cluster around P53 and MYC signaling as well as metabolic pathways (mTORC1, glycolysis, etc.), as well as extracellular matrix remodeling. Interestingly, the p53 pathway is enriched in factor 2 at both the protein and the RNA expression level. We also see that many significantly enriched pathways are shared between protein and RNA levels, while others are distinct for the protein level (e.g. HALLMARK_ESTROGEN_RESPONSE_LATE). Overall, it is interesting to see that both modality-specific and joint processes can be observed in this factor.

Separation of samples

We can check whether the joint modeling of proteomics, RNA, and methylation data improved the separation of samples based on their histological phenotype compared to only the proteomics measurements. Here, we compute the silhouette width of the histological subtypes on the MOFA+ embeddings for the multimodal modeling case and the proteomics-only modeling case, using the proteomics PCA embedding as a reference point.

The silhouette width (\(s\)) is a measure of how similar samples from the same condition are compared to samples from other conditions (separation). It ranges from −1 to +1, where higher values indicate better separation of conditions.

For comparison, we also compute principal components (PCA) and MOFA+ embedding for the single-modality measurements alone, using the same settings:

# Run PCA on every individual modality
sc.pp.pca(methylation, n_comps=mofa_kwargs["n_factors"])
sc.pp.pca(rna, n_comps=mofa_kwargs["n_factors"])
apt.tl.pca(adata, n_comps=mofa_kwargs["n_factors"])

# Run MOFA+
# Runs ca. 6 min on MacOS with M4
muon.tl.mofa(adata, **mofa_kwargs, outfile="mofa-proteomics.model.hdf5")
muon.tl.mofa(rna, **mofa_kwargs, outfile="mofa-rna.model.hdf5")
muon.tl.mofa(methylation, **mofa_kwargs, outfile="mofa-methylation.model.hdf5")

        #########################################################
        ###           __  __  ____  ______                    ### 
        ###          |  \/  |/ __ \|  ____/\    _             ### 
        ###          | \  / | |  | | |__ /  \ _| |_           ### 
        ###          | |\/| | |  | |  __/ /\ \_   _|          ###
        ###          | |  | | |__| | | / ____ \|_|            ###
        ###          |_|  |_|\____/|_|/_/    \_\              ###
        ###                                                   ### 
        ######################################################### 
         


Scaling views to unit variance...

Scaling groups to unit variance...

Loaded view='data' group='group1' with N=141 samples and D=9793 features...


Model options:
- Automatic Relevance Determination prior on the factors: True
- Automatic Relevance Determination prior on the weights: True
- Spike-and-slab prior on the factors: True
- Spike-and-slab prior on the weights: True
Likelihoods:
- View 0 (data): gaussian




######################################
## Training the model with seed 1 ##
######################################



Converged!



#######################
## Training finished ##
#######################


Warning: Output file mofa-proteomics.model.hdf5 already exists, it will be replaced
Saving model in mofa-proteomics.model.hdf5...
Saved MOFA embeddings in .obsm['X_mofa'] slot and their loadings in .varm['LFs'].

        #########################################################
        ###           __  __  ____  ______                    ### 
        ###          |  \/  |/ __ \|  ____/\    _             ### 
        ###          | \  / | |  | | |__ /  \ _| |_           ### 
        ###          | |\/| | |  | |  __/ /\ \_   _|          ###
        ###          | |  | | |__| | | / ____ \|_|            ###
        ###          |_|  |_|\____/|_|/_/    \_\              ###
        ###                                                   ### 
        ######################################################### 
         


Scaling views to unit variance...

Scaling groups to unit variance...

Loaded view='data' group='group1' with N=118 samples and D=14548 features...


Model options:
- Automatic Relevance Determination prior on the factors: True
- Automatic Relevance Determination prior on the weights: True
- Spike-and-slab prior on the factors: True
- Spike-and-slab prior on the weights: True
Likelihoods:
- View 0 (data): gaussian




######################################
## Training the model with seed 1 ##
######################################



Converged!



#######################
## Training finished ##
#######################


Warning: Output file mofa-rna.model.hdf5 already exists, it will be replaced
Saving model in mofa-rna.model.hdf5...
Saved MOFA embeddings in .obsm['X_mofa'] slot and their loadings in .varm['LFs'].

        #########################################################
        ###           __  __  ____  ______                    ### 
        ###          |  \/  |/ __ \|  ____/\    _             ### 
        ###          | \  / | |  | | |__ /  \ _| |_           ### 
        ###          | |\/| | |  | |  __/ /\ \_   _|          ###
        ###          | |  | | |__| | | / ____ \|_|            ###
        ###          |_|  |_|\____/|_|/_/    \_\              ###
        ###                                                   ### 
        ######################################################### 
         


Scaling views to unit variance...

Scaling groups to unit variance...

Loaded view='data' group='group1' with N=113 samples and D=11937 features...


Model options:
- Automatic Relevance Determination prior on the factors: True
- Automatic Relevance Determination prior on the weights: True
- Spike-and-slab prior on the factors: True
- Spike-and-slab prior on the weights: True
Likelihoods:
- View 0 (data): gaussian




######################################
## Training the model with seed 1 ##
######################################



Converged!



#######################
## Training finished ##
#######################


Warning: Output file mofa-methylation.model.hdf5 already exists, it will be replaced
Saving model in mofa-methylation.model.hdf5...
Saved MOFA embeddings in .obsm['X_mofa'] slot and their loadings in .varm['LFs'].

We compute the silhouette width on the top factors / principal components of each embedding. The silhouette width is not well-behaved for very small groups, so we restrict the comparison to histological subtypes with more than 5 samples. This drops LCC (n=4) and SCLC (n=2), retaining AC, SqCC, and LCNEC. The results are qualitatively similar for other cutoffs.

# Minimum number of samples a histological subtype must have to be compared (rationale above)
MIN_GROUP_SIZE = 5

# Define embeddings per model
embeddings = {
    "Methylation (PCA)": (methylation.obsm["X_pca"], methylation.obs["who2015_histology"]),
    "RNA (PCA)": (
        rna.obsm["X_pca"],
        rna.obs["who2015_histology"],
    ),
    "Proteins (PCA)": (adata.obsm["X_pca_obs"], adata.obs["who2015_histology"]),
    "Methylation (MOFA+)": (methylation.obsm["X_mofa"], methylation.obs["who2015_histology"]),
    "RNA (MOFA+)": (
        rna.obsm["X_mofa"],
        rna.obs["who2015_histology"],
    ),
    "Proteins (MOFA+)": (adata.obsm["X_mofa"], adata.obs["who2015_histology"]),
    "Joint (MOFA+)": (mdata.obsm["X_mofa"], mdata.obs["who2015_histology"]),
}

# To make modalities with partially matched samples comparable, we subset to the intersection of samples
intersecting_samples = set.intersection(
    set(adata.obs_names.tolist()), set(rna.obs_names.tolist()), set(methylation.obs_names.tolist())
)

# We subset to groups with a sufficient number of samples
group_sizes = metadata.value_counts("who2015_histology")
kept_groups = group_sizes.index[group_sizes > MIN_GROUP_SIZE]

records = []
for method, (embedding, labels) in embeddings.items():
    keep = (labels.index.isin(intersecting_samples) & labels.isin(kept_groups)).tolist()
    kept_embedding, kept_labels = embedding[keep], labels[keep]

    records.append(
        {
            "method": method,
            "silhouette_width": silhouette_score(kept_embedding, labels=kept_labels, metric="euclidean"),
            "n_groups": kept_labels.nunique(),
        }
    )

silhouette_scores = pd.DataFrame.from_records(records)
silhouette_scores

	method	silhouette_width	n_groups
0	Methylation (PCA)	0.079848	3
1	RNA (PCA)	0.150960	3
2	Proteins (PCA)	0.115391	3
3	Methylation (MOFA+)	0.084754	3
4	RNA (MOFA+)	0.137201	3
5	Proteins (MOFA+)	0.113013	3
6	Joint (MOFA+)	0.154716	3

plot_kwargs = {
    "basis": "X_mofa",
    "color": "who2015_histology",
    "size": 200,
    "palette": apt.pl.BasePalettes().get("qualitative"),
    "show": False,
    "title": "",
}
modalities = {"methylation": methylation.copy(), "rna": rna.copy(), "proteins": adata.copy(), "joint": mdata.copy()}

fig, axm = apt.pl.create_figure(1, 5, figsize=(19, 4))

for idx, (modality, data) in enumerate(modalities.items()):
    muon.pl.embedding(data, ax=axm[idx], **plot_kwargs)
    apt.pl.label_axes(ax=axm[idx], title=modality)


apt.pl.barplot(
    ax=axm[4],
    data=silhouette_scores,
    value_column="silhouette_width",
    grouping_column="method",
    color_dict={
        "Methylation (PCA)": "#eeeeee",
        "RNA (PCA)": "#eeeeee",
        "Proteins (PCA)": "#eeeeee",
        "Methylation (MOFA+)": apt.pl.BaseColors.get("blue"),
        "RNA (MOFA+)": apt.pl.BaseColors.get("blue"),
        "Proteins (MOFA+)": apt.pl.BaseColors.get("blue"),
        "Joint (MOFA+)": apt.pl.BaseColors.get("green"),
    },
)
apt.pl.label_axes(
    ax=axm[4],
    ylabel="Silhouette Width\n($\\longrightarrow$Stronger Separation)",
    title="Subtype separation",
)
plt.xticks(rotation=45, ha="right")
plt.show()

The silhouette width analysis suggests that the joint modeling of the three available modalities improves the separation of histological subtypes moderately compared to the individual modalities.

That the integration of the multimodal data does not necessarily lead to a stronger separation of histological subtypes can have different reasons:

1. Histology is a strong phenotype that is already well separated by each individual modality, so the separation might “top out” at a certain value, independent of the modality.

2. As the original authors noted, there is additional molecular heterogeneity beyond the histological phenotype that cannot be resolved by the histological assessment, i.e. the labels we use do not recapitulate the full latent structure of the cohort.

3. Summary

In this tutorial, we integrated MS-proteomics data with 2 complementary modalities. We jointly modeled them with a multi-factor model in the muon framework and biologically interpreted the resulting factors using decoupler. By representing MS-proteomics data as anndata objects, alphapepttools makes them directly interoperable with the scverse ecosystem and thus with complementary modalities.

Session Info

from session_info2 import session_info

session_info()