I have single cell RNA seq data from two groups. When I use DESeq2 to explore differential regulation between the two groups, I visualize them in a plotMA graph. To my surprise this shows most genes are not around 0 for the y-axe (log2Fold) but -1. I added a link to the plot I made. What am I doing wrong?

https://imgur.com/a/C1ZRSAH

library(DESeq2)

# g1 = cell names group 1

# g2 = cell names group 2

# x = data frame with unique reads (raw) of all cells

x <- x[,c(g1,g2)]
x <- x[!grepl("ERCC",row.names(x))&!grepl("Rn45s", row.names(x)),]

des <- data.frame(row.names=colnames(x),
                  condition= c(rep("con1",length(g1)),
                               rep("con2",length(g2))))
dds <- DESeqDataSetFromMatrix(
  countData = round(x,0),
  colData = des,
  design = ~ condition)

dds <- DESeq(dds)
res <- results(dds)
plotMA(res, ylim = c(min(res$log2FoldChange, na.rm = TRUE), max(res$log2FoldChange, na.rm = TRUE)))

Why do you round the values (line countData = round(x,0))? DESeq2 works on the raw counts, which are integers. If you give it something else, then garbage-in / garbage-out.

Suppose you are feeding it actual read counts, then one possible explanation is that  the default DESeq "normalization" does not do the right thing for these data. You'll probably want to do more EDA to explore the data properties before fitting a complex statistical model. In the course of that EDA, you could try other normalization options, from within DESeq2 and beyond.

On a more general note, the raison d'être of DESeq2 is the shrinkage (of location and scale parameters), which is relevant for designed experiments with small sample sizes, such as your typical bulk RNA-Seq experiment. When comparing hundreds or thousands of cells, there is no need for shrinkage, and almost any old two-groups test should do.

Personally, I'd probably try (ZI)NB-WaVE.

Best wishes
      Wolfgang