I ll recommend you to use Seurat V5. You can use harmony as below before looking your gene of interest expression in vln plot and preforming Findmarkers-
seuratObj <- IntegrateLayers(object = seuratObj, method = HarmonyIntegration, orig.reduction = "pca", new.reduction = "integrated.harmony", verbose = FALSE)
seuratObj <- FindNeighbors(seuratObj, reduction = "integrated.harmony", dims = 1:30)
seuratObj <- FindClusters(seuratObj, resolution = 2, cluster.name = "harmony_clusters")
seuratObj <- RunUMAP(seuratObj, reduction = "integrated.harmony", dims = 1:30, reduction.name = "umap.harmony")
1. To display expression of your gene (for eg. CD8A in this case) of interest in Voiln plot, you can do sth like this-
VlnPlot(seuratObj, features = "CD8A", group.by = "harmony_clusters") + NoLegend() + ggtitle("CD8A - harmony Clusters")
2. You can run FindMarkers in cluster1 as below-
Idents(seuratObj)="harmony_clusters"
cluster1.markers <- FindMarkers(seuratObj, ident.1 = 1, min.pct = 0.25)