_utile4disSuptyper

_utile

#' version 1 (1/31/2018)
#' short description


######### source this ----
rm(list = ls())

#' Get all needed data
LoadMyData <- funtion()
{
  dataFile <- "temp/1-2v3data.RData"
  if( file.exists(dataFile) ){
    load(dataFile, env = globalenv())
  } else {
  #Get w
  load(file="temp/100_NMFsurv.RData", env = globalenv())
  #Get A
  load("results/BRCA15.RData", env = globalenv())
  A <<- Boolcp_3comp$A #patients vs components
  rm(list=ls()[-which(ls() %in% c("w","A"))])

  libs<-c("Packages.R", "permutation.R", "clustering.R", "Plots.R")
  libs<-paste("https://gist.githubusercontent.com/datad/39b9401a53e7f4b44e9bea4d584ac3e8/raw/", libs,sep='')
  sapply(libs, function(u) {source(u)})
  onlyBaseLibs()
  loadls("SNFtool survival ROntoTools org.Hs.eg.db")

  loadlib("plyr", FALSE) #for function count
  loadlib("lmtest", FALSE) #for comparing two models - likelihood
  loadls("survival survAUC", FALSE) #for Cox proportional hazard model
  #For non-negative matrix factorization vvv
  install.packages("digest", repos="http://R-Forge.R-project.org") #required
  loadls("plyr survival lmtest", FALSE)
  loadls("stats registry", FALSE)
  loadls("methods utils pkgmaker registry rngtools cluster", FALSE)
  loadlib("NMF", FALSE) 
  #For non-negative matrix factorization   ^^^^
  loadlib("glmnet") #for lasso regression L1
  
  A <<- A
  
  save.image(file="temp/myData.RData")
  }
}


######### run this ----
LoadMyData()




#' Title
#'
#' @section step 1
#'
#' @examples
foo <- function() {
    
}

#end
sess <- sessionInfo() #save session on variable
save.image(file="temp/x.RData")

clustering.R

#version 2 (1/1/2018)


computeAPN <- function(pathwMolecules, cluster)
{
  measuresAPN <-0
  ## stability validation measures
  for (del in seq_along(pathwMolecules)) {
    molsDel <- pathwMolecules[-del]
    clusterDel <- custerByMols(molsDel)
    stabmeasAPN <- stabilityAPN(cluster, clusterDel)
    measuresAPN <- measuresAPN + stabmeasAPN
  } 
  measuresAPN <- measuresAPN/length(pathwMolecules)
  measuresAPN
}


stabilityAPN <- function(cluster, clusterDel) {

  nc1 <- length(cluster)
  nc2 <- length(clusterDel)
  
  ## measure APN
  ## calculate a ncxnc matrix of proportion of non-overlaps in the two collection of nc clusters
  overlap <- xtabs(~cluster + clusterDel)
  ## measure AD
  ## calculate a ncxnc matrix of average-distance in the two collection of nc clusters
  dij <- matrix(rep(NA,nc1*nc2),nc1,nc2)
  
  
  rs <- matrix(rowSums(overlap),nrow=nrow(overlap),ncol=ncol(overlap),byrow=FALSE)
  #rs <- matrix(colSums(overlap),nrow=nrow(overlap),ncol=ncol(overlap),byrow=TRUE)
  APN1 <- 1-sum(overlap^2/rs)/sum(overlap)
  APN1
}



#an actual pathway
clusterByPathwaySNF2 <- function(pathway.i, surData){
  require(ROntoTools)
  genesOfInterest <- nodes(pathway.i)
  listColVals <- clusterByGenesSNF2(genesOfInterest, surData)
  listColVals
}

#genesOfInterest: are genes and mirs from pathways
clusterByGenesSNF2 <- function(genesOfInterest,surData){
  require(survival)
  listColVals <- list()
  #Number of genes on pathway
  listColVals$nGenes <- length(genesOfInterest)
  genesOfInterest <- geneMIRSetFormat(genesOfInterest)
  
  groups <- getGroupsFromGenesSNF2(genesOfInterest)
  listColVals$groups <- groups
  
  stabilityAPN <- computeAPN(pathwMolecules=genesOfInterest, cluster=groups)
  listColVals$stabilityAPN <- stabilityAPN
  
  #   #begin PREVIOUS METRIC
  #   #Kaplan SNF
  #   fit <- survdiff(Surv(time = Survival, event = Death) ~ groups, data = surData)
  #   listColVals$kaplan <- fit
  #   
  #   #Cox SNF
  #   coxFit <- coxph(Surv(time = Survival, event = Death) ~ groups, 
  #                   data = surData, ties="exact")
  #   resCox <- mySumm(coxFit)
  #   
  #   listColVals$coxPval <- resCox$sctest[3]
  #   listColVals$coxScore <- resCox$sctest[1]
  #   listColVals$resCox <- resCox
  #   #end PREVIOUS METRIC
  
  listColVals
}


#Global = geneEData
#genesOfInterest: are genes from pathways
getGroupsFromGenesSNF2 <- function(genesOfInterest, displayT=FALSE){
  require(SNFtool)
  
  GEofI <- intersect(genesOfInterest,colnames(geneEData))
  MIRofI <- intersect(genesOfInterest,colnames(mirEData))
  
  
  genePath <- geneEData[, GEofI]
  mirPath <- mirEData[, MIRofI]
  
  #2. Calculate the pair-wise distance
  PSMgeneE = dist2(as.matrix(genePath),as.matrix(genePath));
  PSMmirE = dist2(as.matrix(mirPath),as.matrix(mirPath));
  
  #3. construct similarity graphs
  W1 = affinityMatrix(PSMgeneE, K, alpha)
  W2 = affinityMatrix(PSMmirE, K, alpha)
  
  ## These similarity graphs have complementary information about clusters.
  #displayClusters(W1,truelabel);
  #displayClusters(W2,truelabel);
  
  #next, fuse all the graphs
  W = SNF(list(W1,W2), K, T)
  
  ####With this unified graph W of size n x n, you can do either spectral clustering 
  #or Kernel NMF. If you need help with further clustering, please let us know. 
  
  #4. clustering
  groupGE = spectralClustering(W,C);
  
  if(displayT) displayClusters(W,groupGE);
  
  #SNFNMI = Cal_NMI(group, truelabel)
  
  
  ###you can also find the concordance between each individual network and the fused network
  #ConcordanceMatrix = Concordance_Network_NMI(list(W, W1,W2));
  
  
  groups <- factor(groupGE)
  names(groups) <- row.names(geneEData)
  return(groups)
}


#Global var survivalData
thetaSNF2 <- function(x){ 
  require(survival)
  
  indicesOfIPgene <- x[which(x<totalNumberOfGenes)]
  indicesOfIPmir <- x[which(x>totalNumberOfGenes)]
  indicesOfIPmir <- indicesOfIPmir - totalNumberOfGenes
  
  
  genesOfInterest <- allGenes[indicesOfIPgene]
  mirsOfInterest <- allMIRs[indicesOfIPmir]
  molsOfInterest <- c(genesOfInterest,mirsOfInterest)
  groups <- getGroupsFromGenesSNF2(molsOfInterest)
  
  
  ##this was before  
  #   #Kaplan SNF
  #   fit <- survdiff(Surv(time = Survival, event = Death) ~ groups, data = survivalData)
  #   #Cox SNF
  #   coxFit <- coxph(Surv(time = Survival, event = Death) ~ groups, data = survivalData, 
  #                   ties="exact")
  #   resCox <- mySumm(coxFit)
  #   coxPval <- resCox$sctest[3]
  #   coxPval 
  
  stabilityAPN <- computeAPN(molsOfInterest,groups)
  stabilityAPN
}



#Global var survivalData
thetaSNF <- function(x){ 
  require(survival)
  genesOfInterest <- allGenes[x]
  groups <- getGroupsFromGenesSNF(genesOfInterest)
  
  #Kaplan SNF
  fit <- survdiff(Surv(time = Survival, event = Death) ~ groups, data = survivalData)
  #Cox SNF
  coxFit <- coxph(Surv(time = Survival, event = Death) ~ groups, data = survivalData, 
                  ties="exact")
  resCox <- mySumm(coxFit)
  coxPval <- resCox$sctest[3]
  coxPval 
}


#genesOfInterest: are genes from pathways
clusterByGenesSNF <- function(genesOfInterest,surData){
  require(survival)
  listColVals <- list()
  #Number of genes on pathway
  listColVals$nGenes <- length(genesOfInterest)
  genesOfInterest <- geneSetFormat(genesOfInterest)
  groups <- getGroupsFromGenesSNF(genesOfInterest)
  
  #Kaplan SNF
  fit <- survdiff(Surv(time = Survival, event = Death) ~ groups, data = surData)
  listColVals$kaplan <- fit
  
  #Cox SNF
  coxFit <- coxph(Surv(time = Survival, event = Death) ~ groups, 
                  data = surData, ties="exact")
  resCox <- mySumm(coxFit)
  
  listColVals$coxPval <- resCox$sctest[3]
  listColVals$coxScore <- resCox$sctest[1]
  listColVals$resCox <- resCox
  
  listColVals$groups <- groups
  
  listColVals
}

