### R code from vignette source 'vignettes/SeqGSEA/inst/doc/SeqGSEA.Rnw'

###################################################
### code chunk number 1: SeqGSEA
###################################################
  library(SeqGSEA)


###################################################
### code chunk number 2: help (eval = FALSE)
###################################################
## ? SeqGSEA


###################################################
### code chunk number 3: ReadCountSet
###################################################
  rcounts <- cbind(t(sapply(1:10, function(x) {rnbinom(5, size=10, prob=runif(1))} )), 
                   t(sapply(1:10, function(x) {rnbinom(5, size=10, prob=runif(1))} )))
  colnames(rcounts) <- c(paste("S", 1:5, sep=""), paste("C", 1:5, sep="")) 
  geneIDs <- c(rep("G1", 4), rep("G2", 6))
  exonIDs <- c(paste("E", 1:4, sep=""), paste("E", 1:6, sep=""))
  RCS <- newReadCountSet(rcounts, exonIDs, geneIDs)
  RCS 


###################################################
### code chunk number 4: RCS_example
###################################################
  data(RCS_example, package="SeqGSEA")
  RCS_example


###################################################
### code chunk number 5: geneID
###################################################
  length(unique(geneID(RCS_example)))


###################################################
### code chunk number 6: exonTestability
###################################################
  RCS_example <- exonTestability(RCS_example, cutoff = 5)


###################################################
### code chunk number 7: DSanalysis
###################################################
  RCS_example <- estiExonNBstat(RCS_example)
  RCS_example <- estiGeneNBstat(RCS_example)
  head(fData(RCS_example)[, c("exonIDs", "geneIDs", "testable", "NBstat")])


###################################################
### code chunk number 8: DSperm
###################################################
  permuteMat <- genpermuteMat(RCS_example, times=20)
  RCS_example <- DSpermute4GSEA(RCS_example, permuteMat)
  head(RCS_example@permute_NBstat_gene)


###################################################
### code chunk number 9: DSscores
###################################################
  DSscore.normFac <- normFactor(RCS_example@permute_NBstat_gene)
  DSscore <- scoreNormalization(RCS_example@featureData_gene$NBstat, 
                                DSscore.normFac)
  DSscore.perm <- scoreNormalization(RCS_example@permute_NBstat_gene, 
                                     DSscore.normFac)
  DSscore[1:5]
  DSscore.perm[1:5,1:10]


###################################################
### code chunk number 10: DSpval
###################################################
  RCS_example <- DSpermutePval(RCS_example, permuteMat)
  head(DSresultGeneTable(RCS_example))


###################################################
### code chunk number 11: genecount
###################################################
  geneCounts <- getGeneCount(RCS_example)
  dim(geneCounts) # 182  20
  head(geneCounts)


###################################################
### code chunk number 12: DEseq
###################################################
  label <- label(RCS_example)
  DEG <- runDESeq(geneCounts, label)


###################################################
### code chunk number 13: DSNB
###################################################
  DEGres <- DENBStat4GSEA(DEG)
  DEGres[1:5, "NBstat"]


###################################################
### code chunk number 14: DSNBpermut
###################################################
  DEpermNBstat <- DENBStatPermut4GSEA(DEG, permuteMat)
  DEpermNBstat[1:5, 1:10]


###################################################
### code chunk number 15: DEscores
###################################################
  DEscore.normFac <- normFactor(DEpermNBstat)
  DEscore <- scoreNormalization(DEGres$NBstat, DEscore.normFac)
  DEscore.perm <- scoreNormalization(DEpermNBstat, DEscore.normFac)
  DEscore[1:5]
  DEscore.perm[1:5, 1:10]


###################################################
### code chunk number 16: DEpval
###################################################
  DEGres <- DEpermutePval(DEGres, DEpermNBstat) 
  DEGres[1:6, c("NBstat", "perm.pval", "perm.padj")]


