Cook's distance was implemented in DESeq2 for outlier detection. Specifically, the calculation is defined here.

calculateCooksDistance <- function(object, H, modelMatrix) {
  p <- ncol(modelMatrix)
  dispersions <- robustMethodOfMomentsDisp(object, modelMatrix)
  V <- assays(object)[["mu"]] + dispersions * assays(object)[["mu"]]^2
  PearsonResSq <- (counts(object) - assays(object)[["mu"]])^2 / V
  cooks <- PearsonResSq / p  * H / (1 - H)^2
  cooks
}

If I understood correctly, Cook's distance can be calculated using Pearson residuals as such (modified from here):

PearsonResSq / (dispersions * p) * H / (1 - H)^2

which is different from the DESeq2 implementation.

It further complicates things when the dispersion parameter in DESeq2 is defined in robustMethodOfMomentsDisp here as

alpha <- ( v - m ) / m^2

which if I understood correctly, returns the inverse of dispersion according to this source.

This in turns make the calculation of the variance correct, because normally it is expressed as such:

V <- assays(object)[["mu"]] + 1 / dispersions * assays(object)[["mu"]]^2

Therefore, I think the correct formula for calculating cooks should be:

cooks <- PearsonResSq * dispersions / p  * H / (1 - H)^2

However, I'm not the expert in statistics, and all my knowledge came from forum posts. I would love to see what others think. Thanks in advance!

The stackexchange post defines dispersion in a different way.

The Pearson residuals in the DESeq2 code are already weighted for variance V.

returns the inverse of dispersion

Also note that there are two definitions of dispersion, one that scales with variance and one that is the inverse.