Dear Micheal,

Thank you for your answer, as I missed that feature while looking in the documentation. Here is the toy example I made using the normFactors.

For that, I calculated the sizeFactors, then I used them as an initial value that I multiplied by the inverse of my coefficients (to scale up), then, I normalized again the size factors with the geometric mean.

It behaves as I want to, i.e it gives me a statistical value for enrichment, that is modulated by my coefficients. It also corrects for library size. Does it sound statistically sound ? Thank you,

    ```{r}

# Load necessary library
library(DESeq2)
library(dplyr)

# Set seed for reproducibility
set.seed(123)

# Parameters
num_genes <- 1000
num_samples <- 10
conditions <- c("A", "B", "C")
base_count <- 100

# Generate count matrix with Poisson distribution and noise
counts <- matrix(rpois(num_genes * num_samples, lambda = base_count), nrow = num_genes, ncol = num_samples)

noise <- matrix(rnorm(num_genes * num_samples, mean = 0, sd = base_count * 0.05), nrow = num_genes, ncol = num_samples)

counts <- counts + noise
counts[counts < 0] <- 0  # Ensure no negative counts

# Round counts to the nearest integer
counts <- round(counts)


# Assign column names (sample IDs)
colnames(counts) <- paste0("sample", 1:num_samples)

# Generate coldata
coldata <- data.frame(
  sample_id = colnames(counts),
  condition = rep(conditions, length.out = num_samples)
)
rownames(coldata) <- coldata$sample_id


# Add a row to counts with expression 1 in 'A' , 'B' samples, but 1000 in 'C' samples to mimic a true "C" DE gene 
new_row <- rep(1, num_samples)
new_row[coldata$condition == "C"] <- 1000
counts <- rbind(counts, new_row)

# Double the counts for samples where condition is "B" to mimic change in library size 
counts[, coldata$condition == "B"] <- counts[, coldata$condition == "B"] * 2


# Create DESeq2 dataset
dds <- DESeqDataSetFromMatrix(countData = counts, colData = coldata, design = ~ condition)

dds <- estimateSizeFactors(dds)

# Initialize normFactors with size factors
normFactors <- matrix(sizeFactors(dds), nrow = num_genes+1, ncol = num_samples, byrow = TRUE)
print(head(normFactors))


#A has thrice the absolute expression. 
scaling_factor <- 3 
# Multiply columns of condition A by 1/scaling_factor
normFactors[, c(1, 4, 7, 10)] <- normFactors[, c(1, 4, 7, 10)] * (1/scaling_factor)


#Norm factors to 1 geometric mean 
normFactors <- normFactors / exp(rowMeans(log(normFactors)))

print(head(normFactors))

# Assign normFactors. 
normalizationFactors(dds) <- normFactors

#Rerun with these factors 
dds <- DESeq(dds)