I am wondering whether my approach is correct or not as I am getting some different result from Limma.  I have two different stress conditions with three replicates comparing with control. The mapping was done using STAR and count data was generated through HTSeq. Is it also possible to compare within the stressA_vs_stressB? Thank you.

library(DESeq2)
library(apeglm)
sampleFiles <- grep(".counts",list.files(directory),value=TRUE)
sampleName<-c("control1", "control2", "control3", "stressA1", "stressA2", "stressA3", "stressB1", "stressB2", "stressB3")
table<-c("control", "control", "control", "stressA", "stressA", "stressA", "stressB", "stressB", "stressB")
sampleTable<-data.frame(sampleName=sampleName, fileName=sampleFiles, condition=table)
ddsHTSeq<- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable,
                                      directory = directory,
                                      design= ~ condition)

> sampleTable
sampleName           fileName condition
1   control1   1S1_htseq.counts   control
2   control2   1S2_htseq.counts   control
3   control3   1S3_htseq.counts   control
4   stressA1   3S1_htseq.counts   stressA
5   stressA2   3S2_htseq.counts   stressA
6   stressA3   3S3_htseq.counts   stressA
7   stressB1 D20S1_htseq.counts   stressB
8   stressB2 D20S2_htseq.counts   stressB
9   stressB3 D20S3_htseq.counts   stressB

dds<-DESeq(ddsHTSeq)
res<-results(dds)
keep<-rowSums(counts(dds))>=10
ddss<-dds[keep,]
resultsNames(ddss)
[1] "Intercept"                    "condition_stressA_vs_control" "condition_stressB_vs_control"
resLFC<-lfcShrink(ddss,coef="condition_stressA_vs_control", type="apeglm")
ress= subset(resLFC, padj<0.05)
ress <- ress[order(ress$padj),]
resLFC1<-lfcShrink(ddss,coef="condition_stressB_vs_control", type="apeglm")
ress1= subset(resLFC1, padj<0.05)
ress1 <- ress1[order(ress1$padj),]

We have some information about how to do this in the vignette, in the extended section on shrinkage estimators.