1. Read data

 setwd('C:/Users/Xijin.Ge/Downloads')   # Needs to be changed  
 source('iDEP_core_functions.R') 
## 
##  Installing lots of R packages, this may take several hours ... 
## 
##    Each of these packages took years to develop.
##  So be patient while we steal a ton of code.
## 
##    
##  Note 1: Sometimes dependencies needs to be installed manually. 
##    
##  Note 2: If you are using an older version of R, and having trouble with package installation.
##    
##          Sometimes, it is easier to install the new version of R and remove all old packages, and start fresh. 
## 
 inputFile <- 'Downloaded_Converted_Data.csv'  # Expression matrix
 sampleInfoFile <- 'Downloaded_sampleInfoFile.csv' # Experiment design file 
 geneInfoFile <- 'Mouse__mmusculus_gene_ensembl_GeneInfo.csv' #Gene symbols, location etc. 
 geneSetFile <- 'Mouse__mmusculus_gene_ensembl.db'  # pathway database in SQL; can be GMT format 
 STRING10_speciesFile <- 'https://raw.githubusercontent.com/iDEP-SDSU/idep/master/shinyapps/idep/STRING10_species.csv'
 input_missingValue <- 'geneMedian' #Missing values imputation method
 input_dataFileFormat <- 1  #1- read counts, 2 FKPM/RPKM or DNA microarray
 input_minCounts <- 0.5 #Min counts
 input_NminSamples <- 1 #Minimum number of samples 
 input_countsLogStart <- 4  #Pseudo count for log CPM
 input_CountsTransform <- 1 #Methods for data transformation of counts. 1-VST, 2-rlog, 3-EdgeR 
 readData.out <- readData(inputFile) 
 head(readData.out$data)  
##                    p53_mock_1 p53_mock_2 p53_mock_3 p53_mock_4 p53_IR_1
## ENSMUSG00000075014   5.986190  10.245204  11.622005  12.445437 19.58986
## ENSMUSG00000076617  18.174331  18.020764  17.901045  16.892389 18.06544
## ENSMUSG00000015656  15.349228  15.086274  15.392536  14.817324 15.95674
## ENSMUSG00000052305  13.766147  11.436670  10.523028  15.927563 10.57175
## ENSMUSG00000069917  13.511114  11.202619  10.240337  15.915347 10.33765
## ENSMUSG00000075015   4.023733   7.791965   9.161421   9.809036 16.98655
##                    p53_IR_2  p53_IR_3 p53_IR_4 null_mock_1 null_mock_2
## ENSMUSG00000075014 21.07925  8.279954 17.40506    13.23741    12.91427
## ENSMUSG00000076617 18.52019 16.933927 17.68636    17.18909    17.39665
## ENSMUSG00000015656 15.75896 16.233039 16.04712    15.03582    14.98886
## ENSMUSG00000052305 11.22625 10.125088 10.20851    11.81668    17.59129
## ENSMUSG00000069917 11.20676 10.280270 10.32910    11.51218    17.41403
## ENSMUSG00000075015 18.47457  5.767006 14.96939    10.59269    10.25800
##                    null_IR_1 null_IR_2
## ENSMUSG00000075014  6.949322  14.57512
## ENSMUSG00000076617 17.663021  17.21374
## ENSMUSG00000015656 18.603855  17.00075
## ENSMUSG00000052305 10.311098  15.82093
## ENSMUSG00000069917 10.236792  15.73450
## ENSMUSG00000075015  4.742948  12.00540
 readSampleInfo.out <- readSampleInfo(sampleInfoFile) 
 readSampleInfo.out 
##             p53    Treatment
## p53_mock_1  "wt"   "mock"   
## p53_mock_2  "wt"   "mock"   
## p53_mock_3  "wt"   "mock"   
## p53_mock_4  "wt"   "mock"   
## p53_IR_1    "wt"   "IR"     
## p53_IR_2    "wt"   "IR"     
## p53_IR_3    "wt"   "IR"     
## p53_IR_4    "wt"   "IR"     
## null_mock_1 "null" "mock"   
## null_mock_2 "null" "mock"   
## null_IR_1   "null" "IR"     
## null_IR_2   "null" "IR"
 input_selectOrg ="NEW" 
 input_selectGO <- 'GOBP'   #Gene set category 
 input_noIDConversion = TRUE  
 allGeneInfo.out <- geneInfo(geneInfoFile) 
 converted.out = NULL 
 convertedData.out <- convertedData()    
 nGenesFilter()  
## [1] "6882 genes in 12 samples. 6882  genes passed filter (see above). Original gene IDs used."
 convertedCounts.out <- convertedCounts()  # converted counts 

2. Pre-process

 parDefault = par() 
 par(mar=c(12,4,2,2)) 
 # barplot of total read counts 
 barplot( colSums(readData.out$rawCounts)/1e6, 
        col="green",las=3, main="Total read counts (millions)")  

 ### Box plot 
 x = readData.out$data 
 boxplot(x, las = 2, 
    ylab='Transformed expression levels',
    main='Distribution of transformed data') 

 #### Density plot 
 par(parDefault) 
