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

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(glomero.dds))
glomero.dds = 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)

    glomero.dds = 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(glomero.dds), 1, gm_mean)
glomero.dds = estimateSizeFactors(glomero.dds, geoMeans = geoMeans, type = "poscount") # to be honest, I haven't figured out yet what all these lines with gmMeans do

glomero.dds = DESeq(glomero.dds, 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(glomero.dds)
plotPCA(deseq_counts_vst.glom)
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(glomero.dds, 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(glomero.dds, 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:

matrix(resultsNames(glomero.dds))

[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(glomero.dds, 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
ADD COMMENT
0
Entering edit mode
@mikelove
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.

ADD COMMENT

Login before adding your answer.

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

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

Powered by the version 2.3.6