Custom circosPlot function that plots only + or - correlations

Hello. We very much like your mixOmics program. One of our coauthors has requested that we only show positive (red lines) or negative (blue lines) correlations in a single figure.

Is this possible to do in the circosplot?

Thank you in advance

Hi @mixOmics_user,

I created a custom circosPlot2 function to do what you asked for. With the sub-functions it’s greater than word limit for Discourse posts so I break it down into multiple posts (I recommend you save them all in a single R file):

1. Sub-functions and dependencies:

run these first:

library(RColorBrewer)
library(dplyr)
library(tidyr)
library(reshape2)

drawIdeogram = function(R, xc=400, yc=400, cir, W,
                        show.band.labels = FALSE,
                        show.chr.labels = FALSE, chr.labels.R = 0,
                        chrData,
                        size.variables,
                        size.labels,
                        color.blocks,
                        line)
{
    # Draw the main circular plot: segments, bands and labels
    chr.po    = cir 
    chr.po[,1]  = gsub("chr","",chr.po[,1]) 
    chr.num     = nrow(chr.po) 
    
    dat.c     = chrData 
    dat.c[,1] = gsub("chr", "", dat.c[,1]) 
    
    for (chr.i in c(1:chr.num)){
        chr.s  = chr.po[chr.i,1] 
        
        v1 = as.numeric(chr.po[chr.i,2]) 
        v2 = as.numeric(chr.po[chr.i,3]) 
        v3 = as.numeric(chr.po[chr.i,6]) 
        v4 = as.numeric(chr.po[chr.i,7]) 
        
        dat.v = subset(dat.c, dat.c[,1]==chr.s) 
        dat.v = dat.v[order(as.numeric(dat.v[,2])),] 
        for (i in 1:nrow(dat.v)){
            
            #col.v = which(colors()==dat.v[i,5])  #get color index
            #col = colors()[col.v]
            dark.clear = color.blocks#brewer.pal(n = 12, name = 'Paired')
            col.v = which(dark.clear==dat.v[i,5])  #get color index
            col = dark.clear[col.v]
            
            w1 = scale.v(as.numeric(dat.v[i,2]), v1, v2, v3, v4) 
            w2 = scale.v(as.numeric(dat.v[i,3]), v1, v2, v3, v4) 
            
            draw.arc.s(xc, yc, R, w1, w2, col=col, lwd=W) 
            
            if (show.band.labels){
                band.text = as.character(dat.v[i,"name.user"]) 
                
                band.po = ((w1+w2)/2)# - ((w2-w1)/3) #position around the circle
                # print(c(band.po, w1, w2, (w2-w1)/3))
                band.po.in = R-(W/3.0) #position on the band (middle)
                draw.text.rt(xc, yc,band.po.in  , band.po , band.text ,
                             cex = size.variables, segmentWidth = W, side="in" )
            }
        } #End for row
        if (show.chr.labels){
            w.m = (v1+v2)/2 
            chr.t = gsub("chr", "", chr.s)
            if(line == TRUE)
            {
                draw.text.rt(xc, yc, chr.labels.R, w.m, chr.t, cex=size.labels,
                             segmentWidth = W, parallel=TRUE)
            } else {
                #put the labels closer to the circle
                draw.text.rt(xc, xc, chr.labels.R, w.m, chr.t, cex=size.labels,
                             segmentWidth = 75, parallel=TRUE)
            }
        }
    } #End for
}

drawLinks = function(R, xc=400, yc=400, cir, W,
                     mapping=mapping,
                     lineWidth=1, col=rainbow(10, alpha=0.8)[7],  drawIntraChr=FALSE,
                     color.cor = color.cor)
{
    # Draw the links (computed correlation) between features
    chr.po    = cir 
    chr.po[,1]  = gsub("chr","",chr.po[,1]) 
    chr.num     = nrow(chr.po) 
    
    chr.po[,4] = gsub("chr", "", chr.po[,4]) 
    dat.in = mapping 
    dat.in[,1] = gsub("chr", "", dat.in[,1]) 
    dat.in[,4] = gsub("chr", "", dat.in[,4]) 
    
    dat    = dat.in 
    
    for (i in 1:nrow(dat)){
        chr1.s   = dat[i,1] 
        chr2.s   = dat[i,4] 
        po1      = dat[i,2] 
        po2      = dat[i,5] 
        
        chr1     = which(chr.po[,1]==chr1.s) 
        chr2     = which(chr.po[,1]==chr2.s) 
        
        v1 = as.numeric(chr.po[chr1,2]) 
        v2 = as.numeric(chr.po[chr1,3]) 
        v3 = as.numeric(chr.po[chr1,6]) 
        v4 = as.numeric(chr.po[chr1,7]) 
        
        w1 = scale.v(as.numeric(po1), v1, v2, v3, v4) 
        
        v1 = as.numeric(chr.po[chr2,2]) 
        v2 = as.numeric(chr.po[chr2,3]) 
        v3 = as.numeric(chr.po[chr2,6]) 
        v4 = as.numeric(chr.po[chr2,7]) 
        
        w2 = scale.v(as.numeric(po2), v1, v2, v3, v4) 
        # Set the link width depending on the correlation coefficient
        lwd = abs(as.numeric(dat[i,7])) 
        
        # Set link color
        if (as.numeric(dat[i,7]) < 0.0){
            linkCol = color.cor[2]#colors()[128]  #pale red
        } else {
            linkCol = color.cor[1]#colors()[134]  # blue
        }
        linkCol = add.alpha(linkCol, alpha=0.4) 
        
        if (chr1 == chr2){
            if (drawIntraChr == TRUE){
                draw.link(xc, yc, R, w1, w2, col=linkCol, lwd=lineWidth) 
            }
        } else {
            draw.link(xc, yc, R, w1, w2, col=linkCol, lwd=lineWidth) 
        }
    } ### End for
}

