Hi all

I would like to save axes scores for individuals samples from the MDS implemented in edgeR (i.e., the x,y coordinates of each sample that are displayed in the MDS plot). Does anyone have any advice on how to do this?

I'd appreciate any help.

Cheers, Michi

If you look at the help page for plotMDS(), you will see this:

Value:

     A plot is created on the current graphics device.

     An object of class ‘"MDS"’ is invisibly returned. This is a list
     containing the following components:

distance.matrix: numeric matrix of pairwise distances between columns
          of ‘x’

cmdscale.out: output from the function ‘cmdscale’ given the distance
          matrix

dim.plot: dimensions plotted

       x: x-xordinates of plotted points

       y: y-cordinates of plotted points

gene.selection: gene selection method

Whenever you see something that is 'invisibly' returned, this means that you can save the output, but if you don't, it doesn't get dumped onto your screen. Now if we run the example,

> example(plotMDS)

pltMDS> # Simulate gene expression data for 1000 probes and 6 microarrays.
pltMDS> # Samples are in two groups
pltMDS> # First 50 probes are differentially expressed in second group
pltMDS> sd <- 0.3*sqrt(4/rchisq(1000,df=4))

pltMDS> x <- matrix(rnorm(1000*6,sd=sd),1000,6)

pltMDS> rownames(x) <- paste("Gene",1:1000)

pltMDS> x[1:50,4:6] <- x[1:50,4:6] + 2

pltMDS> # without labels, indexes of samples are plotted.
pltMDS> mds <- plotMDS(x,  col=c(rep("black",3), rep("red",3)) )
Hit <Return> to see next plot:

pltMDS> # or labels can be provided, here group indicators:
pltMDS> plotMDS(mds,  col=c(rep("black",3), rep("red",3)), labels= c(rep("Grp1",3), rep("Grp2",3)))
Hit <Return> to see next plot:

You can see that the invisible output for plotMDS is captured in an object called 'mds', which we can then look at:

> mds
An object of class MDS
$dim.plot
[1] 1 2

$distance.matrix
          [,1]      [,2]      [,3]      [,4]      [,5] [,6]
[1,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000    0
[2,] 0.7818712 0.0000000 0.0000000 0.0000000 0.0000000    0
[3,] 0.8015551 0.7439845 0.0000000 0.0000000 0.0000000    0
[4,] 0.9851350 0.9575869 0.9505156 0.0000000 0.0000000    0
[5,] 0.9739221 0.9731381 0.9786097 0.7355927 0.0000000    0
[6,] 0.9515537 0.9769885 0.9814324 0.7991225 0.7772792    0

$cmdscale.out
           [,1]        [,2]
[1,] -0.3506229  0.34937830
[2,] -0.3805315 -0.12195160
[3,] -0.3725708 -0.22330723
[4,]  0.3606187 -0.25868024
[5,]  0.3937108 -0.06231905
[6,]  0.3493958  0.31687982

$top
[1] 500

$gene.selection
[1] "pairwise"

$x
[1] -0.3506229 -0.3805315 -0.3725708  0.3606187  0.3937108  0.3493958

$y
[1]  0.34937830 -0.12195160 -0.22330723 -0.25868024 -0.06231905  0.31687982

To answer the follow-up question from phalange.re:

No, you cannot retrieve the 'stress' value. plotMDS() uses classical MDS implemented in cmdscale() which does not make use of the stress loss function.

If you want to optimize 'stress', you would need to use a different MDS method, see:

  https://stat.ethz.ch/pipermail/r-help/2005-November/082042.html