add 2 files

code/README.md

+# Code
+
+Save command-line scripts and shared R code here.

code/maybe_diffcyt_bug.R

+
+source(here::here("code/my_diffcyt.R"))
+
+library(diffcyt)
+library(CATALYST)
+
+rdsdir <- here::here("output")
+sce <- readRDS(file.path(rdsdir, "sce_tumor_merged.rds"))
+
+
+(ei <- metadata(sce)$experiment_info)
+(mm <- model.matrix(~treatment, data=ei))
+
+da_res <- diffcyt(sce, design = mm, contrast = c(0, -1),
+                  analysis_type = "DA", 
+                  method_DA = "diffcyt-DA-voom",
+                  clustering_to_use = "merging", 
+                  verbose = TRUE)
+
+assay(da_res$d_counts)
+
+# repeat after outlier removal
+outliers <- c("5__untreated","6_7__untreated")
+sce_s <- filterSCE(sce, !(sample_id %in% outliers))
+
+
+(ei <- metadata(sce_s)$experiment_info)
+(mm <- model.matrix(~treatment, data=ei))
+
+da_res <- diffcyt(sce_s, design = mm, contrast = c(0, -1),
+                  analysis_type = "DA", 
+                  method_DA = "diffcyt-DA-voom",
+                  clustering_to_use = "merging", 
+                  verbose = TRUE)
+
+z <- my_diffcyt(sce_s, design = mm, contrast = c(0, -1),
+                  analysis_type = "DA", 
+                  method_DA = "diffcyt-DA-voom",
+                  clustering_to_use = "merging", 
+                  verbose = TRUE)
+
+table(cluster_ids(sce, "merging"), sce$sample_id)
+assay(z)
+table(cluster_ids(sce_s, "merging"), sce_s$sample_id)
+
+sessionInfo()