drawLinePlot = function(mapping=mapping, xc=400, yc=400, col.v=3,
                        R, cir,   W, col='black', scale=FALSE, lineWidth=1,
                        background.lines=FALSE,axis.width=1)
{
    # Generate a linear plot around the main ideogram.
    #
    # fixme: the function writes the same line multiple times
    # it needs to be called only once and process the data/segment
    # separately
    chr.po    = cir 
    chr.po[,1]  = gsub("chr","",chr.po[,1]) 
    chr.num     = nrow(chr.po) 
    
    dat.in   = mapping 
    dat.in[,1] = gsub("chr", "", dat.in[,1]) 
    
    # data set for the chromosome
    for (chr.i in 1:chr.num){
        chr.s = chr.po[chr.i,1] 
        chr.s = gsub("chr","",chr.s) 
        dat   = subset(dat.in, dat.in[,1]==chr.s) 
        dat   = dat[order(as.numeric(dat[,2])),] 
        v1 = as.numeric(chr.po[chr.i,2]) 
        v2 = as.numeric(chr.po[chr.i,3]) 
        v3 = as.numeric(chr.po[chr.i,6]) 
        v4 = as.numeric(chr.po[chr.i,7]) 
        
        # background line
        if (background.lines){
            draw.arc.pg(xc, yc, v1, v2, R, R+W-5, col=colors()[245]) 
        } else {
            draw.arc.s(xc, yc, R, v1, v2, col=colors()[245], lwd=axis.width) 
        }
    }
    
    my.R1 = R + W/5 
    my.R2 = R + W - W/5 
    
    ## for the matrix colors
    num.col = ncol(dat[,col.v:ncol(dat)]) 
    num.row = nrow(dat.in) 
    
    if (length(col) == num.col){
        colors = col 
    } else {
        colors = rainbow(num.col, alpha=0.5) 
    }
    
    for (chr.i in 1:chr.num){
        chr.s = chr.po[chr.i,1] 
        chr.s = gsub("chr","",chr.s) 
        
        dat   = subset(dat.in, dat.in[,1]==chr.s) 
        dat   = dat[order(as.numeric(dat[,2])),] 
        #print(head(dat))
        dat.i   = c(col.v:ncol(dat)) 
        dat.m   = dat.in[,dat.i] 
        dat.m   = as.matrix(dat.m) 
        dat.min = min(as.numeric(dat.m), na.rm=TRUE)
        dat.max = max(as.numeric(dat.m), na.rm=TRUE) 
        
        v1 = as.numeric(chr.po[chr.i,2]) 
        v2 = as.numeric(chr.po[chr.i,3]) 
        v3 = as.numeric(chr.po[chr.i,6]) 
        v4 = as.numeric(chr.po[chr.i,7]) 
        
        col.i = 0 
        
        for (j in col.v:ncol(dat)){
            col.i = col.i + 1 
            col   = colors[col.i] 
            
            my.v      = as.numeric(dat[1,j]) 
            dat.i.old = my.v 
            v.old   = scale.v(my.v, my.R1, my.R2, dat.min, dat.max) 
            
            po      = as.numeric(dat[1,2]) 
            w.from  = scale.v(po, v1, v2, v3, v4) 
            
            
            for (k in 1:nrow(dat)){
                
                dat.i = as.numeric(dat[k,j]) 
                
                if (is.na(dat.i)){
                    next 
                }
                
                v    = scale.v(dat.i, my.R1, my.R2, dat.min, dat.max) 
                w.to = scale.v(as.numeric(dat[k,2]), v1, v2, v3, v4) 
                
                if (w.from > 0){
                    draw.line3(xc, yc, w.from, w.to, v.old, v, col=col,
                               lwd=lineWidth)
                }
                
                dat.i.old = dat.i 
                w.from    = w.to 
                v.old     = v 
            } # end the row
        }   # end the col
        if (scale){
            do.scale(xc, yc, dat.min, dat.max, R+W/5, W-2*(W/5)) 
        }
    }     # end the chr/segment
}