## Warning in par(parDefault): graphical parameter "cin" cannot be set
## Warning in par(parDefault): graphical parameter "cra" cannot be set
## Warning in par(parDefault): graphical parameter "csi" cannot be set
## Warning in par(parDefault): graphical parameter "cxy" cannot be set
## Warning in par(parDefault): graphical parameter "din" cannot be set
## Warning in par(parDefault): graphical parameter "page" cannot be set
 densityPlot()       

 ### Scatter plot of the first two samples 
 plot(x[,1:2],xlab=colnames(x)[1],ylab=colnames(x)[2], 
    main='Scatter plot of first two samples') 

 ####plot gene or gene family
 input_selectOrg ="BestMatch" 
 input_geneSearch <- 'HOXA' #Gene ID for searching 
 genePlot()  

 input_useSD <- 'FALSE' #Use standard deviation instead of standard error in error bar? 
 geneBarPlotError()       

## 3. Heatmap

 # hierarchical clustering tree
 x <- readData.out$data
 maxGene <- apply(x,1,max)
 # remove bottom 25% lowly expressed genes, which inflate the PPC
 x <- x[which(maxGene > quantile(maxGene)[1] ) ,] 
 plot(as.dendrogram(hclust2( dist2(t(x)))), ylab="1 - Pearson C.C.", type = "rectangle") 

 #Correlation matrix
 input_labelPCC <- TRUE #Show correlation coefficient? 
 correlationMatrix() 

 # Parameters for heatmap
 input_nGenes <- 1000   #Top genes for heatmap
 input_geneCentering <- TRUE    #centering genes ?
 input_sampleCentering <- FALSE #Center by sample?
 input_geneNormalize <- FALSE   #Normalize by gene?
 input_sampleNormalize <- FALSE #Normalize by sample?
 input_noSampleClustering <- FALSE  #Use original sample order
 input_heatmapCutoff <- 4   #Remove outliers beyond number of SDs 
 input_distFunctions <- 1   #which distant funciton to use
 input_hclustFunctions <- 1 #Linkage type
 input_heatColors1 <- 1 #Colors
 input_selectFactorsHeatmap <- 'Sample_Name'    #Sample coloring factors 
 png('heatmap.png', width = 10, height = 15, units = 'in', res = 300) 
 staticHeatmap() 
 dev.off()  
## png 
##   2

[heatmap] (heatmap.png)

 heatmapPlotly() # interactive heatmap using Plotly 

## 4. K-means clustering

 input_nGenesKNN <- 2000    #Number of genes fro k-Means
 input_nClusters <- 4   #Number of clusters 
 maxGeneClustering = 12000
 input_kmeansNormalization <- 'geneMean'    #Normalization
 input_KmeansReRun <- 0 #Random seed 

 distributionSD()  #Distribution of standard deviations 

 KmeansNclusters()  #Number of clusters 

 Kmeans.out = Kmeans()   #Running K-means 
 KmeansHeatmap()   #Heatmap for k-Means 

 #Read gene sets for enrichment analysis 
 sqlite  <- dbDriver('SQLite')
 input_selectGO3 <- 'GOBP'  #Gene set category
 input_minSetSize <- 15 #Min gene set size
 input_maxSetSize <- 2000   #Max gene set size 
 GeneSets.out <-readGeneSets( geneSetFile,
    convertedData.out, input_selectGO3,input_selectOrg,
    c(input_minSetSize, input_maxSetSize)  )  
 # Alternatively, users can use their own GMT files by
 #GeneSets.out <- readGMTRobust('somefile.GMT')  
 KmeansGO()  #Enrichment analysis for k-Means clusters
