After performing DESEQ2 on my data, I could able to plot MA, PCA, and EnhancedVolcano plots including a (.xlsx) file consisting of log fold change ratios, base mean values, and p- values adjusted, and normal p-values were obtained. I would like to now generate a heatmap.

This is my R-script until now:

```
countData <- read.table("gene_count_matrix.csv", header = TRUE, sep = ",", row.names = 1)
head(countData)
metaData <- read.table("phenodata.csv", header = TRUE, sep = ",")
head(metaData)
#Deseq2
library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData=countData, colData=metaData, design=~stage_condition)
dds <- DESeq(dds)
dds <- dds[rowSums(counts(dds)) > 0]
res <- results(dds)
head(res)
summary(res)
res <- res[order(res$padj),]
head(res)
resSig <- res[ which(res$pvalue < 0.05)]
res_lfc <- subset(resSig, abs(log2FoldChange) > 2)
head(res_lfc)
#MA plot
plotMA(res, ylim=c(-2,2))
#plot for PCA
log_dds<-rlog(dds)
plotPCAWithSampleNames(log_dds, intgroup="treatment", ntop=40000)
#volcano plot
library("EnhancedVolcano")
EnhancedVolcano(res,lab = rownames(res),x = 'log2FoldChange',y = 'pvalue',labSize = 0, pCutoff = 0.05,FCcutoff = 1,xlim=c(-30,12))
plot(EnhancedVolcano)
```

Please format the code in your posts appropriately in the future using the code formatting button (the one with 1s and 0s) for clarity. I have done it for you this time.

Thanks so much, Sir!