#an actual pathway
clusterByPathwaySNF <- function(pathway.i, surData){
  require(ROntoTools)
  genesOfInterest <- nodes(pathway.i)
  listColVals <- clusterByGenesSNF(genesOfInterest, surData)
  listColVals
}

#before
#  genesOfInterest <- geneSetFormat(genesOfInterest)
#Global = geneEData
#genesOfInterest: are genes from pathways
getGroupsFromGenesSNF <- function(genesOfInterest){
  require(SNFtool)
  genePath <- geneEData[, genesOfInterest]
  ##################
  #2. Calculate the pair-wise distance
  PSMgeneE = dist2(as.matrix(genePath),as.matrix(genePath));
  
  #3. construct similarity graphs
  W1 = affinityMatrix(PSMgeneE, K, alpha)
  
  #Groups with just geneExpression
  groupGE = spectralClustering(W1,C);
  
  groups <- factor(groupGE)
  names(groups) <- row.names(geneEData)
  return(groups)
}


########
#For k-means
########

#before
#  genesOfInterest <- geneSetFormat(genesOfInterest)
#Global = geneEData
#genesOfInterest: are genes from pathways
getGroupsFromGenesK <- function(genesOfInterest){
  genePath <- geneEData[, genesOfInterest]
  
  kresults <- kmeans(genePath, C, iter.max = 10000, nstart = 50)
  groupGE = kresults$cluster
  
  groups <- factor(groupGE)
  return(groups)
}




#Global var survivalData
thetaK <- function(x){ 
  genesOfInterest <- allGenes[x]
  groups <- getGroupsFromGenesK(genesOfInterest)
  
  #Kaplan SNF
  fit <- survdiff(Surv(time = Survival, event = Death) ~ groups, data = survivalData)
  #Cox SNF
  coxFit <- coxph(Surv(time = Survival, event = Death) ~ groups, data = survivalData, 
                  ties="exact")
  resCox <- mySumm(coxFit)
  coxPval <- resCox$sctest[3]
  coxPval 
}



#genesOfInterest: are genes from pathways
clusterByGenesK <- function(genesOfInterest,surData){
  listColVals <- list()
  #Number of genes on pathway
  listColVals$nGenes <- length(genesOfInterest)
  genesOfInterest <- geneSetFormat(genesOfInterest)
  groups <- getGroupsFromGenesK(genesOfInterest)
  
  #Kaplan SNF
  fit <- survdiff(Surv(time = Survival, event = Death) ~ groups, data = surData)
  listColVals$kaplan <- fit
  
  #Cox SNF
  coxFit <- coxph(Surv(time = Survival, event = Death) ~ groups, 
                  data = surData, ties="exact")
  resCox <- mySumm(coxFit)
  
  listColVals$coxPval <- resCox$sctest[3]
  listColVals$coxScore <- resCox$sctest[1]
  listColVals$resCox <- resCox
  
  listColVals$groups <- groups
  
  listColVals
}

#an actual pathway
clusterByPathwayK <- function(pathway.i, surData){
  genesOfInterest <- nodes(pathway.i)
  listColVals <- clusterByGenesK(genesOfInterest, surData)
  listColVals
}

DataFrames.R

#data.frame


# DF=tT
# checkTypeOfDF(DF=controlDat)
# checkTypeOfDF(DF=diseaseDat)
# checkTypeOfDF(DF)
# head(DF)
#   
  
checkTypeOfDF <- function(DF=controlDat) {
  sapply(DF, mode)
  sapply(DF, class)
  #IF WANT TO CHANGE
  #transform(d, char = as.numeric(char))
}

mergeByrowName <- function(controlDat, diseaseDat){
  allSamples <- merge(controlDat, diseaseDat, by="row.names", stringsAsFactors=FALSE )
  allSamples <- allSamples [,-1]
  row.names(allSamples) <- row.names(controlDat)
  allSamples
}


# 
# #create data.frame
# #http://www.dummies.com/how-to/content/how-to-create-a-data-frame-from-scratch-in-r.html
# 
# employee <- c('John Doe','Peter Gynn','Jolie Hope')
# salary <- c(21000, 23400, 26800)
# startdate <- as.Date(c('2010-11-1','2008-3-25','2007-3-14'))
# 
# 
# 
# 
# employ.data <- data.frame(employee, salary, startdate)
# 
# #create empty data frame
# #http://stackoverflow.com/questions/10689055/create-an-empty-data-frame
# df <- data.frame(Date=as.Date(character()),
#                  File=character(), 
#                  User=character(), 
#                  stringsAsFactors=FALSE)
# 
# newrow = c("1/1/2010","thefile","dd")
# df <- rbind(df,newrow)
# 
# #read
# 
# 
# head(employ.data)
# 
# employ.data[,c('employee','salary')]
# 
# 
# employ.data[1:5,c('employee','salary')]
# employ.data[1:length(employ.data[,1]),c('employee','salary')]
# 
# 
# 
# #add a column with a list of two columns
# #in data table
# totalCounts[,pValue:= pVHyper(freqS,freqU) ]
# 
# #in dataframe
# ##http://stackoverflow.com/questions/3651651/adding-a-column-to-a-dataframe-in-r
# #https://stat.ethz.ch/pipermail/r-help/2002-April/020399.html
# 
# DF <- data.frame(start=c(1,3,5,7), end=c(2,6,7,9))
# DF$newcol <- apply(DF,1,function(row) mean( row[1] : row[2]  ))
# data$Z <- data$X + 5 * data$Y
# 
# 
# 
# #Add a Column to a Dataframe From a List of Values
# #http://stackoverflow.com/questions/11130037/add-a-column-to-a-dataframe-from-a-list-of-values
# 
# columnsNames = LETTERS[1:10]
# nRows = 2
# 
# num.iters <- 10
# l <- vector('list', num.iters) 
# for (i in 1:num.iters) {
#   l[[i]] <- rep("NA",nRows)                      # the column data
#   names(l)[i] <- columnsNames[i]  # the column name 
# } 
# do.call(cbind, l)  # ... if your cols are the same datatype and you want a matrix
# data.frame(l)      # otherwise
# 
# 
# #DIVIDE THE DATA FRAME IN TWO ONE WITH COLUMNS EMPTY
# 
# x.sub <- subset(DF, interaction == '')
# 
# #get a vector for data set
# #http://stackoverflow.com/questions/2545228/converting-a-dataframe-to-a-vector-by-rows
# 
# test <- data.frame(x=c(26,21,20),y=c(34,29,28))
# vectorCOlum = as.vector(as.matrix(test))
# vectorRow =as.vector(t(test))
# 
# vectorColum1 = as.vector(as.matrix(test[,'x']))
# vectorRow1 = as.vector(as.matrix(test[1,]))
# 
# 
# #Change from factor to char
# #http://stackoverflow.com/questions/2851015/convert-data-frame-columns-from-factors-to-characters
# 
# employ.data[1,1]
# 
# class(employ.data[1,1])
# 
# employ.data <- data.frame(lapply(employ.data, as.character), stringsAsFactors=FALSE)
# 
# class(employ.data[1,1])
# 
# 
# #change value of column
# #http://stackoverflow.com/questions/12888363/how-to-change-column-values-in-a-data-frame
# employ.data <- within(employ.data, salary <- salary+1000)
# 
# employ.data[,2] = log(employ.data[,2])
# 
# #with a new function declared by us
# #http://stackoverflow.com/questions/11330138/finding-columns-with-all-missing-values-in-r
# interactionsTab[,'interaction'] <- sapply(interactionsTab[,'interaction'], interaction2weight )
# 
# 
# 
# #delete or drop columns
# #http://stackoverflow.com/questions/4605206/drop-columns-r-data-frame
# DF <- data.frame(
#   x=1:10,
#   y=10:1,
#   z=rep(5,10),
#   a=11:20
# )
# drops <- c("x","z")
# DF[,c("a","y")]
# 
# 
# 
# #rename
# #http://www.cookbook-r.com/Manipulating_data/Renaming_columns_in_a_data_frame/
# names(DF) [names(DF)=="a"] <- "two"
# 
# 
# #add a row to the dataframe
# #http://stackoverflow.com/questions/11561856/add-new-row-to-dataframe
# 
# #In a r possition
# 
# existingDF <- as.data.frame(matrix(seq(20),nrow=5,ncol=4))
# r <- 3
# newrow <- seq(4)
# existingDF <- rbind(existingDF[1:r,],newrow,existingDF[-(1:r),])
# 
# 
# #at the end
# #it is the same for adding dataframes at the end
# newrow = c(7:10)
# existingDF = rbind(existingDF,newrow)
# 
# 
# #transpose data frame 
# #http://stackoverflow.com/questions/6778908/r-transposing-a-data-frame
# 
# DF<- as.data.frame(t(DF))
# 
# #delete the name of the rows
# #http://r.789695.n4.nabble.com/Remove-row-names-column-in-dataframe-td856647.html
# 
# row.names(DF) <- NULL 
# 
# #Extract or Replace Parts of a Data Frame
# #http://stat.ethz.ch/R-manual/R-devel/library/base/html/Extract.data.frame.html
# 
# 
# 
# #combine columns
# 
# m <- cbind(1, 1:7) # the '1' (= shorter vector) is recycled
# 
# 
# 
# #merge columns
# #https://stat.ethz.ch/pipermail/r-help/2011-June/280275.html
# 
# 
# prefix <- c("cheap", "budget")
# roots <- c("car insurance", "auto insurance")
# suffix <- c("quote", "quotes")
# 
# df1 <- data.frame(prefix, roots, suffix)
# 
# newThing <- do.call(c, df1)
# 
# 
# 
# 
# class(newThing)
# 
# 
# dat <- paste(df1[,2], df1[,1], sep=" ") 
# head(dat) 
# dat[,c(1,2)] <- NA 
# head(dat) 
# 
# 
# 
# 
# #removing duplicated
# #http://stats.stackexchange.com/questions/6759/removing-duplicated-rows-from-r-data-frame
# 
# DF <- data.frame(
#   x=1:10,
#   y=10:1,
#   z=rep(5,10),
#   a=11:20
# )
# 
# DF = rbind(DF, c(1,2,3,4))
# DF = rbind(DF, c(1,2,3,4))
# DF = rbind(DF, c(1,2,3,4))
# DF = rbind(DF, c(1,2,3,4))
# 
# duplicated(DF)
# 
# DF[duplicated(DF), ]
# 
# DF[!duplicated(DF), ]
# 
# 
# #remove column
# #http://stackoverflow.com/questions/9202413/how-do-you-delete-a-column-in-data-table
# 
# 
# #GET NUMBER OF ROWS LENGHT SIZE IN DATA TABLE
# #http://stats.stackexchange.com/questions/5253/how-do-i-get-the-number-of-rows-of-a-data-frame-in-r
# 
# nrow(dataset)
# 
# dim(dataset)
# 
# 
# 
# #delete or remove a row
# #http://stackoverflow.com/questions/13520515/command-to-remove-row-from-a-data-frame
# 
# eldNew2 <- tt[-1,]
# dim(eldNew)
# 
# 
# 
# 
# #To find whether a column exists in data frame or not
# if("d" %in% colnames(dat))
# {
#   cat("Yep, it's in there!\n");
# }
# 
# 
# #split data frame
# #http://stackoverflow.com/questions/7069076/split-column-at-delimiter-in-data-frame
# 
# df <- data.frame(ID=11:13, FOO=c('a|b','b|c','x|y'))
# foo <- data.frame(do.call('rbind', strsplit(as.character(df$FOO),'|',fixed=TRUE)))
# 
# 

