Hi there, I am doing DE analysis using DESeq2 and would like to know your comments regarding my experimental design and how to extract the contrast I'm interested in. My dataset include two different conditions (case & control) and 3 different regions (1, 2, and 3). each affected status include all 3 regions. I would like to see how the gene expression change in each region. First, I considered the below model:

''' dds <- DESeqDataSetFromMatrix(cluster_counts, colData = cluster_metadata, design = ~ sex+PMI+batchlib++individual+nuclei+group_id) ''' in this model the group_id is the combination of status and region: als_ocu, als_med, als_sc, con_ocu, con_med, and con_sc. I would like to compare case vs control in each region which will be 3 different contrasts: als_ocu vs con_ocu, als_med vs con_med, als_sc vs con_sc. (also:case, con:control, ocu, med, and sc are three different regions) But, I realized that I need to count for interaction of status and region in my model as I expect to see different pattern in different regions. So I did below model: ''' dds <- DESeqDataSetFromMatrix(cluster_counts, colData = cluster_metadata, design = ~ sex+PMI+batchlib+nuclei+status+region+status:region) dds <- DESeq(dds) '''

The problem I have with this new model is that I cannot do the contrasts I would like to do. Below is what resultsNames(dds) returns:

''' [1] "Intercept" "sex_2_vs_1" "PMI_3040_vs_2030" "PMI_7_vs_2030"
[5] "libsize_L_vs_H" "libsize_M_vs_H" "batchlib_2_vs_1" "batchlib_3_vs_1"
[9] "batchlib_4_vs_1" "batchlib_5_vs_1" "batchlib_6_vs_1" "batchlib_7_vs_1"
[13] "status_als_vs_con" "region_med_vs_ocu" "region_sc_vs_ocu"
[16] "statusals.regionmed" "statusals.regionsc"
'''

I was wondering how should I define my contrast for each of the comparison:

1. als_ocu vs con_ocu
2. als_med vs con_med
3. als_sc vs con_sc considering that my reference level for status is con and for region is ocu. I considered below contrast. However, I am not sure for example if in the first one not mentioning about region whether it shows effect of als vs control in specific region, ocu, which is set as base level. Please let me know which resultsNames I need to use for my comparisons.

ocu in alsvs con

res <- results(dds, contrast = c("status","als","con"), alpha = 0.05)

med in als vs con

res <- results(dds, contrast=list(c("status_als_vs_con", "statusals.regionmed")))

resultsNames(dds)

sc in als vs con

res <- results(dds, contrast=list(c("status_als_vs_con", "statusals.regionsl")))

Moreover, I was wondering what does each of "statusals.regionmed" "statusals.regionsc" regarding integration means?

May last step of analysis is that after getting these three sets of DEGs I would like to correlate first set to second and third set of DEGs separately (1 vs 2, 1 vs 3) to see how DEGs in region 1 (als vs con) correlate with other two regions DEGs. Do you think it makes more sense to create a contrast for this? and what would be the best contrast?

sessionInfo( ) R version 4.1.2 (2021-11-01) Platform: x86_64-pc-linux-gnu (64-bit) Running under: Rocky Linux 8.7 (Green Obsidian)

Matrix products: default BLAS/LAPACK: /cvmfs/soft.computecanada.ca/easybuild/software/2020/Core/flexiblas/3.0.4/lib64/libflexiblas.so.3.0

locale: [1] LC_CTYPE=en_CA.UTF-8 LC_NUMERIC=C LC_TIME=en_CA.UTF-8
[4] LC_COLLATE=en_CA.UTF-8 LC_MONETARY=en_CA.UTF-8 LC_MESSAGES=en_CA.UTF-8
[7] LC_PAPER=en_CA.UTF-8 LC_NAME=C LC_ADDRESS=C
[10] LC_TELEPHONE=C LC_MEASUREMENT=en_CA.UTF-8 LC_IDENTIFICATION=C

attached base packages: [1] splines stats4 stats graphics grDevices utils datasets methods base

other attached packages: [1] DESeq2_1.34.0 magrittr_2.0.3 Matrix.utils_0.9.8
[4] muscat_1.8.2 limma_3.50.3 metap_1.8
[7] patchwork_1.1.2 RCurl_1.98-1.12 scales_1.2.1
[10] lubridate_1.9.2 forcats_1.0.0 stringr_1.5.0
[13] purrr_1.0.1 readr_2.1.4 tidyr_1.3.0
[16] tibble_3.2.1 tidyverse_2.0.0 SingleCellExperiment_1.16.0 [19] SummarizedExperiment_1.24.0 GenomicRanges_1.46.1 GenomeInfoDb_1.30.1
[22] IRanges_2.28.0 S4Vectors_0.32.4 MatrixGenerics_1.6.0
[25] matrixStats_0.63.0 RColorBrewer_1.1-3 reshape_0.8.9
[28] cowplot_1.1.1 gplots_3.1.3 pheatmap_1.0.12
[31] monocle_2.22.0 DDRTree_0.1.5 irlba_2.3.5.1
[34] VGAM_1.1-8 ggplot2_3.4.2 Biobase_2.54.0
[37] BiocGenerics_0.40.0 Matrix_1.5-4 dplyr_1.1.1
[40] future_1.32.0 SeuratObject_4.1.3 Seurat_4.3.0