genChr =function (expr, bandWidth = 1.0, color.blocks)
{
    # Generate the segments and calculate the
    # unique positions of the bands
    #
    # Args:
    #   expr : dataframe containing the features expression
    # example: colnames(concatFeatExp) "PAM50"    "Features" "Mean"     "SD"
    # "Dataset"
    #   bandWidth: thickness of each band
    #
    # Return:
    #   a data.frame that can be used with segAnglePo
    
    # expr can contains expression data for multiple diseases
    # here, we only use the Chrom and Dataset column and remove duplicates
    keeps = c("Features","Dataset") 
    expr = expr[keeps] 
    expr = unique(expr) 
    chrLengths = data.frame(table(expr$Dataset))  
    rownames(chrLengths) = chrLengths[,1] 
    chrLengths[,1] = NULL 
    colnames(chrLengths) = c( "Freq") 
    chrLengths[, "Count"] = rep(0, nrow(chrLengths)) 
    
    # Last column contains the bands' color
    # Create a color scheme
    #dark = c("brown3","darkgoldenrod","antiquewhite3","steelblue3")
    #clear = c("brown1","darkgoldenrod1","antiquewhite1","steelblue1")
    dark.clear = color.blocks#brewer.pal(n = 12, name = 'Paired')
    dark = dark.clear[seq(2, 12, by = 2)]
    clear = dark.clear[seq(1, 12, by = 2)]
    chrColScheme = data.frame(dark, clear)
    n_datasets = length(unique(expr$Dataset))
    chrColScheme = chrColScheme[c(1:n_datasets),]
    rownames(chrColScheme) = levels(factor(expr$Dataset))#alphabetical order
    
    seg.out = c() 
    for (i in 1:nrow(expr)){
        chrName    = paste("chr", as.character(expr[i,'Dataset'][[1]]), sep="") 
        dType = as.character(expr[i,'Dataset'][[1]]) 
        pStart = chrLengths[dType,'Count'] * bandWidth 
        pStop = chrLengths[dType,'Count'] * bandWidth + bandWidth  
        chrLengths[dType,'Count'] = chrLengths[dType,'Count'] + 1 
        fName = as.character(as.matrix(expr[i,'Features']))
        # added as.character() by amrit
        # Assign colors
        if (chrLengths[dType,'Count'] %% 2 == 0){
            chrCol = chrColScheme[dType,]$clear 
        } else{
            chrCol = chrColScheme[dType,]$dark 
        }
        seg.out = rbind(seg.out, c(chrName, pStart, pStop,  fName, chrCol)) 
    }
    
    # Use the same names than in omicCircos
    colnames(seg.out) = c("chrom", "chromStart", "chromEnd", "name", "color") 
    seg.out = as.data.frame(seg.out) 
    
    return(seg.out) 
}

genLinks2 = function(chr, corMat, threshold)
{
    
    # to satisfy R CMD check that doesn't recognise x, y and group (in aes)
    Var1=Var2=chrom=NULL
    
    
    # Generates the links corresponding to pairwise correlations
    #
    # Args:
    #   chr: ideogram structure (generated from genChr)
    #   corMat: main correlation matrix
    #   threshold: minimum correlation value
    #
    # Return:
    #   the links data (see omicsCircos doc)
    #
    linkList = c() 
    # Remove matrix diagonal and the the mat in a list
    linkList = subset(melt(corMat), Var1!=Var2 ) 
    # Remove links below the threshold
    if(threshold > 0) {
        linkList = subset(linkList, linkList$value >= threshold) 
    } else {
        linkList = subset(linkList, linkList$value <= threshold) 
    }
    
    
    #First merge
    linkList = dplyr::rename(linkList, feat1=Var1, feat2=Var2)
    # CHANGED BY AMRIT
    linkList = merge(linkList, chr, by.x="feat1", by.y="name") 
    # Set the position in the middle of the band
    linkList$po1 = (as.numeric(linkList$chromStart)
                    + as.numeric(linkList$chromEnd)) / 2.0
    linkList = dplyr::rename(linkList, chr1=chrom)   # CHANGED BY AMRIT
    keeps = c("feat1","feat2","value","chr1","po1") 
    linkList = linkList[keeps] 
    
    #Second merge
    linkList = merge(linkList, chr, by.x="feat2", by.y="name") 
    linkList$po2 = (as.numeric(linkList$chromStart)
                    + as.numeric(linkList$chromEnd)) / 2.0
    linkList = dplyr::rename(linkList, chr2=chrom)   # CHANGED BY AMRIT
    keeps = c("chr1","po1","feat1","chr2","po2","feat2","value") 
    linkList = linkList[keeps] 
    
    return(linkList) 
}

bezierCurve = function(x, y, n=10)  {  
    outx = NULL  
    outy = NULL   
    i = 1	
    for (t in seq(0, 1, length.out=n))		{		
        b = bez(x, y, t)		
        outx[i] = b$x		
        outy[i] = b$y 		
        i = i+1		
    } 	
    return (list(x=outx, y=outy))	
} 

##
bez = function(x, y, t)	{	
    outx = 0	
    outy = 0	
    n = length(x)-1	
    for (i in 0:n)		{		
        outx = outx + choose(n, i)*((1-t)^(n-i))*t^i*x[i+1]		
        outy = outy + choose(n, i)*((1-t)^(n-i))*t^i*y[i+1]		
    } 	
    return (list(x=outx, y=outy))	
}

###########################################
# one value : from a to b
scale.v = function(v, a, b, min.v, max.v) {
    v = v-min.v 
    v = v/(max.v-min.v) 
    v = v*(b-a) 
    v+a
}