###################################################
### code chunk number 17: DEpval2
###################################################
  DEGres <- DENBTest(DEG)
  DEGres <- DEpermutePval(DEGres, DEpermNBstat) 
  DEGres[1:6, c("NBstat", "pval", "padj", "perm.pval", "perm.padj")]


###################################################
### code chunk number 18: genescore_l
###################################################
  gene.score <- geneScore(DEscore, DSscore, method="linear", DEweight = 0.3)
  gene.score.perm <- genePermuteScore(DEscore.perm, DSscore.perm, 
                                      method="linear",  DEweight=0.3)
  plotGeneScore(gene.score, gene.score.perm)


###################################################
### code chunk number 19: genescore_r
###################################################
  gene.score <- geneScore(DEscore, DSscore, method="rank", DEweight = 0.3)
  gene.score.perm <- genePermuteScore(DEscore.perm, DSscore.perm, 
                                      method="rank",  DEweight=0.3)
  plotGeneScore(gene.score, gene.score.perm)


###################################################
### code chunk number 20: genescore_gr
###################################################
  combine <- rankCombine(DEscore, DSscore, DEscore.perm, DSscore.perm, DEweight=0.3) 
  gene.score <- combine$geneScore
  gene.score.perm <- combine$genePermuteScore
  plotGeneScore(gene.score, gene.score.perm)


###################################################
### code chunk number 21: seqgeneset
###################################################
  data(GS_example, package="SeqGSEA") 
  GS_example


###################################################
### code chunk number 22: GSEAmain
###################################################
  GS_example <- GSEnrichAnalyze(GS_example, gene.score, gene.score.perm)
  topGeneSets(GS_example, 5)


###################################################
### code chunk number 23: plotES
###################################################
  plotES(GS_example)


###################################################
### code chunk number 24: plotSig
###################################################
  plotSig(GS_example)


###################################################
### code chunk number 25: plotSigGS
###################################################
  plotSigGeneSet(GS_example, 9, gene.score) # 9th gene set is the most significant one.


###################################################
### code chunk number 26: wrightSigGS
###################################################
  writeSigGeneSet(GS_example, 9, gene.score) # 9th gene set is the most significant one.


###################################################
### code chunk number 27: doParallel1
###################################################
  library(doParallel)
  a <- matrix(1:16, 4, 4)
  b <- t(a)
  foreach(b=iter(b, by='col'), .combine=cbind) %dopar%
    (a %*% b)


###################################################
### code chunk number 28: doParallel2
###################################################
  library(doParallel)
  cl <- makeCluster(2) # specify 2 cores to be used in this computing 
  registerDoParallel(cl)
  getDoParWorkers() # 2
  a <- matrix(1:16, 4, 4)
  b <- t(a)
  foreach(b=iter(b, by='col'), .combine=cbind) %dopar%
    (a %*% b)


###################################################
### code chunk number 29: sysfile
###################################################
system.file("extscripts", package="SeqGSEA", mustWork=TRUE)


###################################################
### code chunk number 30: Initialization
###################################################
  rm(list=ls())
  # input count data files
  data.dir <- system.file("extdata", package="SeqGSEA", mustWork=TRUE)
  case.pattern <- "^SC"  # file name starting with "SC"
  ctrl.pattern <- "^SN"  # file name starting with "SN"
  case.files <- dir(data.dir, pattern=case.pattern, full.names = TRUE)
  control.files <- dir(data.dir, pattern=ctrl.pattern, full.names = TRUE)
  # gene set file
  geneset.file <- system.file("extdata", "gs_symb.txt",  
                              package="SeqGSEA", mustWork=TRUE)
  # output file prefix 
  output.prefix <- "SeqGSEA.test"
  # setup parallel backend 
  library(doParallel)
  cl <- makeCluster(2) # specify 2 cores to be used in computing 
  registerDoParallel(cl) # parallel backend registration
  # setup permutation times
  perm.times <- 20 # change the number to >= 1000 in your analysis 


