# CITE-seq data with CiteFuse and MuData

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

CITE-seq data provide RNA and surface protein counts for the same cells. This tutorial shows how MuData can be integrated into with Bioconductor workflows to analyse CITE-seq data.

## 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(SingleCellExperiment)
    library(MultiAssayExperiment)
    library(CiteFuse)
    library(scater)

    library(rhdf5)

## 0.4 Loading data

We will use CITE-seq data available within `CiteFuse` Bioconductor package.

    data("CITEseq_example", package = "CiteFuse")
    lapply(CITEseq_example, dim)
    #> $RNA
    #> [1] 19521   500
    #>
    #> $ADT
    #> [1]  49 500
    #>
    #> $HTO
    #> [1]   4 500

This dataset contains three matrices — one with `RNA` counts, one with antibody-derived tags (`ADT`) counts and one with hashtag oligonucleotide (`HTO`) counts.

## 0.5 Processing count matrices

While CITE-seq analysis workflows such as CiteFuse should be consulted for more details, below we exemplify simple data transformations in order to demonstrate how their output can be saved to an H5MU file later on.

Following the CiteFuse tutorial, we start with creating a `SingleCellExperiment` object with the three matrices:

    sce_citeseq <- preprocessing(CITEseq_example)
    sce_citeseq
    #> class: SingleCellExperiment
    #> dim: 19521 500
    #> metadata(0):
    #> assays(1): counts
    #> rownames(19521): hg19_AL627309.1 hg19_AL669831.5 ... hg19_MT-ND6
    #>   hg19_MT-CYB
    #> rowData names(0):
    #> colnames(500): AAGCCGCGTTGTCTTT GATCGCGGTTATCGGT ... TTGGCAACACTAGTAC
    #>   GCTGCGAGTTGTGGCC
    #> colData names(0):
    #> reducedDimNames(0):
    #> mainExpName: NULL
    #> altExpNames(2): ADT HTO

We will add a new assay with normalised RNA counts:

    sce_citeseq <- scater::logNormCounts(sce_citeseq)
    sce_citeseq  # new assay: logcounts
    #> class: SingleCellExperiment
    #> dim: 19521 500
    #> metadata(0):
    #> assays(2): counts logcounts
    #> rownames(19521): hg19_AL627309.1 hg19_AL669831.5 ... hg19_MT-ND6
    #>   hg19_MT-CYB
    #> rowData names(0):
    #> colnames(500): AAGCCGCGTTGTCTTT GATCGCGGTTATCGGT ... TTGGCAACACTAGTAC
    #>   GCTGCGAGTTGTGGCC
    #> colData names(1): sizeFactor
    #> reducedDimNames(0):
    #> mainExpName: NULL
    #> altExpNames(2): ADT HTO

To the ADT modality, we will add an assay with normalised counts:

    sce_citeseq <- CiteFuse::normaliseExprs(
      sce_citeseq, altExp_name = "ADT", transform = "log"
    )
    altExp(sce_citeseq, "ADT")  # new assay: logcounts
    #> class: SummarizedExperiment
    #> dim: 49 500
    #> metadata(0):
    #> assays(2): counts logcounts
    #> rownames(49): B220 (CD45R) B7-H1 (PD-L1) ... TCRb TCRg
    #> rowData names(0):
    #> colnames(500): AAGCCGCGTTGTCTTT GATCGCGGTTATCGGT ... TTGGCAACACTAGTAC
    #>   GCTGCGAGTTGTGGCC
    #> colData names(0):

We will also generate reduced dimensions:

    sce_citeseq <- scater::runPCA(
      sce_citeseq, exprs_values = "logcounts", ncomponents = 20
    )

    scater::plotReducedDim(sce_citeseq, dimred = "PCA",
                           by_exprs_values = "logcounts", colour_by = "CD27")

## 0.6 Making a MultiAssayExperiment object

An appropriate structure for multimodal datasets is `MultiAssayExperiment`.

We will make a respective MultiAssayExperiment object from `sce_citeseq`:

    experiments <- list(
      ADT = altExp(sce_citeseq, "ADT"),
      HTO = altExp(sce_citeseq, "HTO")
    )

    # Drop other modalities from sce_citeseq
    altExp(sce_citeseq) <- NULL
    experiments[["RNA"]] <- sce_citeseq

    mae <- MultiAssayExperiment(experiments)

## 0.7 Writing to H5MU

We can write the contents of the MultiAssayExperiment object into an H5MU file:

    writeH5MU(mae, "citefuse_example.h5mu")

We can check that all the modalities were written to the file:

    h5 <- rhdf5::H5Fopen("citefuse_example.h5mu")
    h5ls(H5Gopen(h5, "mod"), recursive = FALSE)
    #>   group name     otype dclass dim
    #> 0     /  ADT H5I_GROUP
    #> 1     /  HTO H5I_GROUP
    #> 2     /  RNA H5I_GROUP

… both assays for ADT — raw counts are stored in `X` and normalised counts are in the corresponding layer:

    h5ls(H5Gopen(h5, "mod/ADT"), FALSE)
    #>   group   name     otype dclass dim
    #> 0     /      X H5I_GROUP
    #> 1     / layers H5I_GROUP
    #> 2     /    obs H5I_GROUP
    #> 3     /    var H5I_GROUP
    h5ls(H5Gopen(h5, "mod/ADT/layers"), FALSE)
    #>   group      name       otype dclass      dim
    #> 0     / logcounts H5I_DATASET  FLOAT 49 x 500

… as well as reduced dimensions (PCA):

    h5ls(H5Gopen(h5, "mod/RNA/obsm"), FALSE)
    #>   group name       otype dclass      dim
    #> 0     /  PCA H5I_DATASET  FLOAT 20 x 500
    # There is an alternative way to access groups:
    # h5&'mod'&'RNA'&'obsm'
    rhdf5::H5close()

