lmweber / diffcyt

R package for differential discovery analyses in high-dimensional cytometry data
MIT License
19 stars 12 forks source link

Create design & contrast matrix for differential state analysis #45

Open Stein-ErikG opened 2 years ago

Stein-ErikG commented 2 years ago

I have a mass cytometry dataset of samples collected from patients with AML (n=30). There is one sample from each patient, and the patients can be divided into two groups based on initial response to experimental therapy (21 responders vs 9 non-responders). I have been using the CATALYST pipeline and the diffcyt() function to do differential state analysis. Using this I find lot if very significant differences between the two groups. But when I use a different script that just do a lot of two sided T-tests (and then FDR adjustment), I find nothing significant.

When I use diffcyt() and my t-test script on a different and paired dataset, they find similar significant differences. So, I surmise that my design and/or contrast is not correctly set up? Can anyone help me? I hope my R output and code snippet below is informative! 😊

design <- createDesign(ei(sce), cols_design = “condition”) (Intercept) conditionRD 1 1 1 2 1 1 3 1 0 4 1 1 5 1 0 6 1 0 7 1 0 8 1 0 9 1 0 10 1 0 11 1 0 12 1 0 13 1 0 14 1 0 15 1 0 16 1 1 17 1 1 18 1 1 19 1 0 20 1 1 21 1 0 22 1 1 23 1 1 24 1 0 25 1 0 26 1 0 27 1 0 28 1 0 29 1 0 30 1 0 attr(,"assign") [1] 0 1 attr(,"contrasts") attr(,"contrasts")$condition [1] "contr.treatment"

contrast <- createContrast(c(0,1)) [,1] [1,] 0 [2,] 1

I am using the diffcyt() (https://rdrr.io/bioc/diffcyt/man/diffcyt.html) function from the diffcyt pipeline, with mostly default conditions. Ive pasted the code im using below:

design <- createDesignMatrix(ei(sce), cols_design = cond) contrast <- createContrast(c(0, 1))

res_DS <- diffcyt(sce, clustering_to_use = clust, analysis_type = "DS", method_DS = "diffcyt-DS-limma", design = design, contrast = contrast, verbose = FALSE) tbl_DS <- rowData(res_DS$res)

plotDiffHeatmap(sce, tbl_DS, all = TRUE, top_n = 30, fdr = 0.1, sort_by = "padj", col_anno = cond, normalize =TRUE)

I assume that “diffcyt-DS-limma” uses limma::lmFit. It looks to me that diffcyt() (https://rdrr.io/bioc/diffcyt/src/R/diffcyt_wrapper.R) uses testDS_limma.R with my input as:

res <- testDS_limma(d_counts, d_medians, design, contrast, block_id, trend, weights, markers_to_test, min_cells, min_samples, plot, path)

Looking further into testDS_limma() it seems like it only usable for paired data:

estimate correlation between paired samples

(note: paired designs only; >2 measures per sample not allowed)

But somehow runs on my data without error? I must be overlooking something…