problem with analyzeComparative() with only one DNA replicate

[robinmeyers]

Hello,

I believe this is a bug when running analyzeComparative() with only a single DNA replicate. It appears as though the DNA design is specified correctly. However for a given enhancer, the model only estimates a single abundance for all barcodes.

Thanks for a great tool!

-Robin

library(MPRAnalyze)
library(magrittr)
library(tidyverse)


mpra_sim <- function(n_enh = 50, n_bcs = 5, reads = 50000, rna_conditions = 2, rna_reps = 2, dna_reps = 1) {

  true_dna_abundance <- rgamma(n_bcs*n_enh, shape = 2, rate = 2) %>%
    set_names(paste0("enhancer_", rep(1:n_enh, each = n_bcs), ".bc_", rep(1:n_bcs, n_enh)))

  true_activity <- rnorm(n_enh)

  dna_counts <- rmultinom(n = dna_reps, size = reads, prob = exp(true_dna_abundance)) %>%
    set_colnames(paste0("dna_", 1:dna_reps)) %>%
    as_tibble(rownames = "element") %>%
    separate(element, into = c("enhancer", "bc"), "\\.") %>%
    pivot_wider(names_from = "bc", values_from = starts_with("dna"), names_glue = "{.value}.{bc}") %>%
    as.data.frame() %>% set_rownames(.$enhancer) %>% select(-enhancer)

  rna_libs <- paste0("rna_", rep(1:rna_conditions, each = rna_reps), ".rep_", rep(1:rna_reps, rna_conditions))

  rna_counts <- rmultinom(n = length(rna_libs), size = reads, prob = exp(true_dna_abundance + rep(true_activity, each = n_bcs))) %>%
    set_colnames(rna_libs) %>%
    as_tibble(rownames = "element") %>%
    separate(element, into = c("enhancer", "bc"), "\\.") %>%
    pivot_wider(names_from = "bc", values_from = starts_with("rna"), names_glue = "{.value}.{bc}") %>%
    as.data.frame() %>% set_rownames(.$enhancer) %>% select(-enhancer)

  dna_annot <- tibble(x = colnames(dna_counts)) %>%
    mutate(barcode = str_extract(x, "bc_\\d+")) %>%
    as.data.frame() %>% set_rownames(.$x) %>% select(-x)

  rna_annot <- tibble(x = colnames(rna_counts)) %>%
    mutate(condition = str_extract(x, "rna_\\d+")) %>%
    mutate(barcode = str_extract(x, "bc_\\d+")) %>%
    as.data.frame() %>% set_rownames(.$x) %>% select(-x)


  mpra_obj <- MpraObject(dnaCounts = as.matrix(dna_counts),
                         rnaCounts = as.matrix(rna_counts),
                         dnaAnnot = dna_annot,
                         rnaAnnot = rna_annot)

  mpra_obj <- analyzeComparative(obj = mpra_obj,
                                 dnaDesign = ~ barcode,
                                 rnaDesign = ~ condition,
                                 reducedDesign = ~ 1,
                                 fit.se = T)
  return(mpra_obj)
}


mpra_obj_1 <- mpra_sim(dna_reps = 1)
#> Fitting model...
#> Fitting reduced model...
#> Analysis Done!
mpra_obj_2 <- mpra_sim(dna_reps = 2)
#> Fitting model...
#> Fitting reduced model...
#> Analysis Done!


mpra_obj_1@designs@dna
#>            (Intercept) barcodebc_2 barcodebc_3 barcodebc_4 barcodebc_5
#> dna_1.bc_1           1           0           0           0           0
#> dna_1.bc_2           1           1           0           0           0
#> dna_1.bc_3           1           0           1           0           0
#> dna_1.bc_4           1           0           0           1           0
#> dna_1.bc_5           1           0           0           0           1
#> attr(,"assign")
#> [1] 0 1 1 1 1
#> attr(,"contrasts")
#> attr(,"contrasts")$barcode
#> [1] "contr.treatment"
mpra_obj_2@designs@dna
#>            (Intercept) barcodebc_2 barcodebc_3 barcodebc_4 barcodebc_5
#> dna_1.bc_1           1           0           0           0           0
#> dna_1.bc_2           1           1           0           0           0
#> dna_1.bc_3           1           0           1           0           0
#> dna_1.bc_4           1           0           0           1           0
#> dna_1.bc_5           1           0           0           0           1
#> dna_2.bc_1           1           0           0           0           0
#> dna_2.bc_2           1           1           0           0           0
#> dna_2.bc_3           1           0           1           0           0
#> dna_2.bc_4           1           0           0           1           0
#> dna_2.bc_5           1           0           0           0           1
#> attr(,"assign")
#> [1] 0 1 1 1 1
#> attr(,"contrasts")
#> attr(,"contrasts")$barcode
#> [1] "contr.treatment"

