Entering edit mode
Hello,
Sorry, I've seen that this is a common problem, but for some reason none of the provided solutions has worked for me.
Basicaly, I am trying to relevel, as explained in the vignete. I have CV vs control, but I want the other way around.
When I get to lfcShrink, i get this error:
Error in lfcShrink(DE_control_CV, coef = resultsNames(DE_control_CV)[2], : 'coef' should specify same coefficient as in results 'res'
because resultsName only has "condition_CV_vs_C", so I cannot choose the correct.
Thanyou in advance for your help,
counts_control_CV <- counts[,c(2,8,10,16,19,23,5,9,12,18,22,24)] # subset the control + CV samples
counts_matrix_control_CV <- data.matrix(counts_control_CV) # matrix
counts_design_control_CV = data.frame(row.names = colnames(counts_control_CV),
condition=factor(unlist(c("C","C","C","C","C","C","CV","CV","CV","CV","CV","CV"))))
dds_control_CV <- DESeqDataSetFromMatrix(countData = counts_matrix_control_CV,
colData = counts_design_control_CV,
design = ~ condition)
dds_control_CV$condition <- relevel(dds_control_CV$condition,ref="C")
DE_control_CV <- DESeq(dds_control_CV) # DE object
vsd_control_CV <- vst(DE_control_CV) # variance stabilizing transformation
res_control_CV <- results(DE_control_CV, contrast = c("condition", "C", "CV")) # results
res_Shrink_control_CV <- lfcShrink(DE_control_CV, coef = resultsNames(DE_control_CV)[2], res = res_control_CV)
sessionInfo( )
R version 4.1.1 (2021-08-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19042)
Matrix products: default
locale:
[1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252 LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C LC_TIME=English_United States.1252
attached base packages:
[1] stats4 grid stats graphics grDevices utils datasets methods base
other attached packages:
[1] cowplot_1.1.1 EnhancedVolcano_1.12.0 ggrepel_0.9.1 DESeq2_1.34.0
[5] SummarizedExperiment_1.24.0 Biobase_2.54.0 MatrixGenerics_1.6.0 matrixStats_0.61.0
[9] GenomicRanges_1.46.0 GenomeInfoDb_1.30.0 IRanges_2.28.0 S4Vectors_0.32.0
[13] BiocGenerics_0.40.0 ComplexHeatmap_2.10.0 forcats_0.5.1 stringr_1.4.0
[17] purrr_0.3.4 readr_2.0.2 tidyr_1.1.4 tibble_3.1.5
[21] tidyverse_1.3.1 reshape2_1.4.4 dplyr_1.0.7 DT_0.19
[25] knitr_1.36 ggplot2_3.3.5
loaded via a namespace (and not attached):
[1] ggbeeswarm_0.6.0 colorspace_2.0-2 rjson_0.2.20 ellipsis_0.3.2 circlize_0.4.13 XVector_0.34.0
[7] GlobalOptions_0.1.2 fs_1.5.0 clue_0.3-60 rstudioapi_0.13 farver_2.1.0 bit64_4.0.5
[13] mvtnorm_1.1-3 apeglm_1.16.0 AnnotationDbi_1.56.1 fansi_0.5.0 lubridate_1.8.0 xml2_1.3.2
[19] codetools_0.2-18 splines_4.1.1 extrafont_0.17 doParallel_1.0.16 cachem_1.0.6 geneplotter_1.72.0
[25] jsonlite_1.7.2 Rttf2pt1_1.3.9 broom_0.7.9 annotate_1.72.0 cluster_2.1.2 dbplyr_2.1.1
[31] png_0.1-7 compiler_4.1.1 httr_1.4.2 backports_1.2.1 assertthat_0.2.1 Matrix_1.3-4
[37] fastmap_1.1.0 cli_3.0.1 htmltools_0.5.2 tools_4.1.1 coda_0.19-4 gtable_0.3.0
[43] glue_1.4.2 GenomeInfoDbData_1.2.7 maps_3.4.0 Rcpp_1.0.7 bbmle_1.0.24 cellranger_1.1.0
[49] vctrs_0.3.8 Biostrings_2.62.0 ggalt_0.4.0 extrafontdb_1.0 iterators_1.0.13 xfun_0.26
[55] rvest_1.0.2 lifecycle_1.0.1 XML_3.99-0.8 zlibbioc_1.40.0 MASS_7.3-54 scales_1.1.1
[61] hms_1.1.1 proj4_1.0-10.1 parallel_4.1.1 RColorBrewer_1.1-2 memoise_2.0.0 ggrastr_0.2.3
[67] emdbook_1.3.12 bdsmatrix_1.3-4 stringi_1.7.5 RSQLite_2.2.8 genefilter_1.76.0 foreach_1.5.1
[73] BiocParallel_1.28.0 shape_1.4.6 rlang_0.4.11 pkgconfig_2.0.3 bitops_1.0-7 lattice_0.20-44
[79] labeling_0.4.2 htmlwidgets_1.5.4 bit_4.0.4 tidyselect_1.1.1 plyr_1.8.6 magrittr_2.0.1
[85] R6_2.5.1 generics_0.1.1 DelayedArray_0.20.0 DBI_1.1.1 pillar_1.6.4 haven_2.4.3
[91] withr_2.4.2 survival_3.2-11 KEGGREST_1.34.0 RCurl_1.98-1.5 ash_1.0-15 modelr_0.1.8
[97] crayon_1.4.2 KernSmooth_2.23-20 utf8_1.2.2 tzdb_0.1.2 GetoptLong_1.0.5 locfit_1.5-9.4
[103] readxl_1.3.1 blob_1.2.2 reprex_2.0.1 digest_0.6.28 xtable_1.8-4 numDeriv_2016.8-1.1
[109] munsell_0.5.0 beeswarm_0.4.0 vipor_0.4.5
If you inverse the conditions A and B in
res_control_CV <- results(DE_control_CV, contrast = c("condition", "C", "CV"))that is to sayres_control_CV <- results(DE_control_CV, contrast = c("condition", "CV", "C"))then it is ok I think, but the first one you wrote is Control vs CVThanyou so much for replying so quickly Dr Love.
Yes, if I use: r
res_control_CV <- results(DE_control_CV, contrast = c("condition", "CV", "C"))It works...
My problem is, and this may be a stupid exercise, I need C vs CV for a sanity check, in order to understand if the difference between C vs CV and CV vc C is just a minus...
Even if i change the above line, and with the relevel command, I always get:
r
> resultsNames(DE_control_CV) [1] "Intercept" "condition_CV_vs_C"And I need r
condition_C_vs_CVI am sorry once again if this is stupid.
And thankyou for this amazing software btw
-jorge
You need to relevel before DESeq() or nbinomWaldTest() as indicated in that section of the vignette.
Sorry, maybe I am missing something, but in the original code provided, the relevel is before DESeq...is the nbiomWaldTest mandatory for this to work?
PS: Could you also confirm that swaping Control vs CV only inverts the signal?
Yes, it only -1x the log2FoldChange.
I see the problem. In R,
referencerefers to the level of the factor you want to compare to.At the beginning you said "I have CV vs control, but I want the other way around". Then you would need to put CV as the reference.
It worked! Thanyou so much for all your patience.
The ref being CV when wanting "C vs CV" was not intuitive for me