### draw.link
draw.link = function(xc, yc, r, w1, w2, col=col, lwd=lwd) {
    # for translocation
    w3  = (w1+w2)/2 
    w1  = w1/360*2*pi 
    w2  = w2/360*2*pi 
    w3  = w3/360*2*pi 
    x0  = xc+r*cos(w1) 
    y0  = yc-r*sin(w1) 
    x1  = xc+r*cos(w2) 
    y1  = yc-r*sin(w2) 
    x = c(x0,xc,xc,x1) 
    y = c(y0,yc,yc,y1) 
    points(bezierCurve(x,y,60), type="l", col=col, lwd=lwd, lend="butt")
}

### draw.link2
draw.link2 = function(xc, yc, r, w1, w2, col=col, lwd=lwd) {
    # for translocation
    w3  = (w1+w2)/2 
    w1  = w1/360*2*pi 
    w2  = w2/360*2*pi 
    w3  = w3/360*2*pi 
    x2  = xc+r/2*cos(w3) 
    y2  = yc-r/2*sin(w3) 
    x0  = xc+r*cos(w1) 
    y0  = yc-r*sin(w1) 
    x1  = xc+r*cos(w2) 
    y1  = yc-r*sin(w2) 
    x = c(x0, x2, x2, x1) 
    y = c(y0, y2, y2, y1) 
    points(bezierCurve(x,y,60), type="l", col=col, lwd=lwd, lend="butt")
}
###

###
draw.link.pg = function(xc, yc, r, w1.1, w1.2, w2.1, w2.2, col=col, lwd=lwd) {
    w1 = w1.1 
    w2 = w2.2 
    w3  = (w1+w2)/2 
    w1  = w1/360*2*pi 
    w2  = w2/360*2*pi 
    w3  = w3/360*2*pi 
    x0  = xc+r*cos(w1) 
    y0  = yc-r*sin(w1) 
    x1  = xc+r*cos(w2) 
    y1  = yc-r*sin(w2) 
    x = c(x0,xc,xc,x1) 
    y = c(y0,yc,yc,y1) 
    bc1 = bezierCurve(x,y,60) 
    
    ang.d = abs(w1.1-w1.2) 
    pix.n = ang.d * 10 
    if (pix.n < 10){
        pix.n = 10 
    }
    
    ang.seq = rev(seq(w1.1,w1.2,length.out=pix.n)) 
    ang.seq = ang.seq/360*2*pi 
    
    fan.1.x = xc + cos(ang.seq) * r 
    fan.1.y = yc - sin(ang.seq) * r 
    
    ######################################################
    w1 = w1.2 
    w2 = w2.1 
    w3  = (w1+w2)/2 
    w1  = w1/360*2*pi 
    w2  = w2/360*2*pi 
    w3  = w3/360*2*pi 
    x0  = xc+r*cos(w1) 
    y0  = yc-r*sin(w1) 
    x1  = xc+r*cos(w2) 
    y1  = yc-r*sin(w2) 
    x = c(x0,xc,xc,x1) 
    y = c(y0,yc,yc,y1) 
    bc2 = bezierCurve(x,y,60) 
    
    ang.d = abs(w2.1-w2.2) 
    pix.n = ang.d * 10 
    if (pix.n < 10){
        pix.n = 10 
    }
    
    ang.seq = rev(seq(w2.1,w2.2,length.out=pix.n)) 
    ang.seq = ang.seq/360*2*pi 
    
    fan.2.x = xc + cos(ang.seq) * r 
    fan.2.y = yc - sin(ang.seq) * r 
    
    polygon(c(bc1$x, fan.2.x, rev(bc2$x), rev(fan.1.x)),
            c(bc1$y, fan.2.y, rev(bc2$y), rev(fan.1.y)),
            fillOddEven=TRUE, border=col, col=col, lwd=lwd) 
}

###
draw.point.w = function(xc, yc, r, w, col=col, cex=cex){
    w = w/360*2*pi 
    x = xc+r*cos(w) 
    y = yc-r*sin(w) 
    points(x, y, pch=20, col=col, cex=cex) 
}

###
draw.text.w = function(xc, yc, r, w, n, col="black", cex=1){
    w = w%%360 
    w = w/360*2*pi 
    x = xc+r*cos(w) 
    y = yc-r*sin(w) 
    text(x,y,labels=n, col=col, cex=cex) 
}

