I am interested in using variance stabilized values of my data (that are the outputs of DESeq2's varianceStabilizingTransformation function) for machine learning and other applications. I ran the following:

dds <- DESeqDataSetFromMatrix(countData = counts, colData = metadata, design = ~Institution + Condition)
vsd <- varianceStabilizingTransformation(dds, blind = FALSE)

Although I tried to include the batch (Institution) in the design argument, I still see batches in my PCA plot. This section of the vignette described a solution for this that I have pasted below:

mat <- assay(vsd)
mat <- limma::removeBatchEffect(mat, vsd$batch)
assay(vsd) <- mat
plotPCA(vsd)

Can the output of removeBatchEffect in this code now be used for all other downstream applications, and not just for visualization?

Thanks very much.

Yes, I have used this kind of workflow to remove unwanted variation from VST data.