##                                                                 Cluster
## RNA processing                                                        A
## Ribonucleoprotein complex biogenesis                                   
## NcRNA metabolic process                                                
## Ribosome biogenesis                                                    
## NcRNA processing                                                       
## Organonitrogen compound metabolic process                              
## Organonitrogen compound biosynthetic process                           
## RRNA processing                                                        
## RRNA metabolic process                                                 
## Translation                                                            
## Cellular response to stress                                           B
## Apoptotic process                                                      
## Cell death                                                             
## Programmed cell death                                                  
## Apoptotic signaling pathway                                            
## Positive regulation of cell death                                      
## Cellular response to DNA damage stimulus                               
## Intrinsic apoptotic signaling pathway                                  
## Positive regulation of programmed cell death                           
## Intrinsic apoptotic signaling pathway in response to DNA damage        
## Immune system process                                                 C
## Positive regulation of response to stimulus                            
## Immune response                                                        
## Establishment of protein localization                                  
## Macromolecular complex assembly                                        
## Regulation of immune system process                                    
## Protein complex assembly                                               
## Protein transport                                                      
## Protein complex biogenesis                                             
## Protein complex subunit organization                                   
## Immune system process1                                                D
## Positive regulation of response to stimulus1                           
## Positive regulation of molecular function                              
## Protein phosphorylation                                                
## Defense response                                                       
## Regulation of intracellular signal transduction                        
## Cell activation                                                        
## Regulation of transport                                                
## Response to external stimulus                                          
## Vesicle-mediated transport                                             
##                                                                      adj.Pval
## RNA processing                                                  2.685004e-144
## Ribonucleoprotein complex biogenesis                            1.992995e-140
## NcRNA metabolic process                                         2.289446e-134
## Ribosome biogenesis                                             3.496906e-117
## NcRNA processing                                                6.343278e-115
## Organonitrogen compound metabolic process                       3.482536e-107
## Organonitrogen compound biosynthetic process                    1.711515e-101
## RRNA processing                                                  6.861087e-93
## RRNA metabolic process                                           2.888202e-92
## Translation                                                      5.045222e-75
## Cellular response to stress                                      2.075602e-25
## Apoptotic process                                                4.559668e-23
## Cell death                                                       4.559668e-23
## Programmed cell death                                            5.845102e-23
## Apoptotic signaling pathway                                      5.845102e-23
## Positive regulation of cell death                                1.232970e-20
## Cellular response to DNA damage stimulus                         1.439213e-20
## Intrinsic apoptotic signaling pathway                            3.577557e-19
## Positive regulation of programmed cell death                     1.011040e-18
## Intrinsic apoptotic signaling pathway in response to DNA damage  1.179594e-18
## Immune system process                                            2.394955e-63
## Positive regulation of response to stimulus                      8.229461e-45
## Immune response                                                  1.180337e-44
## Establishment of protein localization                            1.514473e-42
## Macromolecular complex assembly                                  7.228082e-42
## Regulation of immune system process                              2.287631e-40
## Protein complex assembly                                         2.287631e-40
## Protein transport                                                2.287631e-40
## Protein complex biogenesis                                       2.287631e-40
## Protein complex subunit organization                             2.287631e-40
## Immune system process1                                           9.615041e-94
## Positive regulation of response to stimulus1                     6.158140e-55
## Positive regulation of molecular function                        5.378810e-52
## Protein phosphorylation                                          6.500699e-52
## Defense response                                                 1.857527e-51
## Regulation of intracellular signal transduction                  7.246982e-51
## Cell activation                                                  1.023829e-49
## Regulation of transport                                          2.242061e-47
## Response to external stimulus                                    1.004958e-46
## Vesicle-mediated transport                                       5.131263e-46
##                                                                 Genes
## RNA processing                                                    121
## Ribonucleoprotein complex biogenesis                              106
## NcRNA metabolic process                                           103
## Ribosome biogenesis                                                84
## NcRNA processing                                                   84
## Organonitrogen compound metabolic process                         144
## Organonitrogen compound biosynthetic process                      119
## RRNA processing                                                    65
## RRNA metabolic process                                             65
## Translation                                                        74
## Cellular response to stress                                        36
## Apoptotic process                                                  36
## Cell death                                                         37
## Programmed cell death                                              36
## Apoptotic signaling pathway                                        25
## Positive regulation of cell death                                  24
## Cellular response to DNA damage stimulus                           23
## Intrinsic apoptotic signaling pathway                              18
## Positive regulation of programmed cell death                       22
## Intrinsic apoptotic signaling pathway in response to DNA damage    14
## Immune system process                                             115
## Positive regulation of response to stimulus                        86
## Immune response                                                    76
## Establishment of protein localization                              81
## Macromolecular complex assembly                                    76
## Regulation of immune system process                                69
## Protein complex assembly                                           69
## Protein transport                                                  75
## Protein complex biogenesis                                         69
## Protein complex subunit organization                               72
## Immune system process1                                            147
## Positive regulation of response to stimulus1                       99
## Positive regulation of molecular function                          93
## Protein phosphorylation                                            93
## Defense response                                                   84
## Regulation of intracellular signal transduction                    86
## Cell activation                                                    72
## Regulation of transport                                            88
## Response to external stimulus                                      95
## Vesicle-mediated transport                                         73
##                                                                                                                        Pathways
## RNA processing                                                                                                   RNA processing
## Ribonucleoprotein complex biogenesis                                                       Ribonucleoprotein complex biogenesis
## NcRNA metabolic process                                                                                 NcRNA metabolic process
## Ribosome biogenesis                                                                                         Ribosome biogenesis
## NcRNA processing                                                                                               NcRNA processing
## Organonitrogen compound metabolic process                                             Organonitrogen compound metabolic process
## Organonitrogen compound biosynthetic process                                       Organonitrogen compound biosynthetic process
## RRNA processing                                                                                                 RRNA processing
## RRNA metabolic process                                                                                   RRNA metabolic process
## Translation                                                                                                         Translation
## Cellular response to stress                                                                         Cellular response to stress
## Apoptotic process                                                                                             Apoptotic process
## Cell death                                                                                                           Cell death
## Programmed cell death                                                                                     Programmed cell death
## Apoptotic signaling pathway                                                                         Apoptotic signaling pathway
## Positive regulation of cell death                                                             Positive regulation of cell death
## Cellular response to DNA damage stimulus                                               Cellular response to DNA damage stimulus
## Intrinsic apoptotic signaling pathway                                                     Intrinsic apoptotic signaling pathway
## Positive regulation of programmed cell death                                       Positive regulation of programmed cell death
## Intrinsic apoptotic signaling pathway in response to DNA damage Intrinsic apoptotic signaling pathway in response to DNA damage
## Immune system process                                                                                     Immune system process
## Positive regulation of response to stimulus                                         Positive regulation of response to stimulus
## Immune response                                                                                                 Immune response
## Establishment of protein localization                                                     Establishment of protein localization
## Macromolecular complex assembly                                                                 Macromolecular complex assembly
## Regulation of immune system process                                                         Regulation of immune system process
## Protein complex assembly                                                                               Protein complex assembly
## Protein transport                                                                                             Protein transport
## Protein complex biogenesis                                                                           Protein complex biogenesis
## Protein complex subunit organization                                                       Protein complex subunit organization
## Immune system process1                                                                                    Immune system process
## Positive regulation of response to stimulus1                                        Positive regulation of response to stimulus
## Positive regulation of molecular function                                             Positive regulation of molecular function
## Protein phosphorylation                                                                                 Protein phosphorylation
## Defense response                                                                                               Defense response
## Regulation of intracellular signal transduction                                 Regulation of intracellular signal transduction
## Cell activation                                                                                                 Cell activation
## Regulation of transport                                                                                 Regulation of transport
## Response to external stimulus                                                                     Response to external stimulus
## Vesicle-mediated transport                                                                           Vesicle-mediated transport
 input_seedTSNE <- 0    #Random seed for t-SNE
 input_colorGenes <- TRUE   #Color genes in t-SNE plot? 
 tSNEgenePlot()  #Plot genes using t-SNE 

