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
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) ```
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
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)
}
}
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
Hi,
I've been using your dreamlet package for a number of DE analyses of high-dimension transcriptomic data and came across something that's got me scratching my head in the provided full-analysis scripts examining published data. In many cases, the models used in processAssays()
and subsequently in dreamlet()
are not the same.
A straightforward example comes from https://github.com/GabrielHoffman/dreamlet_analysis/blob/main/COVID_combat/COVID_combat_major.Rmd , wherein sex
and Source
are random variables in the processAssays()
call but subsequently converted to fixed effects. In other examples, variables used in processAssays are simply dropped from the formula passed to dreamlet:
form = ~ Age + **(1|sex) + (1|Source)**
res.proc = processAssays(pb, form,
min.count = 3,
min.cells = 3,
BPPARAM = SnowParam(6))
res.vp = fitVarPart(res.proc, form, BPPARAM=SnowParam(6) )
form = ~ Source + **Age + sex**
fit = dreamlet(res.proc, form, BPPARAM=SnowParam(6))
It's somewhat unclear at the level of the code what is prompting these changes in formula between steps, and I can't find much clarifying why one would want to use different variables in the two steps in either the variancePartition
or the dream
/dreamlet
docs. Some pointers would be much appreciated!
Best regards, and great work with this package! -Bernie Mulvey
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
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, I am trying to use dreamlet
for a dataset with one control and multiple treatment groups per cell type. I wish to compare each treatment group to the same control group. In some other tools, this is achieved by converting the treatment column in the metadata to a factor and setting "Control" as the first level.
Does dreamlet
follow a similar behavior?
Thanks,
Kunal
Hi!
I have 2 genotypes, disease and control; the disease one has 6 timepoints, the control just one (the first of the disease group). In the tutorial, the ctrl group has the same timepoints as the treated one, whereas in my case the disease one has 6 timepoints and the ctrl just 1 (the first of the disease group).
There are two types analysis I am considering, one to compare all timepoints against the wt and one to compare each disease timepoint with the previous one. Plus, I could also model continuous covariates, as the sce object contains eg
n_genes_by_counts # nr of genes with a least 1 count in a cell
total_counts # sum of counts for a gene
total_counts_mito # sum of counts for a mito gene
total_counts_ribo # sum of counts for a ribo gene
# pct_counts_mito # percentage of counts for a mito gene
# pct_counts_ribo # percentage of counts for a ribo gene
Any help will be much apppreciated!
Hi, I’m trying to compare several closely related dreamlet models - one with a simple linear term for experimental timepoint, one that treats timepoint as categorical (I have only three timepoints), and one with linear and quadratic terms (as you nicely described here https://diseaseneurogenomics.github.io/dreamlet/articles/non_lin_eff.html). Would you recommend using AIC scores for this? If so, can you give me guidance on how to calculate this correctly for dreamlet objects? I was planning to look at the median AIC score across all genes for each model. Also open to suggestions if there is a different strategy or metric that might be better. What I tried:
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
Hello, thank you for all your work in providing these great analysis tools.
I am trying to run dreamlet on a merged seurat object which has undergone cca integration. I have aggregated to pseudobulk, splitting by 3 pseudotime bins and am attempting to processAssays/examine variance:
form <- ~ Condition + Donor + nCount_RNA + log1pcounts + mitoPercent + nFeature_RNA
res.proc <- processAssays(pb, form, min.cells = 5,
min.count = 5,
min.samples = 4,
min.prop = 0.4,
verbose = TRUE
)
details(res.proc)
However, I get errors saying certain coefficients not estimable, and its always the last two coefficients in my model form, regardless of what they are. For example:
form <- ~ Condition + Donor + nCount_RNA + log1pcounts + mitoPercent + nFeature_RNA
will show:
Pseudo1...Coefficients not estimable: mitoPercent nFeature_RNA
Warning: Partial NA coefficients for 15792 probe(s)
Warning in voomWithDreamWeights(y[keep, ], formula, data, weights = precWeights, :
The experimental design has no replication. Setting weights to 1.
0.41 secs
OR
form <- ~ Condition + Donor + mitoPercent + nFeature_RNA + nCount_RNA + log1pcounts
will give the error:
Pseudo1...Coefficients not estimable: nCount_RNA + log1pcounts
Warning: Partial NA coefficients for 15792 probe(s)
Warning in voomWithDreamWeights(y[keep, ], formula, data, weights = precWeights, :
The experimental design has no replication. Setting weights to 1.
0.41 secs
I am confused as to whether its simply not listing more things? or if there is some actual error with my model form?
Regardless -- whenever I get these warnings it won't allow me to plot variance for any of these covariates.
Thanks in advance for your help
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
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.
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!
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.