Packages.R

#version 2 (1/1/2018)

#Install libraries automatically
#Clean screen and delete elements
# rm(list=ls())
# cat("\014")
# source("../Utils/utilitary/Packages.R")
# sourceAFolder("../Utils/utilitary/")

#What is this?
#source("http://bioconductor.org/workflows.R")
# workflowInstall("arrays")

loadlib <- function (libname, bioc=TRUE, verbosel = FALSE){
  #newInstruction <- try( require(libname, character.only = TRUE, 
  #                        quietly = !verbosel , warn.conflicts = verbosel), silent=TRUE ) 
  #if( inherits(newInstruction, "try-error") )
  if(!require(libname, character.only = TRUE, quietly = !verbosel , 
              warn.conflicts = verbosel ) ) #previous statement #dd 05/18
  {
    errorMsn <- paste("Could not load package '",libname,"'.")
    if(readline(prompt = paste(errorMsn,
                               "\nDo you like to installed it from Bioconductor? (y/n)")) == 'y')
    {
      source("http://bioconductor.org/biocLite.R")
      biocLite(libname)
      #biocLite(c("GenomicFeatures", "AnnotationDbi"))
      require(libname, character.only = TRUE) || stop(errorMsn)
    } else if(readline(prompt = paste(errorMsn,
                                      "\nDo you like to installed it from CRAN? (y/n)")) == 'y'){
      install.packages(libname)
      require(libname, character.only = TRUE) || stop(errorMsn)
    }
    else{
      stop(errorMsn)
    }
  } 
  if(bioc)
  {
    source("http://bioconductor.org/biocLite.R")
    message(biocValid())
    message("Above is the output of biocValid()")
  }
}

loadls <- function (libs, bioc=TRUE){
  libs <- unlist(strsplit(libs, " "))
  lapply(libs,loadlib,bioc)
}

Packages_Extra.R

#version 2 (1/1/2018)

#Install libraries automatically
#Clean screen and delete elements
# rm(list=ls())
# cat("\014")
# source("../Utils/utilitary/Packages.R")
# sourceAFolder("../Utils/utilitary/")

#What is this?
#source("http://bioconductor.org/workflows.R")
# workflowInstall("arrays")

loadlib <- function (libname, bioc=TRUE){
  if(!require(libname, character.only = TRUE))
  {
    errorMsn <- paste("Could not load package '",libname,"'.")
    if(readline(prompt = paste(errorMsn,
                               "\nDo you like to installed it from Bioconductor? (y/n)")) == 'y')
    {
      source("http://bioconductor.org/biocLite.R")
      biocLite(libname)
      #biocLite(c("GenomicFeatures", "AnnotationDbi"))
      require(libname, character.only = TRUE) || stop(errorMsn)
    } else if(readline(prompt = paste(errorMsn,
                                      "\nDo you like to installed it from CRAN? (y/n)")) == 'y'){
      install.packages(libname)
      require(libname, character.only = TRUE) || stop(errorMsn)
    }
    else{
      stop(errorMsn)
    }
  }
  if(bioc)
  {
    source("http://bioconductor.org/biocLite.R")
    message(biocValid())
    message("Above is the output of biocValid()")
  }
}

loadls <- function (libs, bioc=TRUE){
  libs <- unlist(strsplit(libs, " "))
  lapply(libs,loadlib,bioc)
}

lspack <- function(libname){
  print(ls(paste("package:",libname, sep='')))
}


#
listLoadedLibs <- function()
{
  print((.packages()))
  (.packages())
}


detachlib <- function(libname){
  try(  detach(paste("package:",libname, sep=''), unload=TRUE, character.only=TRUE) , silent = TRUE)
}

detachlibs <- function(libs){
  libs <- unlist(strsplit(libs, " "))
  lapply(libs,detachlib)
}


onlyBaseLibs <- function(){
  baseLibs <- c("stats", "graphics", "grDevices", "utils", "datasets", "methods", "base")
  ldl <- listLoadedLibs()
  detachlibs(ldl[-which(ldl %in% baseLibs)])
}

searchMethods <- function(methodName) {
  loadlib("sos")
  findFn(methodName, maxPages=2, sortby="MaxScore")
}

sourceAFolder <- function(folder) {
  scripts <- dir(folder, pattern = "*.R")
  scripts <- paste(folder,scripts,sep="/")
  sapply(scripts,source)
}

printRVersion <- function(){
  R.Version()
  getRversion()
}

#update bioconductor after instar new R and new Rstudio
updateBiocon <- function(){
  remove.packages("BiocInstaller")
  source("https://bioconductor.org/biocLite.R")
  biocLite("BiocInstaller")
  biocLite()
}

updateAllCRANPackages <- function() {
  update.packages()
  inst <- packageStatus()$inst
  inst[inst$Status != "ok", c("Package", "Version", "Status")]
  loadedPackages <- (.packages())
}

#########
# 
unistallPackages <- function(packageNames){
  #Option B: remove.packages2() (function from reposTools).
  libs <- unlist(strsplit(packageNames, " "))
  lapply(libs,function(x){remove.packages(x)})
}
# 
reinstallPackages <- function(thePackages){
  unistallPackages(thePackages)
  loadls(thePackages)
}
# 
UpgradeR <- function(){
  #http://cran.r-project.org/bin/windows/base/rw-FAQ.html "What's the best way to upgrade?"
  #TODO copy (say) R\win-library\3.0 to R\win-library\3.1
  update.packages(checkBuilt=TRUE, ask=FALSE)
}
# 
# whereLibs <- function(){
#   #to figure out where your packages are located
#   .libPaths() 
# }
# 
# 
# 
# resetAllVariables <- function(){
#   is.connection <- function(x) inherits(x, "connection")
#   
#   get_connections <- function(envir = parent.frame()) 
#   {
#     Filter(is.connection, mget(ls(envir = envir), envir = envir))
#   }
#   
#   close_all_connections <- function() 
#   {
#     lapply(c(sys.frames(), globalenv(), baseenv()), 
#            function(e) lapply(get_connections(e), close))
#   }
#   
#   close_all_connections()
#   
#   reset_options <- function()
#   {
#     is_win <- .Platform$OS.type == "windows"
#     options(
#       add.smooth            = TRUE,
#       browserNLdisabled     = FALSE,
#       CBoundsCheck          = FALSE,
#       check.bounds          = FALSE,
#       continue              = "+ ",
#       contrasts             = c(
#         unordered = "contr.treatment", 
#         ordered   = "contr.poly"
#       ), 
#       defaultPackages       = c(
#         "datasets", 
#         "utils", 
#         "grDevices", 
#         "graphics", 
#         "stats",
#         "methods"
#       ),  
#       demo.ask              = "default",
#       device                = if(is_win) windows else x11,
#       device.ask.default    = FALSE,
#       digits                = 7,
#       echo                  = TRUE,
#       editor                = "internal",
#       encoding              = "native.enc",
#       example.ask           = "default",
#       expressions           = 5000,
#       help.search.types     = c("vignette", "demo", "help"),    
#       help.try.all.packages = FALSE,
#       help_type             = "text",
#       HTTPUserAgent         = with(
#         R.version, 
#         paste0(
#           "R (", 
#           paste(major, minor, sep = "."), 
#           " ", 
#           platform, 
#           " ", 
#           arch, 
#           " ", 
#           os, 
#           ")"
#         )
#       ),
#       internet.info         = 2,
#       keep.source           = TRUE,
#       keep.source.pkgs      = FALSE,
#       locatorBell           = TRUE,
#       mailer                = "mailto",
#       max.print             = 99999,
#       menu.graphics         = TRUE,
#       na.action             = "na.omit",
#       nwarnings             = 50,
#       OutDec                = ".",
#       pager                 = "internal",
#       papersize             = "a4",
#       pdfviewer             = file.path(R.home("bin"), "open.exe"),
#       pkgType               = if(is_win) "win.binary" else "source",
#       prompt                = "> ",
#       repos                 = c(
#         CRAN      = "@CRAN@", 
#         CRANextra = "http://www.stats.ox.ac.uk/pub/RWin"
#       ), 
#       scipen                = 0,
#       show.coef.Pvalues     = TRUE,
#       show.error.messages   = TRUE,
#       show.signif.stars     = TRUE,
#       str                   = list(
#         strict.width = "no",
#         digits.d     = 3,
#         vec.len      = 4
#       ),
#       str.dendrogram.last   = "`",
#       stringsAsFactors      = TRUE,
#       timeout               = 60,
#       ts.eps                = 1e-05,
#       ts.S.compat           = FALSE,
#       unzip                 = "internal",
#       useFancyQuotes        = TRUE,
#       verbose               = FALSE,
#       warn                  = 0,
#       warning.length        = 1000,
#       width                 = 80,
#       windowsTimeouts       = c(100, 500)
#     )
#     )
# 
# 
# 
# source(file.path(R.home("etc"), "Rprofile.site"))
# 
# cleanHistory()
#   }
# reset_options()
# }
# 
# 

