gabrielhoffman / dreamlet Goto Github PK
View Code? Open in Web Editor NEWPerform differential expression analysis on multi-sample single cell datasets using linear mixed models
Home Page: https://gabrielhoffman.github.io/dreamlet
Perform differential expression analysis on multi-sample single cell datasets using linear mixed models
Home Page: https://gabrielhoffman.github.io/dreamlet
Hi,
Thanks for the awesome package!
I'm trying to compare two populations from the same samples using dreamletCompareClusters. The problem is that I'm trying to provide different covariates for different populations within the same sample, for stuff like picard QC metrics generated for population-specific bam files. If I'm not mistaken, dreamlet only reads the shared sample-specific metadata. I can get around the problem by using regular voomWithDreamWeights() and dream(), but it won't use the cell number when estimating the precision weights. Are there any ways to manually feed separate covariates to dreamletCompareClusters()?
Hi,
Thank you for your great package.
I am running into an issue when running processAssays() with certain cell types.
Please see the example code below:
# Load libraries
suppressPackageStartupMessages({
library(SingleCellExperiment)
library(zellkonverter)
library(dreamlet)
library(rhdf5)
library(HDF5Array)})
#read H5AD files
file <- "../datasets/adata_annotated.h5ad"
sce <- readH5AD(file, use_hdf5=TRUE, verbose=FALSE)
#Compute pseudobulk
pb <- aggregateToPseudoBulk(sce,
assay = "counts",
cluster_id = “annotated_celltypes",
sample_id = "sample",
BPPARAM = SnowParam(6, progressbar=TRUE))
#Normalize and apply voom
form <- ~ (1|Batch) + (1|Sex) + (1|Disease) + Age + SV1 + SV2 + SV3 + SV4
res.proc_celltype <- processAssays(pb, form, assays='celltype' min.count=10,
min.prop=0.8, min.cells=5, min.samples=50, BPPARAM = SnowParam(12))
celltype...
Error in weightsMatrix * out$weights: non-conformable arrays
Traceback:
- processAssays(pb, form, assays = “celltype”, min.count = 10,
. min.prop = 0.8, min.cells = 5, min.samples = 50, BPPARAM = SnowParam(12))- lapply(assays, function(k) {
. if (!quiet)
. message(" ", k, "...", appendLF = FALSE)
. startTime <- Sys.time()
. y <- assay(sceObj, k)
. n.cells <- n.cells_full[colnames(y), k, drop = FALSE]
. data <- merge_metadata(data_constant, metadata(sceObj)$aggr_means,
. k, metadata(sceObj)$agg_pars$by)
. res <- processOneAssay(y[, rownames(data), drop = FALSE],
. formula = formula, data = data, n.cells = n.cells[rownames(data),
. , drop = FALSE], min.cells = min.cells, min.count = min.count,
. min.samples = min.samples, min.prop = min.prop, isCounts = isCounts,
. normalize.method = normalize.method, useCountsWeights = useCountsWeights,
. span = span, BPPARAM = BPPARAM, ...)
. if (!quiet)
. message(format(Sys.time() - startTime, digits = 2))
. res
. })- FUN(X[[i]], ...)
- processOneAssay(y[, rownames(data), drop = FALSE], formula = formula,
. data = data, n.cells = n.cells[rownames(data), , drop = FALSE],
. min.cells = min.cells, min.count = min.count, min.samples = min.samples,
. min.prop = min.prop, isCounts = isCounts, normalize.method = normalize.method,
. useCountsWeights = useCountsWeights, span = span, BPPARAM = BPPARAM,
. ...)- voomWithDreamWeights(y[keep, ], formula, data, weights = w_cells,
. BPPARAM = BPPARAM, ..., save.plot = TRUE, quiet = quiet,
. span = span, hideErrorsInBackend = TRUE)
sessionInfo()
R version 4.3.1 (2023-06-16)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)Matrix products: default
BLAS/LAPACK: /conda-envs/dreamlet/lib/libopenblasp-r0.3.24.so; LAPACK version 3.11.0locale:
[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=Ctime zone: Europe/Berlin
tzcode source: system (glibc)attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] baseother attached packages:
[1] readxl_1.4.3 HDF5Array_1.28.1
[3] DelayedArray_0.27.10 SparseArray_1.1.12
[5] S4Arrays_1.1.6 abind_1.4-5
[7] Matrix_1.6-1.1 rhdf5_2.44.0
[9] dreamlet_0.99.28 variancePartition_1.31.20
[11] BiocParallel_1.35.4 limma_3.57.10
[13] ggplot2_3.4.4 zellkonverter_1.11.3
[15] SingleCellExperiment_1.23.0 SummarizedExperiment_1.31.1
[17] Biobase_2.61.0 GenomicRanges_1.53.2
[19] GenomeInfoDb_1.37.6 IRanges_2.35.3
[21] S4Vectors_0.39.3 BiocGenerics_0.47.0
[23] MatrixGenerics_1.12.2 matrixStats_1.0.0loaded via a namespace (and not attached):
[1] splines_4.3.1 pbdZMQ_0.3-10
[3] bitops_1.0-7 filelock_1.0.2
[5] tibble_3.2.1 cellranger_1.1.0
[7] basilisk.utils_1.13.3 graph_1.79.4
[9] XML_3.99-0.14 lifecycle_1.0.3
[11] Rdpack_2.5 mixsqp_0.3-48
[13] edgeR_3.99.2 lattice_0.21-9
[15] MASS_7.3-60 backports_1.4.1
[17] magrittr_2.0.3 reticulate_1.34.0
[19] DBI_1.1.3 minqa_1.2.6
[21] zlibbioc_1.47.0 EnvStats_2.8.1
[23] purrr_1.0.2 rmeta_3.0
[25] msigdbr_7.5.1 RCurl_1.98-1.12
[27] GenomeInfoDbData_1.2.10 ggrepel_0.9.4
[29] pbkrtest_0.5.2 irlba_2.3.5.1
[31] annotate_1.79.0 DelayedMatrixStats_1.23.9
[33] codetools_0.2-19 tidyselect_1.2.0
[35] lme4_1.1-34 base64enc_0.1-3
[37] jsonlite_1.8.7 iterators_1.0.14
[39] tools_4.3.1 progress_1.2.2
[41] snow_0.4-4 Rcpp_1.0.11
[43] zenith_1.3.1 glue_1.6.2
[45] IRdisplay_1.1 dplyr_1.1.3
[47] withr_2.5.1 numDeriv_2016.8-1.1
[49] fastmap_1.1.1 basilisk_1.13.4
[51] boot_1.3-28.1 rhdf5filters_1.12.1
[53] fansi_1.0.5 caTools_1.18.2
[55] digest_0.6.33 truncnorm_1.0-9
[57] R6_2.5.1 colorspace_2.1-0
[59] scattermore_1.2 gtools_3.9.4
[61] RSQLite_2.3.1 RhpcBLASctl_0.23-42
[63] utf8_1.2.3 tidyr_1.3.0
[65] generics_0.1.3 data.table_1.14.8
[67] corpcor_1.6.10 prettyunits_1.2.0
[69] httr_1.4.7 pkgconfig_2.0.3
[71] gtable_0.3.4 blob_1.2.4
[73] XVector_0.41.1 remaCor_0.0.16
[75] htmltools_0.5.6.1 GSEABase_1.63.0
[77] scales_1.2.1 png_0.1-8
[79] fANCOVA_0.6-1 ashr_2.2-63
[81] reshape2_1.4.4 uuid_1.1-1
[83] nlme_3.1-163 nloptr_2.0.3
[85] repr_1.1.6 cachem_1.0.8
[87] stringr_1.5.0 KernSmooth_2.23-22
[89] parallel_4.3.1 RcppZiggurat_0.1.6
[91] AnnotationDbi_1.63.2 pillar_1.9.0
[93] grid_4.3.1 vctrs_0.6.4
[95] gplots_3.1.3 mashr_0.2.69
[97] xtable_1.8-4 Rgraphviz_2.45.0
[99] evaluate_0.22 KEGGgraph_1.61.0
[101] invgamma_1.1 mvtnorm_1.2-3
[103] cli_3.6.1 locfit_1.5-9.8
[105] compiler_4.3.1 rlang_1.1.1
[107] crayon_1.5.2 SQUAREM_2021.1
[109] plyr_1.8.9 stringi_1.7.12
[111] assertthat_0.2.1 babelgene_22.9
[113] lmerTest_3.1-3 munsell_0.5.0
[115] Biostrings_2.69.2 aod_1.3.2
[117] dir.expiry_1.9.0 IRkernel_1.3.2
[119] hms_1.1.3 sparseMatrixStats_1.13.4
[121] bit64_4.0.5 Rhdf5lib_1.22.0
[123] KEGGREST_1.41.4 statmod_1.5.0
[125] rbibutils_2.2.15 Rfast_2.0.8
[127] broom_1.0.5 memoise_2.0.1
[129] bit_4.0.5 EnrichmentBrowser_2.31.5
Thank you for the excellent method development.
I am very interested in your work and am hoping to apply it to my research. However, I have some confusion regarding the formula in this package. The structure of my data is as follows: I have two groups of data. One consists of 5 individuals in the control group, and the other consists of 7 individuals in the therapy-treated group. I aim to identify cell-type-specific DEGs between these two groups while considering the individual as a random effect. How should I design the formula?
form <- ~ treatment_group + (1|individual)
form <- ~ treatment_group + (1|individual) + 0
Which one is correct? (Here, I am not considering other confounding factors)
Dear all,
I understand the dreamlet::fitVarPart()
is internally using the variancePartition package, which uses linear and linear mixed models to quantify the contribution of multiple sources of expression variation at the gene-level.
When calling dreamlet::fitVarPart()
using a mixed-model formula like this ~Sex + Age + (1|Batch)
, where Sex and Batch are categorical and Age is continuous, I encountered the following error:
"Error in run_model_check_mixed(fit, showWarnings, dream, colinearityCutoff, :
Categorical variables modeled as fixed effect: Gender
Must model either all or no categorical variables as random effects here"
From this, I understand the function can only handle mixed formulas with all categorical variables treated as fixed or random factors. Likewise, using the following formulas no longer gives an error: ~Sex+Age+Batch
or ~(1|Sex)+Age+(1|Batch)
.
Here you can find the cross-table showing the number of cells per Batch (columns, which are the sequencing dates) and Sex (rows):
180925 180926 181003 181022 181023 181107 181108 181213 181218
M 454 660 642 196 751 408 245 688 512
F 640 566 533 624 725 705 694 345 612
In the differential expression analysis, we're treating Batch as a random effect and Sex as a fixed effect, but we're only interested in the Sex or Age coefficients/effects, but we still want to control for Batch. Hence, I'm wondering which option would be better (or correct) for the variance partition analysis, in order to be the most similar to the differential expression analysis set-up:
~Sex+Age+Batch
--> all categorical variables as fixed factors~(1|Sex)+Age+(1|Batch)
--> all categorical variables as random factors~Sex+Age
--> in this case, the Batch contribution would be part of the residualsThanks a lot!
Aida
Fantastic package! I have a small request, could you add a function to extract normalized expression values (for boxplots for instance). I wrote code to do that myself but would be nice to have it out of the box.
Hi!
I am currently using the default values on the arguments in the processAssays() function. However, I am missing some important genes after applying this function, which were not filtered out before when I was using another previous QC filtering on the count matrix.
Thus, I would like to ask you if there is any other argument or whether I can tune a bit the function to get the sample/gene statistics for each of the arguments (min.cells, min.count, min.samples, min.prop) used for filtering the samples/genes.
Thanks a lot!
Aida
Thank you for the great package. I have a problem when one contrast doesn't exist for one cell-type (assay). This leads to the below error. I understand that no result can be reported for those cell types. The problem is that when I do multiple group comparisons I have to run them separately for each group. It would be great to fill it either with NAN or skip this element in the results. Currently, I need to rerun the analysis separately.
Error in eval_tidy(expr, data, env): object 'TTT' not found
Traceback:
1. dreamlet(res.proc, form, colData(res.proc), assays = assayNames(res.proc),
. contrasts = c(Responder_vs_Non_responder = "responder_weekResponder.w0 - responder_weekNon_responder.w0"),
. computeResiduals = FALSE, )
2. dreamlet(res.proc, form, colData(res.proc), assays = assayNames(res.proc),
. contrasts = c(Responder_vs_Non_responder = "responder_weekResponder.w0 - responder_weekNon_responder.w0"),
. computeResiduals = FALSE, )
3. lapply(assays, function(k) {
. if (!quiet)
. message(" ", k, "...", appendLF = FALSE)
. startTime = Sys.time()
. geneExpr = assay(x, k)
. ids = intersect(colnames(geneExpr), rownames(data_constant))
. geneExpr = geneExpr[, ids, drop = FALSE]
. data2 = merge_metadata(data_constant[ids, , drop = FALSE],
. metadata(x), k, x@by)
. form_mod = removeConstantTerms(formula, data2)
. form_mod = dropRedundantTerms(form_mod, data2)
. if (length(all.vars(form_mod)) > 0) {
. if (!is.null(contrasts)) {
. L = makeContrastsDream(form_mod, data2, contrasts = contrasts,
. nullOnError = TRUE)
. }
. else {
. L = NULL
. }
. fit = tryCatch({
. dream(geneExpr, form_mod, data2, L = L, BPPARAM = BPPARAM,
. ..., quiet = TRUE)
. }, error = function(e) NULL)
. if (!is.null(fit) && !any(is.na(fit$sigma))) {
. if (!is.null(fit$rdf)) {
. keep = which(fit$rdf >= 1)
. fit = fit[keep, ]
. }
. if (use.eBayes) {
. isCounts = ifelse(is(geneExpr, "EList"), TRUE,
. FALSE)
. fit = eBayes(fit, robust = robust, trend = isCounts)
. }
. }
. else {
. fit = NULL
. }
. }
. else {
. fit = NULL
. }
. if (!quiet)
. message(format(Sys.time() - startTime, digits = 2))
. list(fit = fit, formula = form_mod, n_retain = ncol(geneExpr))
. })
4. FUN(X[[i]], ...)
5. makeContrastsDream(form_mod, data2, contrasts = contrasts, nullOnError = TRUE)
6. lapply(e, eval_tidy_safe, env = L_uni_env, nullOnError = nullOnError)
7. FUN(X[[i]], ...)
8. eval_tidy(expr, data, env) ```
run_mash appears to be set up to run on multiple coefficients, but throws an error when multiple are provided. Specifically the line
if (!coefList %in% coefNames(fit)) {
stop("coef not found in coefNames(fit): ", coefList)
}
Perhaps this should read
for (coef in coefList) {
if (!coef %in% coefNames(fit)) {
stop("coef not found in coefNames(fit): ", coef)
}
}
I've figured out this issue already but just wanted to post in case anyone else has the same issue or in case a more informative error message could be provided to save time for others facing this issue.
I ran a dreamlet script that previously worked and got the following error:
Error in processAssays(pb, form, min.count = 5, min.cells = 5, BPPARAM = SnowParam(num_cores, : All models failed. Consider changing formula
In addition: Warning messages:
1: In voomWithDreamWeights(y[keep, ], formula, data, weights = precWeights : The experimental design has no replication. Setting weights to 1.
In this case, the problem was not the pb object or the model formula. I believe the issue was related to that described here since this is the error I got when I tried to run lme4 separately outside of dreamlet:
Error in initializePtr() :
function 'cholmod_factor_ldetA' not provided by package 'Matrix'
I checked my sessionInfo and saw that I had unintentionally updated the Matrix package to a newer version compared to when I had successfully run the dreamlet script previously. I was able to fix the error using the approach described here: https://stackoverflow.com/questions/77481539/error-in-initializeptr-function-cholmod-factor-ldeta-not-provided-by-pack to reinstall lme4 (and other packages depending on Matrix) from source which fixed the issues with lme4 and dreamlet. Hope this helps someone!
Hi,
I would like to know which default values (min.cells, min.count, min.samples, min.prop) you recommend as I'm not sure the documentation you provide for this specific function is updated documentation. You can see that it says dreamlet 1.1.15, but if you check the general documentation page it says dreamlet 1.1.18.
In my case, using R/4.2.0:
> packageVersion("dreamlet")
[1] '0.99.6'
> ?processAssays
processAssays(
sceObj,
formula,
assays = assayNames(sceObj),
min.cells = 10,
min.count = 10,
min.samples = 4,
min.prop = 0.4,
isCounts = TRUE,
normalize.method = "TMM",
useCountsWeights = TRUE,
quiet = FALSE,
BPPARAM = SerialParam(),
...
)
And, using R/4.3.0:
> packageVersion("dreamlet")
[1] '0.99.16'
> ?processAssays
processAssays(
sceObj,
formula,
assays = assayNames(sceObj),
min.cells = 10,
min.count = 10,
min.samples = 4,
min.prop = 0.4,
isCounts = TRUE,
normalize.method = "TMM",
useCountsWeights = TRUE,
quiet = FALSE,
BPPARAM = SerialParam(),
...
)
Thanks a lot!
Aida
Hi Gabriel,
Thanks for the great package.
How can I retrieve the genes lists from the dreamlet object after DE, without the toptable?
I wish to color specific genes using the plotVolcano function.
Best,
Michal
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.