DESeq2: diff gene expression between species using gene-specific normalization factor

[mra]

Hi.

I'm trying to do differential gene expression analysis between the orthologous genes of two species. Leaving aside the technical and biological complications that this kind of analysis carries, I would like to use DESeq2 to perform the computations.

To do so, I'm thinking in take advantage of the gene-specific normalization factors (NF_i,j) of DESeq2, which are used, alongside the common scale quantity q_i,j, to define the expected value of observed counts in gene i and sample j (mu_i,j):

mu_i,j = q_i,j * NF_i,j

Since in my case besides differences in the library size between samples, I also have to take care of the differences in gene length between species, I'm using a form of NF_i,j defined as:

NF_i,j = s_j * L_i,j

Where L_i,j is the length of gene i on condition j, and s_j has the same meaning that in regular analysis but is computed using length-normalized counts:

s_j = median_i { ( k_i,j / L_i,j ) / [ Prod_{v=1,...,m} k_i,v / L_i,j ]^(1/m) }

I would greatly appreciate any comments on the soundness (or lack of it) of this approach. Better ideas to do the analysis are more than welcome.

Thanks,
mra

[Michael Love]

hi,

This looks correct to me.

So you counted genes using different gene models for each species?

[mra]

Hi Michael.

Yes, I used species-specific annotations files.

[Michael Love]

Then I think your formulation looks correct.

[Michael Love]

This question came up a few more times, so we implemented an argument normMatrix to estimateSizeFactors() (DESeq2 >= 1.6). So normMatrix would be L in your case. The estimateSizeFactors() function will then correct for library size taking into consideration this matrix of multiplicative factors, as you have done above but additionally ensuring that the mean of normalized counts is on the same scale as the mean of unnormalized counts by dividing out the geometric mean. Then it stores a normalizationFactor matrix which combines both normMatrix and the library size correction. The code itself can be seen in DESeq2:::estimateNormFactors (an unexported function).