## 5. PCA and beyond

 input_selectFactors <- 'Sample_Name'   #Factor coded by color
 input_selectFactors2 <- 'Sample_Name'  #Factor coded by shape
 input_tsneSeed2 <- 0   #Random seed for t-SNE 
 #PCA, MDS and t-SNE plots
 PCAplot()  

 MDSplot() 

 tSNEplot()  

 #Read gene sets for pathway analysis using PGSEA on principal components 
 input_selectGO6 <- 'GOBP' 
 GeneSets.out <-readGeneSets( geneSetFile,
    convertedData.out, input_selectGO6,input_selectOrg,
    c(input_minSetSize, input_maxSetSize)  )  
 PCApathway() # Run PGSEA analysis 

 cat( PCA2factor() )   #The correlation between PCs with factors 
## 
##  Correlation between Principal Components (PCs) with factors
## PC1 is correlated with Treatment (p=6.56e-05).

6. DEG1

 input_CountsDEGMethod <- 3 #DESeq2= 3,limma-voom=2,limma-trend=1 
 input_limmaPval <- 0.1 #FDR cutoff
 input_limmaFC <- 2 #Fold-change cutoff
 input_selectModelComprions <- NULL     #Selected comparisons
 input_selectFactorsModel <- NULL   #Selected comparisons
 input_selectInteractions <- NULL   #Selected comparisons
 input_selectBlockFactorsModel <- NULL  #Selected comparisons 
 factorReferenceLevels.out <- NULL 

 limma.out <- limma()
 DEG.data.out <- DEG.data()
 limma.out$comparisons 
## [1] "null_mock-null_IR"  "p53_IR-null_IR"     "p53_mock-null_IR"  
## [4] "p53_IR-null_mock"   "p53_mock-null_mock" "p53_mock-p53_IR"
 input_selectComparisonsVenn = limma.out$comparisons[1:3] # use first three comparisons
 input_UpDownRegulated <- FALSE #Split up and down regulated genes 
 vennPlot() # Venn diagram 

  sigGeneStats() # number of DEGs as figure 

  sigGeneStatsTable() # number of DEGs as table 
##                           Comparisons   Up Down
## null_mock-null_IR   null_mock-null_IR  702 1532
## p53_IR-null_IR         p53_IR-null_IR 1118  577
## p53_mock-null_IR     p53_mock-null_IR  910 1687
## p53_IR-null_mock     p53_IR-null_mock 1741  540
## p53_mock-null_mock p53_mock-null_mock   20   16
## p53_mock-p53_IR       p53_mock-p53_IR  487 1856

7. DEG2

##########################
# 7. DEG2 
##########################
 input_selectContrast <- 'null_mock-null_IR'    #Selected comparisons 
 selectedHeatmap.data.out <- selectedHeatmap.data()
 selectedHeatmap()   # heatmap for DEGs in selected comparison

 # Save gene lists and data into files
 write.csv( selectedHeatmap.data()$genes, 'heatmap.data.csv') 
 write.csv(DEG.data(),'DEG.data.csv' )
 write(AllGeneListsGMT() ,'AllGeneListsGMT.gmt')
 input_selectGO2 <- 'GOBP'  #Gene set category 
 geneListData.out <- geneListData()  
  volcanoPlot()  

  scatterPlot()  

  MAplot()  

  geneListGOTable.out <- geneListGOTable()  
 # Read pathway data again 
 GeneSets.out <-readGeneSets( geneSetFile,
    convertedData.out, input_selectGO2,input_selectOrg,
    c(input_minSetSize, input_maxSetSize)  ) 
 input_removeRedudantSets <- TRUE   #Remove highly redundant gene sets? 
 geneListGO() 
