Establishes some general functions relating to single-cell DGE

[dtm2451]

In a recent DS Working Group meeting, we discussed the utility of adding some standardized functions for performing DGE with various tools. We also laid out a few helper functions -- pseudobulking, gene filtering -- that felt required across tools.

This PR will directly include the helper functions and I'd propose that we use this sc-dge-functions-branch as the base branch that we'll PR all of our tool-specific DGE function builds in to!

Planned functions:

• pseudobulking
  • Initial implementation
    • Uses Seurat::AggregateExpression, grouping by a cell.by metadata and any number of sample.by metadata.
    • Adds a cell count to the metadata & allows trimming pseudobulks based on too few cells (10 by default)
    • Also pulls in additional metadata, as Seurat only grabs the 'group.by' metadata automatically.
    • Allows targeting only certain features or cell identities
    • Allows outputting as either Seurat or a list ("counts" matrix + "metadata" data.frame)
  • add dreamlet and scuttle implementations?
  • add normalization options?
• python pseudobulking
  • Initial implementation
    • Uses scanpy.get.aggregate
    • Allows MuData (Muon) or AnnData (scanpy) objects.
    • Allows outputting as either AnnData or a dict ("counts" array + "feature_names" list + "metadata" DataFrame)
    • Very similar functionality to the R version, but definitely faster!
• gene filtering
• DGE_comparison_assessment: for parsing groups and check if sufficient samples to trust results.

Planned Standardizations:

• function names:
  • run_<method> for the functions that directly run DGE on a set of samples
  • run_<method>_within_cells for the functions that loop across cell types, running the run_<method> within each cell type.
• input names:
  • counts = raw count matrix after feature selection
  • metadata = sample metadata
  • dge_by = column name of metadata holding sample groups
  • case_group = name of the case group to be used as numerator in log2FC calculation
  • reference_group = name of the reference group to be used as denominator in log2FC calculation
  • contrast = vector of contrasts, e.g., c(“varCase - varRef”, “varCase2 - varRef”, etc)
  • fixed_effects = names of metadata columns to be used as fixed effects
  • random_effects = names of metadata columns to be used as random effects
  • min_frac = threshold for the minimum number of samples with CPM > 1 expression to be used in selecting genes to retain for DGE
  • min_cells = threshold for the minimum number of cells that a sample should contain in order to be used.
  • ?cell_by = (?=unconfirmed) column name of metadata holding cell annotation or cluster identities. Only needed in functions that require this!
  • ?cell_targets = (?=unconfirmed) optional string vector holding cell annotation or cluster identities to focus.
  • dge_groups = (=used within internal functions only; allows simpler provision / merging info from multiple inputs above) vector of all groups of dge_by that should be retained for the analysis

Side Note:

I have often found it hard to consistently import python functions across distinct methods of running -- jupyter notebooks, interactive shell, running a script with python -u <script>. The method that has been working for me best recently is:

import sys
sys.path.insert(0, '/path/to/folder/ABOVE/essential_scripts/')
from essential_scripts.python_utils.ts_log import ts_log
from essential_scripts.single_cell.pseudobulk_function import dsco_pseudobulk

[erflynn]

this looks awesome!
in case you want the dreamlet version:

 sce = SingleCellExperiment(list(counts=object@assays$RNA@counts), colData=object@meta.data)
 pb <- aggregateToPseudoBulk(sce,
                              assay = "counts",
                              cluster_id = cell.by,
                              sample_id = sample.by, 
                              verbose = FALSE)

[erflynn]

I think it would be worthwhile to include a pre-pseudobulk filter -- e.g. only pseudobulk a sample/cell type pair if there are at least X cells of that cell type in that sample
I'm looking at the dreamlet::ProcessAssays() again and it does this as well, will make a note on my wrapper.

And possibly downstream a corresponding DEG filter that only pulls a DEG comparison if there are a least N samples per group?

[dtm2451]

I think it would be worthwhile to include a pre-pseudobulk filter -- e.g. only pseudobulk a sample/cell type pair if there are at least X cells of that cell type in that sample

This is included already, here for the R function! It runs after the pseudobulking currently, but could move it to before instead if there's good reason.

[erflynn]

oh awesome! apologies, I should have looked more carefully. I just realized it is not in the dreamlet pseudobulk function, but then is implemented in the processAssays, so was kind of making a note to myself

[erflynn]

I do think the downstream DEG filter to at least min.samples per category though is also useful

[dtm2451]

agreed! got pulled away before posting that half =)

[dtm2451]

And possibly downstream a corresponding DEG filter that only pulls a DEG comparison if there are a least N samples per group?

Hmm Agreed. Perhaps a function that assesses the requested DGE comps per the dge.group.by and related vars, and contrast setups for however we end up setting that up to work instead... Adding a ToDo for this, 'dge_comparison_assement function' 👍, but I'm feeling there are extra bits to scope out before I'd start filling in this one.