permutation.R

# #Example
# results <- resampling(x=1:totalNumberOfGenes,
#                       sampleSize=NgenesIp,
#                       nboot=Nrp, theta)
resampling <- function (x, sampleSize, nboot, theta) 
{
  #with no replacement replace = FALSE
  resampled <- matrix(t(replicate(nboot, sample(x,
                                                size = sampleSize,
                                                replace = FALSE))),
                      nrow = nboot)
  #TESTS
  #expect TRUE:
  #all(apply(bootsam, 1, function(x){ length(unique(x)) == dim(bootsam)[2]})) 
  #expect FALSE:  
  #all(apply(bootsam, 2, function(x){ length(unique(x)) == dim(bootsam)[1]}))
  thetastar <- apply(resampled, 1, theta)
  return(list(thetastar = thetastar, theRandomData = resampled))
}

PlotPathways.R

#version 2 (1/1/2018)
#loadls("igraph")

###########
#This method plots a graph from a graphnel with the nodes on entrez IDs (like before send to ROntoTools)
# layout.auto(graph, dim=2, ...)
# layout.random(graph, params, dim=2)
# layout.circle(graph, params)
# layout.sphere(graph, params)
# layout.fruchterman.reingold(graph, ..., dim=2, params)
# layout.kamada.kawai(graph, ..., dim=2, params)
# layout.spring(graph, ..., params)
# layout.reingold.tilford(graph, ..., params)
# layout.fruchterman.reingold.grid(graph, ..., params)
# layout.lgl(graph, ..., params)
# layout.graphopt(graph, ..., params=list())
# layout.svd(graph, d=shortest.paths(graph), ...)
# layout.norm(layout, xmin = NULL, xmax = NULL, ymin = NULL, ymax = NULL,
#             zmin = NULL, zmax = NULL)
# layoutP = layout.reingold.tilford 
# layoutP = layout.fruchterman.reingold
# layoutP = layout.circle
# layoutP = layout.reingold.tilford (circular =T)
# layoutP = layout.reingold.tilford (pathway.igraph, circular =T)
# layoutP = layout.reingold.tilford (pathway.igraph, niter=10000)
plotPathwayIgraph <- function(pathway.i, name = "A pathway", entrezOrganism = "org.Hs.eg", 
                              layoutP = layout.auto(pathway.i),
                        shape="rectangle", vertex.color="lightblue", edge.color="black"){
  loadlib(paste(entrezOrganism,".db",sep=''))
  entrez2Sym <- get(paste(entrezOrganism,"SYMBOL",sep=''))
  pathway.igraph <- igraph.from.graphNEL(pathway.i)
  entrezID <- V(pathway.igraph)$name
  entrezID <- sub("hsa:","", entrezID)
  symbols <-lapply(entrezID, function(x){get(x, entrez2Sym)} )
  symbols <- unlist(symbols)
  V(pathway.igraph)$name <- symbols
  E(pathway.igraph)$color = edge.color
  
#   layoutP = layout.reingold.tilford (pathway.igraph, 1,)
#   
#   layoutP = layout.norm(layoutP, xmin = 1, xmax = 600, ymin = 1, ymax = 600)
#   
#   layoutP = adjustLayout(layoutP, mindist = 80)
#   
#   layoutP = layout.norm(layoutP, xmin = 1, xmax = 600, ymin = 1, ymax = 600)
#   
  
  plot( pathway.igraph, 
               vertex.shape=shape,layout=layoutP, 
               vertex.color=vertex.color, 
               vertex.size =(nchar(symbols) * 20),
               vertex.size2 = (nchar(symbols) * 5 ),
               vertex.label.dist=0.1,
        edge.arrow.size=0.5, main=name) 
}

#######
# Example: plotPathway(kegg_pathways$"path:hsa04122", name = " ", vertex.color="cadetblue2")
# select the color from http://www.stat.columbia.edu/~tzheng/files/Rcolor.pdf
plotPathway <- function(pathway.i, name = "A pathway", entrezOrganism = "org.Hs.eg", 
                        shape="box", vertex.color="lightblue", edge.color="black", fontsize = 10, showEdgeData = TRUE){
  loadlib(paste(entrezOrganism,".db",sep=''), bioc = FALSE)
  loadlib("ROntoTools", bioc = FALSE)
  par(mgp=c(1.5,0.5,0), #axis title, axis labels and axis line
      mai=c(0.6,0.6,0.02,0.02) #c(bottom, left, top, right)
  )
  
  
  entrez2Sym <- get(paste(entrezOrganism,"SYMBOL",sep=''))
  entrezID <- nodes(pathway.i)
  entrezID <- sub("hsa:","", entrezID)
  symbols <-lapply(entrezID, function(x){get(x, entrez2Sym)} ) 
  symbols <- unlist(symbols)
  nodes(pathway.i) <- symbols
  edgeLabels <- unlist(edgeWeights(pathway.i))
  edgesFromTo <- names(edgeLabels)
  if(showEdgeData){
    edgeLabels <- unlist(edgeData(pathway.i))
  }
  nA = makeNodeAttrs(pathway.i, fixedSize=FALSE, height = "1", 
                     width = "1", fillcolor = vertex.color, 
                     shape = shape, fontsize = fontsize)
  eAttrs <- list()
  edgesFromTo <- gsub("\\.", "~", edgesFromTo) 
  edgesFromTo -> names(edgeLabels)
  eAttrs$label <- edgeLabels
  att = list(graph = list(rankdir = "LR", rank = ""))
  plot(pathway.i, main = name, attrs = att, nodeAttrs = nA, edgeAttrs = eAttrs)
}


#plotTheNelGraph:plot Using graphNEL
#Params
#gR: the graph
#factorsE: the factors 
plotPathway2Colors <- function(pathway.i, subclass, name = "A pathway", entrezOrganism = "org.Hs.eg", 
                              twocolors = c("lightgreen", "black"), twoshapes = c("box", "box"),
                              twofontcolor = c("black", "white") , fontsize = 10, numberofchar = 6, showEdgeData = TRUE){
  loadlib(paste(entrezOrganism,".db",sep=''))
  entrez2Sym <- get(paste(entrezOrganism,"SYMBOL",sep=''))
  entrezID <- nodes(pathway.i)
  entrezID <- sub("hsa:","", entrezID)
  symbols <-lapply(entrezID, function(x){if(exists(x, entrez2Sym)) get(x, entrez2Sym) else x} )
  symbols <- substr(symbols, 1, numberofchar)
  symbols <- unlist(symbols)
  nodes(pathway.i) <- symbols
  edgeLabels <- unlist(edgeWeights(pathway.i))
  edgesFromTo <- names(edgeLabels)
  if(showEdgeData){
    edgeLabels <- unlist(edgeData(pathway.i))
  }
  nodeType <- 1 + (nodes(pathway.i) %in% subclass)
  nA = makeNodeAttrs(pathway.i,fixedSize=FALSE, 
                     height = "1", width = "1",
                     fillcolor = twocolors[nodeType], shape = twoshapes[nodeType], 
                     fontcolor = twofontcolor[nodeType], fontsize = fontsize )
  eAttrs <- list()
  edgesFromTo <- gsub("\\.", "~", edgesFromTo) 
  edgesFromTo -> names(edgeLabels)
  eAttrs$label <- edgeLabels
  att = list(graph = list(rankdir = "LR", rank = ""))
  plot(pathway.i, main = name, attrs = att, nodeAttrs = nA, edgeAttrs = eAttrs)
  legend("bottomright", legend = c("Original", "New"), 
         col = twocolors,
         pch = c(19,19), title = "Box color",
         #text.col = c("black","lightblue"), # terrain.colors(10)[1:2]
         #lwd = 2, 
         #cex = 0.5,
         text.font = 3
  )
}

