Supplementary Protocol 5 – Create Multi dataset Enrichment Map and Theme investigation

#make sure the RCy3 is installed. 
if(!requireNamespace("RCy3", quietly = TRUE)){
  if (!requireNamespace("BiocManager", quietly = TRUE)){
      install.packages("BiocManager")
    }
  BiocManager::install("RCy3")
}
if(!requireNamespace("RCurl", quietly = TRUE)){
  install.packages("RCurl")
}
library(RCy3)
library(RCurl)

0.1 Configurable Parameters

In order to run GSEA automatically through the notebook you will need to download the gsea jar from here. Specify the exact path to the gsea jar in the parameters in order to automatically compute enrichments using GSEA.

working_dir <- file.path(".",params$working_dir)
#path to GSEA jar 
# In order to run GSEA automatically you need to speciry the path to the gsea jar file.
#With the latest release of gsea (4.0.2) they no longer release a bundled jar
# and instead release a scriptted way to launch the gsea client.
# specify the java version as 11 if you are using the later version gsea
# the gsea_jar also needs to be the full path to the GSEA 4.0.2 directory that you
# downloaded from GSEA. for example (/Users/johnsmith/GSEA_4.0.2/gsea-cli.sh)
gsea_jar <- params$gsea_jar
java_version <- params$java_version
#Gsea takes a long time to run.  If you have already run GSEA manually or previously there is no need to re-run GSEA.  Make sure the 
# gsea results are in the current directory and the notebook will be able to find them and use them.
run_gsea = params$run_gsea

0.2 Download the latest pathway definition file

Only Human, Mouse and Rat gene set files are currently available on the baderlab downloads site. If you are working with a species other than human (and it is either rat or mouse) change the gmt_url below to correct species. Check here to see all available species.

gmt_url = "http://download.baderlab.org/EM_Genesets/current_release/Human/symbol/"
#list all the files on the server
filenames = RCurl::getURL(gmt_url)
tc = textConnection(filenames)

Found more than one class "textConnection" in cache; using the first, from namespace 'BiocGenerics'
Also defined by ‘RJSONIO’
Found more than one class "textConnection" in cache; using the first, from namespace 'BiocGenerics'
Also defined by ‘RJSONIO’
Found more than one class "textConnection" in cache; using the first, from namespace 'BiocGenerics'
Also defined by ‘RJSONIO’
Found more than one class "textConnection" in cache; using the first, from namespace 'BiocGenerics'
Also defined by ‘RJSONIO’
Found more than one class "textConnection" in cache; using the first, from namespace 'BiocGenerics'
Also defined by ‘RJSONIO’

contents = readLines(tc)
close(tc)
#get the gmt that has all the pathways and does not include terms inferred from electronic annotations(IEA)
#start with gmt file that has pathways only
rx = gregexpr("(?<=<a href=\")(.*.GOBP_AllPathways_no_GO_iea.*.)(.gmt)(?=\">)",
  contents, perl = TRUE)
gmt_file = unlist(regmatches(contents, rx))
dest_gmt_file <- file.path(multi_em_dir,gmt_file )
download.file(
    paste(gmt_url,gmt_file,sep=""),
    destfile=dest_gmt_file
)

Get all the rank files in the working directory

rank_files <- list.files(working_dir)[grep(pattern = "ranks.rnk",list.files(working_dir))]

Create a directory to hold the pathway analysis folders.

** This is an important step. In order for Enrichment map to create a multi dataset map the assumption is that all the results files will be together in a directory where each sub directory is an individual analysis **

multi_em_dir <- file.path(working_dir,"Multi_dataset_EM_analysis")
if(!dir.exists(multi_em_dir )){
  dir.create(multi_em_dir)
} else {
  multi_em_dir <- file.path(working_dir,paste("Multi_dataset_EM_analysis",Sys.Date(),sep = "_"))
  dir.create(multi_em_dir)
}
# copy the rank files to the new directory - if they aren't already there
if(length(list.files(multi_em_dir,pattern=".rnk")) == 0){
  file.copy(file.path(working_dir, rank_files), file.path(multi_em_dir,rank_files))
}

[1] TRUE TRUE TRUE TRUE TRUE

#after the files are copied over change the working dir to the newly 
# create directory
working_dir <- multi_em_dir

0.3 Run GSEA