head(getModelParameters_DNA(mpra_obj_1))
#>                   disp (Intercept) barcodebc_2 barcodebc_3 barcodebc_4
#> enhancer_1 0.590367627    4.449182          NA          NA          NA
#> enhancer_2 1.297220239    3.501869          NA          NA          NA
#> enhancer_3 2.175638767    2.537478          NA          NA          NA
#> enhancer_4 1.347186336    4.281143          NA          NA          NA
#> enhancer_5 2.785081679    1.734531          NA          NA          NA
#> enhancer_6 0.007270567    5.406162          NA          NA          NA
#>            barcodebc_5
#> enhancer_1          NA
#> enhancer_2          NA
#> enhancer_3          NA
#> enhancer_4          NA
#> enhancer_5          NA
#> enhancer_6          NA
head(getModelParameters_DNA(mpra_obj_2))
#>                disp (Intercept) barcodebc_2 barcodebc_3 barcodebc_4 barcodebc_5
#> enhancer_1 5.492000   0.1508875  -1.6231356  -0.3664045 -0.53633054  -0.6380493
#> enhancer_2 6.241402  -2.1520033   0.2280683   0.4045109  0.63723936   0.6695159
#> enhancer_3 5.633199   0.2703070  -0.7146878  -1.2007639 -1.65773848  -1.3400335
#> enhancer_4 5.860700  -1.2570581   0.5824395  -0.5040491 -0.31031146   0.6695802
#> enhancer_5 6.272792  -1.3330017   0.2277560  -0.1179701 -0.49401180  -0.4425472
#> enhancer_6 5.196522  -0.5280535   1.0842769   1.5746666 -0.09114279  -0.3974384

head(getFits_DNA(mpra_obj_1))
#>            dna_1.bc_1 dna_1.bc_2 dna_1.bc_3 dna_1.bc_4 dna_1.bc_5
#> enhancer_1  154.40052  154.40052  154.40052  154.40052  154.40052
#> enhancer_2  121.39977  121.39977  121.39977  121.39977  121.39977
#> enhancer_3  111.39886  111.39886  111.39886  111.39886  111.39886
#> enhancer_4  278.19692  278.19692  278.19692  278.19692  278.19692
#> enhancer_5   91.80001   91.80001   91.80001   91.80001   91.80001
#> enhancer_6  224.40057  224.40057  224.40057  224.40057  224.40057
head(getFits_DNA(mpra_obj_2))
#>            dna_1.bc_1 dna_1.bc_2 dna_1.bc_3 dna_1.bc_4 dna_1.bc_5 dna_2.bc_1
#> enhancer_1  282.27661   55.68729  195.68045  165.10092  149.13295  282.27661
#> enhancer_2   59.70397   74.99840   89.47053  112.91528  116.61925   59.70397
#> enhancer_3  366.31963  179.25664  110.24910   69.80934   95.91599  366.31963
#> enhancer_4   99.84725  178.76660   60.31570   73.20985  195.04332   99.84725
#> enhancer_5  139.74093  175.48371  124.19093   85.26623   89.76930  139.74093
#> enhancer_6  106.53449  315.05451  514.46904   97.25399   71.59537  106.53449
#>            dna_2.bc_2 dna_2.bc_3 dna_2.bc_4 dna_2.bc_5
#> enhancer_1   55.68729  195.68045  165.10092  149.13295
#> enhancer_2   74.99840   89.47053  112.91528  116.61925
#> enhancer_3  179.25664  110.24910   69.80934   95.91599
#> enhancer_4  178.76660   60.31570   73.20985  195.04332
#> enhancer_5  175.48371  124.19093   85.26623   89.76930
#> enhancer_6  315.05451  514.46904   97.25399   71.59537

^{Created on 2022-04-12 by the reprex package (v2.0.1)}

[taniafabo]

just wanted to follow up on this! running into this issue