Hi everyone,

I am trying to fit a multi factor comparison using EdgeR

I want to analyse expression data from samples collected from different populations at different latitudes (3 populations at both latitudes) which were subjected to a treatment.

I want to know the effect of latitude (lat), treatment (treat) and their interaction on the expression data.

My data looks like this:

sample    pop    lat    treat
1    1    N    1
2    1    N    1
3    1    N    2
4    1    N    2
5    2    N    1
6    2    N    1
7    2    N    2
8    2    N    2
9    3    N    1
10    3    N    1
11    3    N    2
12    3    N    2
13    4    S    1
14    4    S    1
15    4    S    2
16    4    S    2
17    5    S    1
18    5    S    1
19    5    S    2
20    5    S    2
21    6    S    1
22    6    S    1
23    6    S    2
24    6    S    2

When I want to fit a model using this design

design <- model.matrix(~lat+pop+treat+lat:treat)

the matrix is not of full rank because pop is nested in lat.

I have read in another post (DESeq2 paired multifactor test) and in the Edge R user guide (section 3.5) that this can be solved by creating another population variable (pop.nested) like this:

sample    pop    lat    treat    pop.nested
1    1    N    1    1
2    1    N    1    1
3    1    N    2    1
4    1    N    2    1
5    2    N    1    2
6    2    N    1    2
7    2    N    2    2
8    2    N    2    2
9    3    N    1    3
10    3    N    1    3
11    3    N    2    3
12    3    N    2    3
13    4    S    1    1
14    4    S    1    1
15    4    S    2    1
16    4    S    2    1
17    5    S    1    2
18    5    S    1    2
19    5    S    2    2
20    5    S    2    2
21    6    S    1    3
22    6    S    1    3
23    6    S    2    3
24    6    S    2    3

Yet, how does EdgeR handle this data? It seems statistically not sound to me as now population 1 and 4 are assumed to be same in the model?

Could someone comment on this?
Thanks a lot!

The best way to proceed would be something like this:

pop <- factor(rep(1:6, each=4))
lat <- factor(rep(c("N", "S"), each=12))
treatment <- factor(rep(rep(1:2, each=2), 6))
design <- model.matrix(~0 + pop + lat:treatment)
design <- design[,-grep("treatment1", colnames(design))]
colnames(design) <- sub(":", "_", colnames(design))

Now, looking at colnames(design), you get:

[1] "pop1"            "pop2"            "pop3"            "pop4"           
[5] "pop5"            "pop6"            "latN_treatment2" "latS_treatment2"

The first 6 coefficients represent blocking factors for the populations, while the last two represent latitude-specific treatment effects (of treatment 2 over treatment 1). You can drop them separately to test for the effect of treatment at either latitude; or you can compare them to each other, to test for differential treatment effects between latitudes:

con <- makeContrasts(latN_treatment2 - latS_treatment2, levels=design)

... or you can average them, to test for the average effect of treatment across both latitudes:

con <- makeContrasts((latN_treatment2 + latS_treatment2)/2, levels=design)

However, you cannot directly test for the effect of latitude because this is nested in population. Comparisons between populations are not sensible in this model as we've used all the population-level information to estimate the pop* coefficients. If you want to do this, you should consider using voom with duplicateCorrelation, which allows you to compare between levels of a blocking factor.

Thanks for your clear and fast answer!

Hi again,

I actually have the same question regarding DESeq2. How can I implement the above method to account for populations into the design formula of DESeq2? I can not find a way to edit the design matrix as you can not give a matrix as input in DESeq2. At least, I can't find a straightforward way to do that.

Any help is greatly appreciated.

Best,

Janne