###
draw.text.rt = function(xc, yc, r, w, n, col="black", cex=1, side="out",
                        segmentWidth=20, parallel=FALSE){
    w     = w%%360 
    the.o = w 
    
    the.w = 360-w 
    w     = w/360*2*pi 
    x     = xc+r*cos(w) 
    y     = yc-r*sin(w) 
    
    
    num2  = (segmentWidth*2)/2.0 
    b = the.w
    if (side=="out"){
        if (the.w <= 90 ){
            the.pos = 4 
            if (parallel == TRUE){
                the.w = the.w -90 # 180 
            }
        } else if (the.w > 90 & the.w <= 180) {
            if (parallel == TRUE){
                the.w = the.w -90 # 180 
            }
            else {
                the.w = the.w + 180 
            }
            the.pos = 2 
        } else if (the.w > 180 & the.w <= 270){
            the.w = the.w%%180 
            if (parallel == TRUE){
                the.w = the.w -90 # 180 
            }
            the.pos = 2 
        } else if (the.w > 270 & the.w <= 360){
            the.pos = 4 
            if (parallel == TRUE){
                the.w = the.w + 90 
            }
        }
        
        if (the.pos==2){
            x = x+num2 
        }
        if (the.pos==4){
            x = x-num2 
        }
    }
    
    if (side=="in"){
        if (the.w <= 90 ){
            the.pos = 4 
        } else if (the.w > 90 & the.w <= 180) {
            the.w = the.w + 180 
            the.pos = 2 
        } else if (the.w > 180 & the.w <= 270){
            the.w = the.w%%180 
            the.pos = 2 
        } else if (the.w > 270 & the.w <= 360){
            the.pos = 4 
        }
        
        if (the.pos==2){
            x = x+segmentWidth 
        }
        if (the.pos==4){
            x = x-segmentWidth 
        }
    }
    
    
    text(x, y, adj=0, offset=1, labels=n, srt=the.w,
         pos=the.pos, col=col, cex=cex) 
}

###strokeLine2
draw.line = function (xc, yc, w, l1, l2, col=col, lwd=lwd, lend=1) {
    w  = (w/360)*2*pi 
    x1 = xc+l1*cos(w) 
    y1 = yc-l1*sin(w) 
    x2 = xc+l2*cos(w) 
    y2 = yc-l2*sin(w) 
    segments(x1, y1, x2, y2, col=col, lwd=lwd, lend=lend) 
}

###strokeLine3
draw.line2 = function (xc, yc, w, r, l, col=col, lwd=lwd){
    line_w   = l 
    theangle = w 
    l1       = r 
    theangle = (theangle/360)*2*pi 
    x0       = xc+l1*cos(theangle) 
    y0       = yc+l1*sin(theangle) 
    w1       = 45/360*2*pi 
    x1 = xc + sin(w1) * (x0) 
    y1 = yc + cos(w1) * (y0) 
    x2 = xc - sin(w1) * (x0) 
    y2 = yc - cos(w1) * (y0) 
    segments(x1, y1, x2, y2, col=col, lwd=lwd, lend="butt") 
}

###strokeLine by two angles
draw.line3 = function (xc, yc, w1, w2, r1, r2, col=col, lwd=lwd){
    theangle1 = w1 
    theangle2 = w2 
    l1        = r1 
    l2        = r2 
    
    theangle1 = (theangle1/360)*2*pi 
    x1        = xc+l1*cos(theangle1) 
    y1        = yc-l1*sin(theangle1) 
    
    theangle2 = (theangle2/360)*2*pi 
    x2        = xc+l2*cos(theangle2) 
    y2        = yc-l2*sin(theangle2) 
    
    segments(x1, y1, x2, y2, col=col, lwd=lwd, lend="butt") 
}

### plot fan or sector that likes a piece of doughnut (plotFan)
draw.arc.pg = function (xc, yc,
                        w1, w2, r1, r2, col="lightblue", border="lightblue", lwd=0.01
){
    
    ang.d = abs(w1-w2) 
    pix.n = ang.d * 10 
    if (pix.n < 10){
        pix.n = 10 
    }
    
    ang.seq = rev(seq(w1,w2,length.out=pix.n)) 
    ang.seq = ang.seq/360*2*pi 
    
    fan.i.x = xc + cos(ang.seq) * r1 
    fan.i.y = yc - sin(ang.seq) * r1 
    
    
    fan.o.x = xc + cos(ang.seq) * r2 
    fan.o.y = yc - sin(ang.seq) * r2 
    
    polygon(c(rev(fan.i.x), fan.o.x ), c(rev(fan.i.y), fan.o.y),
            fillOddEven=TRUE, border=border, col=col, lwd=lwd, lend=1)
    
}

draw.arc.s = function (xc, yc, r, w1, w2, col="lightblue", lwd=1, lend=1){
    # Draw circular arcs for the main ideogram)
    # s = simple
    # r = radius
    ang.d = abs(w1-w2) 
    pix.n = ang.d * 5 
    if (pix.n < 2){
        pix.n = 2 
    }
    
    ang.seq = rev(seq(w1,w2,length.out=pix.n)) 
    ang.seq = ang.seq/360*2*pi 
    
    fan.i.x = xc + cos(ang.seq) * r 
    fan.i.y = yc - sin(ang.seq) * r 
    ## lend=0(round)  #lend=1(butt)  lend=2(square)
    lines(fan.i.x, fan.i.y, col=col, lwd=lwd, type="l", lend=lend) 
    #points(fan.i.x, fan.i.y, col=col, lwd=lwd, type="l", lend=lend) 
}

#########################################
## segment to angle and position
## segAnglePo
#########################################

