# Getting started with MuData for MultiAssayExperiment

Danila Bredikhin1* and Ilia Kats2**

1European Molecular Biology Laboratory, Heidelberg, Germany
2German Cancer Research Center, Heidelberg, Germany

*danila.bredikhin@embl.de
**i.kats@dkfz-heidelberg.de

#### 2022-11-01

## 0.1 Introduction

Multimodal data format — MuData — has been introduced to address the need for cross-platform standard for sharing large-scale multimodal omics data. Importantly, it develops ideas of and is compatible with AnnData standard for storing raw and derived data for unimodal datasets.

In Bioconductor, multimodal datasets have been stored in MultiAssayExperiment (MAE) objects. This `MuData` package provides functionality to read data from MuData files into MAE objects as well as to save MAE objects into H5MU files.

## 0.2 Installation

The most recent dev build can be installed from GitHub:

    library(remotes)
    remotes::install_github("ilia-kats/MuData")

Stable version of `MuData` will be available in future bioconductor versions.

## 0.3 Loading libraries

    library(MuData)
    library(MultiAssayExperiment)

    library(rhdf5)

## 0.4 Writing H5MU files

We’ll use a simple MAE object from the `MultiAssayExperiment` package that we’ll then save in a H5MU file.

    data(miniACC)
    miniACC
    #> A MultiAssayExperiment object of 5 listed
    #>  experiments with user-defined names and respective classes.
    #>  Containing an ExperimentList class object of length 5:
    #>  [1] RNASeq2GeneNorm: SummarizedExperiment with 198 rows and 79 columns
    #>  [2] gistict: SummarizedExperiment with 198 rows and 90 columns
    #>  [3] RPPAArray: SummarizedExperiment with 33 rows and 46 columns
    #>  [4] Mutations: matrix with 97 rows and 90 columns
    #>  [5] miRNASeqGene: SummarizedExperiment with 471 rows and 80 columns
    #> Functionality:
    #>  experiments() - obtain the ExperimentList instance
    #>  colData() - the primary/phenotype DataFrame
    #>  sampleMap() - the sample coordination DataFrame
    #>  `$`, `[`, `[[` - extract colData columns, subset, or experiment
    #>  *Format() - convert into a long or wide DataFrame
    #>  assays() - convert ExperimentList to a SimpleList of matrices
    #>  exportClass() - save data to flat files

We will now write its contents into an H5MU file with `WriteH5MU`:

    writeH5MU(miniACC, "miniacc.h5mu")

## 0.5 Reading H5MU files

We can manually check the top level of the structure of the file:

    rhdf5::h5ls("miniacc.h5mu", recursive = FALSE)
    #>   group   name     otype dclass dim
    #> 0     /    mod H5I_GROUP
    #> 1     /    obs H5I_GROUP
    #> 2     /   obsm H5I_GROUP
    #> 3     / obsmap H5I_GROUP
    #> 4     /    uns H5I_GROUP
    #> 5     /    var H5I_GROUP

Or dig deeper into the file:

    h5 <- rhdf5::H5Fopen("miniacc.h5mu")
    h5&'mod'
    #> HDF5 GROUP
    #>         name /mod
    #>     filename
    #>
    #>              name     otype dclass dim
    #> 0 Mutations       H5I_GROUP
    #> 1 RNASeq2GeneNorm H5I_GROUP
    #> 2 RPPAArray       H5I_GROUP
    #> 3 gistict         H5I_GROUP
    #> 4 miRNASeqGene    H5I_GROUP
    rhdf5::H5close()

### 0.5.1 Creating MAE objects from H5MU files

