Skip to contents

Goal: demonstrate basic use of the datasets_presence_matrix array.

The presence matrix is a sparse array, indicating which features (var) were present in each dataset. The array has dimensions [n_datasets, n_var], and is stored in the SOMA Measurement varp collection. The first dimension is indexed by the soma_joinid in the census_datasets dataframe. The second is indexed by the soma_joinid in the var dataframe of the measurement.

census <- cellxgene.census::open_soma()
# Grab the experiment containing human data, and the measurement therein with RNA
human <- census$get("census_data")$get("homo_sapiens")
human_rna <- human$ms$get("RNA")

# The census-wide datasets
datasets_df <- as.data.frame(census$get("census_info")$get("datasets")$read())
print(datasets_df)
#> # A tibble: 522 × 8
#>    soma_joinid collection_id     colle…¹ colle…² datas…³ datas…⁴ datas…⁵ datas…⁶
#>          <int> <chr>             <chr>   <chr>   <chr>   <chr>   <chr>     <int>
#>  1           0 43d4bb39-21af-4d… Transc… 10.101… f512b8… Skin    f512b8…   68036
#>  2           1 d36ca85c-3e8b-44… A mole… 10.110… 90d4a6… Fallop… 90d4a6…   60574
#>  3           2 d36ca85c-3e8b-44… A mole… 10.110… d1207c… Ovary … d1207c…   26134
#>  4           3 2b02dff7-e427-4c… Single… 10.101… 36c867… Ileum   36c867…   32458
#>  5           4 e9eec7f5-8519-42… Humora… 10.101… 58b010… A scRN… 58b010…  130908
#>  6           5 a72afd53-ab92-45… Single… 10.103… 456e8b… Single… 456e8b…   44721
#>  7           6 e4c9ed14-e560-49… A mole… 10.103… d8da61… A mole… d8da61…  116313
#>  8           7 4796c91c-9d8f-46… MSK SP… 10.103… 97d923… MSK SP… 97d923…   24025
#>  9           8 4796c91c-9d8f-46… MSK SP… 10.103… e3a7e9… MSK SP… e3a7e9…  221315
#> 10           9 4796c91c-9d8f-46… MSK SP… 10.103… 0caede… MSK SP… 0caede…  166895
#> # … with 512 more rows, and abbreviated variable names ¹​collection_name,
#> #   ²​collection_doi, ³​dataset_id, ⁴​dataset_title, ⁵​dataset_h5ad_path,
#> #   ⁶​dataset_total_cell_count

For convenience, read the entire presence matrix (for Homo sapiens) into a Matrix::sparseMatrix. There is a convenience API providing this capability:

presence_matrix <- cellxgene.census::get_presence_matrix(census, "Homo sapiens", "RNA")
print(dim(presence_matrix))
#> [1]   522 60664

We also need the var dataframe, which is read into an R data frame for convenient manipulation:

var_df <- as.data.frame(human_rna$var$read())
print(var_df)
#> # A tibble: 60,664 × 4
#>    soma_joinid feature_id      feature_name   feature_length
#>          <int> <chr>           <chr>                   <int>
#>  1           0 ENSG00000238009 RP11-34P13.7             3726
#>  2           1 ENSG00000279457 WASH9P                   1397
#>  3           2 ENSG00000228463 AP006222.1               8224
#>  4           3 ENSG00000237094 RP4-669L17.4             6204
#>  5           4 ENSG00000230021 RP11-206L10.17           5495
#>  6           5 ENSG00000237491 LINC01409                8413
#>  7           6 ENSG00000177757 FAM87B                   1947
#>  8           7 ENSG00000225880 LINC00115                1317
#>  9           8 ENSG00000230368 FAM41C                   1971
#> 10           9 ENSG00000230699 RP11-54O7.1              3043
#> # … with 60,654 more rows

Is a feature present in a dataset?

Goal: test if a given feature is present in a given dataset.

Important: the (one-based) indexes in the sparse presence matrix correspond to the (zero-based) soma_joinid + 1. In other words:

  • the first dimension of the presence matrix is (one plus) the dataset’s soma_joinid as stored in the census_datasets dataframe.
  • the second dimension of the presence matrix is (one plus) the feature’s soma_joinid as stored in the var dataframe.
