MarioniLab / scran

Clone of the Bioconductor repository for the scran package.
https://bioconductor.org/packages/devel/bioc/html/scran.html
39 stars 23 forks source link

psuedoBulkDGE #80

Open elachtarajax opened 3 years ago

elachtarajax commented 3 years ago

Hello, Is this function no longer available with scran? Best, Emily

LTLA commented 3 years ago

Should be there. Make sure you're using a reasonably version of the package - see install instruction here.

brunosa3 commented 3 years ago

is it possible that there is a bug in pseudoBulkDGE?

I get the following error when using pseudobulkdge

`samples <- colData(sce) %>% as.data.frame() %>% distinct(sample, .keep_all = TRUE) %>% DataFrame()

design <- model.matrix(~ 0 + health, data = samples) %>% magrittr::set_rownames(samples$sample)

contrast <- limma::makeContrasts(healthad - healthhc, levels = design)

dea <- pseudoBulkDGE(se, sample = se$sample, label = se$celltype, condition = se$health, design = design, contrast = contrast)`

Error in .compute_offsets_by_lfc(design = curdesign, coef = coef, contrast = contrast, : argument "null.lfc" is missing, with no default In addition: Warning message: Error in .compute_offsets_by_lfc(design = curdesign, coef = coef, contrast = contrast, : argument "null.lfc" is missing, with no default

I cannot see where the function assigns the appropriate value to lfc.null

my sessionInfo()

`R version 4.0.3 (2020-10-10) Platform: x86_64-pc-linux-gnu (64-bit) Running under: CentOS Linux 7 (Core)

Matrix products: default BLAS/LAPACK: /usr/prog/OpenBLAS/0.2.20-GCC-6.4.0-2.28/lib/libopenblas_haswellp-r0.2.20.so

locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

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

other attached packages: [1] edgeR_3.32.0 limma_3.46.0
[3] scran_1.18.3 scater_1.18.3
[5] SingleCellExperiment_1.12.0 SummarizedExperiment_1.20.0 [7] Biobase_2.50.0 GenomicRanges_1.42.0
[9] GenomeInfoDb_1.26.1 IRanges_2.24.0
[11] S4Vectors_0.28.0 BiocGenerics_0.36.0
[13] MatrixGenerics_1.2.0 matrixStats_0.57.0
[15] forcats_0.5.0 stringr_1.4.0
[17] dplyr_1.0.2 purrr_0.3.4
[19] readr_1.4.0 tidyr_1.1.2
[21] tibble_3.0.4 ggplot2_3.3.2
[23] tidyverse_1.3.0

loaded via a namespace (and not attached): [1] readxl_1.3.1 backports_1.2.0 plyr_1.8.6
[4] igraph_1.2.6 lazyeval_0.2.2 splines_4.0.3
[7] BiocParallel_1.24.1 listenv_0.8.0 scattermore_0.7
[10] digest_0.6.27 htmltools_0.5.0 viridis_0.5.1
[13] fansi_0.4.1 magrittr_1.5 tensor_1.5
[16] cluster_2.1.0 ROCR_1.0-11 globals_0.14.0
[19] modelr_0.1.8 colorspace_2.0-0 rvest_0.3.6
[22] ggrepel_0.8.2 haven_2.3.1 crayon_1.3.4
[25] RCurl_1.98-1.2 jsonlite_1.7.1 spatstat_1.64-1
[28] spatstat.data_1.7-0 survival_3.2-7 zoo_1.8-8
[31] glue_1.4.2 polyclip_1.10-0 pals_1.6
[34] gtable_0.3.0 zlibbioc_1.36.0 XVector_0.30.0
[37] leiden_0.3.7 DelayedArray_0.16.0 BiocSingular_1.6.0
[40] future.apply_1.7.0 maps_3.3.0 abind_1.4-5
[43] scales_1.1.1 DBI_1.1.0 miniUI_0.1.1.1
[46] Rcpp_1.0.5 viridisLite_0.3.0 xtable_1.8-4
[49] dqrng_0.2.1 reticulate_1.18 rsvd_1.0.3
[52] mapproj_1.2.7 htmlwidgets_1.5.2 httr_1.4.2
[55] RColorBrewer_1.1-2 ellipsis_0.3.1 Seurat_4.0.0
[58] ica_1.0-2 scuttle_1.0.4 pkgconfig_2.0.3
[61] uwot_0.1.10 dbplyr_2.0.0 deldir_0.2-9
[64] locfit_1.5-9.4 tidyselect_1.1.0 rlang_0.4.9
[67] reshape2_1.4.4 later_1.1.0.1 munsell_0.5.0
[70] cellranger_1.1.0 tools_4.0.3 cli_2.2.0
[73] generics_0.1.0 broom_0.7.2 ggridges_0.5.3
[76] fastmap_1.0.1 goftest_1.2-2 fs_1.5.0
[79] fitdistrplus_1.1-3 RANN_2.6.1 sparseMatrixStats_1.2.0
[82] pbapply_1.4-3 future_1.21.0 nlme_3.1-150
[85] mime_0.9 xml2_1.3.2 compiler_4.0.3
[88] rstudioapi_0.13 beeswarm_0.2.3 plotly_4.9.2.1
[91] png_0.1-7 spatstat.utils_2.0-0 reprex_0.3.0
[94] statmod_1.4.35 stringi_1.5.3 bluster_1.0.0
[97] lattice_0.20-41 Matrix_1.2-18 vctrs_0.3.5
[100] pillar_1.4.7 lifecycle_0.2.0 BiocManager_1.30.10
[103] lmtest_0.9-38 BiocNeighbors_1.8.2 RcppAnnoy_0.0.18
[106] data.table_1.13.2 cowplot_1.1.0 bitops_1.0-6
[109] irlba_2.3.3 httpuv_1.5.4 patchwork_1.1.1
[112] R6_2.5.0 promises_1.1.1 KernSmooth_2.23-18
[115] gridExtra_2.3 vipor_0.4.5 parallelly_1.23.0
[118] writexl_1.3.1 codetools_0.2-18 dichromat_2.0-0
[121] MASS_7.3-53 assertthat_0.2.1 withr_2.3.0
[124] SeuratObject_4.0.0 sctransform_0.3.2 GenomeInfoDbData_1.2.4
[127] mgcv_1.8-33 hms_0.5.3 beachmat_2.6.4
[130] grid_4.0.3 rpart_4.1-15 DelayedMatrixStats_1.12.2 [133] Rtsne_0.15 shiny_1.5.0 lubridate_1.7.9.2
[136] ggbeeswarm_0.6.0 `

LTLA commented 3 years ago

There is a bug, but also you should have gotten a deprecation warning along the lines of:

matrix arguments for 'design=' are deprecated. 
Use functions or formulas instead.

Matrices are deprecated because it was really painful to figure out how the rows of the matrix matched up to the subset of data for a particular cell type if that cell type wasn't present in all samples. Then we would have to remove rows of the matrix, and the result might not be of full rank... very tricky. The solution is just to pass:

dea <- pseudoBulkDGE(se, sample = se$sample, label = se$celltype,
    condition = se$health, design = ~0+ health,
    contrast = "healthad - healthhc")`
brunosa3 commented 3 years ago

wow...many thanks you made my day :-D I have overseen this warning - my bad.

It worked 👍

However, I somehow realized that my colleague with the same code (only difference is the adjustment you shared here) has slightly different DGE lists.

Did you change anything in the way you do the statistics when comparing version 1.16 to 1.18 or would you expect the same DGE lists?

LTLA commented 3 years ago

Hm... nothing on my end, at least not from memory. Something may have changed in edgeR itself.

Might be worth installing 1.18.4 from the RELEASE_3_12 branch and seeing what happens with and without the adjustment.

Edit: Hm. There is a difference. Wonder why that happens.

Edit 2: Oh, it's another bug in the old implementation, related to how condition= is handled. The same bug is not present in the new implementation. It is a minor bug that should not affect the validity of the results; I guess I must have just left it in there for back-compatibility, given that I was going to deprecate the old implementation anyway.