The function voom of the limma package returns a EList object with, among others:

• matrix of normalized eset, corrected for library size, zero counts and log2 transformed
• matrix of inverse weights (for each combination of features and samples)

There is an example p.70 of the vignette of the limma package in which the entire object is given to limma, which makes use and combines both items to correct for mean-variance relationship:

v <- voom(counts,design,plot=TRUE,normalize="quantile")
fit <- lmFit(v,design)
fit <- eBayes(fit)

Is there a logical way to obtain the pre-processed data with mean-variance correction, e.g. product of the normalized expression matrix and inverse weights?

Our interest is to visualize the data after mean-variance correction but before differential analysis, e.g. with principal component analysis.

voom doesn't correct the mean:variance relationship, it rather models it. This is then used to generate correct statistics for detecting differential expression. In short, voom cannot provide you the type of data you are look for.

To get what you are after, the edgeR/limma folks suggest calling edgeR::cpm with a prior.count set between 3-5 (and log=TRUE) to get a (more) variance stabilized version of your count data, while the DESeq(2) camp provides the rlog and vst transformations.

I suspect you will find that simply using

plotMDS(v)

will be fine. While this doesn't explicitly adjust for the mean-variance relationship, we have found it to be satisfactory enough.

Usually, as Steve Lianoglou has mentioned, I recommend using cpm() on TMM normalized counts for exploratory purposes, but in your case it is more important that the plot reflects the quantile normalization than to adjust for the mean-variance relationship. The use of quantile normalization means that the various edgeR or DESeq2 solutions will be not quite right for you as they don't reflect the same normalization step.