##                                                      Direction adj.Pval
## Organonitrogen compound metabolic process       Down regulated 6.1e-157
## RNA processing                                                 9.4e-155
## NcRNA metabolic process                                        4.1e-151
## Ribonucleoprotein complex biogenesis                           7.9e-143
## Organonitrogen compound biosynthetic process                   3.3e-136
## NcRNA processing                                               2.7e-120
## Ribosome biogenesis                                            6.8e-120
## Small molecule metabolic process                               4.3e-114
## RRNA processing                                                 6.5e-90
## RRNA metabolic process                                          7.6e-89
## Cellular response to organic substance            Up regulated  3.1e-34
## Immune system process                                           3.0e-32
## Response to oxygen-containing compound                          1.7e-31
## Response to endogenous stimulus                                 2.2e-30
## Regulation of transport                                         1.7e-26
## Response to nitrogen compound                                   1.9e-26
## Lipid metabolic process                                         2.8e-26
## Cellular response to oxygen-containing compound                 3.6e-24
## Response to organonitrogen compound                             5.5e-24
## Transmembrane transport                                         5.5e-24
##                                                 nGenes
## Organonitrogen compound metabolic process          251
## RNA processing                                     159
## NcRNA metabolic process                            135
## Ribonucleoprotein complex biogenesis               131
## Organonitrogen compound biosynthetic process       192
## NcRNA processing                                   104
## Ribosome biogenesis                                102
## Small molecule metabolic process                   202
## RRNA processing                                     76
## RRNA metabolic process                              76
## Cellular response to organic substance              86
## Immune system process                               89
## Response to oxygen-containing compound              72
## Response to endogenous stimulus                     70
## Regulation of transport                             69
## Response to nitrogen compound                       52
## Lipid metabolic process                             57
## Cellular response to oxygen-containing compound     51
## Response to organonitrogen compound                 47
## Transmembrane transport                             57
##                                                                                        Pathways
## Organonitrogen compound metabolic process             Organonitrogen compound metabolic process
## RNA processing                                                                   RNA processing
## NcRNA metabolic process                                                 NcRNA metabolic process
## Ribonucleoprotein complex biogenesis                       Ribonucleoprotein complex biogenesis
## Organonitrogen compound biosynthetic process       Organonitrogen compound biosynthetic process
## NcRNA processing                                                               NcRNA processing
## Ribosome biogenesis                                                         Ribosome biogenesis
## Small molecule metabolic process                               Small molecule metabolic process
## RRNA processing                                                                 RRNA processing
## RRNA metabolic process                                                   RRNA metabolic process
## Cellular response to organic substance                   Cellular response to organic substance
## Immune system process                                                     Immune system process
## Response to oxygen-containing compound                   Response to oxygen-containing compound
## Response to endogenous stimulus                                 Response to endogenous stimulus
## Regulation of transport                                                 Regulation of transport
## Response to nitrogen compound                                     Response to nitrogen compound
## Lipid metabolic process                                                 Lipid metabolic process
## Cellular response to oxygen-containing compound Cellular response to oxygen-containing compound
## Response to organonitrogen compound                         Response to organonitrogen compound
## Transmembrane transport                                                 Transmembrane transport
 # STRING-db API access 
 STRING10_species = read.csv(STRING10_speciesFile)  
 ix = grep('Mus musculus', STRING10_species$official_name )
 findTaxonomyID.out <- STRING10_species[ix,1] # find taxonomyID
 findTaxonomyID.out  
## [1] 10090
 # users can also skip the above and assign NCBI taxonomy id directly by
 # findTaxonomyID.out = 10090 # mouse 10090, human 9606 etc. 
  STRINGdb_geneList.out <- STRINGdb_geneList() #convert gene lists
## Warning:  we couldn't map to STRING 9% of your identifiers
 input_STRINGdbGO <- 'Process'  #'Process', 'Component', 'Function', 'KEGG', 'Pfam', 'InterPro' 
 stringDB_GO_enrichmentData() 
##         Direction adj.Pval nGenes
## 1    Up regulated  7.3e-08     49
## 2                  5.1e-07     41
## 3                  8.4e-07     65
## 4                  8.4e-07     46
## 5                  8.4e-07     53
## 6                  9.5e-07     37
## 7                  1.0e-06     49
## 8                  1.3e-06     60
## 9                  2.3e-06     60
## 10                 9.5e-06     40
## 11                 9.5e-06     36
## 12                 9.5e-06     35
## 13                 5.0e-05     26
## 14                 5.4e-05     46
## 15                 6.6e-05     51
## 16 Down regulated  7.5e-62    107
## 17                 2.8e-55    106
## 18                 2.1e-52     81
## 19                 1.6e-51    136
## 20                 4.2e-50    146
## 21                 2.3e-46     85
## 22                 1.3e-43    181
## 23                 9.8e-41     61
## 24                 7.4e-40     61
## 25                 8.8e-32     84
## 26                 1.9e-29     72
## 27                 1.9e-28    130
## 28                 1.1e-27     72
## 29                 2.2e-27     91
## 30                 8.4e-27    156
##                                           Pathways
## 1                          lipid metabolic process
## 2                    response to nitrogen compound
## 3                          regulation of transport
## 4                          transmembrane transport
## 5           response to oxygen-containing compound
## 6              response to organonitrogen compound
## 7                                    ion transport
## 8        positive regulation of molecular function
## 9                            immune system process
## 10                positive regulation of transport
## 11       positive regulation of hydrolase activity
## 12 cellular response to oxygen-containing compound
## 13          positive regulation of GTPase activity
## 14                regulation of hydrolase activity
## 15               intracellular signal transduction
## 16            ribonucleoprotein complex biogenesis
## 17                         ncRNA metabolic process
## 18                             ribosome biogenesis
## 19                                  RNA processing
## 20    organonitrogen compound biosynthetic process
## 21                                ncRNA processing
## 22       organonitrogen compound metabolic process
## 23                                 rRNA processing
## 24                          rRNA metabolic process
## 25                      amide biosynthetic process
## 26                                     translation
## 27            single-organism biosynthetic process
## 28                    peptide biosynthetic process
## 29                cellular amide metabolic process
## 30                small molecule metabolic process
 # PPI network retrieval and analysis
 input_nGenesPPI <- 100 #Number of top genes for PPI retrieval and analysis 
 stringDB_network1(1) #Show PPI network 

 write(stringDB_network_link(), 'PPI_results.html') # write results to html file 
## Warning:  we couldn't map to STRING 9% of your identifiers
 browseURL('PPI_results.html') # open in browser 

8. Pathway analysis