# get angle if given seg and position
# seg should be ordered by user
segAnglePo = function (seg.dat=seg.dat, seg=seg, angle.start=angle.start,
                       angle.end=angle.end){
    
    if (missing(angle.start)){
        angle.start = 0 
    }
    if (missing(angle.end)){
        angle.end = 360 
    }
    ## check data.frame?
    colnames(seg.dat) = c("seg.name","seg.Start","seg.End","name","gieStain") 
    
    ## get length of the segomosomes
    seg.l   = c() 
    seg.min = c() 
    seg.sum = c() 
    seg.s   = 0 
    seg.num   = length(seg) 
    seg.names = seg 
    
    ########################################################
    ########################################################
    for (i in 1:seg.num){
        seg.n = seg.names[[i]] 
        
        dat.m = subset(seg.dat, seg.dat[,1]==seg.n) 
        seg.full.l = max(as.numeric(dat.m[,"seg.End"])) 
        seg.full.m = min(as.numeric(dat.m[,"seg.Start"])) 
        seg.l      = cbind(seg.l, seg.full.l) 
        seg.min    = cbind(seg.min, seg.full.m) 
        seg.s      = seg.s + seg.full.l 
        seg.sum    = cbind(seg.sum, seg.s) 
    }
    
    ## initial parameters
    gap.angle.size = 2 
    seg.angle.from = angle.start + 270 
    
    seg.full.l  = sum(as.numeric(seg.l)) 
    angle.range = angle.end - angle.start 
    cir.angle.r = (angle.range - seg.num * gap.angle.size)/seg.full.l 
    
    out.s     = c() 
    l.old     = 0 
    gap.angle = 0 
    for (i in 1:seg.num){
        seg.n = seg.names[[i]] 
        dat.m = subset(seg.dat, seg.dat[,1]==seg.n) 
        len   = seg.sum[i] 
        w1    = cir.angle.r*l.old + gap.angle 
        w2    = cir.angle.r*len   + gap.angle 
        out.s     = rbind(out.s, c(seg.n, w1+seg.angle.from, w2+seg.angle.from,
                                   l.old, len, seg.min[i], seg.l[i]))
        gap.angle = gap.angle + gap.angle.size 
        l.old     = len 
    }
    
    colnames(out.s) = c("seg.name","angle.start", "angle.end", "seg.sum.start",
                        "seg.sum.end","seg.start", "seg.end")
    return(out.s) 
}


### start do.scale
do.scale = function(xc=xc, yc=yc, dat.min=dat.min, dat.max=dat.max,
                    R=R, W=W, s.n=1, col="blue"){
    dat.m   = round((dat.min+dat.max)/2, s.n) 
    dat.min = round(dat.min, s.n) 
    dat.max = round(dat.max, s.n) 
    y1      = yc + R  
    y2      = yc + R + W/2 
    y3      = yc + R + W 
    x1      = xc - W/20 
    x2      = x1 - (W/20)*1.2 
    x3      = x1 - (W/20)*3 
    segments(x1, y1, x1, y3, lwd=0.01, col=col) 
    segments(x1, y1, x2, y1, lwd=0.01, col=col) 
    segments(x1, y2, x2, y2, lwd=0.01, col=col) 
    segments(x1, y3, x2, y3, lwd=0.01, col=col) 
    text(x3, y1, dat.min, cex=0.2, col=col) 
    text(x3, y2, dat.m,   cex=0.2, col=col) 
    text(x3, y3, dat.max, cex=0.2, col=col)
}
### end do.scale

## Add an alpha value to a colour
add.alpha = function(col, alpha=1){
    if(missing(col))
        stop("Please provide a vector of colours.")
    apply(sapply(col, col2rgb)/255, 2, 
          function(x) 
              rgb(x[1], x[2], x[3], alpha=alpha))  
}
1 Like

2. Main function

then run this in the same environment