This package provides `ReadH5MU` to create an object with data from an H5MU file. Since H5MU structure has been designed to accommodate more structured information than MAE, only some data will be read. For instance, MAE has no support for loading multimodal embeddings or pairwise graphs.

    acc <- readH5MU("miniacc.h5mu")
    #> Reading as SingleCellExperiment where the original object class is matrix
    acc
    #> A MultiAssayExperiment object of 5 listed
    #>  experiments with user-defined names and respective classes.
    #>  Containing an ExperimentList class object of length 5:
    #>  [1] RNASeq2GeneNorm: SummarizedExperiment with 198 rows and 79 columns
    #>  [2] gistict: SummarizedExperiment with 198 rows and 90 columns
    #>  [3] RPPAArray: SummarizedExperiment with 33 rows and 46 columns
    #>  [4] Mutations: SingleCellExperiment with 97 rows and 90 columns
    #>  [5] miRNASeqGene: SummarizedExperiment with 471 rows and 80 columns
    #> Functionality:
    #>  experiments() - obtain the ExperimentList instance
    #>  colData() - the primary/phenotype DataFrame
    #>  sampleMap() - the sample coordination DataFrame
    #>  `$`, `[`, `[[` - extract colData columns, subset, or experiment
    #>  *Format() - convert into a long or wide DataFrame
    #>  assays() - convert ExperimentList to a SimpleList of matrices
    #>  exportClass() - save data to flat files

Importantly, we recover the information from the original MAE object:

    head(colData(miniACC)[,1:4])
    #> DataFrame with 6 rows and 4 columns
    #>                 patientID years_to_birth vital_status days_to_death
    #>               <character>      <integer>    <integer>     <integer>
    #> TCGA-OR-A5J1 TCGA-OR-A5J1             58            1          1355
    #> TCGA-OR-A5J2 TCGA-OR-A5J2             44            1          1677
    #> TCGA-OR-A5J3 TCGA-OR-A5J3             23            0            NA
    #> TCGA-OR-A5J4 TCGA-OR-A5J4             23            1           423
    #> TCGA-OR-A5J5 TCGA-OR-A5J5             30            1           365
    #> TCGA-OR-A5J6 TCGA-OR-A5J6             29            0            NA
    head(colData(acc)[,1:4])
    #> DataFrame with 6 rows and 4 columns
    #>                 patientID years_to_birth vital_status days_to_death
    #>               <character>      <integer>    <integer>     <integer>
    #> TCGA-OR-A5J1 TCGA-OR-A5J1             58            1          1355
    #> TCGA-OR-A5J2 TCGA-OR-A5J2             44            1          1677
    #> TCGA-OR-A5J3 TCGA-OR-A5J3             23            0            NA
    #> TCGA-OR-A5J4 TCGA-OR-A5J4             23            1           423
    #> TCGA-OR-A5J5 TCGA-OR-A5J5             30            1           365
    #> TCGA-OR-A5J6 TCGA-OR-A5J6             29            0            NA

Features metadata is also recovered:

    head(rowData(miniACC[["gistict"]]))
    #> DataFrame with 6 rows and 3 columns
    #>        Gene.Symbol    Locus.ID    Cytoband
    #>        <character> <character> <character>
    #> DIRAS3      DIRAS3        9077      1p31.3
    #> MAPK14      MAPK14        1432     6p21.31
    #> YAP1          YAP1       10413     11q22.1
    #> CDKN1B      CDKN1B        1027     12p13.1
    #> ERBB2        ERBB2        2064       17q12
    #> G6PD          G6PD        2539        Xq28
    head(rowData(acc[["gistict"]]))
    #> DataFrame with 6 rows and 3 columns
    #>        Gene.Symbol    Locus.ID    Cytoband
    #>        <character> <character> <character>
    #> DIRAS3      DIRAS3        9077      1p31.3
    #> MAPK14      MAPK14        1432     6p21.31
    #> YAP1          YAP1       10413     11q22.1
    #> CDKN1B      CDKN1B        1027     12p13.1
    #> ERBB2        ERBB2        2064       17q12
    #> G6PD          G6PD        2539        Xq28

#### 0.5.1.1 Backed objects

It is possible to read H5MU files while keeping matrices (both `.X` and `.layers`) on disk.

    acc_b <- readH5MU("miniacc.h5mu", backed = TRUE)
    #> Reading as SingleCellExperiment where the original object class is matrix
    assay(acc_b, "RNASeq2GeneNorm")[1:5,1:3]
    #> <5 x 3> matrix of class DelayedMatrix and type "double":
    #>        TCGA-OR-A5J1-01A-11R-A29S-07 TCGA-OR-A5J2-01A-11R-A29S-07
    #> DIRAS3                    1487.0317                       9.6631
    #> MAPK14                     778.5783                    2823.6469
    #> YAP1                      1009.6061                    2305.0590
    #> CDKN1B                    2101.3449                    4248.9584
    #> ERBB2                      651.2968                     246.4098
    #>        TCGA-OR-A5J3-01A-11R-A29S-07
    #> DIRAS3                      18.9602
    #> MAPK14                    1061.7686
    #> YAP1                      1561.2502
    #> CDKN1B                    1348.5410
    #> ERBB2                       90.0607

