Skip to contents

# Census datasets presence

Source: `vignettes/census_dataset_presence.Rmd`
`census_dataset_presence.Rmd`

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