Skip to content

get.exp.stat() returns all-NA min.lfc when pairwise logFC contains partial NAs #128

Description

@Zzzweiguang

Hi,

I am trying to follow the official BayesPrism workflow for bulk RNA-seq deconvolution using a scRNA-seq reference.

My workflow is:

  1. Use raw scRNA-seq counts as reference
  2. Filter genes with cleanup.genes() using default parameters
  3. Filter protein-coding genes with select.gene.type()
  4. Run get.exp.stat()
  5. Run select.marker()

However, in my dataset, get.exp.stat() produces all-NA min.lfc, and then select.marker() returns 0 signature genes.

After debugging, I found that the NA values seem to appear during the pairwise marker calculation step, likely from scran::pairwiseTTests() / combineMarkers(). Some pairwise logFC.* columns become NA for comparisons involving specific cell types, especially Neutrophil and Plasma cell.

In my case, after gene filtering, some cells from these cell types have zero library size across the retained genes. For example:

  • Downsample 300 cells per cell type: no zero-library cells, get.exp.stat() succeeds
  • Downsample 400 cells per cell type: Neutrophil has 1 zero-library cell and Plasma cell has 2 zero-library cells, min.lfc becomes all NA
  • Full dataset: Neutrophil has 7 zero-library cells and Plasma cell has 76 zero-library cells after filtering

It looks like partial NA values in the pairwise logFC matrix are propagated when calculating min.lfc, possibly because row-wise min/max is calculated without na.rm = TRUE.

My questions are:

  1. Is this expected behavior?
  2. Should users remove cells with zero library size after gene filtering before running get.exp.stat()?
  3. Would it be appropriate for get.exp.stat() to ignore partial NA values when calculating min.lfc, while keeping genes as NA only when all pairwise comparisons are NA?
  4. Is there a recommended official workaround for this situation?

A simplified description of the issue is:

exp.stat <- get.exp.stat(
  input = sc.counts,
  cell.type.labels = cell.type.labels,
  cell.state.labels = cell.state.labels
)

# min.lfc becomes all NA
# select.marker(exp.stat, ...) returns 0 genes

I am happy to provide more debugging output, package versions, or a small reproducible example if helpful.

Thank you!

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions