when I calculate coefficient vectors for each group the results are not the same as when the coefficients estimated by DEseq · Issue #9 · tavareshugo/tutorial_DESeq2_contrasts (original) (raw)

Dear @tavareshugo,
I am currently conducting single-nuclei RNA sequencing and have been following your tutorial closely. However, I've encountered a discrepancy in my results when comparing the extraction of coefficients from the model versus calculating the coefficient vector for each group.
Here is the design model I've been using:
design <- ~ sex + PMI + batchlib + group_id + region + group_id:region

#pre-filtering keep <- rowSums(counts(dds)>=10)>= 0.9*37
dds <- dds[keep,] #assign ref level dds$group_id <- relevel(dds$group_id, ref="con") dds$region <- relevel(dds$region, ref = "ocu") #running DESeq dds <- DESeq(dds, fitType='local') resultsNames(dds) [1] "Intercept" "sex_2_vs_1" "PMI" "batchlib_2_vs_1"
[5] "batchlib_3_vs_1" "batchlib_4_vs_1" "batchlib_5_vs_1" "batchlib_6_vs_1"
[9] "batchlib_7_vs_1" "group_id_als_vs_con" "region_med_vs_ocu" "region_sc_vs_ocu"
[13] "region_sl_vs_ocu" "group_idals.regionmed" "group_idals.regionsc" "group_idals.regionsl"

med_als <- colMeans(mod_mat[dds$region == "med" & dds$group_id == "als", ]) med_con <- colMeans(mod_mat[dds$region == "med" & dds$group_id == "con", ])

res1 <- results(dds, contrast = med_als - med_con) res2 <- results(dds, contrast = list("region_med_vs_ocu"))

res1 and res2 are not the same. I am not sure what point I did not consider to get consistent results, with res1 I have way more number of DEGs compared to res2.
Another question; if I am interested in the effect of "als" in "ocu" region and I extract "group_id_als_vs_con" considering that I put the reference level for region as "ocu" I expect to get the effect of "als" in ocu. However, I am not sure if this gives me the effect in batch 1 as it is the ref level. I need it to regress batch.
Thank you so much for your input!
Paria