plotPathway2ColorsIGrpah <- function(pathway.i, subclass, name = "A pathway", entrezOrganism = "org.Hs.eg", layoutP = layout.auto, shape="rectangle", 
                               vertex.color="lightblue", vertex.color.subClass="red", edge.color="black"){
  loadlib(paste(entrezOrganism,".db",sep=''))
  entrez2Sym <- get(paste(entrezOrganism,"SYMBOL",sep=''))
  pathway.igraph <- igraph.from.graphNEL(pathway.i)
  entrezID <- V(pathway.igraph)$name
  entrezID <- sub("hsa:","", entrezID)
  symbols <-lapply(entrezID, function(x){if(exists(x, entrez2Sym)) get(x, entrez2Sym) else x} )
  symbols <- unlist(symbols)
  V(pathway.igraph)$name <- symbols
  vertexColors <- sapply(symbols %in% subclass, function(x) {if(x) vertex.color else vertex.color.subClass})
  E(pathway.igraph)$color = edge.color
  plot( pathway.igraph, vertex.shape=shape,layout=layoutP, 
        vertex.color=vertexColors, 
        vertex.size = (nchar(symbols) * 20),
        vertex.size2 = (nchar(symbols) * 5 ),
        vertex.label.dist=0.1, edge.arrow.size=0.5, main=name)   
}

#tkplot
#http://stackoverflow.com/questions/17627052/change-background-color-of-tkplot-igraph-r
tkplotplotPathway <- function(pathway.igraph){
  loadlib(paste(entrezOrganism,".db",sep=''))
  entrez2Sym <- get(paste(entrezOrganism,"SYMBOL",sep=''))
  pathway.igraph <- igraph.from.graphNEL(pathway.i)
  entrezID <- V(pathway.igraph)$name
  entrezID <- sub("hsa:","", entrezID)
  symbols <-lapply(entrezID, function(x){get(x, entrez2Sym)} )
  symbols <- unlist(symbols)
  V(pathway.igraph)$name <- symbols
  E(pathway.igraph)$color = edge.color
  
  #   layoutP = layout.reingold.tilford (pathway.igraph, 1,)
  #   
  #   layoutP = layout.norm(layoutP, xmin = 1, xmax = 600, ymin = 1, ymax = 600)
  #   
  #   layoutP = adjustLayout(layoutP, mindist = 80)
  #   
  #   layoutP = layout.norm(layoutP, xmin = 1, xmax = 600, ymin = 1, ymax = 600)
  #   
  
  temporalT <- tkplot( pathway.igraph, 
        vertex.shape=shape,layout=layoutP, 
        vertex.color=vertex.color, 
        vertex.size =(nchar(symbols) * 20),
        vertex.size2 = (nchar(symbols) * 5 ),
        vertex.label.dist=0.1,
        edge.arrow.size=0.5, main=name) 
  

  tkconfigure(igraph:::.tkplot.get(temporalT)$canvas, "bg"="white")
  
  
  
  
  tkplot( pathway.igraph, 
               vertex.label = V(pathway.igraph)$name,
               vertex.label.color = vertex.color,
               vertex.frame.color = vertex.color,
               #vertex.shape=shape, 
               vertex.size =(nchar(symbols) * 5),
               #vertex.size2 = (nchar(symbols) * 5 ),
               vertex.label.dist=0.1, 
               edge.arrow.size=0.5, 
               edge.width = 0.2,
               main=name, 
               edge.curved = T,
               layout=layout.reingold.tilford) 
}


rglplot <- function(g){
  rglplot(g)
}



###############
#GRAPHNEL
#More  Examples http://www.bioconductor.org/packages/release/bioc/vignettes/graph/inst/doc/graphAttributes.R
simpleTestGRAPHNEL <- function(){
  library(graph)
  set.seed(123)
  V <- LETTERS[1:4]
  edL <- vector("list", length=4)
  names(edL) <- V
  for(i in 1:4)
    edL[[i]] <- list(edges=5-i, weights=runif(1))
  gR <- new("graphNEL", nodes=V, edgeL=edL)
  edges(gR)
  edgeWeights(gR)
  head(edL)
  
  edL[[1]] <- list(edges=2)
  edL$'C' <- list(edges='D')
  gR <- new("graphNEL", nodes=V, edgeL=edL, edgemode="directed")
  plot(gR)
  
  
  V <- c(V,"z")
  V
  plot(gR)
  gX <- addEdge("A", "C", gR, 1)
  
  plot(gX)
  
  
  #count adjacent nodes
  
  nodes(gR)
  adj(gR, "A")
  
  
  set.seed(123)
  g1 = randomEGraph(LETTERS[1:15], edges=100)
  g1
  plot(g1)
  
}

#Rgraphviz
#other graphs
#http://students.washington.edu/mclarkso/documents/gplot%20Ver2.pdf
#labels
#http://mae.ucdavis.edu/dsouza/Lectures/Rgraphviz.pdf
simpleTestRgraphviz <- function(){
  library("Rgraphviz")
  
  defAttrs <- getDefaultAttrs()
  set.seed(123)
  V <- letters[1:10]
  M <- 1:4
  g1 <- randomGraph(V, M, 0.2)
  plot(g1, attrs = list(node = list(label = "foo", fillcolor = "lightgreen"),
                        edge = list(color = "cyan"), graph = list(rankdir = "LR")))
  
  nAttrs <- list()
  eAttrs <- list()
  
  z <- strsplit(packageDescription("Rgraphviz")$Description, " ")[[1]]
  z <- z[1:numNodes(g1)]
  names(z) = nodes(g1)
  nAttrs$label <- z
  eAttrs$label <- c("a~h" = "Label 1", "c~h" = "Label 2")
  attrs <- list(node = list(shape = "ellipse", fixedsize = FALSE))
  plot(g1, nodeAttrs = nAttrs, edgeAttrs = eAttrs, attrs = attrs)
  
  eAttrs <- list()
  ew <- as.character(unlist(edgeWeights(g1)))
  ew <- ew[setdiff(seq(along = ew), removedEdges(g1))]
  names(ew) <- edgeNames(g1)
  eAttrs$label <- ew
  plot(g1, nodeAttrs = nAttrs, edgeAttrs = eAttrs, attrs = attrs)
  
  
  
  nAttrs$color <- c(a = "red", b = "red", g = "green", d = "blue")
  eAttrs$color <- c("a~d" = "blue", "c~h" = "purple")
  nAttrs$fillcolor <- c(j = "yellow")
  nAttrs$fontcolor <- c(e = "green", f = "red")
  eAttrs$fontcolor <- c("a~h" = "green", "c~h" = "brown")
  nAttrs
  
  plot(g1, nodeAttrs = nAttrs, attrs = attrs)
  
  attrs$node$shape <- "ellipse"
  nAttrs$shape <- c(g = "box", f = "circle", j = "box", a = "plaintext")
  
  
  nodes <- buildNodeList(g1)
  edges <- buildEdgeList(g1)
  
  nodes[[1]]
  nodes <- buildNodeList(g1, nodeAttrs = nAttrs, defAttrs = defAttrs$node)
  edges <- buildEdgeList(g1, edgeAttrs = eAttrs, defAttrs = defAttrs$edge)
  nodes[[1]]
  for (j in c("a~e", "a~h")) edges[[j]]@attrs$arrowhead <- "open"
  
  vv <- agopen(name = "foo", nodes = nodes, edges = edges, attrs = attrs,
               edgeMode = "undirected")
  plot(vv)
}

plotsExamples <-function(){
  
  data(graphExamples)
  z <- graphExamples[[8]]
  nNodes <- length(nodes(z))
  nA <- list()
  nA$fixedSize <- rep(FALSE, nNodes)
  nA$height <- nA$width <- rep("1", nNodes)
  nA$label <- rep("Z", nNodes)
  nA$color <- rep("green", nNodes)
  nA$fillcolor <- rep("orange", nNodes)
  nA$shape <- rep("box", nNodes)
  nA$fontcolor <- rep("blue", nNodes)
  nA$fontsize <- rep(40, nNodes)
  nA <- lapply(nA, function(x) {
    names(x) <- nodes(z)
    x
  })
  plot(z, nodeAttrs = nA)
}


bipartitegraphs<-function(){
  set.seed(123)
  nodes1 <- paste(0:7)
  nodes2 <- letters[1:10]
  ft <- cbind(sample(nodes1, 24, replace = TRUE), sample(nodes2,24, replace = TRUE))
  ft <- ft[!duplicated(apply(ft, 1, paste, collapse = "")), ]
  g <- ftM2graphNEL(ft, edgemode = "directed")
  g
  
  
  twocolors <- c("#D9EF8B", "#E0F3F8")
  nodeType <- 1 + (nodes(g) %in% nodes1)
  nA = makeNodeAttrs(g, fillcolor = twocolors[nodeType])
  sg1 = subGraph(nodes1, g)
  sgL = list(list(graph = sg1, cluster = FALSE, attrs = c(rank = "sink")))
  att = list(graph = list(rankdir = "LR", rank = ""))
  
  
  plot(g, attrs = att, nodeAttrs = nA, subGList = sgL)
  
}