###################################################
### code chunk number 31: DS_analysis
###################################################
  # load exon read count data
  RCS <- loadExonCountData(case.files, control.files)
  # remove genes with low exprssion 
  RCS <- exonTestability(RCS, cutoff=5)
  geneTestable <- geneTestability(RCS)
  RCS <- subsetByGenes(RCS, unique(geneID(RCS))[ geneTestable ])
  # get gene IDs, which will be used in initialization of gene set
  geneIDs <- unique(geneID(RCS))
  # calculate DS NB statistics
  RCS <- estiExonNBstat(RCS)
  RCS <- estiGeneNBstat(RCS)
  # calculate DS NB statistics on the permutation data sets
  permuteMat <- genpermuteMat(RCS, times=perm.times)
  RCS <- DSpermute4GSEA(RCS, permuteMat)


###################################################
### code chunk number 32: DE_analysis
###################################################
  # get gene read counts
  geneCounts <- getGeneCount(RCS)
  # calculate DE NB statistics 
  label <- label(RCS)
  DEG <-runDESeq(geneCounts, label)
  DEGres <- DENBStat4GSEA(DEG)
  # calculate DE NB statistics on the permutation data sets
  DEpermNBstat <- DENBStatPermut4GSEA(DEG, permuteMat) # permutation


###################################################
### code chunk number 33: score_int
###################################################
  # DE score normalization
  DEscore.normFac <- normFactor(DEpermNBstat)
  DEscore <- scoreNormalization(DEGres$NBstat, DEscore.normFac)
  DEscore.perm <- scoreNormalization(DEpermNBstat, DEscore.normFac)
  # DS score normalization
  DSscore.normFac <- normFactor(RCS@permute_NBstat_gene)
  DSscore <- scoreNormalization(RCS@featureData_gene$NBstat, DSscore.normFac)
  DSscore.perm <- scoreNormalization(RCS@permute_NBstat_gene, DSscore.normFac)
  # score integration
  gene.score <- geneScore(DEscore, DSscore, DEweight=0.5)
  gene.score.perm <- genePermuteScore(DEscore.perm, DSscore.perm, DEweight=0.5)
  # visilization of scores 
  # NOT run in the example; users to uncomment the following 6 lines to run 
  #plotGeneScore(DEscore, DEscore.perm, pdf=paste(output.prefix,".DEScore.pdf",sep=""),
  #              main="Expression")
  #plotGeneScore(DSscore, DSscore.perm, pdf=paste(output.prefix,".DSScore.pdf",sep=""), 
  #              main="Splicing")
  #plotGeneScore(gene.score, gene.score.perm, 
  #              pdf=paste(output.prefix,".GeneScore.pdf",sep=""))


###################################################
### code chunk number 34: main_gsea
###################################################
  # load gene set data
  gene.set <- loadGenesets(geneset.file, geneIDs, geneID.type="ensembl",
                             genesetsize.min = 5, genesetsize.max = 1000)
  # enrichment analysis
  gene.set <- GSEnrichAnalyze(gene.set, gene.score, gene.score.perm, weighted.type=1)
  # format enrichment analysis results
  GSEAres <- GSEAresultTable(gene.set, TRUE)
  # output results 
  # NOT run in the example; users to uncomment the following 4 lines to run 
  #write.table(GSEAres, paste(output.prefix,".GSEA.result.txt",sep=""), 
  #            quote=FALSE, sep="\t", row.names=FALSE)
  #plotES(gene.set, pdf=paste(output.prefix,".GSEA.ES.pdf",sep=""))
  #plotSig(gene.set, pdf=paste(output.prefix,".GSEA.FDR.pdf",sep=""))


###################################################
### code chunk number 35: sessionInfo
###################################################
  sessionInfo()


###################################################
### code chunk number 36: <closeConnetions
###################################################
allCon <- showConnections()
socketCon <- as.integer(rownames(allCon)[allCon[, "class"] == "sockconn"])
sapply(socketCon, function(ii) close.connection(getConnection(ii)) )