var_joinid <- var_df$soma_joinid[var_df$feature_id == "ENSG00000286096"]
dataset_joinid <- datasets_df$soma_joinid[datasets_df$dataset_id == "97a17473-e2b1-4f31-a544-44a60773e2dd"]
is_present <- presence_matrix[dataset_joinid + 1, var_joinid + 1]
cat(paste("Feature is", if (is_present) "present." else "not present."))
#> Feature is present.

What datasets contain a feature?

Goal: look up all datasets that have a feature_id present.

# Grab the feature's soma_joinid from the var dataframe
var_joinid <- var_df$soma_joinid[var_df$feature_id == "ENSG00000286096"]

# The presence matrix is indexed by the joinids of the dataset and var dataframes,
# so slice out the feature of interest by its joinid.
dataset_joinids <- datasets_df$soma_joinid[presence_matrix[, var_joinid + 1] != 0]

print(datasets_df[dataset_joinids + 1, ])
#> # A tibble: 24 × 8
#>    soma_joinid collection_id     colle…¹ colle…² datas…³ datas…⁴ datas…⁵ datas…⁶
#>          <int> <chr>             <chr>   <chr>   <chr>   <chr>   <chr>     <int>
#>  1          89 283d65eb-dd53-49… Transc… 10.110… 07b1d7… Dissec… 07b1d7…   27210
#>  2         102 283d65eb-dd53-49… Transc… 10.110… 7c1c3d… Dissec… 7c1c3d…   49512
#>  3         103 283d65eb-dd53-49… Transc… 10.110… 9372df… Dissec… 9372df…   33794
#>  4         131 283d65eb-dd53-49… Transc… 10.110… dd03ce… Dissec… dd03ce…   23732
#>  5         145 283d65eb-dd53-49… Transc… 10.110… 7a0a88… Superc… 7a0a88…   74979
#>  6         147 283d65eb-dd53-49… Transc… 10.110… d2b5ef… Dissec… d2b5ef…   36886
#>  7         151 283d65eb-dd53-49… Transc… 10.110… f8dda9… Dissec… f8dda9…   31899
#>  8         154 283d65eb-dd53-49… Transc… 10.110… 3a7f3a… Superc… 3a7f3a…  291833
#>  9         156 283d65eb-dd53-49… Transc… 10.110… bdb26a… Superc… bdb26a…  227671
#> 10         158 283d65eb-dd53-49… Transc… 10.110… 5e5ab9… Dissec… 5e5ab9…   32306
#> # … with 14 more rows, and abbreviated variable names ¹​collection_name,
#> #   ²​collection_doi, ³​dataset_id, ⁴​dataset_title, ⁵​dataset_h5ad_path,
#> #   ⁶​dataset_total_cell_count

What features are in a dataset?

Goal: lookup the features present in a given dataset.

This example also demonstrates the ability to do the query on multiple datasets.

# Slice the dataset(s) of interest, and get the joinid(s)
dataset_joinids <- datasets_df$soma_joinid[datasets_df$collection_id == "17481d16-ee44-49e5-bcf0-28c0780d8c4a"]

# Slice the presence matrix by the first dimension, i.e., by dataset
var_joinids <- var_df$soma_joinid[which(Matrix::colSums(presence_matrix[dataset_joinids + 1, ]) > 0)]

print(var_df[var_joinids + 1, ])
#> # A tibble: 27,211 × 4
#>    soma_joinid feature_id      feature_name   feature_length
#>          <int> <chr>           <chr>                   <int>
#>  1           0 ENSG00000238009 RP11-34P13.7             3726
#>  2           1 ENSG00000279457 WASH9P                   1397
#>  3           2 ENSG00000228463 AP006222.1               8224
#>  4           3 ENSG00000237094 RP4-669L17.4             6204
#>  5           4 ENSG00000230021 RP11-206L10.17           5495
#>  6           5 ENSG00000237491 LINC01409                8413
#>  7           6 ENSG00000177757 FAM87B                   1947
#>  8           7 ENSG00000225880 LINC00115                1317
#>  9           8 ENSG00000230368 FAM41C                   1971
#> 10           9 ENSG00000230699 RP11-54O7.1              3043
#> # … with 27,201 more rows