(GSEA)[http://software.broadinstitute.org/gsea/index.jsp] is a stand alone java program with many customizable options. It can be easily run through its integrated user interface. To make this a seemless pipeline we can run GSEA from the command line with a set of options. Any of the supplied options can be customized and there are many additional options that can be specified. For more details see (here)[http://software.broadinstitute.org/gsea/doc/GSEAUserGuideTEXT.htm#_Running_GSEA_from]

In the below command the following options have been specified:

• rnk - path to the rank file
• gmx - path to the gene set definition (gmt) file
• collapse - true/false indicates whether the expression/rnk file needs to be collapsed from probes to gene symbols
• nperm - number of permutations
• permute - permute gene sets or phentoypes. For GSEA preranked you can only permute genesets.
• scoring_scheme -
• rpt_label - name of the directory with output
• num - number of results to plot output file for
• rnd_seed - random seed to use
• set_max - maximum size for individual gene sets. In GSEA interface this is set to 500 but we prefer to use a more stringent setting of 200.
• set_min - minimum size for individual gene sets
• zip_report - true/false to zip output directory
• out - directory where to place the result directory.
• gui - true/false. When running GSEA from the commandline this needs to be false.

If you are using GSEA 4.0.2 you will need to update your system to use Java 11. The GSEA command line interface has changed slightly in this version of GSEA. Instead of launching GSEA using java you need to use one of the scripts supplied by GSEA. (gsea-cli.bat for Windows and gsea-cli.sh for Mac or linux systems).

In the GSEA 4.0.2 command the following options can been specified:

• rnk - path to the rank file
• gmx - path to the gene set definition (gmt) file
• collapse - true/false indicates whether the expression/rnk file needs to be collapsed from probes to gene symbols
• nperm - number of permutations
• scoring_scheme -
• rpt_label - name of the directory with output
• rnd_seed - random seed to use
• set_max - maximum size for individual gene sets. In GSEA interface this is set to 500 but we prefer to use a more stringent setting of 200.
• set_min - minimum size for individual gene sets
• zip_report - true/false to zip output directory
• out - directory where to place the result directory.

If you have encounter an out of memory error you will need to open you gsea-cli script and manually update your xmx parameter (more information on this can be found in the GSEA_4.0.2/README file)

#run GSEA for each of the rank files.
for(i in 1:length(rank_files)){
    current_rank_file <- rank_files[i]
    
    analysis_name <- unlist(strsplit(rank_files[i],split = "_"))[1]
    #if you are using GSEA 4.0.2 then you need to use the command script
    # as opposed to launching GSEA through java. 
    # in the later version of GSEA command line implementation the following 
    # parameters are no longer valid: -permute gene_set,  -num 100, -gui false
    # no longer need to specify the whole path to the GseaPreranked package
    if(run_gsea && java_version == "11"){
      command <- paste("",gsea_jar,  "GSEAPreRanked -gmx ", dest_gmt_file, "-rnk" ,file.path(working_dir,current_rank_file ), "-collapse false -nperm 1000 -scoring_scheme weighted -rpt_label ",analysis_name,"  -plot_top_x 20 -rnd_seed 12345  -set_max 200 -set_min 15 -zip_report false -out" ,working_dir, " > gsea_output.txt",sep=" ")
      system(command)
    } else {
      command <- paste("java  -Xmx1G -cp",gsea_jar,  "xtools.gsea.GseaPreranked -gmx", dest_gmt_file, "-rnk" ,file.path(working_dir,current_rank_file ), "-collapse false -nperm 1000 -permute gene_set -scoring_scheme weighted -rpt_label ",analysis_name,"  -num 100 -plot_top_x 20 -rnd_seed 12345  -set_max 200 -set_min 15 -zip_report false -out" ,working_dir, "-gui false > gsea_output.txt",sep=" ")
      system(command)
    }
}

WARNING: package com.sun.java.swing.plaf.windows not in java.desktop
WARNING: package sun.awt.windows not in java.desktop
java.io.FileNotFoundException: File not found: /Users/ruthisserlin/Dropbox (Bader Lab)/Ruth Isserlin's files/Sourcecode/Cytoscape_workflows/EnrichmentMapPipeline/./data/Multi_dataset_EM_analysis_2020-07-24/Human_GOBP_AllPathways_no_GO_iea_July_01_2020_symbol.gmt
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.ParserFactory.createInputStream(ParserFactory.java:1056)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.ParserFactory.read(ParserFactory.java:799)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.ParserFactory.read(ParserFactory.java:787)
    at org.gsea_msigdb.gsea/xtools.api.param.GeneSetMatrixMultiChooserParam._getObjects(GeneSetMatrixMultiChooserParam.java:123)
    at org.gsea_msigdb.gsea/xtools.api.param.GeneSetMatrixMultiChooserParam._getGeneSets(GeneSetMatrixMultiChooserParam.java:167)
    at org.gsea_msigdb.gsea/xtools.api.param.GeneSetMatrixMultiChooserParam.getGeneSetMatrixCombo(GeneSetMatrixMultiChooserParam.java:62)
    at org.gsea_msigdb.gsea/xtools.gsea.GseaPreranked.execute(GseaPreranked.java:101)
    at org.gsea_msigdb.gsea/xtools.api.AbstractTool.module_main(AbstractTool.java:434)
    at org.gsea_msigdb.gsea/org.genepattern.modules.GseaPrerankedWrapper.main(GseaPrerankedWrapper.java:228)
    at org.gsea_msigdb.gsea/xapps.gsea.CLI.main(CLI.java:31)
WARNING: package com.sun.java.swing.plaf.windows not in java.desktop
WARNING: package sun.awt.windows not in java.desktop
java.io.FileNotFoundException: File not found: /Users/ruthisserlin/Dropbox (Bader Lab)/Ruth Isserlin's files/Sourcecode/Cytoscape_workflows/EnrichmentMapPipeline/./data/Multi_dataset_EM_analysis_2020-07-24/Human_GOBP_AllPathways_no_GO_iea_July_01_2020_symbol.gmt
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.ParserFactory.createInputStream(ParserFactory.java:1056)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.ParserFactory.read(ParserFactory.java:799)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.ParserFactory.read(ParserFactory.java:787)
    at org.gsea_msigdb.gsea/xtools.api.param.GeneSetMatrixMultiChooserParam._getObjects(GeneSetMatrixMultiChooserParam.java:123)
    at org.gsea_msigdb.gsea/xtools.api.param.GeneSetMatrixMultiChooserParam._getGeneSets(GeneSetMatrixMultiChooserParam.java:167)
    at org.gsea_msigdb.gsea/xtools.api.param.GeneSetMatrixMultiChooserParam.getGeneSetMatrixCombo(GeneSetMatrixMultiChooserParam.java:62)
    at org.gsea_msigdb.gsea/xtools.gsea.GseaPreranked.execute(GseaPreranked.java:101)
    at org.gsea_msigdb.gsea/xtools.api.AbstractTool.module_main(AbstractTool.java:434)
    at org.gsea_msigdb.gsea/org.genepattern.modules.GseaPrerankedWrapper.main(GseaPrerankedWrapper.java:228)
    at org.gsea_msigdb.gsea/xapps.gsea.CLI.main(CLI.java:31)
WARNING: package com.sun.java.swing.plaf.windows not in java.desktop
WARNING: package sun.awt.windows not in java.desktop
java.io.FileNotFoundException: File not found: /Users/ruthisserlin/Dropbox (Bader Lab)/Ruth Isserlin's files/Sourcecode/Cytoscape_workflows/EnrichmentMapPipeline/./data/Multi_dataset_EM_analysis_2020-07-24/Human_GOBP_AllPathways_no_GO_iea_July_01_2020_symbol.gmt
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.ParserFactory.createInputStream(ParserFactory.java:1056)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.ParserFactory.read(ParserFactory.java:799)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.ParserFactory.read(ParserFactory.java:787)
    at org.gsea_msigdb.gsea/xtools.api.param.GeneSetMatrixMultiChooserParam._getObjects(GeneSetMatrixMultiChooserParam.java:123)
    at org.gsea_msigdb.gsea/xtools.api.param.GeneSetMatrixMultiChooserParam._getGeneSets(GeneSetMatrixMultiChooserParam.java:167)
    at org.gsea_msigdb.gsea/xtools.api.param.GeneSetMatrixMultiChooserParam.getGeneSetMatrixCombo(GeneSetMatrixMultiChooserParam.java:62)
    at org.gsea_msigdb.gsea/xtools.gsea.GseaPreranked.execute(GseaPreranked.java:101)
    at org.gsea_msigdb.gsea/xtools.api.AbstractTool.module_main(AbstractTool.java:434)
    at org.gsea_msigdb.gsea/org.genepattern.modules.GseaPrerankedWrapper.main(GseaPrerankedWrapper.java:228)
    at org.gsea_msigdb.gsea/xapps.gsea.CLI.main(CLI.java:31)
WARNING: package com.sun.java.swing.plaf.windows not in java.desktop
WARNING: package sun.awt.windows not in java.desktop
java.io.FileNotFoundException: File not found: /Users/ruthisserlin/Dropbox (Bader Lab)/Ruth Isserlin's files/Sourcecode/Cytoscape_workflows/EnrichmentMapPipeline/./data/Multi_dataset_EM_analysis_2020-07-24/Human_GOBP_AllPathways_no_GO_iea_July_01_2020_symbol.gmt
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.ParserFactory.createInputStream(ParserFactory.java:1056)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.ParserFactory.read(ParserFactory.java:799)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.ParserFactory.read(ParserFactory.java:787)
    at org.gsea_msigdb.gsea/xtools.api.param.GeneSetMatrixMultiChooserParam._getObjects(GeneSetMatrixMultiChooserParam.java:123)
    at org.gsea_msigdb.gsea/xtools.api.param.GeneSetMatrixMultiChooserParam._getGeneSets(GeneSetMatrixMultiChooserParam.java:167)
    at org.gsea_msigdb.gsea/xtools.api.param.GeneSetMatrixMultiChooserParam.getGeneSetMatrixCombo(GeneSetMatrixMultiChooserParam.java:62)
    at org.gsea_msigdb.gsea/xtools.gsea.GseaPreranked.execute(GseaPreranked.java:101)
    at org.gsea_msigdb.gsea/xtools.api.AbstractTool.module_main(AbstractTool.java:434)
    at org.gsea_msigdb.gsea/org.genepattern.modules.GseaPrerankedWrapper.main(GseaPrerankedWrapper.java:228)
    at org.gsea_msigdb.gsea/xapps.gsea.CLI.main(CLI.java:31)

0.4 Launch Cytoscape

Create EM through Cyrest interface - make sure you open cytoscape with a -R 1234 (to enable rest functionality) and allow R to talk directly to cytoscape.

Launch Cytoscape (by default cytoscape will automatically enable rest so as long as cytoscape 3.3 or higher is open R should be able to communicate with it)

0.5 Set up connection from R to cytoscape

   cytoscapePing ()

[1] "You are connected to Cytoscape!"

    cytoscapeVersionInfo ()

      apiVersion cytoscapeVersion 
            "v1"          "3.8.0" 

0.6 Create a multi dataset EM

#defined threshold for GSEA enrichments (need to be strings for cyrest call)
pvalue_gsea_threshold <- params$pval_thresh
qvalue_gsea_threshold <- params$fdr_thresh
similarity_threshold <- "0.375"
similarity_metric = "COMBINED"
cur_model_name <- "Multi_dataset_EM"
gsea_results_path <- multi_em_dir
#although there is a gmt file in the gsea edb results directory it have been filtered to 
#contain only genes represented in the expression set.  If you use this fltered file you 
#will get different pathway connectivity depending on the dataset being used.  We recommend 
#using original gmt file used for the gsea analysis and not the filtered one in the results directory.
gmt_gsea_file <- file.path(getwd(),dest_gmt_file)
gsea_root_dir <- file.path(getwd(),multi_em_dir)
#######################################
#create EM
current_network_name <- paste(cur_model_name,pvalue_gsea_threshold,qvalue_gsea_threshold,sep="_")
em_command = paste('enrichmentmap mastermap commonGMTFile=',gmt_gsea_file,
                   'pvalue=',pvalue_gsea_threshold, 'qvalue=',qvalue_gsea_threshold,
                   'similaritycutoff=',similarity_threshold,
                   'coefficients=',similarity_metric,'rootFolder=', 
                   gsea_root_dir,
                   'filterByExpressions=false',
                   'commonExpressionFile=',file.path(getwd(),params$expression_file),
                   sep=" ")
#enrichment map command will return the suid of newly created network.
response <- RCy3::commandsGET(em_command)

Error in curl::curl_fetch_memory(url, handle = handle) : 
  Operation was aborted by an application callback

0.7 Change the coloring on the network

By default the network will colour each section of the pie using the NES value (GSEA analysis) and FDR value for any other analysis. This doesn't let you see clearly which pathways are significant in each dataset. By changing the colouring type to "DATA_SET" each node will only be coloured if the pathway is significant in the given dataset (by the user defined thresholds).

** The colours are set by default. There is no progrommatic way to change the colours. It can be done manually though.**

response <- RCy3::commandsGET('enrichmentmap chart data="DATA_SET"')

0.8 Define themes

Annotate the Enrichment map in order to calculate the themes

response <- RCy3::commandsGET(paste('autoannotate annotate-clusterBoosted network=',current_network_suid,sep=""))

0.9 Create a summary network

theme_network <- as.numeric(RCy3::commandsGET(paste('autoannotate summary network=', current_network_suid,sep="")))

0.10 Get all the theme information

Depending on the data type of the node attributes when the theme network is created they will be collapsed differently. You can set this behaviour manually in cytoscape from Edit -> Preference -> Group Preferences. For more info on this see here

  default_node_table <- RCy3::getTableColumns(table= "node",network = theme_network)

0.11 Output all Themes found in all the datasets

num_datasets <- length(rank_files)
dataset_chart <- t(sapply(default_node_table$`EnrichmentMap::Dataset_Chart`, function (x) { return (x)}))
rownames(dataset_chart) <- default_node_table$name
#calculate the number of datasets for each theme
dataset_per_theme <- rowSums(dataset_chart)
#output the themes that are found in the most datasets
rownames(dataset_chart)[which(dataset_per_theme == max(dataset_per_theme))]

[1] "axonemal dynein microtubule"

0.12 Output the pathways specific to each dataset

#Go through each index of the dataset chart and output the pathways specific for each group
specific_to_theme_summary <- list()
for(i in 1:dim(dataset_chart)[2]){
  current_dataset <- i
  
  specific_to_theme_summary[[i]] <- rownames(dataset_chart)[intersect(which(dataset_per_theme == 1) , which(dataset_chart[,i]==1))]
}
#unfortunately there is no way currently to get which id corresponds to which dataset.  Guess which according to 
# the order of the columns in the returned table.
names(specific_to_theme_summary) <- unlist(lapply(colnames(default_node_table)[grep(colnames(default_node_table),pattern="pvalue")],FUN=function(x){unlist(strsplit(x, split =" "))[2]}))
specific_to_theme_summary

$`(DiffvsRest.GseaPreranked)`
[1] "signal kinetochores amplification" "chromosome segregation sister"     "pid aurora signaling"             

$`(ImmunovsRest.GseaPreranked)`
 [1] "camera-type eye sensory"           "central nervous neuron"            "limb morphogenesis appendage"     
 [4] "ncam signaling neurite"            "cardiac septum morphogenesis"      "bmp tgf-beta pathway"             
 [7] "positive interleukin-4 production" "nephron epithelium kidney"         "digestive tract development"      
[10] "negative regulation smoothened"    "involved wnt carcinoma"            "growth factor stimulus"           
[13] "lymphocyte activation selection"  

$`(MesenvsRest.GseaPreranked)`
 [1] "gene silencing rna"                              "trna modification processing"                   
 [3] "mitochondrial translational termination"         "glycosaminoglycan chondroitin sulfate"          
 [5] "substrate adhesion-dependent spreading"          "mature transcript cytoplasm"                    
 [7] "coupled electron nadh"                           "negative blood vessel"                          
 [9] "pid avb3 integrin"                               "o-glycosylation domain-containing glycosylation"
[11] "flagellum-dependent motility cilium-dependent"   "assembly collagen formation"                    
[13] "blood vessel morphogenesis"                      "positive vasculature endothelial"               
[15] "protein-dna nucleosome centromere"               "spliceosomal snrnp ribonucleoprotein"           
[17] "pre-mrna splicing u12"                           "regulation chondrocyte cartilage"               
[19] "smooth muscle proliferation"                     "ribosomal subunit biogenesis"                   
[21] "replication dna dna-dependent"                   "regulation igf activity"                        
[23] "nucleotide excision repair"                      "beta3 basement membranes"                       
[25] "primary germ endoderm"                          

$`(ProlifvsRest.GseaPreranked)`
 [1] "erk1 erk2 cascade"                      "complement maturation protein"         
 [3] "platelet activation coagulation"        "blood coagulation hemostasis"          
 [5] "pid il4 2pathway"                       "osteoclast differentiation myeloid"    
 [7] "leukocyte homeostasis number"           "regulation secretion negative"         
 [9] "pid il8 cxcr2"                          "peptidyl-tyrosine phosphorylation stat"
[11] "pid gmcsf pathway"                      "cellular lipoprotein particle"         
[13] "beta2 integrin cell"                    "pattern receptor toll-like"