diff --git a/code/README.md b/code/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..4c2402fd9d755278177f69df0c97cfb4fe8ccbf8
--- /dev/null
+++ b/code/README.md
@@ -0,0 +1,3 @@
+# Code
+
+Save command-line scripts and shared R code here.
diff --git a/code/maybe_diffcyt_bug.R b/code/maybe_diffcyt_bug.R
new file mode 100644
index 0000000000000000000000000000000000000000..34911efd4fd65cb95832b73c7b1d13f1c93b0186
--- /dev/null
+++ b/code/maybe_diffcyt_bug.R
@@ -0,0 +1,46 @@
+
+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()