Hi,
I am analyzing microRNAs across different taxonomic levels: Suborder, family, and species. I found that many microRNAs have low expression levels in certain groups, and I am trying to correct with shrinkage using lfcShrink. This was straightforward with Suborder, only two levels, I did the Wald test followed by lfcShrink, and as suggested I used LRT for the other level, family: 3 levels, species: 4, and then I did lfcShrink.
I am interested in doing pairwise tests (Wald) with DESEq2 and use lfcShrink, but I get an error. I created the dataset and ran DESeq().
miRNA.cet.Jan.15 <- DESeqDataSetFromMatrix(countData = subset.miRNA.cetacean.Jan15,
colData = subset.metadata,
design = ~ sp)
miRNA.cet.Jan.15<- DESeq(miRNA.cet.Jan.15)
I read that the approach is: 1. Define contrasts, 2. extract results table, and then 3. shrink the log2 fold changes. And then I get an error:
contrast<- c("sp","Mnov","Ddel")
res.MnovDde<- results(miRNA.cet.Jan.15, contrast=contrast)
res.MnovDde<- lfcShrink(miRNA.cet.Jan.15, contrast=contrast, res=res.MnovDde)
Error in lfcShrink(miRNA.cet.Jan.15, contrast = contrast, res = res.MnovDde) :
type='apeglm' shrinkage only for use with 'coef'
According to the DESEq2 guide, one"can use coef with either the name or numeric index from resultsNames(dds) to specify which coefficient to shrink". Then I ran it and got:
resultsNames(miRNA.cet.Jan.15)
[1] "Intercept" "sp_Mnov_vs_Ddel" "sp_Ttru_vs_Ddel" "sp_Zcav_vs_Ddel"
And there are only three comparison, and my issue is that I have four species.
levels(miRNA.cet.Jan.15$sp)
[1] "Ddel" "Mnov" "Ttru" "Zcav"
Does DESeq only do 3 comparisons? I read that it sets the first level as reference (alphabetic order), and then compare the others against it. In my case, I have four sps: Ttru, Ddel, Zcav, Mnov. I want to compare Zcav vs all others, and also Mnov_vs_Ddel. As shown, there are only 2 that I am interested: Mnov_vs_Ddel, and Zcav_vs_Ddel, but the others are missing.
I tried it as here: https://www.biostars.org/p/325009/#325036 and it failed:
MnoDdel.cet<- results(miRNA.cet.Jan.15,contrast = c("sp","Mnov","Ddel"))
MnoDdel.cet<- lfcShrink(miRNA.cet.Jan.15,contrast = c("sp","Mnov","Ddel"),res=MnoDdel.cet)
Error in lfcShrink(miRNA.cet.Jan.15, contrast = c("sp", "Mnov", "Ddel"), :
type='apeglm' shrinkage only for use with 'coef'
Is there a way to do all the pairwise comparisons with lfcShrink without having a reference, or do I have to run it with different references to get all the comparisons I want?
sessionInfo( )
R version 4.3.2 (2023-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 11 x64 (build 26100)
Matrix products: default
locale:
[1] LC_COLLATE=English_United States.utf8 LC_CTYPE=English_United States.utf8
[3] LC_MONETARY=English_United States.utf8 LC_NUMERIC=C
[5] LC_TIME=English_United States.utf8
time zone: America/New_York
tzcode source: internal
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods base
other attached packages:
[1] DESeq2_1.42.1 SummarizedExperiment_1.32.0 Biobase_2.62.0
[4] MatrixGenerics_1.14.0 matrixStats_1.5.0 GenomicRanges_1.54.1
[7] GenomeInfoDb_1.38.8 IRanges_2.36.0 S4Vectors_0.40.2
[10] BiocGenerics_0.48.1
loaded via a namespace (and not attached):
[1] gtable_0.3.6 ggplot2_4.0.0 lattice_0.22-7 numDeriv_2016.8-1.1
[5] vctrs_0.6.5 tools_4.3.2 bitops_1.0-9 generics_0.1.4
[9] parallel_4.3.2 tibble_3.2.1 pkgconfig_2.0.3 Matrix_1.6-4
[13] SQUAREM_2021.1 RColorBrewer_1.1-3 S7_0.2.0 lifecycle_1.0.4
[17] GenomeInfoDbData_1.2.11 truncnorm_1.0-9 compiler_4.3.2 farver_2.1.2
[21] codetools_0.2-20 RCurl_1.98-1.17 pillar_1.11.1 crayon_1.5.3
[25] MASS_7.3-60 BiocParallel_1.36.0 DelayedArray_0.28.0 emdbook_1.3.14
[29] abind_1.4-8 tidyselect_1.2.1 locfit_1.5-9.12 bdsmatrix_1.3-7
[33] mvtnorm_1.3-3 dplyr_1.1.4 ashr_2.2-63 grid_4.3.2
[37] colorspace_2.1-2 cli_3.6.4 invgamma_1.2 SparseArray_1.2.4
[41] magrittr_2.0.4 S4Arrays_1.2.1 dichromat_2.0-0.1 scales_1.4.0
[45] XVector_0.42.0 coda_0.19-4.1 irlba_2.3.5.1 bbmle_1.0.25.1
[49] rlang_1.1.5 Rcpp_1.1.0 mixsqp_0.3-54 glue_1.8.0
[53] pkgload_1.4.1 apeglm_1.24.0 rstudioapi_0.17.1 R6_2.6.1
[57] plyr_1.8.9 zlibbioc_1.48.2
