R/plotSubSEScurve.R

In psSubpathway: Flexible Identification of Phenotype-Specific Subpathways

##'plotSubSEScurve
##'
##'
##' @title Plot subtype set sample enrichment score curve graph
##' @description Draw a sample enrichment score curve graph of a single or all subtypes.
##' @param inputdata A list of result data generated by function `SubSEA` or `DCSA`.
##' @param spwid The subpathway id which the user wants to plot.
##' @param phenotype The `phenotype`` specifies which phenotypic phenotype set enrichment curve is drawn for subpathway.
##' When `phenotype="all"` (default), the phenotype set enrichment score curve graph of subpathway under all phenotypes is drawn.
##' @details
##' Plot a phenotype set enrichment score curve graph of a subpathway under all phenotypes or specified phenotype, including the
##' location of the maximum enrichment score (ES) and the leading-edge subset. This function can only be used for the results of
##' the `SubSEA` function.
##' @return a plot
##' @importFrom graphics layout
##' @importFrom graphics plot
##' @importFrom graphics abline
##' @importFrom graphics lines
##' @importFrom grDevices colors
##' @importFrom graphics text
##' @return an enrichment score curve graph
##' @author Xudong Han,
##' Junwei Han,
##' Qingfei Kong
##' @examples
##' # get the results of the SubSEA function for breast cancer subtypes.
##' Subspwresult<-get("Subspwresult")
##' # plot enrichment score curve of the subpathway 00120_9 in all breast cancer subtypes.
##' plotSubSEScurve(Subspwresult,spwid="00120_9",phenotype="all")
##' # plot enrichment score curve of the subpathway 00120_9 in the basal breast cancer subtypes.
##' plotSubSEScurve(Subspwresult,spwid="00120_9",phenotype="Basal")
##' @export

