gabrielhoffman / dreamlet Goto Github PK

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:

1. ~Sex+Age+Batch --> all categorical variables as fixed factors
2. ~(1|Sex)+Age+(1|Batch) --> all categorical variables as random factors
3. ~Sex+Age --> in this case, the Batch contribution would be part of the residuals

Thanks 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:

• Using limma::selectModel with the design matrix from my res.dl object. However, I was not sure if the calculation is correct for dreamlet objects (I don’t think it includes random effects parameters?).
• Calculating manually using the formula for AIC and the logLik slot in the res.dl object. I was just wondering exactly how to calculate the number of parameters - is this just equal to the number of fixed and random effects parameters or are there additional parameters that should be added related to the limma-voom method? Is there a way to get the number of parameters from the res.dl object (in case some variables were dropped for certain cell types, for example)?
  Thank you!!

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:

1. processAssays(pb, form, assays = “celltype”, min.count = 10,
   . min.prop = 0.8, min.cells = 5, min.samples = 50, BPPARAM = SnowParam(12))
2. 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
   . })
3. FUN(X[[i]], ...)
4. 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,
   . ...)
5. 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.0

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

time zone: Europe/Berlin
tzcode source: system (glibc)

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

other 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.0

loaded 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!