circosPlot2 = function(object,
                       comp = 1 : min(object$ncomp),
                       cutoff,
                       color.Y,
                       color.blocks,
                       color.cor,
                       var.names = NULL,
                       showIntraLinks = FALSE,
                       line=TRUE,
                       size.legend=0.8,
                       ncol.legend=1,
                       size.variables = 0.25,
                       size.labels=1,
                       legend = TRUE)
{
    # to satisfy R CMD check that doesn't recognise x, y and group (in aes)
    Features = Exp = Dataset = Mean = linkColors = chrom = po = NULL
    
    
    options(stringsAsFactors = FALSE)
    figSize = 800
    segmentWidth = 25
    linePlotWidth = 90
    
    ##############################
    ###   networkDiagram_core.R
    ###
    ###   Authors: Michael Vacher (minor changes by Amrit :)
    ###
    ###   Parts of this code has been modified from the original OmicCircos
    ###     package obtained from:
    ###   Ying Hu Chunhua Yan <yanch@mail.nih.gov> (2015). OmicCircos:
    ###     High-quality circular visualization of omics data. v1.6.0.
    ##############################
    
    # check input object
    if (!is(object, "block.splsda"))
        stop("circosPlot is only available for 'block.splsda' objects")
    
    if (length(object$X) < 2)
        stop("This function is only available when there are more than 3 blocks
    (2 in object$X + an outcome object$Y)") # so 2 blocks in X + the outcome Y
    
    if (missing(cutoff))
        stop("'cutoff' is missing", call.=FALSE) # so 2 blocks in X + the outcome Y
    
    if(missing(color.Y))
    {
        color.Y = color.mixo(1:nlevels(object$Y))
    } else {
        if(length(color.Y) != nlevels(object$Y))
            stop("'color.Y' must be of length ", nlevels(object$Y))
        
    }
    
    if(missing(color.blocks))
    {
        color.blocks = brewer.pal(n = 12, name = 'Paired') #why 12??
    } else {
        if(length(color.blocks) != length(object$X))
            stop("'color.blocks' must be of length ", length(object$X))
        
        color.blocks.adj = adjustcolor(color.blocks, alpha.f = 0.5)
        #to get two shades of the same color per block
        
        color.blocks = c(rbind(color.blocks, color.blocks.adj))
        # to put the color next to its shaded color
    }
    
    if(missing(color.cor))
    {
        color.cor = c(colors()[134],  # blue, negative correlation
                      colors()[128])  # pale red, positive correlation
    } else {
        if(length(color.cor) != 2)
            stop("'color.cor' must be of length 2")
    }
    
    X = object$X
    Y = object$Y
    
    #need to reorder variates and loadings to put 'Y' in last
    indY = object$indY
    object$variates = c(object$variates[-indY], object$variates[indY])
    object$loadings = c(object$loadings[-indY], object$loadings[indY])
    object$ncomp = c(object$ncomp[-indY], object$ncomp[indY])
    
    #check var.names
    sample.X = lapply(object$loadings[-length(object$loadings)],
                      function(x){1 : nrow(x)})
    if (is.null(var.names))
    {
        var.names.list = unlist(sapply(object$loadings[
            -length(object$loadings)], rownames))
    } else if (is.list(var.names)) {
        if (length(var.names) != length(object$loadings[
            -length(object$loadings)]))
            stop.message('var.names', sample.X)
        
        if(sum(sapply(1 : length(var.names), function(x){
            length(var.names[[x]]) == length(sample.X[[x]])})) !=
           length(var.names))
            stop.message('var.names', sample.X)
        
        var.names.list = var.names
    } else {
        stop.message('var.names', sample.X)
    }
    
    
    if(any(comp > min(object$ncomp)))
    {
        warning("Limitation to ",min(object$ncomp),
                " components, as determined by min(object$ncomp)")
        comp[which(comp > min(object$ncomp))] = min(object$ncomp)
    }
    comp = unique(sort(comp))
    
    
    keepA = lapply(object$loadings, function(i)
        apply(abs(i)[, comp, drop = FALSE], 1, sum) > 0)
    cord = mapply(function(x, y, keep){
        cor(x[, keep], y[, comp], use = "pairwise")
    }, x=object$X, y=object$variates[-length(object$variates)],
    keep = keepA[-length(keepA)],SIMPLIFY = FALSE)
    
    simMatList = vector("list", length(X))
    for(i in 1:length(cord))
    {
        for(j in 1:length(cord))
        {
            simMatList[[i]][[j]] = cord[[i]] %*% t(cord[[j]])
        }
    }
    corMat = do.call(rbind, lapply(simMatList, function(i) do.call(cbind, i)))
    
    ## Expression levels
    Xdat = as.data.frame(do.call(cbind, X)[, colnames(corMat)])
    
    AvgFeatExp0 = Xdat %>% mutate(Y = Y) %>% gather(Features, Exp, -Y) %>%
        group_by(Y, Features) %>% dplyr::summarise(Mean = mean(Exp), SD = sd(Exp))
    AvgFeatExp0$Dataset = factor(rep(names(X), unlist(lapply(cord, nrow))),
                                 levels = names(X))[match(AvgFeatExp0$Features,colnames(Xdat))]
    # to match Xdat that is reordered in AvgFeatExp0
    featExp = AvgFeatExp0 %>% group_by(Dataset, Y) %>% arrange(Mean)
    #Generate a circular plot (circos like) from a correlation matrix (pairwise)
    #
    # Args:
    #   corMat: the main correlation matrix.
    #         -> colnames == rownames (pairwise)  values = correlations
    #   featExp: data.frame holding the expression data.
    #   cutoff: minimum value for correlations (<threshold will be ignored)
    #   figSize: figure size
    #   segmentWidth: thickness of the segment (main circle)
    #   linePlotWidth: thickness of the line plot (showing expression data)
    #   showIntraLinks = display links intra segments
    
    
    
    # 1) Generate karyotype data
    chr = genChr(featExp, color.blocks = color.blocks)
    chr.names = unique(chr$chrom) # paste("chr", 1:seg.num, sep="") 
    # Calculate angles and band positions
    db = segAnglePo(chr, seg=chr.names)
    db = data.frame(db)
    
    # 2) Generate Links
    links = genLinks2(chr, corMat, threshold=cutoff)
    if (nrow(links) < 1)
        warning("Choose a lower correlation threshold to highlight
    links between datasets")
    
    # 3) Plot
    # Calculate parameters
    circleR = (figSize / 2.0) -  segmentWidth - linePlotWidth
    linksR = circleR - segmentWidth
    linePlotR = circleR + segmentWidth
    chrLabelsR = (figSize / 2.0)
    
    # replace chr$name by the ones in var.names (matching)
    # matching var.names.list with object$loadings
    ind.match = match(chr$name, unlist(sapply(object$loadings[
        -length(object$loadings)],rownames)))
    chr$name.user = unlist(var.names.list)[ind.match]
    
    opar1=par("mar")
    par(mar=c(2, 2, 2, 2))
    
    plot(c(1,figSize), c(1,figSize), type="n", axes=FALSE, xlab="",
         ylab="", main="")
    
    #save(list=ls(),file="temp.Rdata")
    # Plot ideogram
    drawIdeogram(R=circleR, cir=db, W=segmentWidth,  show.band.labels=TRUE,
                 show.chr.labels=TRUE, chr.labels.R= chrLabelsR, chrData=chr,
                 size.variables = size.variables, size.labels=size.labels,
                 color.blocks = color.blocks, line = line)
    # Plot links
    if(nrow(links)>0)
        drawLinks(R=linksR, cir=db,   mapping=links,   col=linkColors,
                  drawIntraChr=showIntraLinks, color.cor = color.cor)
    
    # Plot expression values
    cTypes = levels(Y)
    #unique(featExp[,1]) #Get the different disease/cancer types (lines)
    #lineCols = rainbow(nrow(cTypes), alpha=0.5)
    lineCols = color.Y
    #color.mixo(1:nlevels(Y))#color.mixo(match(levels(Y), levels(Y)))
    
    # Fixme: remove this loop and send the whole expr dframe to drawLinePlot
    if(line==TRUE)
    {
        for (i in 1:length(chr.names)){
            seg.name = gsub("chr","",chr.names[i])
            #Get data for each segment
            expr = subset(featExp,featExp$Dataset==seg.name)
            
            expr = dcast(expr, formula = Features ~ Y, value.var="Mean")
            ## changed PAM50 to Y
            expr = merge(expr, chr, by.x="Features", by.y="name")
            expr$po = (as.numeric(expr$chromStart) +
                           as.numeric(expr$chromEnd)) / 2.0
            expr = dplyr::rename(expr, seg.name = chrom, seg.po = po)
            
            # Reorder columns
            cOrder = c(c(grep("seg.name", colnames(expr)),
                         grep("seg.po", colnames(expr))),c(1:length(cTypes)+1))
            expr = expr[, cOrder]
            
            # Plot data on each sub segment
            subChr = subset(db, db$seg.name == chr.names[i] )
            drawLinePlot(R=linePlotR, cir=subChr,   W=linePlotWidth,
                         lineWidth=1, mapping=expr, col=lineCols, scale=FALSE)
        }
    }
    opar=par("xpd")
    par(xpd=TRUE) # to authorise the legend to be written outside the margin,
    #       otherwise it's too small
    # Plot legend
    if(legend == TRUE)
    {
        # First legeng bottom left corner
        legend(x=5, y = (circleR/4), title="Correlations",
               c("Positive Correlation", "Negative Correlation"),
               col = color.cor, pch = 19, cex=size.legend, bty = "n")
        # Second legend bottom righ corner
        if(line==TRUE)
            legend(x=figSize-(circleR/3), y = (circleR/3), title="Expression",
                   legend=levels(Y),  ## changed PAM50 to Y
                   col = lineCols, pch = 19, cex=size.legend, bty = "n",ncol=ncol.legend)
        # third legend top left corner
        legend(x=figSize-(circleR/2), y = figSize, title="Correlation cut-off",
               legend=paste("r", cutoff, sep = "="),
               col = "black", cex=size.legend, bty = "n")
        
        legend(x=-circleR/4, y = figSize, legend=paste("Comp",
                                                       paste(comp,collapse="-")),
               col = "black", cex=size.legend, bty = "n")
    }
    par(xpd=opar,mar=opar1)# put the previous defaut parameter for xpd
    return(invisible(corMat))
}

3. Use examples

library(mixOmics)
data(nutrimouse)
Y = nutrimouse$diet
data = list(gene = nutrimouse$gene, lipid = nutrimouse$lipid)
design = matrix(c(0,1,1,1,0,1,1,1,0), ncol = 3, nrow = 3, byrow = TRUE)


nutrimouse.sgccda <- wrapper.sgccda(X=data,
                                    Y = Y,
                                    design = design,
                                    keepX = list(gene=c(10,10), lipid=c(15,15)),
                                    ncomp = 2,
                                    scheme = "horst")

## cor < -0.7
circosPlot2(nutrimouse.sgccda, cutoff = -0.7, ncol.legend = 2, size.legend = 1.1)
## cor > 0.7 - excluding legend
circosPlot2(nutrimouse.sgccda, cutoff = 0.7, ncol.legend = 2, size.legend = 1.1, legend = FALSE)

## side by side - you can adjust sizes in markdown chunks or RStudio's 'Zoom'
par(mfrow = c(1,2))

circosPlot2(nutrimouse.sgccda, cutoff = -0.7, ncol.legend = 2, size.legend = 1.1, legend = FALSE)

circosPlot2(nutrimouse.sgccda, cutoff = 0.7, ncol.legend = 2, size.legend = 1.1, legend = FALSE)

Hope it helps

Changing cutoff to -0.7 or 0.7 doesn’t work - both positive and negative correlations are still shown on mixOmics circosplot.

Are you sure you’ve copied all the source code from these posts and initialised them in your environment? Also make sure you’re calling circosPlot2() and not circosPlot()