Modifying R pairs() function

Question

[shania90.lk]

Hi to the experts,

I'm attempting to generate a pairs() plot using the plotMEpairs() function in WGCNA. My current code looks very simple as following;
plotMEpairs(MEs, traitsDat$blade_Na, cex.labels=12)

My MEs data frame looks like this

head(MEs)

  darkgoldenrod4 darkseagreen3  brown3 palevioletred3 darkolivegreen
1      0.1547894    0.15358785  0.00863432     0.21155555     0.14404473
2      0.2276100    0.18018804  0.10890920     0.27717332     0.14444103
3      0.1427315    0.14522530 -0.38375679    -0.11940695    -0.23971664
4      0.1709225    0.15518469  0.13983453     0.40029139    -0.07366996
5     -0.1231839    0.03224154  0.28591815    -0.04757507     0.06252794
6     -0.2057145   -0.05495319 -0.46779093    -0.38743274    -0.10632115

My traitsDat data frame looks like this

head(traitsDat)

blade_Na sheath_Na   blade_K sheath_K blade_K_Na sheath_K_Na SBK_SBNa
L1A-A-R 2351.8181      2887  823.9930      578  0.3503643   0.2002216     0.57
L2A-A-R 2573.5617      3309  823.6208      450  0.3200315   0.1361210     0.43
L3A-A-R 2339.7940      3101  659.0464      467  0.2816686   0.1506410     0.53
L4A-A-R 2242.5629      3080  685.5601      496  0.3057038   0.1609653     0.53
L1B-B-R 1315.1276      1313  932.4009     1210  0.7089814   0.9212283     1.30
L2B-B-R  627.4885       563 1098.9011     1439  1.7512690   2.5588874     1.46

I can come up with an output like this

Now, what I actually want is to

• match each of the diagonal histogram colour to the corresponding colnames(MEs) and,
• have the scatter plots in the upper triangle to be coloured as following; the first 4 values for a plot in one colour, 5-8 in another, 9-12 in another ... 20-24 in another (NROW(MEs) is 24).

I'm quite sure that this is possible through manipulating the panel.hist and/or panel.smooth. But I'm not able to figure it out after 2 days of struggling. Any help would be very much appreciated.

Thanks a lot in advance

Shani

[shania90.lk]

Hi Kevin,

Thank you for this. However, In the panel.hist() however an error occurs : Error: '...' used in an incorrect context.

Yes you are right pairs() is a very powerful function to master. Any help with this one for me? Thank you.

Answer 1

[Kevin Blighe]

To do this, I edited the base pairs function code and have hosted this on my GitHub page: https://github.com/kevinblighe/Biostars/blob/master/pairsKB.R

It was too much text to put into a single answer here.

In the code on GitHub, I have added comments at parts that I have edited to help you see what I've done. Just for a quick glimpse, I've also put the edited parts here (the first edited part is not necessary - I just don't like the ticks on the pairs plot):

pairsKB <- function (x, labels, panel=points, ..., horInd=1:nc, verInd=1:nc, lower.panel=panel, upper.panel=panel, diag.panel=NULL, text.panel=textPanel, label.pos=0.5+has.diag/3, line.main=3, cex.labels=NULL, font.labels=1, row1attop=TRUE, gap=1, log="") 
{
    ...
        if (side%%2L == 1L)
            ################################################
            #Edited next 3 lines, 27th November 2017 by Kevin Blighe
            #Added 'labels=FALSE, tick=FALSE,'
            ################################################
            Axis(x, side = side, xpd = xpd, labels=FALSE, tick=FALSE, ...)
        else
            Axis(y, side = side, xpd = xpd, labels=FALSE, tick=FALSE, ...)
    }

    ...
                    if (has.diag)
                        ################################################
                        #Edited next line, 27th November 2017 by Kevin Blighe
                        #Now passing column index, 'i', to diagonal panel function
                        ################################################
                        localDiagPanel(as.vector(x[, i]), i, ...)
    ...
                else if (i < j)
                    ################################################
                    #Edited next line, 27th November 2017 by Kevin Blighe
                    #Add col and pch arguments
                    ################################################
                    localLowerPanel(as.vector(x[, j]), as.vector(x[, i]), col=c(rep("royalblue", 2), rep("black", 2), rep("gold", 2)), pch=16, ...)

        else
            ################################################
            #Edited next line, 27th November 2017 by Kevin Blighe
            #Add col and pch arguments
            ################################################
            localUpperPanel(as.vector(x[, j]), as.vector(x[, i]), col=c(rep("royalblue", 2), rep("black", 2), rep("gold", 2)), pch=16, ...)
    ...
}

The localLowerPanel and localUpperPanel is where you can specify the colours of the data points in the scatter plots. Here, with your sample data, I've only got 6 values, so, I just create a colour vector of 6 with c(rep("royalblue", 2), rep("black", 2), rep("gold", 2)). Edit this to suit your own colour designations.

