Hi,

Does anyone know how to get the loading scores for MDS plot in limma?

So, for example in the following example code in the plotMDS in limma, I'd like to know the contribution of each of the simulated gene in x for each dimension and which gene contribute the most.

# Simulate gene expression data for 1000 probes and 6 microarrays.
# Samples are in two groups
# First 50 probes are differentially expressed in second group
sd <- 0.3*sqrt(4/rchisq(1000,df=4))
x <- matrix(rnorm(1000*6,sd=sd),1000,6)
rownames(x) <- paste("Gene",1:1000)
x[1:50,4:6] <- x[1:50,4:6] + 2
# without labels, indexes of samples are plotted.
mds <- plotMDS(x,  col=c(rep("black",3), rep("red",3)) )
# or labels can be provided, here group indicators:
plotMDS(mds,  col=c(rep("black",3), rep("red",3)), labels= c(rep("Grp1",3), rep("Grp2",3)))

I found this discussion LIMMA: plotMDS but I cannot find a clear answer on how to do this.

Thank you in advance!

plotMDS doesn't use gene loading scores, so there is no measure of the contribution of each gene. Indeed, plotMDS may use different genes to distinguish each pair of samples, so the relative importance of each gene may vary from one sample pair to another.

I said the same thing 11 years ago in the answer that you link to in your question. Note also that the plotMDS code has changed in the past 11 years (the function no longer uses cmdscale) so the code shown in the 11-year-old thread is out-of-date.

The limma documentation describes what each function does. You can type

names(mds)

to see all the components that are returned by the function. If you type help(plotMDS), the help page will give a description of each component of mds.

If you simply want to find genes are correlated with the x or y-axes of the MDS plot, you can do that by:

design <- cbind(Intercept=1,xaxis=mds$x,yaxis=mds$y)
fit <- lmFit(x, design)
fit <- eBayes(fit)
topTable(fit, coef="xaxis")
topTable(fit, coef="yaxis")

Using the plotMDS data example:

> example(plotMDS)
> design <- cbind(Intercept=1,xaxis=mds$x,yaxis=mds$y)
> fit <- lmFit(x, design)
> fit <- eBayes(fit)
> topTable(fit, coef="xaxis")
           logFC   AveExpr        t      P.Value    adj.P.Val        B
Gene 43 3.138563 1.0178138 11.51615 3.493344e-06 0.0007955732 5.127485
Gene 20 2.870733 1.0552325 11.13657 4.480185e-06 0.0007955732 4.896124
Gene 39 2.983242 0.8659482 11.13176 4.494566e-06 0.0007955732 4.893129
Gene 44 3.209721 1.4578645 10.85037 5.432145e-06 0.0007955732 4.715350
Gene 49 2.883946 0.9043585 10.71216 5.971502e-06 0.0007955732 4.626032
Gene 10 2.994383 0.9223217 10.70109 6.017229e-06 0.0007955732 4.618820
Gene 8  2.733811 1.0181440 10.65553 6.209599e-06 0.0007955732 4.589048
Gene 12 2.934456 0.9709351 10.61996 6.364586e-06 0.0007955732 4.565701
Gene 38 2.684030 1.0137147 10.43894 7.223665e-06 0.0008026294 4.445457
Gene 2  2.707727 0.7085983 10.29002 8.028297e-06 0.0008028297 4.344738
> topTable(fit, coef="yaxis")
             logFC     AveExpr         t      P.Value adj.P.Val           B
Gene 346 -4.009426 -0.01718160 -6.570801 0.0001922147 0.1506713  0.64234269
Gene 331 -2.736757  0.18273787 -5.942688 0.0003749138 0.1506713  0.18865372
Gene 412  5.630653 -0.23139630  5.774588 0.0004520139 0.1506713  0.05680294
Gene 698  3.527671  0.34003146  4.454978 0.0022419641 0.3540862 -1.14999378
Gene 541 -2.293136 -0.10443727 -4.427329 0.0023246646 0.3540862 -1.17873556
Gene 979  1.809885 -0.42864149  4.249250 0.0029434847 0.3540862 -1.36735877
Gene 895  1.630157  0.04463720  4.202459 0.0031342481 0.3540862 -1.41792657
Gene 121 -2.069048  0.27158813 -4.195380 0.0031642528 0.3540862 -1.42561264
Gene 477 -2.389437 -0.02027189 -4.083645 0.0036814134 0.3540862 -1.54820289
Gene 165 -1.901643 -0.04190741 -4.070279 0.0037491460 0.3540862 -1.56302551
> summary(decideTests(fit))
       Intercept xaxis yaxis
Down           0     0     0
NotSig       952   951  1000
Up            48    49     0

The results tell you that about 50 genes contribute to the x-axis while no genes contribute significantly to the y-axis, and that is a correct reflection of the simulation.

> sessionInfo()
R version 4.3.1 (2023-06-16 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)

Matrix products: default


locale:
[1] LC_COLLATE=English_Australia.utf8  LC_CTYPE=English_Australia.utf8   
[3] LC_MONETARY=English_Australia.utf8 LC_NUMERIC=C                      
[5] LC_TIME=English_Australia.utf8    

time zone: Australia/Melbourne
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] limma_3.57.6

loaded via a namespace (and not attached):
[1] compiler_4.3.1 tools_4.3.1