R/permutationsOut.r

In adricaba/hrunbiased: Survival analysis for gene signatures

## Function permutationsOut and plotPermutationsOut
permutationsOut <- function(nda, sigSize, perms, nite = 1000, type,
   mc.cores = 1, GS = NULL, executation.info = TRUE,
   GSassociation = TRUE, adj.var.GS = "GSadj", signature.method="zscore"){


  ln <- length(sigSize)
  new.coefs.nda.ori <- new.coefs.nda.perm <- list()
  new.coefs.nda.perm.GS <- new.coefs.nda.ori.GS <- NULL
  corrGS <- list()
  if(type == "norm01")
  {
          exprs(nda) <- t(apply(exprs(nda),1, sample))
          type <- "none"
  }
  if(is.null(GS)) GS <- apply(exprs(nda),2,mean)

  perms.a <- as.matrix(seq_len(dim(nda)[2]))
  for(L in seq_len(ln)){
     sigsize <- sigSize[L]
     randsign <- vapply(seq_len(nite), function(o)
                       sample(rownames(exprs(nda)), sigsize), character(sigsize))
     new.coefs.nda.perm[[L]] <- coxme.three.options(nda, perms = perms,
            randsign = randsign, type = type,  mc.cores =  mc.cores, GS = GS,
            executation.info = executation.info, adj.var.GS = adj.var.GS,
            signature.method = signature.method)
     new.coefs.nda.ori[[L]] <- coxme.three.options(nda, perms = perms.a,
            randsign = randsign, type = type,  mc.cores =  mc.cores, GS = GS,
            executation.info = executation.info, adj.var.GS = adj.var.GS,
            signature.method = signature.method)
     if(sigsize==1){
         randsign <- t(as.matrix(randsign))
         corrGS[[L]] <- apply(randsign,2,function(o) cor(exprs(nda)[o,],GS))
     }
     else
         corrGS[[L]] <- apply(randsign,2,function(o) cor(apply(exprs(nda)[o,],2,mean),GS))
     }
     if(GSassociation & type!="Gsignature"){
         new.coefs.nda.perm.GS <- coxme.three.options(nda, perms = perms,
            randsign = as.matrix(rownames(nda)), type = type, adj.var.GS = adj.var.GS,
            mc.cores =  mc.cores, GS = GS, executation.info = executation.info)
         new.coefs.nda.ori.GS <- coxme.three.options(nda, perms = perms.a,
            randsign = as.matrix(rownames(nda)), type = type, adj.var.GS = adj.var.GS,
            mc.cores =  mc.cores, GS = GS, executation.info = executation.info)
     }

  return(list(new.coefs.nda.perm = new.coefs.nda.perm,
              new.coefs.nda.ori = new.coefs.nda.ori,
              new.coefs.nda.perm.GS = new.coefs.nda.perm.GS,
              new.coefs.nda.ori.GS = new.coefs.nda.ori.GS, corrGS = corrGS))
}

