Hi Michael Michael Love , thanks for your reply. We have already checked for other possible batch effects by looking at all the available variables. We are aware that there are some variables affecting gene expression (such as sex), but these were used for the matching, so by using the Pair variable they should have been already taken into account. In addition, if there were other batch effects that would cause this high dispersion I guess we would be able to see a similar high dispersion behavior with design = Batch + Case/Control. Here's the PCA with samples labeled as Case = 1 or Control = 0.

We also looked at plotCounts for a gene with high dispersion (left) and a gene with normal dispersion (right). I find it hard to see any difference between them.

We currently believe that the problem might be related to the number of pairs, but we don't have a full understanding of why/how this happens. We tried to reproduce our problem with simulated data. We created a dataset with 296 samples and ran DESeq. Here dispersion and volcano plots look normal.

If we add 2 more samples to the dataset, dispersion and volcano plots become weird and similar to what we obtained.

library(DESeq2)

# Generate data with 296 samples
#-------------------------------------------------------------------------------
nSamples <- 296
nGenes <- 500
set.seed(1234)
# Create count matrix
dds <- makeExampleDESeqDataSet(n=nGenes, m=nSamples)
# Create colData
colData <- data.frame(pair=as.factor(rep(1:(nSamples/2), each=2)),
                      condition=as.factor(rep(c(0,1), nSamples/2)))
# Create dds object
dds <- DESeqDataSetFromMatrix(assay(dds), colData, ~ pair+condition)
# Run DESeq
dds <- DESeq(dds)
res.dds <- results(dds)
# Dispersion plot
plotDispEsts(dds)
# Volcano plot
plot(res.dds$log2FoldChange, -1*log(res.dds$pvalue, 10))

# Both plots look "normal"
#-------------------------------------------------------------------------------

# Generate and add 2 more samples
#-------------------------------------------------------------------------------
set.seed(5678)
# Create count matrix for 2 other samples
assay.new <- assay(makeExampleDESeqDataSet(n=nGenes, m=2))
colnames(assay.new) <- c("sample297", "sample298")
# Create colData for other 2 samples and combine it with the previous one
colData.1 <- rbind(colData, data.frame(pair=as.factor(rep(nSamples/2+1 , 2)), condition=as.factor(c(0,1))))
# Create dds object
dds.1 <- DESeqDataSetFromMatrix(cbind(assay(dds), assay.new) , colData.1, ~ pair+condition)
# Run DESeq
dds.1 <- DESeq(dds.1)
# Results
res.dds.1 <- results(dds.1)
# Dispersion plot
plotDispEsts(dds.1)
# Volcano plot
plot(res.dds.1$log2FoldChange, -1*log(res.dds.1$pvalue, 10))

# Plots look weird
#-------------------------------------------------------------------------------