##########################
# 8. Pathway analysis 
##########################
 input_selectContrast1 <- 'null_mock-null_IR'   #select Comparison 
 #input_selectContrast1 = limma.out$comparisons[3] # manually set
 input_selectGO <- 'GOBP'   #Gene set category 
 #input_selectGO='custom' # if custom gmt file
 input_minSetSize <- 15 #Min size for gene set
 input_maxSetSize <- 2000   #Max size for gene set 
 # Read pathway data again 
 GeneSets.out <-readGeneSets( geneSetFile,
    convertedData.out, input_selectGO,input_selectOrg,
    c(input_minSetSize, input_maxSetSize)  ) 
 input_pathwayPvalCutoff <- 0.2 #FDR cutoff
 input_nPathwayShow <- 30   #Top pathways to show
 input_absoluteFold <- FALSE    #Use absolute values of fold-change?
 input_GenePvalCutoff <- 1  #FDR to remove genes 

 input_pathwayMethod = 1  # 1  GAGE
 gagePathwayData.out <- gagePathwayData()  # pathway analysis using GAGE  
 gagePathwayData.out  
##                                                Direction
## NcRNA metabolic process                             Down
## Ribosome biogenesis                                     
## NcRNA processing                                        
## Ribonucleoprotein complex biogenesis                    
## RRNA processing                                         
## RRNA metabolic process                                  
## RNA processing                                          
## RNA modification                                        
## TRNA metabolic process                                  
## Ribosomal large subunit biogenesis                      
## TRNA processing                                         
## TRNA modification                                       
## Maturation of SSU-rRNA                                  
## Protein folding                                         
## Translation                                             
## Organonitrogen compound biosynthetic process            
## Peptide biosynthetic process                            
## Ribonucleoprotein complex subunit organization          
## Ribonucleoprotein complex assembly                      
## Amide biosynthetic process                              
## Cellular amino acid metabolic process                   
## RNA localization                                        
## Ribosomal small subunit biogenesis                      
## MRNA metabolic process                                  
## Translational initiation                                
## DNA metabolic process                                   
## Nucleoside metabolic process                            
## Ribonucleoside metabolic process                        
## Peptide metabolic process                               
## MRNA processing                                         
##                                                          GAGE  analysis: null_mock vs null_IR
## NcRNA metabolic process                                               NcRNA metabolic process
## Ribosome biogenesis                                                       Ribosome biogenesis
## NcRNA processing                                                             NcRNA processing
## Ribonucleoprotein complex biogenesis                     Ribonucleoprotein complex biogenesis
## RRNA processing                                                               RRNA processing
## RRNA metabolic process                                                 RRNA metabolic process
## RNA processing                                                                 RNA processing
## RNA modification                                                             RNA modification
## TRNA metabolic process                                                 TRNA metabolic process
## Ribosomal large subunit biogenesis                         Ribosomal large subunit biogenesis
## TRNA processing                                                               TRNA processing
## TRNA modification                                                           TRNA modification
## Maturation of SSU-rRNA                                                 Maturation of SSU-rRNA
## Protein folding                                                               Protein folding
## Translation                                                                       Translation
## Organonitrogen compound biosynthetic process     Organonitrogen compound biosynthetic process
## Peptide biosynthetic process                                     Peptide biosynthetic process
## Ribonucleoprotein complex subunit organization Ribonucleoprotein complex subunit organization
## Ribonucleoprotein complex assembly                         Ribonucleoprotein complex assembly
## Amide biosynthetic process                                         Amide biosynthetic process
## Cellular amino acid metabolic process                   Cellular amino acid metabolic process
## RNA localization                                                             RNA localization
## Ribosomal small subunit biogenesis                         Ribosomal small subunit biogenesis
## MRNA metabolic process                                                 MRNA metabolic process
## Translational initiation                                             Translational initiation
## DNA metabolic process                                                   DNA metabolic process
## Nucleoside metabolic process                                     Nucleoside metabolic process
## Ribonucleoside metabolic process                             Ribonucleoside metabolic process
## Peptide metabolic process                                           Peptide metabolic process
## MRNA processing                                                               MRNA processing
##                                                statistic Genes adj.Pval
## NcRNA metabolic process                         -11.8689   154  3.7e-23
## Ribosome biogenesis                             -10.8929   110  1.1e-18
## NcRNA processing                                -10.5686   114  2.0e-18
## Ribonucleoprotein complex biogenesis             -10.562   154  2.9e-19
## RRNA processing                                  -9.7918    81  1.3e-14
## RRNA metabolic process                           -9.6836    82  1.6e-14
## RNA processing                                   -9.4752   216  5.3e-17
## RNA modification                                 -7.8231    42  3.5e-08
## TRNA metabolic process                           -7.5376    51  1.0e-08
## Ribosomal large subunit biogenesis               -6.7767    31  5.7e-06
## TRNA processing                                  -6.1338    30  1.2e-05
## TRNA modification                                -5.7857    20  4.8e-04
## Maturation of SSU-rRNA                           -5.6703    21  4.8e-04
## Protein folding                                  -5.6113    72  1.1e-05
## Translation                                      -5.6082   172  5.6e-06
## Organonitrogen compound biosynthetic process     -5.4954   397  5.7e-06
## Peptide biosynthetic process                     -5.1789   178  3.3e-05
## Ribonucleoprotein complex subunit organization   -5.1249    70  7.8e-05
## Ribonucleoprotein complex assembly               -4.9952    66  1.3e-04
## Amide biosynthetic process                       -4.9881   203  7.4e-05
## Cellular amino acid metabolic process            -4.9853    70  1.3e-04
## RNA localization                                 -4.9661    67  1.4e-04
## Ribosomal small subunit biogenesis               -4.8412    27  7.3e-04
## MRNA metabolic process                           -4.6528   120  3.4e-04
## Translational initiation                         -4.6098    33  1.3e-03
## DNA metabolic process                            -4.4818   288  4.9e-04
## Nucleoside metabolic process                     -4.4411   132  6.9e-04
## Ribonucleoside metabolic process                 -4.2843   127  1.2e-03
## Peptide metabolic process                         -4.247   227  1.2e-03
## MRNA processing                                  -4.2087    94  1.7e-03
  pathwayListData.out = pathwayListData() 
 enrichmentPlot(pathwayListData.out, 25  ) 

  enrichmentNetwork(pathwayListData.out )  

  enrichmentNetworkPlotly(pathwayListData.out) 

 input_pathwayMethod = 3  # 1  fgsea 
 fgseaPathwayData.out <- fgseaPathwayData() #Pathway analysis using fgsea
 fgseaPathwayData.out 
