Without knowing what platform gave rise to your expression data, I'd like to start by pointing you to the documentation for some popular methods of DE analysis.

From the edgeR manual (https://bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf):

edgeR is primarily intended for use with data including biological replication. Nevertheless, RNA-Seq and ChIP-Seq are still expensive technologies, so it sometimes happens that only one library can be created for each treatment condition. In these cases there are no replicate libraries from which to estimate biological variability. In this situation, the data analyst is faced with the following choices, none of which are ideal. We do not recommend any of these choices as a satisfactory alternative for biological replication.

From the DEseq2 manual (https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#can-i-use-deseq2-to-analyze-a-dataset-without-replicates):

Can I use DESeq2 to analyze a dataset without replicates? No. This analysis is not possible in DESeq2.

You can check out limma-voom (https://genomebiology.biomedcentral.com/articles/10.1186/gb-2014-15-2-r29) as well, but it also relies on varition for its models and I think it would be very hard to draw any conclusions.

I would treat this as a pilot study and just explore the differences in expression between your groups of interest, for example by plotting. Maybe only run DE analysis comparing types or treatments for Stage3 for now, which does have replicates?