Hi everybody. I'm working with RNASeq data and I am relatively new in that field. I would like to know if my approach is correct or not (or how I can improve it). Assuming that I have 3 different condition or treatment (high, low and mut) what I did is:

dds <- DESeqDataSetFromMatrix(
  countData = count.table,
  colData = data.frame(treatment=treatment),
  design = ~ treatment)

vsd <- varianceStabilizingTransformation(dds, blind=TRUE)
vst <- assay(vsd)
colnames(vst) <- colnames(count.table)
# To avoid lengthy unnecessary explanations, we simply shift all values so that non expressed genes have a vst
#value of 0
vst <- vst - min(vst)
dds1 <- DESeq(dds)

res.highCNA.low_CNA <- results(dds1, contrast=c("treatment","High_CNA", "low_CNA"))
res.highCNA.mut <- results(dds1, contrast=c("treatment","High_CNA", "mut"))
res.lowCNA.mut <- results(dds1, contrast=c("treatment", "low_CNA", "mut"))

#to select the 500 most expressed genes 

mat.high.low <- assay(vsd)[head(order(res.highCNA.low_CNA$padj), 500), ]
pheatmap(mat.high.low) 

mat.high.mut <- assay(vsd)[head(order(res.highCNA.mut$padj), 500), ]
pheatmap(mat.high.mut) 

mat.low.mut <- assay(vsd)[head(order(res.lowCNA.mut$padj), 500), ]
pheatmap(mat.low.mut)

Thanks in advance!!!