#http://svitsrv25.epfl.ch/R-doc/library/graph/html/inEdges.html
getedgesincoming<-function(){
  V <- LETTERS[1:4]
  edL3 <- vector("list", length=4)
  for(i in 1:4)
    edL3[[i]] <- list(edges=(i%%4)+1, weights=i)
  names(edL3) <- V
  gR3 <- new("graphNEL", nodes=V, edgeL=edL3, "directed")
  inEdges(c("A", "B"), gR3)
  
  plot(gR3)
  inEdges("A", gR3)
}


#FOR PAINTING NICE GRAPHS USE IGRAPH
#http://igraph.sourceforge.net/doc/R/graphNEL.html
simpleTestGRAPHNEL <- function(){
  
  library (igraph)
  
  
  plotPathway <- function(pathway.i){
    pathway.i <- igraph.from.graphNEL(pathway.i)
    shapes1 <- rep("square", length(E(pathway.i)))
    shapes2 <- rep("cicle", length(E(pathway.i)))
    theColors <- c("red", "black", "green")
    twoshapes <- c("circle", "square")
    nodeType <- 1 + (nodes(gR) %in% factorsE[,"factor"])
    edgeType <- 2 + (E(pathway.i)$weight )
    shapes = sapply(nodeType,function(x) twoshapes[x])
    colors = sapply(nodeType,function(x) theColors[x])
    E(pathway.i)$color = sapply(edgeType,function(x) theColors[x])
    E(pathway.i)$weight = sapply(E(pathway.i)$weight,function(x) {if(x==-1 |  x==1) 2 else 1 })
    if(tkplotV){
      temporalT <- tkplot(pathway.i, vertex.shape=shapes,layout=layoutP,
                          vertex.color=colors, vertex.size=4, edge.width=E(pathway.i)$weight, 
                          edge.arrow.size=0.5)
      tkconfigure(igraph:::.tkplot.get(temporalT)$canvas, "bg"="white")
    }
    else{
      plot( pathway.i, vertex.shape=shapes,layout=layoutP, 
            vertex.color=colors, vertex.size=4, edge.width=E(pathway.i)$weight, 
            vertex.label.dist=0.3, edge.arrow.size=0.5) #edge.label=E(pathway.i)$weight,
    }
  }
  
  
  #SOME PAINTS
  #http://igraph.sourceforge.net/doc/R/plot.common.html
  g <- graph.ring(10)
  
  # plotting a random graph, set the parameters in the command arguments
  g <- barabasi.game(100)
  plot(g, layout=layout.fruchterman.reingold, vertex.size=4,
       vertex.label.dist=0.5, vertex.color="red", edge.arrow.size=0.5)
  
  # plot a random graph, different color for each component
  g <- erdos.renyi.game(100, 1/100)
  comps <- clusters(g)$membership
  colbar <- rainbow(max(comps)+1)
  V(g)$color <- colbar[comps+1]
  plot(g, layout=layout.fruchterman.reingold, vertex.size=5, vertex.label=NA)
  
  # plot communities in a graph
  g <- graph.full(5) %du% graph.full(5) %du% graph.full(5)
  g <- add.edges(g, c(0,5, 0,10, 5,10))
  com <- spinglass.community(g, spins=5)
  V(g)$color <- com$membership+1
  g <- set.graph.attribute(g, "layout", layout.kamada.kawai(g))
  plot(g, vertex.label.dist=1.5)
  
  # draw a bunch of trees, fix layout
  igraph.options(plot.layout=layout.reingold.tilford)
  plot(graph.tree(20, 2))
  plot(graph.tree(50, 3), vertex.size=3, vertex.label=NA)
  tkplot(graph.tree(50, 2, mode="undirected"), vertex.size=10,
         vertex.color="green")
  
  #draw the edges
  #https://sites.google.com/site/daishizuka/toolkits/sna/weighted-edges
  plot.igraph(net,vertex.label=V(net)$name,layout=layout.fruchterman.reingold, edge.color="black",edge.width=E(net)$weight)
  
  #name of the weigth
  #http://cran.r-project.org/web/packages/igraph/igraph.pdf
  plot( g3, layout=layout.circle, edge.label=E(g3)$weight)
  
  
  #bipartite
  g <- graph.ring(10)
  bipartite.mapping(g)
  ## A star is fine, too
  g2 <- graph.star(10)
  bipartite.mapping(g2)
  ## A graph containing a triangle is not fine
  g3 <- graph.ring(10)
  g3 <- add.edges(g3, c(1,3))
  bipartite.mapping(g3)
  plot(g)
  
  
  
  ##vertex shape
  #http://igraph.sourceforge.net/doc/R/igraph.vertex.shapes.html
  
  # all vertex shapes, minus "raster", that might not be available
  #http://igraph.sourceforge.net/screenshots2.html
  shapes <- setdiff(vertex.shapes(), "")
  g <- graph.ring(length(shapes))
  set.seed(42)
  plot(g, vertex.shape=shapes, vertex.label=shapes, vertex.label.dist=1,
       vertex.size=15, vertex.size2=15,
       vertex.pie=lapply(shapes, function(x) if (x=="pie") 2:6 else 0),
       vertex.pie.color=list(heat.colors(5)), edge.color="black",edge.width=2)
  
  # add new vertex shape, plot nothing with no clipping
  add.vertex.shape("nil")
  plot(g, vertex.shape="nil")
  
  
  #R get edge weigth value
  #http://igraph.sourceforge.net/doc/R/attributes.html
  
  g <- graph.ring(10)
  g <- set.graph.attribute(g, "name", "RING")
  # It is the same as
  g$name <- "RING"
  g$name
  
  g <- set.vertex.attribute(g, "color", value=c("red", "green"))
  get.vertex.attribute(g, "color")
  
  g <- set.edge.attribute(g, "weight", value=runif(ecount(g)))
  get.edge.attribute(g, "weight")
  
  E(g)
  
  list.edge.attributes(g)
  
  get.edge.attribute(g, "weight", index=1)
  
  # The following notation is more convenient
  
  g <- graph.star(LETTERS[1:9])
  
  V(g)$color <- c("red", "green")
  V(g)$color
  
  E(g)$weight <- runif(ecount(g))
  E(g)$weight
  
  print(g, g=TRUE, v=TRUE, e=TRUE)
  
  plot(g2)
  
  #neighbors incoming
  neighbors(g2, "B", mode="in")
  E(g2)
  V(g2)
  
  incident(g, 1)
  theName <- V(g2)[1]
  are.connected(g2, theName, "B")
  
  
  
  
  
  are.connected(g, 2, 4)
  get.edges(g, 1)
  
  get.edge(g, 4)
  
  get.vertex(2)
  
  
  
  ei <- get.edge.ids(g, c(1,2, 4,5))
  
  
  
  #
  
  library(igraph)
  camp <- graph.formula(Harry:Steve:Don:Bert - Harry:Steve:Don:Bert,
                        Pam:Brazey:Carol:Pat - Pam:Brazey:Carol:Pat,
                        Holly   - Carol:Pat:Pam:Jennie:Bill,
                        Bill    - Pauline:Michael:Lee:Holly,
                        Pauline - Bill:Jennie:Ann,
                        Jennie  - Holly:Michael:Lee:Ann:Pauline,
                        Michael - Bill:Jennie:Ann:Lee:John,
                        Ann     - Michael:Jennie:Pauline,
                        Lee     - Michael:Bill:Jennie,
                        Gery    - Pat:Steve:Russ:John,
                        Russ    - Steve:Bert:Gery:John,
                        John    - Gery:Russ:Michael)
  
  V(camp)$label <- V(camp)$name
  set.seed(42)   ## to make this reproducable
  co <- layout.auto(camp)
  
  plot(0, type="n", ann=FALSE, axes=FALSE, xlim=extendrange(co[,1]), 
       ylim=extendrange(co[,2]))
  plot(camp, layout=co, rescale=FALSE, add=TRUE,
       vertex.shape="rectangle",
       vertex.size=(strwidth(V(camp)$label) + strwidth("oo")) * 100,
       vertex.size2=strheight("I") * 2 * 100)
  
  
  #nice plots
  #http://rulesofreason.wordpress.com/tag/tkplot/
  #http://cran.r-project.org/web/packages/igraph/igraph.pdf
  
  
  
  library(igraph)
  L3 <- LETTERS[1:8]
  d <- data.frame(start = sample(L3, 16, replace = T), end = sample(L3, 16, replace = T),
                  weight = c(20,40,20,30,50,60,20,30,20,40,20,30,50,60,20,30))
  
  
  g <- graph.data.frame(d, directed = T)
  
  V(g)$name 
  E(g)$weight
  
  ideg <- degree(g, mode = "in", loops = F)
  
  col=rainbow(12) # For edge colors
  
  plot.igraph(g, 
              vertex.label = V(g)$name, vertex.label.color = "gray20",
              vertex.size = ideg*25 + 40, vertex.size2 = 30,
              vertex.color = "gray90", vertex.frame.color = "gray20",
              vertex.shape = "rectangle",
              edge.arrow.size=0.5, edge.color=col, edge.width = E(g)$weight / 10,
              edge.curved = T, 
              layout = layout.reingold.tilford)
  
  
  
  ###############################
  #http://stackoverflow.com/questions/12088473/spacing-vertices-evenly-in-igraph-in-r
  
  set.seed(1410)
  df<-data.frame(
    "site.x"=c(rep("a",3),rep("b",3),rep("c",3),rep("d",3)),
    "site.y"=c(rep(c("e","f","g"),4)),
    "bond.strength"=sample(1:100,12, replace=TRUE))
  
  library(igraph)
  df<-graph.data.frame(df)
  V(df)$names <- c("a","b","c","d","e","f","g")
  
  layOUT<-data.frame(x=c(rep(1,4),rep(2,3)),y=c(4:1,3:1))
  layOUT[1,2] <- 3
  
  E(df)[ bond.strength < 101 ]$color <- "red"
  E(df)[ bond.strength < 67 ]$color <- "yellow"
  E(df)[ bond.strength < 34 ]$color <- "green"
  V(df)$color <- "white"
  l<-as.matrix(layOUT)
  plot(df,layout=l,vertex.size=10,vertex.label=V(df)$names,
       edge.arrow.size=0.01,vertex.label.color = "black")
  
  
  
  ################################
  #Other
  #http://stackoverflow.com/questions/7521381/draw-network-in-r-control-edge-thickness-plus-non-overlapping-edges
  
  Edges <- data.frame(
    from = rep(c("asdsf","iorejoriejo","asfasfas","asdfsdf","fasdfsawwwwww"),each=5),
    to = rep(c("asdsf","iorejoriejo","asfasfas","asdfsdf","fasdfsawwwwww"),times=5),
    thickness = abs(rnorm(25)))
  
  Edges <- subset(Edges,from!=to)
  
  qgraph(Edges,esize=5,gray=TRUE)
}