plotPermutationsOut <- function(perm.out,
            plot.type = c("perm.violin", "perm.GSvsEvents", "perms.corr.GS"),
            perms.to.show = c(1:10), legend.out= TRUE, perms = NULL, GS = NULL,
            ev.info = NULL, sigSize, ...){

    ln <- length(sigSize)
    perms.to.show <-  perms.to.show[perms.to.show <= dim(perms)[2]]
    lp <- length(perms.to.show)

    if("perm.violin" %in% plot.type){
       ylims <- max(c(vapply(perm.out$new.coefs.nda.ori, function(x)
         max(abs(x[1,,4])), double(1)),vapply(perm.out$new.coefs.nda.perm, function(x)
           max(abs(x[perms.to.show,,4])),double(1))))
       ylims <- c(ylims + ylims/20)
       if(legend.out) par(mar=c(5.1, 5.1, 4.1, 9.1), xpd=TRUE)
       colss <- grey.colors(ln)[ln:1]
       plot(1, 1, xlim = c(0,1), ylim =c(-ylims, ylims), type = 'n',
            #xlab = 'permuted outcome instance', ylab = 'lHR.stat',
            xaxt="n", ...)
       for(L in seq_len(ln))
       {
         for(i in seq_len(lp))
            vioplot(perm.out$new.coefs.nda.perm[[L]][perms.to.show[i],,4],
                  at = i/(lp+1) -  1/(2*(lp+1)), drawRect= FALSE, add = TRUE,
                  col = rgb(1,1,1,alpha=0.3), border=colss[L], lty=L, wex = .1, lwd=2,
                  colMed = colss[L])
            vioplot(perm.out$new.coefs.nda.ori[[L]][1,,4], at=1 -1/(2*(lp+1)),
               add = TRUE,  col = rgb(1,1,1,alpha=0.3), drawRect= FALSE,
               border=brewer.pal(ln, "Greens")[L], lty=L, wex = .1, lwd=2,
               colMed = brewer.pal(ln, "Greens")[ln+1-L])
           }
       if(!is.null(perm.out$new.coefs.nda.perm.GS)){
          text(seq_len(lp)/(lp+1) -  1/(2*(lp+1)),
               perm.out$new.coefs.nda.perm.GS[perms.to.show,,4], "GS",
               col = 2)
          text(1 -1/(2*(lp+1)), perm.out$new.coefs.nda.ori.GS[1,,4], "GS",
                col = 2)
       }
       text(1 -1/(2*(lp+1)), -ylims, "original")
       coord <- par("usr")
       lines(c(coord[1],coord[2]),c(0,0),col=5, lty=2)
    #   abline(h=0,col=5, lty=2)
       if(legend.out)
         legend(x = coord[2] * 1.05, y = coord[4], legend=sigSize, col = colss,
                lty = seq_len(ln), title = "signature size", cex = 1.2)
   }

   if(any(c("perm.GSvsEvents","perms.corr.GS") %in% plot.type))
      shuff.ev <- apply(perms, 1, function(o) ev.info[o])
   if("perm.GSvsEvents" %in% plot.type){
    if(legend.out) par(mar=c(5.1, 5.1, 4.1, 9.1), xpd=TRUE)
    ylims <- c(min(GS) - (max(GS)-min(GS))/20, max(GS) + (max(GS)-min(GS))/20)
    plot(1, 1, xlim = c(0,1), ylim =ylims, type = 'n',
      #   xlab = 'permuted outcome instance', ylab = 'Global signature',
         xaxt = "n", ...)
    for(i in seq_len(lp)){
      boxplot(GS[which(shuff.ev[i,] ==0)],
          at = ((i-1)*2 +1)/((lp+1)*2) -  1/(4*(lp+1)), add = TRUE,
          col = "green", border="darkgreen", lty=1, boxwex = 0.05,
          xaxt ="n", yaxt= "n")

      boxplot(GS[which(shuff.ev[i,] ==1)],
          at = ((i-1)*2 +2)/((lp+1)*2) -  1/(4*(lp+1)),
          drawRect = FALSE, add = TRUE, col = "blue", border="darkblue", lty=1,
          wex = .1, lwd=2, colMed = "black", boxwex = 0.05, xaxt ="n",
          yaxt= "n")
      lines(c(((i-1)*2 +1)/((lp+1)*2) -  1/(4*(lp+1)),((i-1)*2 +2)/((lp+1)*2) -
            1/(4*(lp+1))), c(median(GS[which(shuff.ev[i,] ==0)]),
            median(GS[which(shuff.ev[i,] ==1)])), col=2, lty=1, lwd=2)
    }

    boxplot(GS[which(ev.info == 0)], at = 1 - 3/(4*(lp+1)) , drawRect= FALSE,
         add = TRUE, col = "green", border = "darkgreen", lty=1, wex = .1,
         lwd = 2, colMed = "black", boxwex = 0.05, xaxt ="n", yaxt= "n")
    boxplot(GS[which(ev.info  == 1)], at = 1 -   1/(4*(lp+1)), drawRect= FALSE,
         add = TRUE, col = "blue", border="darkblue", lty=1, wex = .1, lwd=2,
         colMed = "black", boxwex = 0.05, xaxt ="n", yaxt= "n")
    lines(c( 1 - 3/(4*(lp+1)),  1 -  1/(4*(lp+1))),
          c(median(GS[which(ev.info == 0)]), median(GS[which(ev.info == 1)])),
          col = 2, lty = 1, lwd = 2)
    text(c(2 - 3/(4*(lp+1)) -  1/(4*(lp+1)))/2, ylims[1], "original")
    coord <- par("usr")
    if(legend.out)
     legend(x = coord[2] * 1.05, y = coord[4], legend=c("no event","event"),
            col = c("green", "blue"), pch = c(15,15), cex = 1.2)
   }
   if("perms.corr.GS" %in% plot.type){
     if(legend.out) par(mar=c(5.1, 5.1, 4.1, 9.1), xpd=TRUE)
     cors <- array(0,dim=c(dim(perms)[2],2,ln))
     for(L in seq_len(ln)){
       cors[,,L] <- t(sapply(seq_len(dim(perms)[2]), function(i){
                c(mean(t(perm.out$new.coefs.nda.perm[[L]][,,4])[,i]),
                c(mean(GS[which(shuff.ev[i,] ==1)]) -
                   mean(GS[which(shuff.ev[i,] ==0)])))
            }))
       }
     colss <- grey.colors(ln)[ln:1]
     L <- 2
     plot(cors[,1,1], cors[,2,1],
#          xlab = "mean lHR.stat", ylab = "diff GS (mean.ev - mean.noev)",
          col = colss[1], xlim = c(min(cors[,1,]), max(cors[,1,])),
           ylim = c(min(cors[,2,]), max(cors[,2,])), pch = 15+1, ...)
     while(L <=ln){
         points(cors[,1,L], cors[,2,L], col = colss[L], pch = 15+L)
         L <- L+1
     }
     coord <- par("usr")
     if(legend.out)
      legend(x = coord[2] * 1.05, y = coord[4],
        legend = paste0(sigSize," / ", round(apply(cors,3,function(x)
           cor(x)[1,2]),2)), col = colss, pch = seq_len(ln)+15,
        title = "sign size / cor", cex = 1.2)
    }
#   par(mfrow=c(1,1))

}