Question: Seeking confirmation on manually setting contrasts for a 3-factor design in DESeq2 with multiple interaction terms
1
3.3 years ago by
Australia
piensaglobalmente10 wrote:

Hello,

I am seeking confirmation on manually setting contrasts for a 3-factor design in DESeq2 with multiple interaction terms

I have a 3-factor experiment, with each factor containing 2 levels:

Type: Juvenile vs Sediment

Shore: Inshore vs Offshore

Sector: North vs Central

I am interested in the following contrasts:

1. Average main effects of Type, Shore and Sector across all the levels.
2. Interactive effect of Type:Shore, specifically:
1. Juvenile: Inshore vs Offshore
2. Sediment: Inshore vs Offshore
3. Interactive effect of Type: Sector, specifically:
1. Juvenile: North vs Central
2. Sediment: North vs Central
4. Interactive effect of Type vs Shore:Sector, specifically:
1. Juvenile: Inshore North vs. Inshore Central
2. Juvenile: Inshore North vs. Offshore Central
3. Juvenile: Inshore North vs. Offshore North
4. Juvenile: Offshore Central vs. Inshore Central
5. Juvenile: Offshore North vs. offshore Central
6. Juvenile: Offshore North vs. Inshore Central

5. 1-6: Sediment: with same 6 comparisons above

My experimental design formula is as follows:

~Type*Shore*Sector

Base levels for all factors are: Juvenile, Inshore, Central

And my Model Matrix column names outputted with resultsNames function:
resultsNames(varstab_geomLove_T.S.S)

[1] "Intercept"

[2] "Type_Sediment_vs_Juvenile"

[3] "Shore_Offshore_vs_Inshore"

[4] "Sector_North_vs_Central"

[5] "TypeSediment.ShoreOffshore"

[6]"TypeSediment.SectorNorth"

[7] "ShoreOffshore.SectorNorth"

[8]"TypeSediment.ShoreOffshore.SectorNorth"
 Effect Function Ok? Main effect of Type (aka Juvenile vs Sediments) results(varstab2, contrast=c("Type", "Juvenile", "Sediment")) yes Main effect of Shore (aka Inshore vs Offshore) results(varstab2, contrast=c("Shore", "Inshore", "Offshore")) yes Main effect of Sector (aka North or Central) results(varstab2, contrast=c("Sector", "North", "Central")) yes Juvenile: Inshore vs Offshore Subtract column with all inshore vs offshore comparison from sediment: inshore vs offshore   [3] "Shore_Offshore_vs_Inshore"     1,0,1,0,0,0,0,0 [5] "TypeSediment.ShoreOffshore"  1,0,0,0,1,0,0,0 Subtract the two to get the following: results(varstab2, contrast=c(0,0,1,0,-1,0,0,0)) Not sure Sediment: Inshore vs Offshore results(varstab2, name="TypeSediment.ShoreOffshore") yes Juvenile: North vs Central Subtract column with all north vs central comparison from sediment: north vs central   [4] "Sector_North_vs_Central"     1,0,0,1,0,0,0,0 [6]"TypeSediment.SectorNorth"  1,0,0,0,0,1,0,0 Subtract the two to get the following: results(varstab2, contrast=c(0,0,0,1,0,-1,0,0)) Not sure Sediment: North vs Central results(varstab2, name=”TypeSediment.SectorNorth”) yes Juvenile: Inshore North vs. Inshore Central [1] "Intercept"   1,0,0,0,0,0,0,0  which is juvenile, inshore, central [4] "Sector_North_vs_Central"   0,0,0,1,0,0,0,0   Subtract the two to get the following: results(varstab2, contrast=c(1,0,0,-1,0,0,0,0)) Not sure Juvenile: Inshore North vs. Offshore Central ? Juvenile: Inshore North vs. Offshore North ? Juvenile: Offshore Central vs. Inshore Central ? Juvenile: Offshore North vs. offshore Central [7] "ShoreOffshore.SectorNorth"    results(varstab2, name="ShoreOffshore.SectorNorth") Ok Juvenile: Offshore North vs. Inshore Central ? Sediments: Inshore North vs. Inshore Central ? Sediments: Inshore North vs. Offshore Central ? Sediments: Inshore North vs. Offshore North ? Sediments: Offshore Central vs. Inshore Central ? Sediments: Offshore North vs. offshore Central [8]"TypeSediment.ShoreOffshore.SectorNorth" results(varstab2, name="TypeSediment.ShoreOffshore.SectorNorth") Ok Sediments: Offshore North vs. Inshore Central ?

I feel confident in the comparisons with the a "Yes" in the last column. However, could someone please check that I have specified the contrasts in the rows with "Not sure?" Especially the contrasts "Juvenile: Inshore North vs. Inshore Central." If I can get a confirmation that that contrast is set correctly, then I think I can manage the rest.

Thank you for any insights!

cheers,

k

modified 3.3 years ago by Michael Love24k • written 3.3 years ago by piensaglobalmente10
1
3.3 years ago by
Michael Love24k
United States
Michael Love24k wrote:

hi,

This is a good question for a statistical consultant or a local statistician who you could collaborate with.

There is nothing specific about DESeq2 here, you would use the same contrasts for a normal linear model. So it's not really a DESeq2 software question.

Thank you Michael. I can keep looking through the glmm literature now that you confirm that contrasts in DESeq2 and other negative binomial models are specified the same. I just wasn't sure if there was some unique DESeq2 functionality I wasn't finding that would draw out these contrasts.

cheers,

Hello Michael:

I would like to propose a possible solution to such contrast questions and run it by you.

To estimate higher order contrasts:

1) Construct a design matrix in GLM coding. An example is given here:

2) Each categorical column (“A1”, “AB21”, etc) in the design can be seen as a separate binary factor. I can also pretend that every column is a quantitative factor (please let me know if it makes a difference for DESeq2). Either way, an additive model that includes all of such factors is fed to DESeq2. Since interactions are actually present, shrinkage of regression coefficients is disabled (betaPrior = FALSE).

3) Obtain the contrast vector as a difference between the two corresponding lsmeans vectors. E.g. for “(A1, B1) vs (A3, B1)” contrast the vector is LSM(AB11) – LSM(AB31). How to obtain lsmeans is described here:

https://support.sas.com/documentation/cdl/en/statug/63962/HTML/default/viewer.htm#statug_glm_a0000000871.htm

4) Feed the contrast vector to results() function from DESeq2. Here I hope to rely on the option: “a numeric contrast vector with one element for each element in resultsNames(object) (most general case)”.

Please tell me if DESeq2 can digest an overparameterized design “as is” and whether the results() will work the way I expect.

You can pass in any full rank design matrix to the 'full' argument, and use numeric contrasts to compare coefficients using results(). The current release of DESeq2 has betaPrior=FALSE by default (see NEWS).

If the user can only feed full rank designs to DESeq2, do you know if there is an R package that could help figure out how to specify contrast coefficients for non-trivial cases? If not, it's very hard to use the design produced by model.matrix() if the goal is to estimate non-trivial contrasts.

The contrast package works with 'lm' like objects. So it doesn't work with DESeq2 but you could extract the linear contrast from the output of the contrast() function from that package, for example running on fake data for 'y' but using the same model.matrix. See the example of the contrast package here.

For non-trivial contrasts I find it just as hard to specify the desired contrast using that package, as working it out myself.