#this looks nice http://www.vesnam.com/Rblog/page/2/
#http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3702256/

Plots.R

#version 2 (1/1/2018)

myDensity <-function(fileNam, ptitle)
{
  #myDensityDes(fileNam, ptitle)
  myDensityHis(fileNam, ptitle)
}

myDensityHis <-function(fileNam,ptitle)
{
  font_size_times = 2.5 #big plot
  pdf(file = fileNam ,  onefile = TRUE, pagecentre = FALSE, compress = FALSE)
  
  par(mgp=c(2.5,0.8,0),
      mai=c(1.2,1.5,1.0,0.5) ) #c(bottom, left, top, right)
  hist(hx, freq=FALSE,  breaks = 50, main="", col = "darkseagreen1" ,
       lty = 1, #line type solid
       cex.lab=font_size_times, #font size       
       xlab = "Cox p-value",
       cex.sub = font_size_times,
       cex.axis = font_size_times
  )
  abline(v=FinalPCP,col="red", lwd=3 )  
  legend("topright",
         paste("p-value=",FinalPCP),
         #text.font=2, #2 bold
         cex=font_size_times-0.5,
         box.lty=0
  )
  
  title(ptitle, line = 0.5, 
        adj=0.4, #centers the titles (0 left - 1 right)
        cex.main=(font_size_times+0.3))
  
  dev.off()
}


myDensityDes <-function(fileNam, ptitle)
{
  pdf(file = fileNam ,  onefile = TRUE)
  denshx <- density(
    hx, kernel = "gaussian", from = 0, to = 1, bw = "SJ" )
  plot(denshx, main = ptitle)
  abline(v = FinalPCP,col = 4)
  abline(v = crcPthpval,col = 4)
  legend("topright", paste( "p-val=",round(FinalPCP,5),
                            "\n Cox p-val=", round(crcPthpval,5) ))
  dev.off()
}


########
#plot.survfit manual
#https://stat.ethz.ch/R-manual/R-devel/library/survival/html/plot.survfit.html
# @param groups is a factor
# @param mainTitle is a string
# @param survivalData is data.frame with headers "PatientID" "Survival"  "Death" 
#> head(survivalGBM)
#                     PatientID Survival Death
#1 TCGA-02-0001-01C-01T-0179-07      358     1
#2 TCGA-02-0003-01A-01T-0179-07      144     1
#3 TCGA-02-0007-01A-01T-0179-07      705     1
#
# Example 1:
# pdf( file=paste(resultsFolder, "SNF_all_kaplan.pdf", sep="")  ,  onefile=TRUE)  
#      plot.SurvivalK3 (groups, mainTitle, survivalData, "bottomright")
# dev.off()
#
#
# Example 2:
# mainTitle <- "b) Survival curve, SNF, selected genes"
# pdf(
#   file = paste(resultsFolder, "SNF_union_kaplan.pdf", sep = ""),  onefile =
#     TRUE, width = 9, height = 7
# )
# plot.SurvivalK3 (groups, mainTitle, survivalData,"bottomright",
#                  colorsP = c("red","blue","black"))
# dev.off()
#
#
# Example 3:
# labelsGroups <- c(paste("group 1 (",table(groups)[2], ")", sep="" ),
#                   paste("group 2 (",table(groups)[1], ")", sep="" ),
#                   paste("group 3 (",table(groups)[3], ")", sep="" ) )
# colorsLines <- c("black","red","blue")
# colorsLabels <- c("red","black","blue")
#  pdf( file = pdfName,  onefile = TRUE, width = 9, height = 7)
#  plot.SurvivalK3 (groups, mainTitle, survivalData,"bottomright",
#                   labelClu = labelsGroups,
#                   colorsP = colorsLines,
#                   colorsL = colorsLabels,
#                   centerT = centerT )
#  dev.off()
#
#
#TRY THIS: http://www.r-statistics.com/2013/07/creating-good-looking-survival-curves-the-ggsurv-function/
plot.SurvivalK3 <- function(groups, mainTitle, survData,llocation, 
                            colorsP = c("red","black","blue"),
                            colorsL = c("red","black","blue"),
                            labelClu = c("group1","group2", "group3"),
                            centerT = 0.4){
  
  font_size_times = 1.5 
  require("survival")
  par( mgp=c(1.8,0.3,0), #axis title, axis labels and axis line
       mai=c(2.0,2.0,1.0,1.0) ) #c(bottom, left, top, right)
  
  coxFit <- coxph(Surv(time = Survival, event = Death) ~ groups, 
                  data = survData, ties="exact")
  resCox <- mySumm(coxFit)
  mfit <- survfit(Surv(time = Survival, event = Death) ~ groups, data = survData)
  pVal<-format.pval(resCox$sctest[3], digits = 3)
  plot(mfit, col=colorsP,
       lty = 1, #line type solid
       lwd=5, #line tightness
       #main = mainTitle,
       cex.lab=font_size_times, #font size       
       #xaxt="n", #do not plot x axe
       xlab = "Years",
       conf.int=F,
       #cex.main=(font_size_times+0.1),
       cex.sub =font_size_times,
       cex.axis =font_size_times,
       mark=3, #type of sensored  mark http://www.statmethods.net/advgraphs/parameters.html
       cex = 2, #size of mark
       xscale = 365.25, #values in years
       ylab="Survival probability",
       xaxs="i", #survival curve should touch the y-axis
       tck=0.01 #length of tick marks
  )
  title(mainTitle, line = 0.5, 
        adj=centerT, #centers the titles (0 left - 1 right)
        cex.main=(font_size_times+0.3))
  
  
  #   axis(1, at = seq(0, max(survData$Survival), length.out = 15),
  #       labels=round(seq(0, 4.1, length.out = 15), digits=1), las=1) #customize x axe
  legend(llocation, 
         labelClu,   
         col=colorsL, 
         text.col=colorsL,
         text.font=2, #bold
         fill=colorsL, #boxes
         cex=font_size_times,
         title = paste(" Cox p-value = ", pVal), 
         title.col="black")
}


#or use briss
#delete margins on plots
plot.NoMargins <- function(){
  par(mgp=c(1.5,0.5,0), #axis title, axis labels and axis line
      mai=c(0.6,0.6,0.02,0.02) #c(bottom, left, top, right)
  )
}






# Example:
# pdf( file=paste(resultsFolder, "SNF_all_kaplan.pdf", sep="")  ,  onefile=TRUE)  
# plot.SurvivalK3 (groups, mainTitle, survivalData, "bottomright")
# dev.off()
#TRY THIS: http://www.r-statistics.com/2013/07/creating-good-looking-survival-curves-the-ggsurv-function/
bigplot.SurvivalK3 <- function(groups, mainTitle, survData,llocation){
  require("survival")
  par(mgp=c(2.0,0.4,0), #axis title, axis labels and axis line
      mai=c(1.0,1.0,0.02,0.02) #c(bottom, left, top, right)
  )
  
  colorsP <- c("red","black","blue")
  labelClu <- c("group1","group2", "group3")
  coxFit <- coxph(Surv(time = Survival, event = Death) ~ groups, 
                  data = survData, ties="exact")
  resCox <- mySumm(coxFit)
  mfit <- survfit(Surv(time = Survival, event = Death) ~ groups, data = survData)
  pVal<-format.pval(resCox$sctest[3], digits = 4)
  plot(mfit, col=colorsP, lty = 1, #line type solid
       lwd=15, #line tightness
       main = mainTitle,
       cex.lab=2.5, #font size       
       xaxt="n", #do not plot x axe
       xlab = "Years",conf.int=F, 
       ylab="Survival probability"
  )
  axis(1, at = seq(0, max(survData$Survival), length.out = 15),
       labels=round(seq(0, 4.1, length.out = 15), digits=1), las=1) #customize x axe
  legend(llocation, 
         labelClu,   col=colorsP, 
         text.col=colorsP,
         text.font=2, #bold
         fill=colorsP, #boxes
         cex=2,
         title = paste(" Cox p-value = ", pVal), 
         title.col="black")
}


