GSEA enrichment plot

#' Calculate Enrichment score for GSEA plot

#' 

#' This function is slightly modified from GSEA.1.0.R downloaded from msigDB.\cr

#' Specific site address is below.\cr

#' http://software.broadinstitute.org/gsea/msigdb/download_file.jsp?filePath=/resources/software/GSEA-P-R.1.0.zip\cr

#' Original functions name are 'GSEA.EnrichmentScore' and 'GSEA.GeneRanking'\cr

#' I combined two functions into single function calculate enrichment score vector.\cr

#' It won't produce any statics including p-values and normalized enrichment score. Parameters are default.

#' @param exp.mat(default=NULL) : Matrix/data.frame of gene expression values (rows are genes, columns are samples) 

#' @param cls(default=NULL) : List type of class\cr

#' (ex. list(GPC1=c('A','B','C','D'),GPC2=c('E','F','G','H')))

#' @param geneSet(default=NULL) :  GeneSet list\cr

#' (ex. list(GeneSetA=c('geneA','geneB','geneC','geneD'),GeneSetB=c('geneE','geneF','geneG')))

#' @keywords enrichment score, GSEA, ES

#' @export

#' @examples

#' exp.mat=matrix(rnorm(10000),1000,10)

#' rownames(exp.mat)=paste0('gene',1:nrow(exp.mat))

#' colnames(exp.mat)=paste0('sample',1:ncol(exp.mat))

#' genes=rownames(exp.mat)

#' geneSet=list('A'=sample(genes,100),'B'=sample(genes,50),'C'=sample(genes,150))

#' samples=sample(colnames(exp.mat),10)

#' cls=list('groupA'=samples[1:5],'groupB'=samples[6:10])

#' draw.GSEAplot(geneSet=geneSet,exp.mat=exp.mat,cls=cls,color=c('red','pink','skyblue','blue'),line.color='green')

#===================================

# Drawing GSEA enrichment plot

# This script is not calculate p-value for GSEA

# It only provide plot for selected genesets

#===================================

draw.GSEAplot=function(geneSet=NULL,exp.mat=NULL,cls=NULL,color=c('red','pink','skyblue','blue'),line.color='green'){

  exp.mat=as.matrix(exp.mat[,unlist(cls)])

  # Calculate signal2Noise across all genes (Group1.mean-Group2.mean)/(Group1.std+Group2.std)

    message('Calculate Signal2Noise and sort table')

    EScore=calc.GSEA_EScore(geneSet = geneSet,exp.mat = exp.mat,cls = cls)

  

  # Drawing part

    # Set layout

      layout.mat=as.matrix(c(1:4))

      layout(layout.mat,heights = c(3,1,0.5,2))

      l_ply(1:length(EScore$ES),.progress='text',function(k){

        print(k)

        ES=EScore$ES[[k]]

        # Mount plot of Enrichment score

          par(mar=c(0,6,3,1))

          plot(ES$EnrichmentScore,ylab='Enrichment score',type='l',col=line.color,lwd=5,ann=T,xaxt='n',

               main=names(EScore$ES)[[k]])

          

          grid(lwd=2)

        # tag plot

          par(mar=c(0,6,0,1))

          plot(NULL,xlim=c(0,length(ES[[1]])),ylim=c(0,1),ann=F,xaxt='n',yaxt='n')

          abline(v=ES$Location.gene.tag)

        # Gradient plot

          col.vec=1:14

          mat=t(rbind(col.vec,col.vec,col.vec))

          image(mat,col=colorRampPalette(color)(14),xaxt='n',yaxt='n')

        # Signal2Noise plot

          par(mar=c(3,6,0,1))

          plot(1:length(EScore$S2N),EScore$S2N,ylab='Signal2Noise',xaxt='n',type='h')

          axis(side=1,at = seq(0,to=length(EScore$S2N),by=1000))

          grid(lwd=2)

      })

}