## 0.8 References

- Muon: multimodal omics analysis framework preprint

- mudata (Python) documentation

- muon documentation and web page

- Kim HJ, Lin Y, Geddes TA, Yang P, Yang JYH (2020). “CiteFuse enables multi-modal analysis of CITE-seq data.” Bioinformatics, 36(14), 4137–4143. https://doi.org/10.1093/bioinformatics/btaa282.

- Ramos M, Schiffer L, Re A, Azhar R, Basunia A, Cabrera CR, Chan T, Chapman P, Davis S, Gomez-Cabrero D, Culhane AC, Haibe-Kains B, Hansen K, Kodali H, Louis MS, Mer AS, Reister M, Morgan M, Carey V, Waldron L (2017). “Software For The Integration Of Multi-Omics Experiments In Bioconductor.” Cancer Research, 77(21); e39-42.

## 0.9 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] scater_1.26.0               ggplot2_3.3.6
    #>  [3] scuttle_1.8.0               CiteFuse_1.10.0
    #>  [5] MultiAssayExperiment_1.24.0 SingleCellExperiment_1.20.0
    #>  [7] SummarizedExperiment_1.28.0 Biobase_2.58.0
    #>  [9] GenomicRanges_1.50.0        GenomeInfoDb_1.34.0
    #> [11] IRanges_2.32.0              MatrixGenerics_1.10.0
    #> [13] matrixStats_0.62.0          MuData_1.2.0
    #> [15] rhdf5_2.42.0                S4Vectors_0.36.0
    #> [17] BiocGenerics_0.44.0         Matrix_1.5-1
    #> [19] BiocStyle_2.26.0
    #>
    #> loaded via a namespace (and not attached):
    #>   [1] ggbeeswarm_0.6.0          Rtsne_0.16
    #>   [3] colorspace_2.0-3          ggridges_0.5.4
    #>   [5] bluster_1.8.0             XVector_0.38.0
    #>   [7] BiocNeighbors_1.16.0      farver_2.1.1
    #>   [9] graphlayouts_0.8.3        ggrepel_0.9.1
    #>  [11] fansi_1.0.3               codetools_0.2-18
    #>  [13] splines_4.2.1             sparseMatrixStats_1.10.0
    #>  [15] cachem_1.0.6              knitr_1.40
    #>  [17] polyclip_1.10-4           jsonlite_1.8.3
    #>  [19] cluster_2.1.4             kernlab_0.9-31
    #>  [21] uwot_0.1.14               pheatmap_1.0.12
    #>  [23] ggforce_0.4.1             BiocManager_1.30.19
    #>  [25] compiler_4.2.1            dqrng_0.3.0
    #>  [27] assertthat_0.2.1          fastmap_1.1.0
    #>  [29] limma_3.54.0              cli_3.4.1
    #>  [31] tweenr_2.0.2              BiocSingular_1.14.0
    #>  [33] htmltools_0.5.3           tools_4.2.1
    #>  [35] rsvd_1.0.5                igraph_1.3.5
    #>  [37] gtable_0.3.1              glue_1.6.2
    #>  [39] GenomeInfoDbData_1.2.9    reshape2_1.4.4
    #>  [41] dplyr_1.0.10              Rcpp_1.0.9
    #>  [43] jquerylib_0.1.4           vctrs_0.5.0
    #>  [45] rhdf5filters_1.10.0       nlme_3.1-160
    #>  [47] DelayedMatrixStats_1.20.0 ggraph_2.1.0
    #>  [49] xfun_0.34                 stringr_1.4.1
    #>  [51] beachmat_2.14.0           lifecycle_1.0.3
    #>  [53] irlba_2.3.5.1             statmod_1.4.37
    #>  [55] edgeR_3.40.0              zlibbioc_1.44.0
    #>  [57] MASS_7.3-58.1             scales_1.2.1
    #>  [59] tidygraph_1.2.2           parallel_4.2.1
    #>  [61] RColorBrewer_1.1-3        yaml_2.3.6
    #>  [63] gridExtra_2.3             sass_0.4.2
    #>  [65] segmented_1.6-0           stringi_1.7.8
    #>  [67] highr_0.9                 ScaledMatrix_1.6.0
    #>  [69] randomForest_4.7-1.1      scran_1.26.0
    #>  [71] BiocParallel_1.32.0       rlang_1.0.6
    #>  [73] pkgconfig_2.0.3           bitops_1.0-7
    #>  [75] evaluate_0.17             lattice_0.20-45
    #>  [77] purrr_0.3.5               Rhdf5lib_1.20.0
    #>  [79] labeling_0.4.2            cowplot_1.1.1
    #>  [81] tidyselect_1.2.0          plyr_1.8.7
    #>  [83] magrittr_2.0.3            bookdown_0.29
    #>  [85] R6_2.5.1                  magick_2.7.3
    #>  [87] generics_0.1.3            metapod_1.6.0
    #>  [89] DelayedArray_0.24.0       DBI_1.1.3
    #>  [91] pillar_1.8.1              withr_2.5.0
    #>  [93] mixtools_1.2.0            survival_3.4-0
    #>  [95] RCurl_1.98-1.9            tibble_3.1.8
    #>  [97] utf8_1.2.2                rmarkdown_2.17
    #>  [99] viridis_0.6.2             locfit_1.5-9.6
    #> [101] grid_4.2.1                propr_4.2.6
    #> [103] digest_0.6.30             tidyr_1.2.1
    #> [105] dbscan_1.1-11             munsell_0.5.0
    #> [107] beeswarm_0.4.0            viridisLite_0.4.1
    #> [109] vipor_0.4.5               bslib_0.4.0