Analyzing scRNA seq data and want to analyze it's DEG with DESeq2

Question

[kayah]

I have four dataset which is young_male, young_female, old_male, old_female. And I want to analyze it with DESeq2. The question is which data should I use for DESeq2. If I analyze it with bulk RNA seq then I'm gonna use feature counts file but if I want to analyze sc RNA seq then what file should I use? for further explanation of my data, I already did 'RunHarmony' integration process. and I upload my code for more understanding! ps. I want to transform my data assay to SCT but I failed...

library(Seurat)
library(SeuratObject)
library(dplyr)
library(tidyverse)
library(patchwork)

#data loading
WAT_M_Y.data <- Read10X(data.dir = "~/Desktop/GSE137869/2ndtry/WAT-M-Y/") 

colnames(WAT_M_Y.data) <- gsub("_", "-", colnames(WAT_M_Y.data))

WAT_M_Y <- CreateSeuratObject(counts = WAT_M_Y.data, 
                              project = "WAT_M_Y",
                              min.cells = 3,
                              min.features = 200)

WAT_M_Y [["percent.mt"]]<-PercentageFeatureSet(WAT_M_Y, pattern =  "^Mt-")


WAT_M_O.data <- Read10X(data.dir = 
                          "~/Desktop/GSE137869/GSE137869_aging/WAT-M-O/") 
colnames(WAT_M_O.data) <- gsub("_", "-", colnames(WAT_M_O.data))
WAT_M_O <- CreateSeuratObject(counts = WAT_M_O.data, 
                              project = "WAT_M_O",
                              min.cells = 3,
                              min.features = 200)

WAT_M_O [["percent.mt"]]<-PercentageFeatureSet(WAT_M_O, pattern =  "^Mt-")

WAT_F_Y.data <- Read10X(data.dir = 
                          "~/Desktop/GSE137869/GSE137869_aging/WAT-F-Y/") 
colnames(WAT_F_Y.data) <- gsub("_", "-", colnames(WAT_F_Y.data))
WAT_F_Y <- CreateSeuratObject(counts = WAT_F_Y.data, 
                              project = "WAT_F_Y",
                              min.cells = 3,
                              min.features = 200)

WAT_F_Y [["percent.mt"]]<-PercentageFeatureSet(WAT_F_Y, pattern =  "^Mt-")


WAT_F_O.data <- Read10X(data.dir = 
                          "~/Desktop/GSE137869/GSE137869_aging/WAT-F-O/") 
colnames(WAT_F_O.data) <- gsub("_", "-", colnames(WAT_F_O.data))
WAT_F_O <- CreateSeuratObject(counts = WAT_F_O.data, 
                              project = "WAT_F_O",
                              min.cells = 3,
                              min.features = 200)

WAT_F_O [["percent.mt"]]<-PercentageFeatureSet(WAT_F_O, pattern =  "^Mt-")

young_1 <- subset(WAT_M_Y, subset  = nFeature_RNA > 200 & nFeature_RNA < 2500 & nCount_RNA<15000 & percent.mt < 10)
young_2 <- subset(WAT_F_Y, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & nCount_RNA<15000 & percent.mt < 10)
old_1 <- subset(WAT_M_O, subset  = nFeature_RNA > 200 & nFeature_RNA < 2500 & nCount_RNA<15000 & percent.mt < 10)
old_2 <- subset(WAT_F_O, subset  = nFeature_RNA > 200 & nFeature_RNA < 2500 & nCount_RNA<15000 & percent.mt < 10)

saveRDS(young_1, "seurat_object_young1.rds")
saveRDS(young_2, "seurat_object_young2.rds")
saveRDS(old_1, "seurat_object_old1.rds")
saveRDS(old_2, "seurat_object_old2.rds")

library(devtools) library(dplyr) library(ggplot2) library(Seurat) library(SeuratObject) library(sp) library(patchwork) library(harmony) library(Rcpp) ?library(SeuratWrappers) library(clustree) library(cowplot) library(RColorBrewer) library(pheatmap) install.packages("FindvariableFeatures")

young_1 <- readRDS('~/Desktop/GSE137869/2ndtry/seurat_object_young1.rds')
young_2 <- readRDS('~/Desktop/GSE137869/2ndtry/seurat_object_young2.rds')
old_1 <- readRDS('~/Desktop/GSE137869/2ndtry/seurat_object_old1.rds')
old_2 <- readRDS('~/Desktop/GSE137869/2ndtry/seurat_object_old2.rds')

young_1 <- NormalizeData(young_1)
young_2 <- NormalizeData(young_2)
old_1 <- NormalizeData(old_1)
old_2 <- NormalizeData(old_2)

library(glmGamPoi)

WAT <- merge(young_1, y = c(young_2, old_1, old_2), 
             add.cell.ids = c("young1", "young2", "old1", "old2"),
             project = "WAT")
WAT <- NormalizeData(WAT) 
WAT <- FindVariableFeatures(WAT,selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(WAT)
WAT <-ScaleData(WAT, features = all.genes)
View([e-mail])
#WAT <- SCTransform(WAT, vars.to.regress = c("percent.mt"))
WAT <- RunPCA(WAT, npcs = 50)
View([e-mail])
DimPlot(WAT, reduction = "pca")

table(WAT$type)

harmonized_WAT <- RunHarmony(WAT,
                             group.by.vars = "orig.ident",
                             reduction ="pca",
                             reduction.save = "harmony")

harmonized_WAT <- RunUMAP(harmonized_WAT, reduction = "harmony",dims = 1:20)

harmonized_WAT <- FindNeighbors(harmonized_WAT, reduction = "harmony",dims = 1:20)

harmonized_WAT_0.6 <- FindClusters(harmonized_WAT, resolution = c(0.6))
DimPlot(harmonized_WAT_0.6, 
        label = TRUE, 
        split.by = "orig.ident")

And my data file looks like this

[kayah]

I do this for finding differently expressed genes between young and old group. My final destination is senescence marker(gene). I also want to draw Heatmap. Thank you so much