############
#data2plot for example gene expression data
#colorscales Greens, Greys https://plot.ly/r/heatmaps/
#save with preview
potheatmap <- function(data2plot, colorscale="Hot"){
  require("plotly")
  
  f <- list(
    family = "Calibri, monospace",
    size = 30,
    color= "#000000"
  )
  x <- list(
    tickfont = f
  )
  plot_ly(z = data2plot, cl.lim = c(0,1),colorscale=colorscale,
          type = "heatmap")
  layout(xaxis = x, yaxis = x)
}


plotLineSimple <- function (line1, line2, line3){
  loadls("Rcpp ggplot2")
  minY = min(line1, line2, line3)
  maxY = max(line1, line2, line3)
  plot (df1, xlim=c(minX, maxX), ylim=c(minY, maxY))
  plot(line1,type="l", col = "red", ylim=c(minY, maxY))
  lines(line2, col = "blue")
  lines(line3, col = "green")
}


plotLines <- function (line1, line2, line3,lab1, lab2, lab3,xlab, ylab ){
  loadls("Rcpp ggplot2")
  df <- data.frame(x = rep(c(1:length(line1)), 3), y = c(line1, line2, line3), 
                   colors = rep(c(lab1, lab2, lab3), each=length(line1)) )
  plines <- ggplot(data = df, aes(x=x, y=y, col=colors ) )+ geom_line() + xlab(xlab) + ylab(ylab)
  plines + guides(fill=guide_legend(title=NULL))
  plines + theme(legend.title=element_blank())
  plines
}



# Example


# full documentation http://www.inside-r.org/packages/cran/VennDiagram/docs/venn.diagram
plot.venn <- function (){
  require("VennDiagram")
  
  
  # sample four-set Venn Diagram
  A <- sample(1:1000, 400, replace = FALSE);
  B <- sample(1:1000, 600, replace = FALSE);
  C <- sample(1:1000, 350, replace = FALSE);
  D <- sample(1:1000, 550, replace = FALSE);
  E <- sample(1:1000, 375, replace = FALSE);
  
  venn.plot <- venn.diagram(
    x = list(
      A = A,
      D = D,
      B = B,
      C = C
    ),
    filename = "Venn_4set_pretty.tiff",
    col = "transparent",
    fill = c("cornflowerblue", "green", "yellow", "darkorchid1"),
    alpha = 0.50,
#     label.col = c("orange", "white", "darkorchid4", "white", 
#                   "white", "white", "white", "white", "darkblue", "white", 
#                   "white", "white", "white", "darkgreen", "white"),
    cex = 0,
    #fontfamily = "serif",
    #fontface = "bold",
    cat.col = c("darkblue", "darkgreen", "orange", "darkorchid4"),
    cat.cex = 0,
    cat.pos = 0,
    cat.dist = 0.07,
    #cat.fontfamily = "serif",
   # rotation.degree = 270,
    margin = 0.2
  );
  plot(venn.plot)
  
  # sample five-set Venn Diagram
  venn.plot <- venn.diagram(
    x = list(
      A = A,
      B = B,
      C = C,
      D = D,
      E = E
    ),
    filename = "Venn_5set_pretty.tiff",
    col = "black",
    fill = c("dodgerblue", "goldenrod1", "darkorange1", "seagreen3", "orchid3"),
    alpha = 0.50,
    cex = c(1.5, 1.5, 1.5, 1.5, 1.5, 1, 0.8, 1, 0.8, 1, 0.8, 1, 0.8,
            1, 0.8, 1, 0.55, 1, 0.55, 1, 0.55, 1, 0.55, 1, 0.55, 1, 1, 1, 1, 1, 1.5),
    cat.col = c("dodgerblue", "goldenrod1", "darkorange1", "seagreen3", "orchid3"),
    cat.cex = 1.5,
    cat.fontface = "bold",
    margin = 0.05
  );
  
  
  
  plot(venn.plot)
}


#save plot pdf nice
#pdf( file=paste(resultsFolder, "plots2.pdf", sep="")  ,  onefile=TRUE)  
# layout(matrix(c(1,1,2,3),ncol=2,byrow=TRUE),widths=c(1,1),heights=c(3,2))
# par(mar=c(0,0,5,0))
# plot.new(); plot.window(xlim=c(0,1),ylim=c(0,1))
# title("Title",cex.main=1.5)
# text(0.4,0.5,adj=c(0,0),lab=
#        "> print(myList)
# $a
# [1] 1
# 
# $b
# [1] 2")
# par(mar=c(5,4,1,1))
# boxplot(1:10)
# hist(1:10)
# dev.off()

survival_analysis.R

geneSetFormat <- function(genesOfInterest){
  genesOfInterest <- gsub("hsa:","",genesOfInterest)
  genesOfInterest <- lapply(genesOfInterest, find_sym )
  genesOfInterest <- unlist(genesOfInterest[!is.na(genesOfInterest)])
  #Number of mRNA on pathway
  genesOfInterest <- intersect(colnames(geneEData), as.character(genesOfInterest))
  genesOfInterest
}


#get symbols
#x <- org.Hs.egSYMBOL
# Get the gene symbol that are mapped to an entrez gene identifiers
# mapped_genes <- mappedkeys(x)
#example
find_sym <- function(g){
  require(org.Hs.eg.db)
  SYMB <- org.Hs.egSYMBOL[[g]]
  if(is.null(SYMB) ){
    return(NA)
  }else{
    return(SYMB)
  }
}



#function copied from survival:::summary.coxph
mySumm <- function (object, conf.int = 0.95, scale = 1, ...) 
{
  cox <- object
  beta <- cox$coefficients
  if (is.null(cox$coefficients)) {
    return(object)
  }
  nabeta <- !(is.na(beta))
  beta2 <- beta[nabeta]
  if (is.null(beta) | is.null(cox$var)) 
    stop("Input is not valid")
  se <- sqrt(diag(cox$var))
  if (!is.null(cox$naive.var)) 
    nse <- sqrt(diag(cox$naive.var))
  rval <- list(call = cox$call, fail = cox$fail, na.action = cox$na.action, 
               n = cox$n, loglik = cox$loglik)
  if (!is.null(cox$nevent)) 
    rval$nevent <- cox$nevent
  if (is.null(cox$naive.var)) {
    tmp <- cbind(beta, exp(beta), se, beta/se, 1 - pchisq((beta/se)^2, 
                                                          1))
    dimnames(tmp) <- list(names(beta), c("coef", "exp(coef)", 
                                         "se(coef)", "z", "Pr(>|z|)"))
  }
  else {
    tmp <- cbind(beta, exp(beta), nse, se, beta/se, 1 - pchisq((beta/se)^2, 
                                                               1))
    dimnames(tmp) <- list(names(beta), c("coef", "exp(coef)", 
                                         "se(coef)", "robust se", "z", "Pr(>|z|)"))
  }
  rval$coefficients <- tmp
  if (conf.int) {
    z <- qnorm((1 + conf.int)/2, 0, 1)
    beta <- beta * scale
    se <- se * scale
    tmp <- cbind(exp(beta), exp(-beta), exp(beta - z * se), 
                 exp(beta + z * se))
    dimnames(tmp) <- list(names(beta), c("exp(coef)", "exp(-coef)", 
                                         paste("lower .", round(100 * conf.int, 2), sep = ""), 
                                         paste("upper .", round(100 * conf.int, 2), sep = "")))
    rval$conf.int <- tmp
  }
  df <- length(beta2)
  logtest <- -2 * (cox$loglik[1] - cox$loglik[2])
  rval$logtest <- c(test = logtest, df = df, pvalue = 1 - pchisq(logtest, 
                                                                 df))
  rval$sctest <- c(test = cox$score, df = df, pvalue = 1 - 
                     pchisq(cox$score, df))
  rval$rsq <- c(rsq = 1 - exp(-logtest/cox$n), maxrsq = 1 - 
                  exp(2 * cox$loglik[1]/cox$n))
  rval$waldtest <- c(test = as.vector(round(cox$wald.test, 
                                            2)), df = df, pvalue = 1 - pchisq(as.vector(cox$wald.test), 
                                                                              df))
  if (!is.null(cox$rscore)) 
    rval$robscore <- c(test = cox$rscore, df = df, pvalue = 1 - 
                         pchisq(cox$rscore, df))
  rval$used.robust <- !is.null(cox$naive.var)
  if (!is.null(cox$concordance)) {
    if (is.matrix(cox$concordance)) 
      temp <- colSums(cox$concordance)
    else temp <- cox$concordance
    rval$concordance <- c(concordance = (temp[1] + temp[3]/2)/sum(temp[1:3]), 
                          se = temp[5]/(2 * sum(temp[1:3])))
  }
  rval
}