Hi,

I would like to normalise my bulk data using DESeq2. Since the samples are normalised against an averaged reference, 50% of the genes must be constant across samples [1]. I was wondering why aren't the samples normalised against one of the samples. This would only require the weaker assumption that between any pair of samples (more specifically, the pairs including the reference), less than 50% of the genes are differential (see [1], section "Clustering to weaken the non-DE assumption"). This approach is followed for normalising single cells across clusters in the scran package and for normalising single cells across batches in the multiBatchNorm function of the batchelor package.

References: [1] Lun, Aaron TL, Karsten Bach, and John C. Marioni. "Pooling across cells to normalize single-cell RNA sequencing data with many zero counts." Genome biology 17.1 (2016): 75.

You can use whichever size factor estimation method you like in DESeq2:

sf <- ...
sizeFactors(dds) <- sf
dds <- DESeq(dds)

DESeq2 uses the median ratio method by default because it works well on most bulk datasets, where the samples are typically from the same tissue but the organisms or cells differ by treatment, etc.

It's not a requirement that 50% of the genes be constant for median ratio to work. For example, all of the genes could be DE and median ratio could still provide effectize size factor estimation. It's more about the scale and balance of DE. In this three line simulation all of the genes are DE, the second sample is sequenced at 2x depth, and median ratio correctly recovers the size factors:

set.seed(1)
cts <- matrix(rpois(800,1:400),ncol=2)
cts[,2] <- cts[,2] * 2^rnorm(400,0,.5) # DE
cts[,2] <- cts[,2] * 2 # seq depth
DESeq2::estimateSizeFactorsForMatrix(cts)

[1] 0.7123276 1.4038484

From the scran perspective, it's just easier to use an averaged reference, because the single-cell profiles are too variable to use for pairwise normalization on the sample level.

As for the 50% requirement: it is easy to provide an example of what happens just above that:

set.seed(0)

# 0% DE, no difference in the size factors is expected.
cts <- matrix(rpois(800,1:400),ncol=2)
DESeq2::estimateSizeFactorsForMatrix(cts)
## [1] 1 1

# 40% DE, begins to diverge but is still reasonable.
cts2 <- cts
cts2[1:160,1] <- cts2[1:160,1] * 10
DESeq2::estimateSizeFactorsForMatrix(cts2)
## [1] 1.0431312 0.9586522

# >50% DE, woah.
cts2 <- cts
cts2[1:201,1] <- cts[1:201,1]*100
DESeq2::estimateSizeFactorsForMatrix(cts2)
## [1] 5.5470020 0.1802776

As Mike says, it is possible to retain accuracy up to 100% DEGs if the number and effect size of the DEGs are the same in both directions. But this is a pretty strong assumption, especially for single-cell data. Besides, most normalization methods will do pretty well in such balanced cases - I could even use good ol' library size normalization, which has a theoretical breakdown point of 0% that can be reported as 100% with perfect balance.