##      Direction               GSEA  analysis: null_mock vs null_IR     NES
## 1247      Down                            NcRNA metabolic process  -2.636
## 1234                                             NcRNA processing -2.6069
## 1339                                          Ribosome biogenesis -2.5705
## 923                          Ribonucleoprotein complex biogenesis -2.5522
## 289                                               RRNA processing  -2.511
## 804                                        RRNA metabolic process -2.5083
## 294                                        TRNA metabolic process  -2.438
## 525                                               TRNA processing -2.4074
## 292                                                RNA processing -2.3922
## 619                                              RNA modification -2.3841
## 1343                           Ribosomal large subunit biogenesis -2.2866
## 295                                             TRNA modification -2.2736
## 322                         Cellular amino acid metabolic process -2.2331
## 304                                               Protein folding -2.2127
## 2352           RNA phosphodiester bond hydrolysis endonucleolytic -2.1758
## 992                                        Maturation of SSU-rRNA -2.1705
## 300                                                   Translation -2.1502
## 2262               Ribonucleoprotein complex subunit organization -2.1476
## 925                            Ribonucleoprotein complex assembly -2.1474
## 32                                           Telomere maintenance -2.1332
## 1099                                        Telomere organization -2.1332
## 301                                      Translational initiation  -2.106
## 303                   TRNA aminoacylation for protein translation -2.0919
## 1402                                 Peptide biosynthetic process -2.0795
## 298                                              RNA localization -2.0793
## 731                 Telomere maintenance via telomere lengthening -2.0743
## 1344                           Ribosomal small subunit biogenesis -2.0724
## 1400                                        Amino acid activation -2.0694
## 1401                                          TRNA aminoacylation -2.0694
## 2124                           Chaperone-mediated protein folding -2.0685
##      Genes adj.Pval
## 1247   154  3.6e-03
## 1234   114  3.6e-03
## 1339   110  3.6e-03
## 923    154  3.6e-03
## 289     81  3.6e-03
## 804     82  3.6e-03
## 294     51  3.6e-03
## 525     30  3.6e-03
## 292    216  3.6e-03
## 619     42  3.6e-03
## 1343    31  3.6e-03
## 295     20  3.6e-03
## 322     70  3.6e-03
## 304     72  3.6e-03
## 2352    19  3.6e-03
## 992     21  3.6e-03
## 300    172  3.6e-03
## 2262    70  3.6e-03
## 925     66  3.6e-03
## 32      35  3.6e-03
## 1099    35  3.6e-03
## 301     33  3.6e-03
## 303     19  3.6e-03
## 1402   178  3.6e-03
## 298     67  3.6e-03
## 731     24  3.6e-03
## 1344    27  3.6e-03
## 1400    20  3.6e-03
## 1401    20  3.6e-03
## 2124    29  3.6e-03
  pathwayListData.out = pathwayListData() 
 enrichmentPlot(pathwayListData.out, 25  ) 

  enrichmentNetwork(pathwayListData.out )  

  enrichmentNetworkPlotly(pathwayListData.out) 

   PGSEAplot() # pathway analysis using PGSEA 
## 
## Computing P values using ANOVA

9. Chromosome

 input_selectContrast2 <- 'null_mock-null_IR'   #select Comparison 
 #input_selectContrast2 = limma.out$comparisons[3] # manually set
 input_limmaPvalViz <- 0.1  #FDR to filter genes
 input_limmaFCViz <- 2  #FDR to filter genes 
 genomePlotly() # shows fold-changes on the genome 
## Warning in eval(quote(list(...)), env): NAs introduced by coercion
## Warning in genomePlotly(): NAs introduced by coercion

## 10. Biclustering

 input_nGenesBiclust <- 1000    #Top genes for biclustering
 input_biclustMethod <- 'BCCC()'    #Method: 'BCCC', 'QUBIC', 'runibic' ... 
 biclustering.out = biclustering()  # run analysis

 input_selectBicluster <- 1 #select a cluster 
 biclustHeatmap()   # heatmap for selected cluster 

 input_selectGO4 <- 'GOBP'  #Gene set 
 # Read pathway data again 
 GeneSets.out <-readGeneSets( geneSetFile,
    convertedData.out, input_selectGO4,input_selectOrg,
    c(input_minSetSize, input_maxSetSize)  )  
 geneListBclustGO()  # enrichment analysis 