plotSubSEScurve<-function(inputdata,spwid="",phenotype="all"){

  spwmatrix<-inputdata$spwmatrix
  class.labels <-colnames(spwmatrix)
  class.phen <- names(table(class.labels))


  spwmatrix.r<-length(spwmatrix[,1])
  spwmatrix.c<-length(spwmatrix[1,])
  spw.rname<-row.names(spwmatrix)

  ESdata<-list()
  for(j in 1:length(class.phen)){
    inphen<-class.phen[j]
    ident<-class.labels
    ident[ident!=inphen]<-0
    ident[ident==inphen]<-1
    ident<-as.numeric(ident)

    Obs.ES <- vector(length = spwmatrix.r, mode = "numeric")
    Obs.arg.ES <- vector(length = spwmatrix.r, mode = "numeric")
    Obs.indicator <- matrix(nrow= spwmatrix.r, ncol=spwmatrix.c)
    Obs.RES <- matrix(nrow= spwmatrix.r, ncol=spwmatrix.c)
    rowve<-vector(length = spwmatrix.c,mode="numeric")
    for(e in 1:spwmatrix.r){
      rowve<-spwmatrix[e,]
      ordindex<-order(rowve,decreasing = T)
      correl<-rowve[ordindex]
      sampleorder<-ident[ordindex]
      o<-SEAscore(sampleorder,correl.vector =correl)
      Obs.ES[e] <- o$ES
      Obs.arg.ES[e] <- o$arg.ES
      Obs.RES[e ,] <- o$RES
      Obs.indicator[e ,] <- o$indicator
    }
    ESdata[[j]]<-list(Obs.ES=Obs.ES,Obs.arg.ES=Obs.arg.ES,Obs.RES=Obs.RES,Obs.indicator=Obs.indicator)
  }
  names(ESdata)<-class.phen


  spwmatrix.r <- length(spwmatrix[,1])
  spwmatrix.c <- length(spwmatrix[1,])
  spwmatrix.rnames<-row.names(spwmatrix)
  spwmatrix.cnames<-colnames(spwmatrix)
  pl<-length(ESdata)
  gsindex<-which(spwmatrix.rnames==spwid)

  c<-apply(spwmatrix,1,scale)
  c<-t(apply(c,2,sort,decreasing = T))
  colnames(c)<-spwmatrix.cnames
  zhindex<-which(c[gsindex,]<0)[1]


  if(phenotype=="all"){
    layout(matrix(c(1:pl),nrow = 1,ncol = pl,byrow=T))


    for(i in 1:pl){
      topindex<-ESdata[[i]][["Obs.arg.ES"]][gsindex]
      ind <- 1:spwmatrix.c
      obs.s2n<-as.numeric(c[gsindex,])
      RE.max<-max(ESdata[[i]][["Obs.RES"]][gsindex,])
      RE.min<-min(ESdata[[i]][["Obs.RES"]][gsindex,])
      tag.sj<-RE.min-0.1
      tag.xj<-RE.min-0.6
      maxs2n<-max(obs.s2n)
      Obs.correl.vector.norm<-obs.s2n-(maxs2n-RE.min+0.7)
      s2n.sj<-max(Obs.correl.vector.norm)
      s2n.zx<- -(maxs2n-RE.min+0.7)
      s2n.xj<-min(Obs.correl.vector.norm)

      labnumber<-length(which(ESdata[[i]][["Obs.indicator"]][1,]==1))

      sub.string <- paste("Number of samples: ", spwmatrix.c, " (in list), ", labnumber, " (in phenotype set)", sep = "", collapse="")
      main.string <- paste("phenotype Set ",  ":", names(ESdata)[i])
      plot(ind, ESdata[[i]][["Obs.RES"]][gsindex,], main = main.string, sub = sub.string, xlab = spwmatrix.rnames[gsindex], ylab = "Running Enrichment Score (RES)", ylim=c(s2n.xj-0.3,RE.max),type = "o", lwd = 2, cex = 1, col = (i+1))
      abline(h=0,lty = 1)
      abline(h=s2n.zx)
      abline(v=zhindex,lty = 2)
      abline(v=topindex,lty = 6)
      for (j in 1:spwmatrix.c) {
        if (ESdata[[i]][["Obs.indicator"]][gsindex, j] == 1) {
          lines(c(j, j), c(tag.sj, tag.xj), lwd = 1, lty = 1, cex = 1, col = 1)
        }
      }

      for (j in seq(1, spwmatrix.c, 1)) {
        lines(c(j, j), c(s2n.zx, Obs.correl.vector.norm[j]), lwd = 1, cex = 1, col = colors()[12])
      }
      leg.txt <- paste("\"", "High", "\" ", sep="", collapse="")
      text(x=1, y=s2n.xj-0.3, adj = c(0, 0), labels=leg.txt, cex = 1.0)

      leg.txt <- paste("\"", "Low", "\" ", sep="", collapse="")
      text(x=spwmatrix.c, y=s2n.xj-0.3, adj = c(1, 0), labels=leg.txt, cex = 1.0)
    }



  }
  else{
    topindex<-ESdata[[phenotype]][["Obs.arg.ES"]][gsindex]
    ind <- 1:spwmatrix.c
    phenindex<-which(names(ESdata)==phenotype)
    obs.s2n<-as.numeric(sort(spwmatrix[gsindex,]))
    RE.max<-max(ESdata[[phenindex]][["Obs.RES"]][gsindex,])
    RE.min<-min(ESdata[[phenindex]][["Obs.RES"]][gsindex,])
    tag.sj<-RE.min-0.1
    tag.xj<-RE.min-0.6
    maxs2n<-max(obs.s2n)
    Obs.correl.vector.norm<-obs.s2n-(maxs2n-RE.min+0.7)
    s2n.sj<-max(Obs.correl.vector.norm)
    s2n.zx<- -(maxs2n-RE.min+0.7)
    s2n.xj<-min(Obs.correl.vector.norm)

    labnumber<-length(which(ESdata[[phenindex]][["Obs.indicator"]][1,]==1))
    col <- which(names(table(names(ESdata)))==phenotype)
    sub.string <- paste("Number of samples: ", spwmatrix.c, " (in list), ", labnumber, " (in phenotype set)", sep = "", collapse="")
    main.string <- paste("phenotype Set ",  ":", phenotype)
    plot(ind, ESdata[[phenindex]][["Obs.RES"]][gsindex,], main = main.string, sub = sub.string, xlab = spwmatrix.rnames[gsindex], ylab = "Running Enrichment Score (RES)", ylim=c(s2n.xj-0.3,RE.max),type = "o", lwd = 2, cex = 1, col = (col+1))
    abline(h=0,lty = 1)
    abline(h=s2n.zx)
    abline(v=zhindex,lty = 2)
    abline(v=topindex,lty = 6)
    for (j in 1:spwmatrix.c) {
      if (ESdata[[phenindex]][["Obs.indicator"]][gsindex, j] == 1) {
        lines(c(j, j), c(tag.sj, tag.xj), lwd = 1, lty = 1, cex = 1, col = 1)  # enrichment tags
      }
    }

    for (j in seq(1, spwmatrix.c, 1)) {
      lines(c(j, j), c(s2n.zx, Obs.correl.vector.norm[j]), lwd = 2, cex = 1, col = colors()[12]) # shading of correlation plot
    }
    leg.txt <- paste("\"", "High", "\" ", sep="", collapse="")
    text(x=1, y=s2n.xj-0.3, adj = c(0, 0), labels=leg.txt, cex = 1.0)

    leg.txt <- paste("\"", "Low", "\" ", sep="", collapse="")
    text(x=spwmatrix.c, y=s2n.xj-0.3, adj = c(1, 0), labels=leg.txt, cex = 1.0)
  }
}