An introduction to the MoleculeExperiment Class
Bárbara Zita Peters Couto
School of Mathematics and Statistics, The University of Sydney, Sydney, NSW, 2006, AustraliaCharles Perkins Centre, The University of Sydney, Sydney, NSW, 2006, AustraliaSydney Precision Data Science Centre, The University of Sydney, Sydney, NSW, 2006, AustraliaNick Robertson
School of Mathematics and Statistics, The University of Sydney, Sydney, NSW, 2006, AustraliaCharles Perkins Centre, The University of Sydney, Sydney, NSW, 2006, AustraliaSydney Precision Data Science Centre, The University of Sydney, Sydney, NSW, 2006, AustraliaEllis Patrick
School of Mathematics and Statistics, The University of Sydney, Sydney, NSW, 2006, AustraliaCharles Perkins Centre, The University of Sydney, Sydney, NSW, 2006, AustraliaSydney Precision Data Science Centre, The University of Sydney, Sydney, NSW, 2006, AustraliaWestmead Institute for Medical Research, University of Sydney, AustraliaShila Ghazanfar
School of Mathematics and Statistics, The University of Sydney, Sydney, NSW, 2006, AustraliaCharles Perkins Centre, The University of Sydney, Sydney, NSW, 2006, AustraliaSydney Precision Data Science Centre, The University of Sydney, Sydney, NSW, 2006, Australia30 October 2023
Source:vignettes/MoleculeExperiment.Rmd
MoleculeExperiment.RmdMoleculeExperiment
The R package MoleculeExperiment contains functions to create and work with objects from the new MoleculeExperiment class. We introduce this class for analysing molecule-based spatial transcriptomics data (e.g., Xenium by 10X, CosMx SMI by Nanostring, and Merscope by Vizgen, among others).
Why the MoleculeExperiment class?
The goal of the MoleculeExperiment class is to: 1. Enable analysis of
spatial transcriptomics (ST) data at the molecule level, independent of
aggregating to the cell or tissue level. 2. Standardise molecule-based
ST data across vendors, to hopefully facilitate comparison of different
data sources and common analytical and visualisation workflows. 3.
Enable aggregation to a SpatialExperiment object given
combinations of molecules and segmentation boundaries.
Installation
The latest release of MoleculeExperiment can be installed using:
if (!require("BiocManager", quietly = TRUE)) {
install.packages("BiocManager")
}
BiocManager::install("MoleculeExperiment")Minimal example
- Load required libraries.
- Create MoleculeExperiment object with example Xenium data, taken over a small patch.
repoDir <- system.file("extdata", package = "MoleculeExperiment")
repoDir <- paste0(repoDir, "/xenium_V1_FF_Mouse_Brain")
me <- readXenium(repoDir, keepCols = "essential")
me
#> MoleculeExperiment class
#>
#> molecules slot (1): detected
#> - detected:
#> samples (2): sample1 sample2
#> -- sample1:
#> ---- features (137): 2010300C02Rik Acsbg1 ... Zfp536 Zfpm2
#> ---- molecules (962)
#> ---- location range: [4900,4919.98] x [6400.02,6420]
#> -- sample2:
#> ---- features (143): 2010300C02Rik Acsbg1 ... Zfp536 Zfpm2
#> ---- molecules (777)
#> ---- location range: [4900.01,4919.98] x [6400.16,6419.97]
#>
#>
#> boundaries slot (1): cell
#> - cell:
#> samples (2): sample1 sample2
#> -- sample1:
#> ---- segments (5): 67500 67512 67515 67521 67527
#> -- sample2:
#> ---- segments (9): 65043 65044 ... 65070 65071- Use standardised data in ME object for molecule-level analyses. For example, plot a simple digital in-situ, with cell boundaries overlaid.
ggplot_me() +
geom_polygon_me(me, assayName = "cell", fill = "grey") +
geom_point_me(me) +
# zoom in to selected patch area
coord_cartesian(
xlim = c(4900, 4919.98),
ylim = c(6400.02, 6420)
)
- Finally, it is also possible to go from a MoleculeExperiment object to a SpatialExperiment object. This enables the transition from a molecule-level analysis to a cell-level analysis with already existing tools.
# transform ME to SPE object
spe <- countMolecules(me)
spe
#> class: SpatialExperiment
#> dim: 178 14
#> metadata(0):
#> assays(1): counts
#> rownames(178): 2010300C02Rik Acsbg1 ... Zfp536 Zfpm2
#> rowData names(0):
#> colnames(14): sample1.67500 sample1.67512 ... sample2.65070
#> sample2.65071
#> colData names(4): sample_id cell_id x_location y_location
#> reducedDimNames(1): spatial
#> mainExpName: NULL
#> altExpNames(0):
#> spatialCoords names(2) : x_location y_location
#> imgData names(0):The ME object in detail
Constructing an ME object
Use case 1: from dataframes to ME object
Here we demonstrate how to work with an ME object from toy data, representing a scenario where both the detected transcripts information and the boundary information have already been read into R. This requires the standardisation of the data with the dataframeToMEList() function.
The flexibility of the arguments in dataframeToMEList() enable the creation of a standard ME object across dataframes comming from different vendors of molecule-based spatial transcriptomics technologies.
- Generate a toy transcripts data.frame.
moleculesDf <- data.frame(
sample_id = rep(c("sample1", "sample2"), times = c(30, 20)),
features = rep(c("gene1", "gene2"), times = c(20, 30)),
x_coords = runif(50),
y_coords = runif(50)
)
head(moleculesDf)
#> sample_id features x_coords y_coords
#> 1 sample1 gene1 0.6041390 0.74815289
#> 2 sample1 gene1 0.3198386 0.73074275
#> 3 sample1 gene1 0.9764606 0.12925932
#> 4 sample1 gene1 0.2047070 0.22801227
#> 5 sample1 gene1 0.8933710 0.78820939
#> 6 sample1 gene1 0.4270442 0.07544415- Generate a toy boundaries data.frame.
boundariesDf <- data.frame(
sample_id = rep(c("sample1", "sample2"), times = c(16, 6)),
cell_id = rep(
c(
"cell1", "cell2", "cell3", "cell4",
"cell1", "cell2"
),
times = c(4, 4, 4, 4, 3, 3)
),
vertex_x = c(
0, 0.5, 0.5, 0,
0.5, 1, 1, 0.5,
0, 0.5, 0.5, 0,
0.5, 1, 1, 0.5,
0, 1, 0, 0, 1, 1
),
vertex_y = c(
0, 0, 0.5, 0.5,
0, 0, 0.5, 0.5,
0.5, 0.5, 1, 1,
0.5, 0.5, 1, 1,
0, 1, 1, 0, 0, 1
)
)
head(boundariesDf)
#> sample_id cell_id vertex_x vertex_y
#> 1 sample1 cell1 0.0 0.0
#> 2 sample1 cell1 0.5 0.0
#> 3 sample1 cell1 0.5 0.5
#> 4 sample1 cell1 0.0 0.5
#> 5 sample1 cell2 0.5 0.0
#> 6 sample1 cell2 1.0 0.0- Standardise transcripts dataframe to ME list format.
moleculesMEList <- dataframeToMEList(moleculesDf,
dfType = "molecules",
assayName = "detected",
sampleCol = "sample_id",
factorCol = "features",
xCol = "x_coords",
yCol = "y_coords"
)
str(moleculesMEList, max.level = 3)
#> List of 1
#> $ detected:List of 2
#> ..$ sample1:List of 2
#> .. ..$ gene1: tibble [20 × 2] (S3: tbl_df/tbl/data.frame)
#> .. ..$ gene2: tibble [10 × 2] (S3: tbl_df/tbl/data.frame)
#> ..$ sample2:List of 1
#> .. ..$ gene2: tibble [20 × 2] (S3: tbl_df/tbl/data.frame)- Standardise boundaries dataframe to ME list format.
boundariesMEList <- dataframeToMEList(boundariesDf,
dfType = "boundaries",
assayName = "cell",
sampleCol = "sample_id",
factorCol = "cell_id",
xCol = "vertex_x",
yCol = "vertex_y"
)
str(boundariesMEList, 3)
#> List of 1
#> $ cell:List of 2
#> ..$ sample1:List of 4
#> .. ..$ cell1: tibble [4 × 2] (S3: tbl_df/tbl/data.frame)
#> .. ..$ cell2: tibble [4 × 2] (S3: tbl_df/tbl/data.frame)
#> .. ..$ cell3: tibble [4 × 2] (S3: tbl_df/tbl/data.frame)
#> .. ..$ cell4: tibble [4 × 2] (S3: tbl_df/tbl/data.frame)
#> ..$ sample2:List of 2
#> .. ..$ cell1: tibble [3 × 2] (S3: tbl_df/tbl/data.frame)
#> .. ..$ cell2: tibble [3 × 2] (S3: tbl_df/tbl/data.frame)- Create an ME object by using the MoleculeExperiment object constructor.
toyME <- MoleculeExperiment(
molecules = moleculesMEList,
boundaries = boundariesMEList
)
toyME
#> MoleculeExperiment class
#>
#> molecules slot (1): detected
#> - detected:
#> samples (2): sample1 sample2
#> -- sample1:
#> ---- features (2): gene1 gene2
#> ---- molecules (30)
#> ---- location range: [0.05,0.99] x [0.05,0.96]
#> -- sample2:
#> ---- features (1): gene2
#> ---- molecules (20)
#> ---- location range: [0.06,0.96] x [0.01,0.91]
#>
#>
#> boundaries slot (1): cell
#> - cell:
#> samples (2): sample1 sample2
#> -- sample1:
#> ---- segments (4): cell1 cell2 cell3 cell4
#> -- sample2:
#> ---- segments (2): cell1 cell2- Add boundaries from an external segmentation algorithm.
In this example, we use the extent of the molecules of generated for
toyME to allign the boundaries with the molecules. In
general, the extent of the segmentation is required for this
alignment.
repoDir <- system.file("extdata", package = "MoleculeExperiment")
segMask <- paste0(repoDir, "/BIDcell_segmask.tif")
boundaries(toyME, "BIDcell_segmentation") <- readSegMask(
# use the molecule extent to define the boundary extent
extent(toyME, assayName = "detected"),
path = segMask, assayName = "BIDcell_segmentation",
sample_id = "sample1", background_value = 0
)
toyME
#> MoleculeExperiment class
#>
#> molecules slot (1): detected
#> - detected:
#> samples (2): sample1 sample2
#> -- sample1:
#> ---- features (2): gene1 gene2
#> ---- molecules (30)
#> ---- location range: [0.05,0.99] x [0.05,0.96]
#> -- sample2:
#> ---- features (1): gene2
#> ---- molecules (20)
#> ---- location range: [0.06,0.96] x [0.01,0.91]
#>
#>
#> boundaries slot (2): cell BIDcell_segmentation
#> - cell:
#> samples (2): sample1 sample2
#> -- sample1:
#> ---- segments (4): cell1 cell2 cell3 cell4
#> -- sample2:
#> ---- segments (2): cell1 cell2
#> - BIDcell_segmentation:
#> samples (1): sample1
#> -- sample1:
#> ---- segments (89): 1 2 ... 88 89Use case 2: from machine’s output directory to ME object
The MoleculeExperiment package also provides functions to directly work with the directories containing output files of commonly used technologies. This is especially useful to work with data from multiple samples. Here we provide convenience functions to read in data from Xenium (10X Genomics), CosMx (Nanostring) and Merscope (Vizgen).
repoDir <- system.file("extdata", package = "MoleculeExperiment")
repoDir <- paste0(repoDir, "/xenium_V1_FF_Mouse_Brain")
me <- readXenium(repoDir, keepCols = "essential", addBoundaries = "cell")
me
#> MoleculeExperiment class
#>
#> molecules slot (1): detected
#> - detected:
#> samples (2): sample1 sample2
#> -- sample1:
#> ---- features (137): 2010300C02Rik Acsbg1 ... Zfp536 Zfpm2
#> ---- molecules (962)
#> ---- location range: [4900,4919.98] x [6400.02,6420]
#> -- sample2:
#> ---- features (143): 2010300C02Rik Acsbg1 ... Zfp536 Zfpm2
#> ---- molecules (777)
#> ---- location range: [4900.01,4919.98] x [6400.16,6419.97]
#>
#>
#> boundaries slot (1): cell
#> - cell:
#> samples (2): sample1 sample2
#> -- sample1:
#> ---- segments (5): 67500 67512 67515 67521 67527
#> -- sample2:
#> ---- segments (9): 65043 65044 ... 65070 65071readXenium() standardises the transcript and boundary information such that the column names are consistent across technologies when handling ME objects.
In addition, readXenium() enables the user to decide if they want to keep all data that is vendor-specific (e.g., column with qv score), some columns of interest, or only the essential columns. The latter refers to feature names and locations of the detected transcripts, and segment ids and boundary locations of the segmentation results.
For CosMx and Merscope data we provide convenience functions that standardise the raw transcripts data into a MoleculeExperiment object and additionally read the boundaries included in the standard data releases.
repoDir <- system.file("extdata", package = "MoleculeExperiment")
repoDir <- paste0(repoDir, "/nanostring_Lung9_Rep1")
meCosmx <- readCosmx(repoDir, keepCols = "essential", addBoundaries = "cell")
meCosmx
#> MoleculeExperiment class
#>
#> molecules slot (1): detected
#> - detected:
#> samples (2): sample_1 sample_2
#> -- sample_1:
#> ---- features (969): AATK ABL1 ... YES1 ZFP36
#> ---- molecules (23844)
#> ---- location range: [924.01,1031.94] x [26290,26398]
#> -- sample_2:
#> ---- features (943): AATK ABL1 ... YBX3 ZFP36
#> ---- molecules (7155)
#> ---- location range: [2894.03,3002] x [26290.05,26398]
#>
#>
#> boundaries slot (1): cell
#> - cell:
#> samples (2): sample_1 sample_2
#> -- sample_1:
#> ---- segments (113): 1 2 ... 112 113
#> -- sample_2:
#> ---- segments (83): 114 115 ... 195 196
repoDir <- system.file("extdata", package = "MoleculeExperiment")
repoDir <- paste0(repoDir, "/vizgen_HumanOvarianCancerPatient2Slice2")
meMerscope <- readMerscope(repoDir,
keepCols = "essential",
addBoundaries = "cell"
)
meMerscope
#> MoleculeExperiment class
#>
#> molecules slot (1): detected
#> - detected:
#> samples (1): vizgen_HumanOvarianCancerPatient2Slice2
#> -- vizgen_HumanOvarianCancerPatient2Slice2:
#> ---- features (486): ACKR3 ACTA2 ... ZBED2 ZEB1
#> ---- molecules (15160)
#> ---- location range: [10219.02,10386.83] x [8350.93,8395.3]
#>
#>
#> boundaries slot (1): cell
#> - cell:
#> samples (1): vizgen_HumanOvarianCancerPatient2Slice2
#> -- vizgen_HumanOvarianCancerPatient2Slice2:
#> ---- segments (24): 45862 45865 ... 54562 54564ME object structure
A MoleculeExperiment object contains a @molecules slot and an optional @boundaries slot.
Both slots have a hierarchical list structure that consists of a nested list, ultimately ending in a data.frame/tibble. Traditional rectangular data structures, like dataframes, redundantly store gene names and sample IDs for the millions of transcripts. In contrast, data in a list enables us to avoid this redundancy and work with objects of smaller size.
molecules slot
The @molecules slot contains molecule-level information. The essential data it contains is the feature name (e.g., gene names) and x and y locations of the detected molecules (e.g., transcripts), in each sample. Nevertheless, the user can also decide to keep all molecule metadata (e.g., subcellular location: nucleus/cytoplasm).
The nested list in the molecules slot has the following hierarchical structure: “assay name” > “sample ID” > “feature name” > dataframe/tibble with X and Y locations (and other additional columns of interest).
showMolecules(me)
#> List of 1
#> $ detected:List of 2
#> ..$ sample1:List of 137
#> .. ..$ 2010300C02Rik : tibble [11 × 2] (S3: tbl_df/tbl/data.frame)
#> .. ..$ Acsbg1 : tibble [6 × 2] (S3: tbl_df/tbl/data.frame)
#> .. .. [list output truncated]
#> ..$ sample2:List of 143
#> .. ..$ 2010300C02Rik: tibble [9 × 2] (S3: tbl_df/tbl/data.frame)
#> .. ..$ Acsbg1 : tibble [10 × 2] (S3: tbl_df/tbl/data.frame)
#> .. .. [list output truncated]boundaries slot
The @boundaries slot contains information from segmentation analyses (e.g., cell boundaries, or nucleus boundaries).
The nested list in the boundaries slot has the following hierarchical structure: “assay name” > “sample ID” > “segment ID” > dataframe/tibble with the vertex coordinates defining the boundaries for each segment. For example, if the boundary information is for cells, the assay name can be set to “cell”; or “nucleus” if one is using nucleus boundaries.
showBoundaries(me)
#> List of 1
#> $ cell:List of 2
#> ..$ sample1:List of 5
#> .. ..$ 67500: tibble [13 × 2] (S3: tbl_df/tbl/data.frame)
#> .. ..$ 67512: tibble [13 × 2] (S3: tbl_df/tbl/data.frame)
#> .. .. [list output truncated]
#> ..$ sample2:List of 9
#> .. ..$ 65043: tibble [13 × 2] (S3: tbl_df/tbl/data.frame)
#> .. ..$ 65044: tibble [13 × 2] (S3: tbl_df/tbl/data.frame)
#> .. .. [list output truncated]Methods
Here we introduce basic methods to access and manipulate data in an ME object, i.e., getters and setters, respectively.
Getters
The main getters are molecules() and boundaries(). NOTE: the output of these methods is the ME nested list, which can be very large on screen. Thus, these getters should be used when wanting to work with the data. To quickly view the slot contents, use showMolecules() and showBoundaries() instead.
# NOTE: output not shown as it is too large
# access molecules slot (note: assay name needs to be specified)
MoleculeExperiment::molecules(me, assayName = "detected")
# access cell boundary information in boundaries slot
boundaries(me, assayName = "cell")For ease of use, these getters have arguments that enable the transformation of the data from a nested ME list format to a data.frame format.
MoleculeExperiment::molecules(me, assayName = "detected", flatten = TRUE)
#> # A tibble: 1,739 × 4
#> x_location y_location feature_id sample_id
#> * <dbl> <dbl> <chr> <chr>
#> 1 4918. 6411. 2010300C02Rik sample1
#> 2 4901. 6417. 2010300C02Rik sample1
#> 3 4901. 6417. 2010300C02Rik sample1
#> 4 4910. 6417. 2010300C02Rik sample1
#> 5 4908. 6413. 2010300C02Rik sample1
#> 6 4911. 6407. 2010300C02Rik sample1
#> 7 4915. 6411. 2010300C02Rik sample1
#> 8 4916. 6412. 2010300C02Rik sample1
#> 9 4901. 6415. 2010300C02Rik sample1
#> 10 4906. 6417. 2010300C02Rik sample1
#> # ℹ 1,729 more rows
boundaries(me, assayName = "cell", flatten = TRUE)
#> # A tibble: 182 × 4
#> x_location y_location segment_id sample_id
#> * <dbl> <dbl> <chr> <chr>
#> 1 4905. 6400. 67500 sample1
#> 2 4899. 6401. 67500 sample1
#> 3 4894. 6408. 67500 sample1
#> 4 4890. 6418. 67500 sample1
#> 5 4887. 6423. 67500 sample1
#> 6 4887. 6425. 67500 sample1
#> 7 4890. 6427. 67500 sample1
#> 8 4891. 6427. 67500 sample1
#> 9 4894. 6426. 67500 sample1
#> 10 4908. 6421. 67500 sample1
#> # ℹ 172 more rowsOther getters include: features() and segmentIDs().
# get features in sample 1
features(me, assayName = "detected")[[1]]
#> [1] "2010300C02Rik" "Acsbg1" "Adamts2"
#> [4] "Adamtsl1" "Angpt1" "Aqp4"
#> [7] "Arc" "Arhgap12" "Arhgef28"
#> [10] "BLANK_0022" "BLANK_0414" "BLANK_0424"
#> [13] "Bdnf" "Bhlhe40" "Btbd11"
#> [16] "Cabp7" "Calb1" "Calb2"
#> [19] "Car4" "Cbln4" "Ccn2"
#> [22] "Cdh20" "Cdh4" "Cdh6"
#> [25] "Cdh9" "Chrm2" "Clmn"
#> [28] "Cntn6" "Cntnap5b" "Cplx3"
#> [31] "Cpne4" "Cpne6" "Crh"
#> [34] "Cux2" "Dkk3" "Dner"
#> [37] "Dpy19l1" "Epha4" "Fhod3"
#> [40] "Fos" "Gad1" "Galnt14"
#> [43] "Garnl3" "Gfra2" "Gjc3"
#> [46] "Gli3" "Gng12" "Gsg1l"
#> [49] "Gucy1a1" "Hat1" "Hpcal1"
#> [52] "Id2" "Igf2" "Igfbp4"
#> [55] "Igfbp5" "Igfbp6" "Igsf21"
#> [58] "Inpp4b" "Kcnh5" "Kctd12"
#> [61] "Lamp5" "Lypd6" "Lyz2"
#> [64] "Mapk4" "Meis2" "Myl4"
#> [67] "Myo16" "Ndst3" "Necab1"
#> [70] "NegControlProbe_00031" "NegControlProbe_00062" "Nell1"
#> [73] "Neto2" "Neurod6" "Npnt"
#> [76] "Nrep" "Nrn1" "Nrp2"
#> [79] "Ntsr2" "Nwd2" "Opalin"
#> [82] "Orai2" "Pde11a" "Pde7b"
#> [85] "Pdzd2" "Pdzrn3" "Penk"
#> [88] "Pip5k1b" "Plch1" "Plcxd2"
#> [91] "Plekha2" "Pou3f1" "Ppp1r1b"
#> [94] "Prdm8" "Prox1" "Prr16"
#> [97] "Pvalb" "Rasgrf2" "Rasl10a"
#> [100] "Rbp4" "Rfx4" "Rims3"
#> [103] "Rmst" "Ror1" "Rorb"
#> [106] "Rprm" "Rspo1" "Satb2"
#> [109] "Sema3a" "Sema6a" "Shisa6"
#> [112] "Sipa1l3" "Slc17a6" "Slc17a7"
#> [115] "Slc39a12" "Slit2" "Sncg"
#> [118] "Sorcs3" "Sox10" "Sox11"
#> [121] "Spi1" "Strip2" "Syndig1"
#> [124] "Syt2" "Tanc1" "Thsd7a"
#> [127] "Tle4" "Tmem132d" "Tox"
#> [130] "Trpc4" "Unc13c" "Vat1l"
#> [133] "Vip" "Vwc2l" "Wfs1"
#> [136] "Zfp536" "Zfpm2"
segmentIDs(me, "cell")
#> $sample1
#> [1] "67500" "67512" "67515" "67521" "67527"
#>
#> $sample2
#> [1] "65043" "65044" "65051" "65055" "65063" "65064" "65067" "65070" "65071"Setters
Main setters include molecules<- and
boundaries<-. For example, with
boundaries<- one can add new segmentation assay
information to the boundaries slot. Here we demonstrate this with the
nucleus boundaries.
repoDir <- system.file("extdata", package = "MoleculeExperiment")
repoDir <- paste0(repoDir, "/xenium_V1_FF_Mouse_Brain")
nucleiMEList <- readBoundaries(
dataDir = repoDir,
pattern = "nucleus_boundaries.csv",
segmentIDCol = "cell_id",
xCol = "vertex_x",
yCol = "vertex_y",
keepCols = "essential",
boundariesAssay = "nucleus",
scaleFactorVector = 1
)
boundaries(me, "nucleus") <- nucleiMEList
me # note the addition of the nucleus boundaries to the boundaries slot
#> MoleculeExperiment class
#>
#> molecules slot (1): detected
#> - detected:
#> samples (2): sample1 sample2
#> -- sample1:
#> ---- features (137): 2010300C02Rik Acsbg1 ... Zfp536 Zfpm2
#> ---- molecules (962)
#> ---- location range: [4900,4919.98] x [6400.02,6420]
#> -- sample2:
#> ---- features (143): 2010300C02Rik Acsbg1 ... Zfp536 Zfpm2
#> ---- molecules (777)
#> ---- location range: [4900.01,4919.98] x [6400.16,6419.97]
#>
#>
#> boundaries slot (2): cell nucleus
#> - cell:
#> samples (2): sample1 sample2
#> -- sample1:
#> ---- segments (5): 67500 67512 67515 67521 67527
#> -- sample2:
#> ---- segments (9): 65043 65044 ... 65070 65071
#> - nucleus:
#> samples (2): sample1 sample2
#> -- sample1:
#> ---- segments (5): 67500 67512 67515 67521 67527
#> -- sample2:
#> ---- segments (5): 65044 65051 65055 65063 65064The additional boundaries can be accessed, e.g. for visualisation.
ggplot_me() +
# add cell segments and colour by cell id
geom_polygon_me(me, byFill = "segment_id", colour = "black", alpha = 0.1) +
# add molecule points and colour by feature name
geom_point_me(me, byColour = "feature_id", size = 0.1) +
# add nuclei segments and colour the border with red
geom_polygon_me(me, assayName = "nucleus", fill = NA, colour = "red") +
# zoom in to selected patch area
coord_cartesian(xlim = c(4900, 4919.98), ylim = c(6400.02, 6420))
From MoleculeExperiment to SpatialExperiment
If one is interested in continuing downstream analysis at the cell-level, the MoleculeExperiment package also provides a convenience function, countMolecules(), that enables the transition from a MoleculeExperiment object to a SpatialExperiment object. With this functionality, it is possible to use already existing methods for cell-level data analysis.
spe <- countMolecules(me, boundariesAssay = "nucleus")
spe
#> class: SpatialExperiment
#> dim: 178 10
#> metadata(0):
#> assays(1): counts
#> rownames(178): 2010300C02Rik Acsbg1 ... Zfp536 Zfpm2
#> rowData names(0):
#> colnames(10): sample1.67500 sample1.67512 ... sample2.65063
#> sample2.65064
#> colData names(4): sample_id cell_id x_location y_location
#> reducedDimNames(1): spatial
#> mainExpName: NULL
#> altExpNames(0):
#> spatialCoords names(2) : x_location y_location
#> imgData names(0):Case Study: MoleculeExperiment and napari
Load the demonstration data, which includes molecules for 2 genes.
data(small_me)Read in virtual dissection CSV file, exported from napari (screenshot), of the morphology image.
bds_colours <- setNames(
c("#aa0000ff", "#ffaaffff"),
c("Region 1", "Region 2")
)
file_path <- system.file("extdata/tiny_brain_shape2.csv", package = "MoleculeExperiment")
bds_shape_raw <- read.csv(file = file_path, header = TRUE)
bds_shape_raw$sample_id <- "xenium_tiny_subset"
bds_shape_raw$regionName <- names(bds_colours)[bds_shape_raw$index + 1]
bds_shape <- dataframeToMEList(bds_shape_raw,
dfType = "boundaries",
assayName = "virtualDissection",
sampleCol = "sample_id",
factorCol = "regionName",
xCol = "axis.1",
yCol = "axis.0",
scaleFactor = 0.2125
)
boundaries(small_me, "virtualDissection") <- bds_shapeWe can plot the resulting MoleculeExperiment using the following code.
g <- ggplot() +
geom_point_me(
small_me,
assayName = "detected", byColour = "feature_id", size = 0.2
) +
geom_polygon_me(
small_me,
assayName = "cell", fill = NA, colour = "grey50", size = 0.1
) +
geom_polygon_me(
small_me,
assayName = "nucleus", fill = NA, colour = "black", size = 0.1
) +
geom_polygon_me(
small_me,
assayName = "virtualDissection", byFill = "segment_id", alpha = 0.3
) +
scale_y_reverse() +
theme_classic() +
theme(axis.text = element_blank()) +
theme(axis.ticks = element_blank()) +
coord_fixed() +
scale_fill_manual(values = bds_colours) +
NULL
g
Now that we have added the virtual dissection boundaries, we can use countMolecules to generate psuedobulk expressions over these regions.
spe <- countMolecules(
small_me, boundariesAssay = "virtualDissection")
speSessionInfo
sessionInfo()
#> R version 4.3.1 (2023-06-16)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: EndeavourOS
#>
#> Matrix products: default
#> BLAS: /usr/lib/libblas.so.3.11.0
#> LAPACK: /usr/lib/liblapack.so.3.11.0
#>
#> locale:
#> [1] LC_CTYPE=en_AU.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_AU.UTF-8 LC_COLLATE=en_AU.UTF-8
#> [5] LC_MONETARY=en_AU.UTF-8 LC_MESSAGES=en_AU.UTF-8
#> [7] LC_PAPER=en_AU.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_AU.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: Australia/Sydney
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] SpatialExperiment_1.12.0 GenomicRanges_1.54.0
#> [3] Matrix_1.6-1.1 SummarizedExperiment_1.32.0
#> [5] SingleCellExperiment_1.24.0 IRanges_2.36.0
#> [7] S4Vectors_0.40.1 ggplot2_3.4.4
#> [9] MoleculeExperiment_1.1.4 BiocStyle_2.30.0
#>
#> loaded via a namespace (and not attached):
#> [1] tidyselect_1.2.0 EBImage_4.44.0 farver_2.1.1
#> [4] dplyr_1.1.3 bitops_1.0-7 fastmap_1.1.1
#> [7] RCurl_1.98-1.12 digest_0.6.33 lifecycle_1.0.3
#> [10] terra_1.7-55 magrittr_2.0.3 compiler_4.3.1
#> [13] rlang_1.1.1 sass_0.4.7 tools_4.3.1
#> [16] utf8_1.2.4 yaml_2.3.7 data.table_1.14.8
#> [19] knitr_1.44 labeling_0.4.3 S4Arrays_1.2.0
#> [22] htmlwidgets_1.6.2 bit_4.0.5 DelayedArray_0.28.0
#> [25] abind_1.4-5 BiocParallel_1.36.0 withr_2.5.1
#> [28] purrr_1.0.2 BiocGenerics_0.48.0 desc_1.4.2
#> [31] grid_4.3.1 stats4_4.3.1 fansi_1.0.5
#> [34] colorspace_2.1-0 Rhdf5lib_1.24.0 scales_1.2.1
#> [37] cli_3.6.1 rmarkdown_2.25 crayon_1.5.2
#> [40] ragg_1.2.6 generics_0.1.3 rjson_0.2.21
#> [43] rhdf5_2.46.0 cachem_1.0.8 stringr_1.5.0
#> [46] zlibbioc_1.48.0 parallel_4.3.1 BiocManager_1.30.22
#> [49] XVector_0.42.0 tiff_0.1-11 matrixStats_1.0.0
#> [52] vctrs_0.6.4 jsonlite_1.8.7 bookdown_0.36
#> [55] fftwtools_0.9-11 bit64_4.0.5 systemfonts_1.0.5
#> [58] jpeg_0.1-10 magick_2.8.1 locfit_1.5-9.8
#> [61] jquerylib_0.1.4 glue_1.6.2 pkgdown_2.0.7
#> [64] codetools_0.2-19 stringi_1.7.12 gtable_0.3.4
#> [67] GenomeInfoDb_1.38.0 munsell_0.5.0 tibble_3.2.1
#> [70] pillar_1.9.0 rhdf5filters_1.14.0 htmltools_0.5.6.1
#> [73] GenomeInfoDbData_1.2.11 R6_2.5.1 textshaping_0.3.7
#> [76] rprojroot_2.0.3 evaluate_0.22 lattice_0.22-5
#> [79] Biobase_2.62.0 png_0.1-8 memoise_2.0.1
#> [82] bslib_0.5.1 Rcpp_1.0.11 SparseArray_1.2.0
#> [85] xfun_0.40 fs_1.6.3 MatrixGenerics_1.14.0
#> [88] pkgconfig_2.0.3