Perfect, thanks for bearing with me. I'll summarized what we discussed here for anyone who will read this post later.

In my pseudo-bulk DEG analysis, I have samples from four tissue regions of both cases and controls. The variables I am including in my model are as follows: design = ~ sex + PMI + batchlib + group_id + region + group_id:region

I need to test different hypothesis:

First, I need to measure the effect of disease in each region separately: null hypothesis: Disease has no effect on different regions. The below analysis are done after aggregation of data to the sample level.

#Create DESeq2 object
dds <- DESeqDataSetFromMatrix(cluster_counts, 
                              colData = cluster_metadata, 
                              design = ~ sex + PMI + batchlib + group_id + region + group_id:region)

#pre-filtering (37 is sample size)
keep <- rowSums(counts(dds)>=10)>= 0.9*37   
dds <- dds[keep,]

#determining ref level
#to measure the effect of disease in s=each region change the ref level to that region.
dds$group_id <- relevel(dds$group_id, ref="con")
dds$region <- relevel(dds$region, ref = "ocu")

#Running DESeq2
dds <- DESeq(dds, fitType='local') 
#then, I extract the coefficient "group_id_als_vs_con".

Second, I need to obtain region specific differences which means how much the effect of disease is different in regions compared to ref region (which is ocu here). null hypothesis: There is no region specific difference in regions (med, sc, sl) compared to ocu

#Create DESeq2 object excluding region
dds <- DESeqDataSetFromMatrix(cluster_counts, 
                              colData = cluster_metadata, 
                              design = ~ sex + PMI + batchlib + group_id + group_id:region)

#pre-filtering (37 is sample size)
keep <- rowSums(counts(dds)>=10)>= 0.9*37   
dds <- dds[keep,]

#determining ref level
#to measure the effect of disease in s=each region change the ref level to that region.
dds$group_id <- relevel(dds$group_id, ref="con")
dds$region <- relevel(dds$region, ref = "ocu")

#Running DESeq2
dds <- DESeq(dds, fitType='local') 
#then extracting group_idals.regionsl, or group_idals.regionsc, or group_idals.regionmed

Third, I would like to know how the impact of the disease is different in all the regions together compared to the ref region. The difference with the second approach is that here we look for a patterns, like genes that show different pattern in all regions and not looking at regions two by two. I used LRT which tests multiple terms at once, for example testing 3 or more levels of a factor at once, or all interactions between two variables. null hypothesis: there is no genes whose expression is different in all regions.

#Create DESeq2 object with full model
dds <- DESeqDataSetFromMatrix(cluster_counts, 
                              colData = cluster_metadata, 
                              design = ~ sex + PMI + batchlib + group_id + region + group_id:region)

# pre-filtering
keep <- rowSums(counts(dds)>=10)>= 0.9*37   
dds <- dds[keep,]

dds$group_id <- relevel(dds$group_id, ref="con")
dds$region <- relevel(dds$region, ref = "ocu")


#Running DESeq2 with reduced model
dds <- DESeq(dds, test="LRT", reduced = ~ sex + PMI + batchlib + group_id + region)

#then extracting any coefficient will output the effect of the variable which was reduced

Fourth, to measure the overall effect of the disease I removed interaction term; this way, status effect represents the overall effect controlling for differences due to region.

#Create DESeq2 object with full model
dds <- DESeqDataSetFromMatrix(cluster_counts, 
                              colData = cluster_metadata, 
                               design = ~ sex + PMI + batchlib + group_id + region)
# pre-filtering
keep <- rowSums(counts(dds)>=10)>= 0.9*37   
dds <- dds[keep,]

dds$group_id <- relevel(dds$group_id, ref="con")

dds <- DESeq(dds, fitType='local')
#then, extracting the "group_id_als_vs_con" coefficient.

If there is anything wrong please notify me dear Michael Love

Thanks,

Paria