The data in the assay is a `DelayedMatrix` object:

    class(assay(acc_b, "RNASeq2GeneNorm"))
    #> [1] "DelayedMatrix"
    #> attr(,"package")
    #> [1] "DelayedArray"

This is in contrast to the `acc` object that has matrices in memory:

    assay(acc, "RNASeq2GeneNorm")[1:5,1:3]
    #>        TCGA-OR-A5J1-01A-11R-A29S-07 TCGA-OR-A5J2-01A-11R-A29S-07
    #> DIRAS3                    1487.0317                       9.6631
    #> MAPK14                     778.5783                    2823.6469
    #> YAP1                      1009.6061                    2305.0590
    #> CDKN1B                    2101.3449                    4248.9584
    #> ERBB2                      651.2968                     246.4098
    #>        TCGA-OR-A5J3-01A-11R-A29S-07
    #> DIRAS3                      18.9602
    #> MAPK14                    1061.7686
    #> YAP1                      1561.2502
    #> CDKN1B                    1348.5410
    #> ERBB2                       90.0607
    class(assay(acc, "RNASeq2GeneNorm"))
    #> [1] "matrix" "array"

## 0.6 References

- Muon: multimodal omics analysis framework preprint

- mudata (Python) documentation

- muon documentation and web page

## 0.7 Session Info

    sessionInfo()
    #> R version 4.2.1 (2022-06-23 ucrt)
    #> Platform: x86_64-w64-mingw32/x64 (64-bit)
    #> Running under: Windows Server x64 (build 20348)
    #>
    #> Matrix products: default
    #>
    #> locale:
    #> [1] LC_COLLATE=C
    #> [2] LC_CTYPE=English_United States.utf8
    #> [3] LC_MONETARY=English_United States.utf8
    #> [4] LC_NUMERIC=C
    #> [5] LC_TIME=English_United States.utf8
    #>
    #> attached base packages:
    #> [1] stats4    stats     graphics  grDevices utils     datasets  methods
    #> [8] base
    #>
    #> other attached packages:
    #>  [1] SingleCellMultiModal_1.9.1  scater_1.26.0
    #>  [3] ggplot2_3.3.6               scuttle_1.8.0
    #>  [5] CiteFuse_1.10.0             MultiAssayExperiment_1.24.0
    #>  [7] SingleCellExperiment_1.20.0 SummarizedExperiment_1.28.0
    #>  [9] Biobase_2.58.0              GenomicRanges_1.50.0
    #> [11] GenomeInfoDb_1.34.0         IRanges_2.32.0
    #> [13] MatrixGenerics_1.10.0       matrixStats_0.62.0
    #> [15] MuData_1.2.0                rhdf5_2.42.0
    #> [17] S4Vectors_0.36.0            BiocGenerics_0.44.0
    #> [19] Matrix_1.5-1                BiocStyle_2.26.0
    #>
    #> loaded via a namespace (and not attached):
    #>   [1] utf8_1.2.2                    R.utils_2.12.1
    #>   [3] tidyselect_1.2.0              RSQLite_2.2.18
    #>   [5] AnnotationDbi_1.60.0          grid_4.2.1
    #>   [7] BiocParallel_1.32.0           Rtsne_0.16
    #>   [9] DropletUtils_1.18.0           munsell_0.5.0
    #>  [11] ScaledMatrix_1.6.0            codetools_0.2-18
    #>  [13] statmod_1.4.37                scran_1.26.0
    #>  [15] withr_2.5.0                   colorspace_2.0-3
    #>  [17] filelock_1.0.2                highr_0.9
    #>  [19] knitr_1.40                    labeling_0.4.2
    #>  [21] GenomeInfoDbData_1.2.9        polyclip_1.10-4
    #>  [23] bit64_4.0.5                   farver_2.1.1
    #>  [25] pheatmap_1.0.12               vctrs_0.5.0
    #>  [27] generics_0.1.3                xfun_0.34
    #>  [29] BiocFileCache_2.6.0           randomForest_4.7-1.1
    #>  [31] R6_2.5.1                      ggbeeswarm_0.6.0
    #>  [33] graphlayouts_0.8.3            rsvd_1.0.5
    #>  [35] locfit_1.5-9.6                bitops_1.0-7
    #>  [37] rhdf5filters_1.10.0           cachem_1.0.6
    #>  [39] DelayedArray_0.24.0           assertthat_0.2.1
    #>  [41] promises_1.2.0.1              scales_1.2.1
    #>  [43] ggraph_2.1.0                  beeswarm_0.4.0
    #>  [45] gtable_0.3.1                  beachmat_2.14.0
    #>  [47] tidygraph_1.2.2               rlang_1.0.6
    #>  [49] splines_4.2.1                 BiocManager_1.30.19
    #>  [51] yaml_2.3.6                    reshape2_1.4.4
    #>  [53] httpuv_1.6.6                  tools_4.2.1
    #>  [55] bookdown_0.29                 SpatialExperiment_1.8.0
    #>  [57] ellipsis_0.3.2                jquerylib_0.1.4
    #>  [59] RColorBrewer_1.1-3            ggridges_0.5.4
    #>  [61] Rcpp_1.0.9                    plyr_1.8.7
    #>  [63] sparseMatrixStats_1.10.0      zlibbioc_1.44.0
    #>  [65] purrr_0.3.5                   RCurl_1.98-1.9
    #>  [67] dbscan_1.1-11                 viridis_0.6.2
    #>  [69] cowplot_1.1.1                 ggrepel_0.9.1
    #>  [71] cluster_2.1.4                 magrittr_2.0.3
    #>  [73] magick_2.7.3                  mime_0.12
    #>  [75] evaluate_0.17                 xtable_1.8-4
    #>  [77] gridExtra_2.3                 compiler_4.2.1
    #>  [79] tibble_3.1.8                  crayon_1.5.2
    #>  [81] R.oo_1.25.0                   htmltools_0.5.3
    #>  [83] segmented_1.6-0               later_1.3.0
    #>  [85] propr_4.2.6                   tidyr_1.2.1
    #>  [87] DBI_1.1.3                     tweenr_2.0.2
    #>  [89] ExperimentHub_2.6.0           dbplyr_2.2.1
    #>  [91] MASS_7.3-58.1                 rappdirs_0.3.3
    #>  [93] cli_3.4.1                     R.methodsS3_1.8.2
    #>  [95] parallel_4.2.1                metapod_1.6.0
    #>  [97] igraph_1.3.5                  pkgconfig_2.0.3
    #>  [99] vipor_0.4.5                   bslib_0.4.0
    #> [101] dqrng_0.3.0                   XVector_0.38.0
    #> [103] stringr_1.4.1                 digest_0.6.30
    #> [105] Biostrings_2.66.0             rmarkdown_2.17
    #> [107] uwot_0.1.14                   edgeR_3.40.0
    #> [109] DelayedMatrixStats_1.20.0     curl_4.3.3
    #> [111] kernlab_0.9-31                shiny_1.7.3
    #> [113] rjson_0.2.21                  lifecycle_1.0.3
    #> [115] nlme_3.1-160                  jsonlite_1.8.3
    #> [117] Rhdf5lib_1.20.0               BiocNeighbors_1.16.0
    #> [119] viridisLite_0.4.1             limma_3.54.0
    #> [121] fansi_1.0.3                   pillar_1.8.1
    #> [123] lattice_0.20-45               KEGGREST_1.38.0
    #> [125] fastmap_1.1.0                 httr_1.4.4
    #> [127] survival_3.4-0                interactiveDisplayBase_1.36.0
    #> [129] glue_1.6.2                    png_0.1-7
    #> [131] bluster_1.8.0                 BiocVersion_3.16.0
    #> [133] bit_4.0.4                     BiocBaseUtils_1.0.0
    #> [135] HDF5Array_1.26.0              ggforce_0.4.1
    #> [137] stringi_1.7.8                 sass_0.4.2
    #> [139] mixtools_1.2.0                blob_1.2.3
    #> [141] BiocSingular_1.14.0           AnnotationHub_3.6.0
    #> [143] memoise_2.0.1                 dplyr_1.0.10
    #> [145] irlba_2.3.5.1