FindMarkers with DESeq2 test not converting counts to integer after AggregateExpression pseudo-bulk analysis

[mblanche]

I'm trying to do pseudobulk analysis on GEO dataset from Lake et al.. I did many similar analysis without a problem. However, has a sanity check, following the Seurat tutorial work just fine:

Expected:

Seurat tutorial on the ifnb dataset

library(Seurat)
library(SeuratData)
ifnb <- LoadData("ifnb")

# load the inferred sample IDs of each cell
ctrl <- read.table(url("https://raw.githubusercontent.com/yelabucsf/demuxlet_paper_code/master/fig3/ye1.ctrl.8.10.sm.best"), head = T, stringsAsFactors = F)
stim <- read.table(url("https://raw.githubusercontent.com/yelabucsf/demuxlet_paper_code/master/fig3/ye2.stim.8.10.sm.best"), head = T, stringsAsFactors = F)
info <- rbind(ctrl, stim)

# rename the cell IDs by substituting the '-' into '.'
info$BARCODE <- gsub(pattern = "\\-", replacement = "\\.", info$BARCODE)

# only keep the cells with high-confidence sample ID
info <- info[grep(pattern = "SNG", x = info$BEST), ]

# remove cells with duplicated IDs in both ctrl and stim groups
info <- info[!duplicated(info$BARCODE) & !duplicated(info$BARCODE, fromLast = T), ]

# now add the sample IDs to ifnb
rownames(info) <- info$BARCODE
info <- info[, c("BEST"), drop = F]
names(info) <- c("donor_id")
ifnb <- AddMetaData(ifnb, metadata = info)

# remove cells without donor IDs
ifnb$donor_id[is.na(ifnb$donor_id)] <- "unknown"
ifnb <- subset(ifnb, subset = donor_id != "unknown")

# pseudobulk the counts based on donor-condition-celltype
pseudo_ifnb <- AggregateExpression(
  ifnb,
  assays = "RNA",
  return.seurat = T,
  group.by = c("stim", "donor_id", "seurat_annotations")
)

# each 'cell' is a donor-condition-celltype pseudobulk profile
pseudo_ifnb$celltype.stim <- paste(pseudo_ifnb$seurat_annotations, pseudo_ifnb$stim, sep = "_")
Idents(pseudo_ifnb) <- "celltype.stim"

## Do DESeq2 pseudobulk DGE
bulk.mono.de <- FindMarkers(
  object = pseudo_ifnb,
  ident.1 = "CD14 Mono_STIM",
  ident.2 = "CD14 Mono_CTRL",
  test.use = "DESeq2"
)
head(bulk.mono.de, n = 15)

When running on the tutorial data, FindMarkers() does the conversation to integer automatically

FindMarkers messages

r$> bulk.mono.de <- FindMarkers(
      object = pseudo_ifnb,
      ident.1 = "CD14 Mono_STIM",
      ident.2 = "CD14 Mono_CTRL",
      test.use = "DESeq2"
    )
converting counts to integer mode
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates

Problem:

However, doing the same on the Lake et al. dataset return the error Error in DESeqDataSet(se, design = design, ignoreRank) : some values in assay are not integers.

What am I missing? I did DESeq2 pseudobulk analysis in the past with other dataset from GEO with no issues...

library(Seurat)
library(SeuratDisk)
options(timeout = 60 * 60)

seuratURL <- "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE183276&format=file&file=GSE183276%5FKidney%5FHealthy%2DInjury%5FCell%5FAtlas%5FscCv3%5FSeurat%5F03282022%2Eh5Seurat"
seuratFile <- "GSE183276_Kidney_Healthy-Injury_Cell_Atlas_scCv3_Seurat_03282022.h5Seurat"
download.file(seuratURL, seuratFile)
seurat <- LoadH5Seurat(seuratFile)


bulk <- AggregateExpression(
  seurat,
  assays = "RNA",
  return.seurat = TRUE,
  group.by = c("condition.l1", "specimen", "subclass.l1")
)

bulk$cellType_disease <- factor(paste(bulk$subclass.l1, bulk$condition.l1, sep = "_"))
Idents(bulk) <- "cellType_disease"

de <- FindMarkers(
  object = bulk,
  ident.1 = "DCT_AKI",
  ident.2 = "DCT_Ref",
  test.use = "DESeq2"
)

I could always convert the counts fromt bulk object to integer manually then run DESeq2 manually but it seems like re-inventing what's already coded in FindMarkers()

Thanks for any guidance here

[ejscience]

I am going to guess you used SoupX to remove ambient RNA. On their vignette they remind us that "by default adjustCounts will return non-integer counts. If you require integers for downstream processing either pass out through round or set roundToInt=TRUE when running adjustCounts." So I wold just run

bulk@assays$RNA$counts <- round(bulk@assays$RNA$counts)

That fixed it for me.

[FionaMoon]

Hi,
The integration will also cause the counts not integers.

[samuel-marsh]

@FionaMoon integration does not alter the counts. In previous versions v3/4 an integrated count matrix was returned but it should not be used for analysis.