loaded via a namespace (and not attached): [1] scattermore_0.8 bit64_4.0.5 multcomp_1.4-23
[4] DelayedArray_0.20.0 data.table_1.14.8 KEGGREST_1.34.0
[7] doParallel_1.0.17 generics_0.1.3 ScaledMatrix_1.2.0
[10] RhpcBLASctl_0.23-42 TH.data_1.1-2 RSQLite_2.3.1
[13] RANN_2.6.1 combinat_0.0-8 bit_4.0.5
[16] tzdb_0.3.0 mutoss_0.1-13 spatstat.data_3.0-1
[19] httpuv_1.6.9 viridis_0.6.2 hms_1.1.3
[22] promises_1.2.0.1 fansi_1.0.4 progress_1.2.2
[25] caTools_1.18.2 igraph_1.4.2 DBI_1.1.3
[28] geneplotter_1.72.0 htmlwidgets_1.6.2 sparsesvd_0.2-2
[31] spatstat.geom_3.1-0 ellipsis_0.3.2 backports_1.4.1
[34] annotate_1.72.0 aod_1.3.2 deldir_1.0-6
[37] sparseMatrixStats_1.6.0 vctrs_0.6.1 ROCR_1.0-11
[40] abind_1.4-5 cachem_1.0.7 withr_2.5.0
[43] grr_0.9.5 progressr_0.13.0 sctransform_0.3.5
[46] prettyunits_1.1.1 goftest_1.2-3 mnormt_2.1.1
[49] cluster_2.1.2 lazyeval_0.2.2 crayon_1.5.2
[52] genefilter_1.76.0 spatstat.explore_3.1-0 labeling_0.4.2
[55] edgeR_3.36.0 pkgconfig_2.0.3 slam_0.1-50
[58] qqconf_1.3.2 vipor_0.4.5 nlme_3.1-153
[61] blme_1.0-5 rlang_1.1.0 globals_0.16.2
[64] lifecycle_1.0.3 miniUI_0.1.1.1 sandwich_3.0-2
[67] mathjaxr_1.6-0 rsvd_1.0.5 polyclip_1.10-4
[70] lmtest_0.9-40 boot_1.3-28 zoo_1.8-12
[73] beeswarm_0.4.0 ggridges_0.5.4 GlobalOptions_0.1.2
[76] png_0.1-8 viridisLite_0.4.1 rjson_0.2.21
[79] bitops_1.0-7 KernSmooth_2.23-20 Biostrings_2.62.0
[82] blob_1.2.4 DelayedMatrixStats_1.16.0 shape_1.4.6
[85] parallelly_1.35.0 spatstat.random_3.1-4 beachmat_2.10.0
[88] memoise_2.0.1 plyr_1.8.8 ica_1.0-3
[91] zlibbioc_1.40.0 compiler_4.1.2 HSMMSingleCell_1.14.0
[94] plotrix_3.8-2 clue_0.3-64 lme4_1.1-32
[97] fitdistrplus_1.1-8 cli_3.6.1 XVector_0.34.0
[100] lmerTest_3.1-3 listenv_0.9.0 pbapply_1.7-0
[103] TMB_1.9.4 MASS_7.3-54 tidyselect_1.2.0
[106] stringi_1.7.12 densityClust_0.3.2 BiocSingular_1.10.0
[109] locfit_1.5-9.7 ggrepel_0.9.3 grid_4.1.2
[112] tools_4.1.2 timechange_0.2.0 future.apply_1.10.0
[115] parallel_4.1.2 circlize_0.4.15 foreach_1.5.2
[118] gridExtra_2.3 farver_2.1.1 Rtsne_0.16
[121] digest_0.6.31 FNN_1.1.3.2 shiny_1.7.4
[124] qlcMatrix_0.9.7 Rcpp_1.0.10 broom_1.0.4
[127] scuttle_1.4.0 later_1.3.0 RcppAnnoy_0.0.20
[130] httr_1.4.5 AnnotationDbi_1.56.2 ComplexHeatmap_2.10.0
[133] Rdpack_2.4 colorspace_2.1-0 XML_3.99-0.14
[136] tensor_1.5 reticulate_1.28 uwot_0.1.14
[139] sn_2.1.1 spatstat.utils_3.0-2 scater_1.22.0
[142] sp_1.6-0 multtest_2.50.0 plotly_4.10.1
[145] xtable_1.8-4 jsonlite_1.8.4 nloptr_2.0.3
[148] R6_2.5.1 TFisher_0.2.0 pillar_1.9.0
[151] htmltools_0.5.5 mime_0.12 glue_1.6.2
[154] fastmap_1.1.1 minqa_1.2.5 BiocParallel_1.28.3
[157] BiocNeighbors_1.12.0 codetools_0.2-18 mvtnorm_1.1-3
[160] utf8_1.2.3 lattice_0.20-45 spatstat.sparse_3.0-1
[163] pbkrtest_0.5.2 numDeriv_2016.8-1.1 ggbeeswarm_0.7.1
[166] leiden_0.4.3 gtools_3.9.4 survival_3.5-5
[169] glmmTMB_1.1.7 docopt_0.7.1 fastICA_1.2-3
[172] munsell_0.5.0 GetoptLong_1.0.5 GenomeInfoDbData_1.2.7
[175] iterators_1.0.14 variancePartition_1.24.1 reshape2_1.4.4
[178] gtable_0.3.3 rbibutils_2.2.13

Thanks, Paria ```