----------------------------------------

For the diagonal histogram colouring per WGCNA module, you need to edit the standard histogram function too. Here is the code:

################################################
#Edited next line, 27th November 2017 by Kevin Blighe
#Modified standard panel.hist function to additionally accept column index being plotted
################################################
panel.hist <- function(x, i, ...)
{
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(usr[1:2], 0, 1.5) )
    h <- hist(x, plot = FALSE)
    breaks <- h$breaks; nB <- length(breaks)
    y <- h$counts; y <- y/max(y)

    ################################################
    #Edited next line, 27th November 2017 by Kevin Blighe
    #Now colouring each histogram based on WGCNA module names
    ################################################
    rect(breaks[-nB], 0, breaks[-1], y, col=colnames(MEs)[i], ...)
}

This will be invoked by localDiagPanel(as.vector(x[, i]), i, ...) in the main pairsKB function. It's necessary to pass the 'i' variable, which contains the column index being plotted.

For this to work, column ordering in MEs must match that of traitsDat.

--------------------------------------

To then execute this:

pairsKB(traitsDat, diag.panel=panel.hist)

[Kevin Blighe]

PS - if you want to add the loess regression fitted red line (as in the figure to which you've linked), just let me know and I'll do it. plotMEpairs is just a light-weight wrapper of pairs and is not doing much extra

[shania90.lk]

I added a reply that supposed to be a comment for this post from you. Can you please have a look Kevin? Thank you.

[Kevin Blighe]

I have now been able to do what I believe you were aiming to do originally.

Firstly, you'll still require a further modified pairsKB function, again on my GitHub page. Only the following lines are now edited:

• 16-22 - to remove axis ticks and numbers (you may want these, in which case set labels=TRUE, tick=TRUE
• 147-151 - to pass an index value to the diagonal function

-----------------------

I also further edited your other functions (not all):

panel.smooth.ordered.categorical <- function (x, y, col = par("col"), bg = NA, pch = par("pch"), cex = 1, col.smooth = "red", span = 2/3, iter = 3,point.size.rescale = 1.0, ...)
{ 
    ################################################
    #Edited next line, 27th November 2017 by Kevin Blighe
    #A vector to be used for colouring the dots
    ################################################
    dotColours <- c(rep("royalblue", 2), rep("black", 2), rep("forestgreen", 2))

    require(reshape2)

    z <- merge(data.frame(x,y), melt(table(x ,y)),sort =F)$value

    z <- point.size.rescale*z/ (length(x))

    symbols( x, y, circles = z,
        #rep(0.1, length(x)),
        #sample(1:2, length(x), replace = T),

        ################################################
        #Edited next line, 27th November 2017 by Kevin Blighe
        #Supply the dotColours vector, and setting dot border to black
        ################################################
        inches=F, bg=dotColours, fg="black", add = T)

    #points(x, y, pch = pch, col = col, bg = bg, cex = cex)

    ok <- is.finite(x) & is.finite(y)

    if (any(ok))
        lines(stats::lowess(x[ok], y[ok], f = span, iter = iter), col = col.smooth, ...)
}

To suit the colour designations for the points, you'll have to set the dotColours vector manually, eg: dotColours <- c(rep("royalblue", 4), rep("black", 4), rep("forestgreen", 4), ...), for first 4 values in one colour, 5-8 in another, et cetera.

-----------------------

################################################
#Edited next line, 27th November 2017 by Kevin Blighe
#Modified standard panel.hist function to additionally accept column index being plotted
################################################
panel.hist <- function(x, i, ...)
{
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(usr[1:2], 0, 1.5) )
    h <- hist(x, plot = FALSE)
    breaks <- h$breaks; nB <- length(breaks)
    y <- h$counts; y <- y/max(y)

    ################################################
    #Edited next 4 lines, 27th November 2017 by Kevin Blighe
    #Now colouring each histogram based on WGCNA module names
    ################################################
    if ( is.na(colnames(MEs)[i]) )
        rect(breaks[-nB], 0, breaks[-1], y, col="white", ...)
    else
        rect(breaks[-nB], 0, breaks[-1], y, col=colnames(MEs)[i], ...)
}

Set the value of the histograms where there's no colour to be used to "white" (you can change this to anything), otherwise, set the colour as the name of the WGCNA module

-----------------------

pairs.ordered.categorical <- function(xx,...) 
{
    ################################################
    #Edited next line, 27th November 2017 by Kevin Blighe
    #Now invoking pairsKB, not pairs
    ################################################
    pairsKB(xx, diag.panel=panel.hist, lower.panel=panel.smooth.ordered.categorical, upper.panel=panel.cor.ordered.categorical, cex.labels=0.8,...)
}

Invokes my pairsKB function to get the histogram colouring correct,

-----------------------

traits_needed <-c("blade_Na","sheath_Na","blade_K","sheath_K","blade_K_Na")
needed <- c("darkgoldenrod4","darkseagreen3","brown3","palevioletred3","darkolivegreen")
pairs.ordered.categorical(cbind(MEs[, needed], traitsDat[, traits_needed]))

The colouring of the histograms will only work if MEs are positioned first in your object, but I think that's the way that you want.

There are undoubtedly more ways of doing this. The pairs function is quite powerful if one can grasp how it works. I devoted a full day to understanding its code, going line by line.

Kevin

[shania90.lk]

Hi Kevin,

Thank you for this. However, In the panel.hist() however an error occurs : Error: '...' used in an incorrect context.

Yes you are right pairs() is a very powerful function to master. Any help with this one for me? Thank you.

[Kevin Blighe]

Ah, I do apologise (I'm British - we apologise for everything), it was actually a formatting issue with how the system on Biostars interprets code. There's a glitch whereby a parenthesis can 'vanish'.

I have edited the code again. It was just the line if ( is.na(colnames(MEs)[i]) ) in the hist function.

I would also just be wary of the colour vector in panel.smooth.ordered.categorical. I'm just trying to double confirm that the symbols function does actually colour the dots in the order they appear in the data (one would expect that they do).

[Kevin Blighe]

It does appear to colour them according to how they appear in the data, which is good. It's just reversed because the data moves from positive to negative.

Answer 2

[shania90.lk]

What a fantastic effort. Greatly appreciate it.
However I would finally like to see a plot like this but with the histogram bars for the diagonal panel colored according to the mentioned colour. If there's no colour to be mentioned the histogram bars can be wither black or null. I'm extremely sorry for not properly formulating my question earlier to make this clearer. I'll add the code I've used to generate the plot I've attached in the image below:
panel.cor.ordered.categorical <- function(x, y, digits=2, prefix="", cex.cor) { usr <- par("usr"); on.exit(par(usr)) par(usr = c(0, 1, 0, 1))

r <- abs(cor(x, y, method = "pearson")) # notie we use pearson, non parametric correlation here r.no.abs <- cor(x, y, method = "pearson")

txt <- format(c(r.no.abs , 0.123456789), digits=digits)[1] txt <- paste(prefix, txt, sep="") if(missing(cex.cor)) cex <- 0.8/strwidth(txt)

test <- cor.test(x,y, method = "pearson") # borrowed from printCoefmat Signif <- symnum(test$p.value, corr = FALSE, na = FALSE, cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1), symbols = c("***", "**", "*", ".", " ")) text(0.5, 0.5, txt, cex = cex * r)
text(.8, .8, Signif, cex=cex, col=2)
}

panel.smooth.ordered.categorical <- function (x, y, col = par("col"), bg = NA, pch = par("pch"), cex = 1, col.smooth = "red", span = 2/3, iter = 3,point.size.rescale = 1.0, ...)
{
require(reshape2)
z <- merge(data.frame(x,y), melt(table(x ,y)),sort =F)$value
z <- point.size.rescale*z/ (length(x))
symbols( x, y, circles = z,#rep(0.1, length(x)), #sample(1:2, length(x), replace = T) , inches=F, bg= "grey", fg = bg, add = T)
# points(x, y, pch = pch, col = col, bg = bg, cex = cex)
ok <- is.finite(x) & is.finite(y)
if (any(ok))
lines(stats::lowess(x[ok], y[ok], f = span, iter = iter), col = col.smooth, ...)
}

panel.hist <- function(x, ...)
{
usr <- par("usr"); on.exit(par(usr))
par(usr = c(usr[1:2], 0, 1.5) )
h <- hist(x, plot = FALSE)
breaks <- h$breaks; nB <- length(breaks)
y <- h$counts; y <- y/max(y)
rect(breaks[-nB], 0, breaks[-1], y, col=NULL, ...)
}

pairs.ordered.categorical <- function(xx,...)
{
pairs(xx ,
diag.panel = panel.hist ,
lower.panel=panel.smooth.ordered.categorical,
upper.panel=panel.cor.ordered.categorical,
cex.labels = 0.8,...)
}

Ran the following to get the output
traits_needed <-c("SB_Na","SB_K", "blade_K_Na", "sheath_K_Na")
needed <- c("darkgoldenrod4", "paleturquoise4", "coral4", "blueviolet")
pairs.ordered.categorical(cbind(MEs[, needed], traitsDat[, traits_needed]))

This is no means my own work. I adopted this from https://www.r-statistics.com/2010/04/correlation-scatter-plot-matrix-for-ordered-categorical-data/

Will you please be able too help me edit this code so it'll fit what I want to do (i.e. color the histogram only when there's a color present and color the scatter plot as you've done in your example)

Thank you very very much
Shani.

[Kevin Blighe]

Okay, let me take a look in a while!

[Kevin Blighe]

Hi Shani, I have responded above.