Nested, multifactorial design formular for targeted contrast of a factor interaction
Entering edit mode
Misemi9 • 0
Last seen 10 days ago

Hi there, I'm trying to run a differential abundance analyses for my microbiome sequencing data from a multifactorial field trial which looks like this:

  • My data is replicated in 'blocks' (bl) that I used as random factor in mixed-effect models (I assume equivalent to 'batch effects' in RNA-seq).
  • Within these blocks, I have my nested treatment factors: fertiliser (ft) is nested within variety (vr) which is nested in crop protection (cp).
  • based on this experimental design, I recorded the following design formular:
mm = model.matrix(~bl+cp*vr*ft+cp:vr+cp:ft+vr:ft+cp:vr:ft, colData( = phyloseq_to_deseq2(ps.filter.glom, design = mm) #ps.filter.glom is a phyloseq object with an otu_table that contains 0s. Depending if I add a pseudo-count or not, I get different results (e.g. if I add 1, it won't run DESeq-function because no overdisersion) = phyloseq_to_deseq2(ps.filter.glom, design = mm)
gm_mean = function(x, na.rm=TRUE){
  exp(sum(log(x[x > 0]), na.rm=na.rm) / length(x))
geoMeans = apply(counts(, 1, gm_mean) = estimateSizeFactors(, geoMeans = geoMeans, type = "poscount") # to be honest, I haven't figured out yet what all these lines with gmMeans do = DESeq(, test="Wald", fitType="mean")

I ran a PERMANOVA with the adonis function (vegan package) based on a euclidean distance matrix:

deseq_counts_vst.glom <- varianceStabilizingTransformation(
vst_trans_count_tab_full.glom <- assay(deseq_counts_vst.glom)
samp.dist.glom <- dist(t(vst_trans_count_tab_full.glom), method = "euclidian")

adonis(data = meta_field, samp.dist.glom~cp*vr*ft, method = "bray",  strata = meta_field$bl) # meta_field is the metadata file from my phyloseq object, strata adds 'block' (bl) as random factor

Permanova results

and found out that I have a significant different microbial composition (p = 0.014) with the interaction of variety x crop protection, and specifically, in the comparison of the treatments CCPASZ vs. CCPSKY (CCP means conventional crop protection, ASZ and SKY are the two varieties) which I found out with pairwise.adonis()

pairwise.adonis2(samp.dist.glom ~ cpvr, data = meta_field, strata = 'bl')

Pairwise comparison results

Now my problem: the results-function in DESeq2 does not show me any differences that I found with adonis. But when I change the formular, I get different results. However, I cannot specify the contrasts in my result()-function for my cp x vr interaction. When I enter

res.A.B = results(, cooksCutoff = FALSE, contrast = list("cpvr","CCPASZ","CCPSKYZ")) # cpvr is the factor that I used in my meta-data file as a fused factor of cp and vr

I get 'Error in results(, cooksCutoff = FALSE, contrast = c("cpvr", : only list- and numeric-type contrasts are supported for user-supplied model matrices'

When I try to use the interaction from resultsNames(), I get a list with all the contrast that DESeq produces:


[1,] "Intercept"
[2,] "bl2"
[3,] "bl3"
[4,] "bl4"
[5,] "cpOCP"
[6,] "vrSKY"
[7,] "ftCM"
[8,] "ftMN"
[9,] "ftZE"
[10,] "cpOCP.vrSKY"
[11,] "cpOCP.ftCM"
[12,] "cpOCP.ftMN"
[13,] "cpOCP.ftZE"
[14,] "cpCCP.vrSKY.ftCM" [15,] "cpOCP.vrSKY.ftCM" [16,] "cpCCP.vrSKY.ftMN" [17,] "cpOCP.vrSKY.ftMN" [18,] "cpCCP.vrSKY.ftZE" [19,] "cpOCP.vrSKY.ftZE"

Number 10: [10,] "cpOCP.vrSKY" would be the right contrast of treatments (cp:vr), but somehow I cannot use this contrast, I either get the same error as before or other error messages.

sigtab.A.B = res.res.A.B = results(, cooksCutoff = FALSE, contrast = c("cpOCP.vrSKY","cpCCP.vrASZ","cpCCP.vrSKY"))

Where is my mistake? Perhpas there are simply no differences between treatment groups according to DESeq2, but I have the feeling that there is sth wrong with my design formular. And if someone could let me know how to use contrast in this specific case, I would be very pleased. I'm very sorry this is far aways from a reproducable example, I don't know how to share the whole phyloseq-object here, but I'm happy to do so!

Thanks a lot!

sessionInfo( )

sessionInfo() R version 4.1.0 (2021-05-18) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 10 x64 (build 19041)

attached base packages: 1 parallel stats4 stats graphics grDevices utils datasets methods base

other attached packages: 1 pairwiseAdonis_0.4 cluster_2.1.2 plyr_1.8.6
[4] RColorBrewer_1.1-2 seqinr_4.2-8 scales_1.1.1
[7] extrafont_0.17 vegan_2.5-7 lattice_0.20-44
[10] permute_0.9-5 dendextend_1.15.1 kableExtra_1.3.4
[13] biomformat_1.20.0 phangorn_2.7.0 ape_5.5
[16] DECIPHER_2.20.0 RSQLite_2.2.7 Biostrings_2.60.1
[19] XVector_0.32.0 mctoolsr_0.1.1.2 ggpubr_0.4.0
[22] decontam_1.12.0 DESeq2_1.32.0 SummarizedExperiment_1.22.0 [25] Biobase_2.52.0 MatrixGenerics_1.4.0 matrixStats_0.59.0
[28] GenomicRanges_1.44.0 GenomeInfoDb_1.28.0 IRanges_2.26.0
[31] S4Vectors_0.30.0 BiocGenerics_0.38.0 plotly_4.9.4.1
[34] forcats_0.5.1 stringr_1.4.0 dplyr_1.0.7
[37] purrr_0.3.4 readr_1.4.0 tidyr_1.1.3
[40] tibble_3.1.2 ggplot2_3.3.5 tidyverse_1.3.1
[43] phyloseq_1.36.0 devtools_2.4.2 usethis_2.0.1
[46] ANCOMBC_1.2.2 Maaslin2_1.6.0 dada2_1.20.0
[49] Rcpp_1.0.6

DESeq2 • 116 views
Entering edit mode
Last seen 9 hours ago
United States

Perhaps there are simply no differences between treatment groups according to DESeq2, but I have the feeling that there is something wrong with my design formula. And if someone could let me know how to use contrast in this specific case, I would be very pleased.

For statistical consultation on the appropriate design and contrasts for various hypotheses, I recommend contacting a local statistician at your institute to talk through these considerations.


Login before adding your answer.

Traffic: 339 users visited in the last hour
Help About
Access RSS

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6