##                                           adj.Pval Genes
## Immune system process                     3.1e-124   207
## Organonitrogen compound metabolic process  1.5e-89   157
## Positive regulation of molecular function  4.7e-86   148
## Cell death                                 5.4e-86   155
## Homeostatic process                        2.0e-83   139
## Macromolecular complex assembly            3.2e-81   136
## Programmed cell death                      3.5e-80   146
## Apoptotic process                          2.9e-79   144
## Regulation of cell death                   1.1e-75   132
## Cell-cell adhesion                         2.2e-75   117
##                                                                            Pathways
## Immune system process                                         Immune system process
## Organonitrogen compound metabolic process Organonitrogen compound metabolic process
## Positive regulation of molecular function Positive regulation of molecular function
## Cell death                                                               Cell death
## Homeostatic process                                             Homeostatic process
## Macromolecular complex assembly                     Macromolecular complex assembly
## Programmed cell death                                         Programmed cell death
## Apoptotic process                                                 Apoptotic process
## Regulation of cell death                                   Regulation of cell death
## Cell-cell adhesion                                               Cell-cell adhesion

11. Co-expression network

##########################
# 11. Co-expression network 
##########################
 input_mySoftPower <- 5 #SoftPower to cutoff
 input_nGenesNetwork <- 1000    #Number of top genes
 input_minModuleSize <- 20  #Module size minimum 
 wgcna.out = wgcna()   # run WGCNA
## ==========================================================================
## *
## *  Package WGCNA 1.63 loaded.
## *
## *    Important note: It appears that your system supports multi-threading,
## *    but it is not enabled within WGCNA in R. 
## *    To allow multi-threading within WGCNA with all available cores, use 
## *
## *          allowWGCNAThreads()
## *
## *    within R. Use disableWGCNAThreads() to disable threading if necessary.
## *    Alternatively, set the following environment variable on your system:
## *
## *          ALLOW_WGCNA_THREADS=<number_of_processors>
## *
## *    for example 
## *
## *          ALLOW_WGCNA_THREADS=8
## *
## *    To set the environment variable in linux bash shell, type 
## *
## *           export ALLOW_WGCNA_THREADS=8
## *
## *     before running R. Other operating systems or shells will
## *     have a similar command to achieve the same aim.
## *
## ==========================================================================
## Warning: executing %dopar% sequentially: no parallel backend registered
##    Power SFT.R.sq   slope truncated.R.sq mean.k. median.k. max.k.
## 1      1   0.9720  2.7600          0.966   572.0     612.0    711
## 2      2   0.9490  1.2200          0.956   412.0     449.0    578
## 3      3   0.9050  0.7070          0.923   325.0     353.0    495
## 4      4   0.7480  0.4190          0.797   269.0     287.0    439
## 5      5   0.5350  0.2490          0.615   230.0     239.0    398
## 6      6   0.1820  0.1090          0.255   200.0     204.0    367
## 7      7   0.0040  0.0131          0.185   178.0     175.0    342
## 8      8   0.0965 -0.0738          0.144   160.0     152.0    321
## 9      9   0.2500 -0.1240          0.238   145.0     133.0    303
## 10    10   0.4570 -0.1930          0.431   132.0     118.0    287
## 11    12   0.5670 -0.2940          0.514   113.0      98.9    263
## 12    14   0.6340 -0.3830          0.587    97.8      84.7    243
## 13    16   0.6890 -0.4420          0.655    86.2      73.5    227
## 14    18   0.7120 -0.4930          0.673    76.9      63.8    213
## 15    20   0.7290 -0.5420          0.683    69.3      55.7    201
## TOM calculation: adjacency..
## ..will not use multithreading.
##  Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
 softPower()  # soft power curve 

  modulePlot()  # plot modules  

  listWGCNA.Modules.out = listWGCNA.Modules() #modules
 input_selectGO5 <- 'GOBP'  #Gene set 
 # Read pathway data again 
 GeneSets.out <-readGeneSets( geneSetFile,
    convertedData.out, input_selectGO5,input_selectOrg,
    c(input_minSetSize, input_maxSetSize)  ) 
 input_selectWGCNA.Module <- '1. turquoise (495 genes)' #Select a module
 input_topGenesNetwork <- 10    #SoftPower to cutoff
 input_edgeThreshold <- 0.4 #Number of top genes 
 moduleNetwork()    # show network of top genes in selected module
##  softConnectivity: FYI: connecitivty of genes with less than 4 valid samples will be returned as NA.
##  ..calculating connectivities..

 input_removeRedudantSets <- TRUE   #Remove redundant gene sets 
 networkModuleGO()  # Enrichment analysis of selected module 
##                                              adj.Pval Genes
## Organonitrogen compound metabolic process     6.2e-68   106
## Ribonucleoprotein complex biogenesis          3.0e-62    61
## RNA processing                                6.9e-62    68
## Organonitrogen compound biosynthetic process  4.0e-61    84
## NcRNA metabolic process                       2.4e-56    57
## Immune system process                         4.2e-53   101
## Macromolecular complex assembly               2.0e-49    81
## Ribosome biogenesis                           3.0e-48    46
## Establishment of protein localization         3.4e-48    84
## Cellular amide metabolic process              4.7e-48    63
##                                                                                  Pathways
## Organonitrogen compound metabolic process       Organonitrogen compound metabolic process
## Ribonucleoprotein complex biogenesis                 Ribonucleoprotein complex biogenesis
## RNA processing                                                             RNA processing
## Organonitrogen compound biosynthetic process Organonitrogen compound biosynthetic process
## NcRNA metabolic process                                           NcRNA metabolic process
## Immune system process                                               Immune system process
## Macromolecular complex assembly                           Macromolecular complex assembly
## Ribosome biogenesis                                                   Ribosome biogenesis
## Establishment of protein localization               Establishment of protein localization
## Cellular amide metabolic process                         Cellular amide metabolic process