Pathway and Network Analysis of Omics Data Integrated Assignment

This work is licensed under a Creative Commons Attribution-ShareAlike 3.0 Unported License. This means that you are able to copy, share and modify the work, as long as the result is distributed under the same license.

Pathway and Network Analysis Integrated Assignment

By Veronique Voisin

Goal

Familiarize yourself with g:Profiler, GSEA , EnrichmentMap using the Esophageal adenocarcinoma gene expression data (DATASET 1). Familiarize yourself with ReactomeFI and GeneMANIA using mutation data (DATASET 2).

NOTE: Network layouts are flexible and can be rearranged. What you see when you perform these exercises may not be identical to what you see in the tutorial, or what you have seen other times that you have performed the exercises. Exact layouts and predictions can also be affected by updates to the networks database that GeneMANIA uses. However it is expected that the network weights and predicted genes will be similar to those shown here.

DATASET 1

Background

Gene expression data from Esophageal adenocarcinoma (EAC) is used for this first part of the integrated assignment. Esophageal adenocarcinoma (EAC) has a rising incidence and a 5-year survival of only 15%. The single major risk factor for development of EAC is chronic heartburn, which eventually leads to a change in the lining of the esophagus called Barrett’s Esophagus (BE).

Specimens were collected from patients with normal esophagus (NE) and Barrett’s esophagus (BE). RNA was extracted from these samples and expression profiling was assessed using Affymetrix HG-U133A microarray PMID:24714516. Differentially expressed genes between BE and NE were determined.

Data processing

The Affymetrix data are stored in the Gene Expression Omnibus (GEO) repository under the accession number GSE39491 PMID:24714516. The RMA (Robust Multichip Average) normalized data were downloaded from GEO and further processed using the Bioconductor package limma to estimate differential expression between the groups. The results of the limma t-tests were corrected for multiple hypothesis testing using the Benjamini-HochBerg method (FDR).

For g:Profiler, genes with a FDR equal or less than 0.0001 and a logFC of 2 were retrieved and stored in a text file. For GSEA, a rank file has been created by ranking the genes from the highest t statistics value (up-regulated in BE compared to NE) to the lowest t values (down-regulated in BE compared to NE). The code used to process the data is available from the course page. Please feel free to adapt it and use it with your own data.

PART 1: run g:Profiler

1. Open g:Profiler

2. In Options, check Significantly only, No electronic GO annotations

3. Choose the databases: GO BP and Reactome.

4. Set minimum size of functional category to 3 and maximum to 500. Set 2 for the Size of query/ term intersection.

5. Set Benjamini-Hochberg in the ‘Significance threshold’ box.

6. Run analysis of the genes differentially altered between BE and normal: copy and paste the gene list into the g:Profiler input window BEonly_genelist.txt.

7. Question: What is the most significant GO:term? What is the p-value for this GO:term?

8. Question: Is this p-value already corrected for multiple hypothesis testing? What type of correction is used for your current analysis?

PART 2: run g:Profiler p-value

1. Re-run the analysis with User p-value threshold set to 0.0001. Question: What has been changed?

PART 3: run g:Profiler ordered query

An important feature of g:Profiler is an ability to work with sorted or ranked lists. The top of such a list is given more weight in determining the functional connections to GO:terms and/or pathways (minimum hypergeometric test). Our gene list was initially sorted based on significance (FDR) in an increasing order (most significant genes (smallest FDR) at the top of the list).

1. Open g:Profiler in a new window and paste genes from BEonly_genelist.txt .

2. Use same options as part 2 (User p-value threshold set to 0.0001) and select also Ordered query.

3. Question: Do you seen any changes in the output in comparison to the analysis of the unordered gene list (PART 2)?

PART 4: save as Generic Enrichment Map output (BE)

Now we have to generate an output from the enrichment analysis and save it in appropriate format for EnrichmentMap. Please, change the output type to Generic Enrichment Map (TAB).

Run it using options used in PART 3. Download data in Generic EnrichmentMap (GEM) format and save it on your computer. We will need this file to create an enrichment map.

PART 5: save as Generic Enrichment Map output (NE)

Generate and save the Generic EnrichmentMap for genes in NConly_genelist.txt. It contains the genes specific of the normal tissue samples. Run g:Profiler with this list using same options as in PART 4 selecting Generic Enrichment Map (GEM) format as output type. We will need this file for EnrichmentMap.

PART 6: create an enrichment map

Create an enrichment map to visualize the outputs from g:Profiler. Let’s create an EnrichmentMap for the pathways that were enriched by the genes specific of the BE samples and one for the genes specific of the NE samples.

1. Make sure to rename your g:Profiler results so you know which one is BE and which one is NE. Create a directory and copy the results of your 2 g:Profiler runs ( 1 for BE and 1 fo NE) in this folder. Let’s name this folder ‘IA_gprofiler_results’. This folder should contains only these 2 files.

2. Open Cytoscape

3. Go: Apps and click on EnrichmentMap. A ‘Create Enrichment Map’ dialog box appears.

4. Drag and Drop your ‘IA_gprofiler_results’ folder in the ‘Data Sets:’ window. It populated automatically two data sets on for the BE results and one for the NE results. Make sure that for the 2 datasets the ‘Analysis Type’ is set to ‘Generic/gProfiler’ and that the g:Profiler result file has been correctly uploaded in the ‘Enrichments’ field.

5. As ‘Network Name’, uncheck ‘Use Default’ and write “BEvsNE” or the name of your choice for this network.

6. Set the ‘FDR q-value cutoff’ to 0.05.

7. ‘Build’ the map.

8. If successful, you will see a network where each node represents a pathway and edges connect pathways with shared genes. Blue edges connect nodes from dataset1 (BE in my case) and Green edges connect nodes from dataset2 (NE in my case).

9. In Control Panel and in the ‘EnrichmentMap’ tab, go to ‘Style’ and change the ‘Chart Data:’ to ‘Color by Data Set’. Now the nodes are colored in blue for dataset1 and in green for dataset2.

10. Try different layouts if you’d like. Zoom in and move nodes around to be able to read the labels.

11. Try to separate blue nodes and green nodes. In ‘Control Panel’ located at the left style, choose the ‘Select’ tab. Click on the ‘+’ plus sign to add a column to filter your data. Choose “NodeL:EM1_Colouring” corresponding to the NE datasets. Do not touch the slide bar and click on “Apply” located at the bottom of the tab. Now it has selected all the NE nodes. Click on one NE nodes and drag them all to the one side of the window, apart from the BE nodes. Click once on a blank space of the window to unselect the nodes.

12. Select a node of your choice. When the node is highlighted, the ‘EM Heat Map’ in ‘Table Panel’ will display the genes in this pathway that are overlapping with your input gene list. A gray square means that the gene is absent in the dataset. Note: you also could create and upload an expression file when you build the enrichment map, and the expression values for each gene in the pathways will be displayed here in the ‘EM Heat Map’.

13. Click on any edge (the line between nodes). In the ‘Table panel’ (‘EM Heat Map’) you should see a heatmap of all genes both gene-sets connected by this edge have in common.

14. Select several nodes and edges. EM Heat map will show the union of all genes (Genes: All) or genes in common (Genes:Common) in the selected gene sets.

15. In Control Panel, go to the EnrichmentMap tab. Change Q-value as well as Edge (Similarity) cutoffs and see how the network changes. Redo the layout. Save the file.

Question What conclusions can you make based on these networks?

BE map Q-value 0.01

NC map Q-value 0.01

Optional: when you build the enrichment map using the g:Profiler results, you have the option to upload a .gmt file (downloaded from the g:Profiler website) hsapiens.pathways.NAME.gmt and an expression file BE_vs_NE_expression.txt. The .gmt file provides the advantage to show not only the genes that are in your gene list and the enriched pathway (selected node) but all the genes in the original pathway. The expression file will be displayed in the EM Heat Map when a node is selected and it will display expression values for each gene in the selected pathway (selected node). Try it and adapt it to your needs.

PART 7: GSEA (run and create an enrichment map)

1. Launch GSEA.

2. Run GSEA using the rank file that has been created from the differential expression test comparing BE vs NE BEvsNE_ranks.rnk and the pathway file Human_GOBP_AllPathways_no_GO_iea_May_24_2016_symbol.gmt.

 * open GSEA and first import the files using the “Load data” window: upload the .rnk and .gmt files.

• Go the ‘Run GSEAPreranked’ window and select the correct gmt file and the rank file
• /!\ Use 100 permutations for the lab exercise /!\ but use 1000 for your own data analysis.
• Choose a name for your analysis, a destination folder and run GSEA.

1. Create an enrichment map:
   • Open Cytoscape and the EnrichmentMap app. The enrichment results are 2 excel files called gsea_report_for_na_neg and gsea_report_for_na_pos within the GSEA folder saved on your computer but you should be able to drag and drop the whole GSEA folder and that will populate the required fields automatically.

• use an FDR q-value cutoff of 0.05. Upload the expression file BE_vs_NE_expression.txt(right click, save link as).

1. Examine the results as you did for the g:Profiler map (e.g move nodes around, use the slide bar to adjust q value to 0.01 and redo the layout, separate blue and red nodes). Save the file. Save an image. Keep your session open for Part 8.

Optional: Autoannotate your map (see below screenshot for results)

PART 8: iRegulon

1. Export extracellular matrix genes.
   • Using your GSEA map at q-value 0.01, select all nodes from the extracellular matrix/collagen cluster. Go to Table Panel, and click on the Settings icon (located on the right, shaped like a ship wheel) and click on ‘Export as TXT’. Save the text file under the name ‘extracellularmatrix_cluster_genes.txt’ or use this file extracellularmatrix_cluster_genes.txt.
2. Import the extracellular matrix genes as a network.
   • In Cytoscape, go to the menu bar and select, File, Import, Network, File…
   • Browse your computer and select the ‘extracellularmatrix_cluster_genes.txt’ file and click on open.
   • An ‘Import Network From Table’ opens and in the table preview, make sure that the ‘Gene’column is the source node (green dot). Click on ‘OK’. A ‘Confirmation’ dialog box saying that ‘No edges will be created in the network’ opens. Click on ‘Yes’.
3. Select nodes and run iRegulon.
   • If successful, you should see a grid of gray nodes. Zoom out to see all the nodes and select them all using the mouse.
   • Go the Cytoscape menu and select ‘Apps’, ‘iRegulon’, ‘Predict regulators and targets’.
   • Click on ‘Submit’.
   • Observe the iRegulon results in the Results Panel.
4. Add the first 2 hits (TCF12 and BATF) to the network.
   • Go to the “Transcription Factors’ tab and click on the first hit (TCF12) to select it.
   • Add it to the network using the green ‘+’ button .
   • Execute the same steps for the second hit (BATF).
   • If successful, you should see targets of TCF12 and BATF linked to these 2 genes by edges (lines).
5. Create a subnetwork with all nodes connected to TCF12 and BATF.
   • using the mouse select TCF12 and all edges around this node and pressing the shift key, select also BATF and all the edges around this node. All selected edges should now be highlighted in red and the 2 transcription factors in yellow.
   • In the Cytoscape menu bar, go to Select, Nodes, Nodes connected by selected edges. More nodes should be selected now and the edges still highlighted in red.
   • Select the subnetwork icon (‘New Network from Selection (all edges)’)from the Cytoscape toolbar. If successful, you should have created a subnetwork containing only the targeted genes and the two transcription factors.
6. Arrange the network such that we can distinguish genes linked to TCF12 only , linked to BATF only or linked to both transcription factors.
   • go to the Cytoscape menu, Layout, Circular Layout, all Nodes. Feel free to use your own strategy.
7. optional. Import the .rnk file that we use for GSEA as attribute and color the nodes accordingly to the score values.
   • in the menu bar, select, File, Import, File, select the rank file and click on ‘Open’. A dialog box (‘Import Columns From Table’) opens. Uncheck ‘Use first line as column names’ and click ‘OK’. Rename Column2 as ‘myscore’. Click ‘OK’.
   • in Control Panel, go to Style and in the Node tab, expand the ‘Fill Color’ tab. Retrieve and select the ‘myscore’ column in the ‘column’ fied. Make sure that the ‘Mapping type’ is set to ‘Continuous Mapping’. The score should ranged from -13.16 to 13.16. Adjust the color if necessary.

Screenshot of resulting network:

DATASET 2

Stomach cancer or gastric cancer is a cancer developing from the lining of the stomach. The most common cause is infection by the bacteria Helicopter pylori, which accounts for more than 60% of cases. Certain types of ‘H. pylori’ have greater risks than others. Other common causes include eating pickled vegetables and smoking.

MutSig - is a mutation signal processing tool created by the Broad Institute. It estimates the significance of the gene mutation rate based on abundances of the mutations, clustering of the mutations in hotspots and conservation of the mutated positions.

The gene list for this assignment is the output from MutSig run based on Stomach Adenocarcinoma somatic mutations found in ~300 samples. It is publicly available through TCGA portal.

File provided: STAD_MutSig.txt

Goal: familiarize yourself with ReactomeFI and GeneMANIA.

PART 1: ReactomeFI

Create a network using ReactomeFI.

1. Open Cytoscape.
2. Choose App -> Reactome FI -> Gene set/mutation analysis
3. Upload STAD_MutSig.txt and built a network without linkers:

/!\ : the latest interface contains 2016 version of the ReactomeFI underlying database but this exercise has been prepared with the 2015 version. Please use the 2015 version to complete this exercise but the most up to date version when you work with your own data.

• tip: upload your file or copy and paste the gene names in the gene set field.

1. Fetch FI annotation (Hint: right click on a blank space, Reactome FI, Fetch FI Annotations). You should see arrows (directed graph).

2. Run Pathway enrichment (Hint: right click on the network panel and follow the path showed below). Question What is the pathway with the lowest (best) FDR?

3. Do a subnetwork of ErbB signaling pathway (K)

Hint:select the pathway in the table, that should highlight the genes in yellow. Use the subnetwork icon on the Cytoscape tool bar to create it (“New network from selection”).

1. Cluster the network and perform pathway enrichment on the network. Questions How many clusters did the analysis retrieve? What is the FDR value of the most significant pathway of module 0

PART 2: GeneMANIA

Use the same mutation data STAD_MutSig.txt to create a network using GeneMANIA in order to visualize which genes are known to physically interact with each other.

Install the latest version of the Cytoscape GeneMANIA app (Version 3.5.0)

• Go the Cytoscape menu bar and select Apps, App Manager… . In the App Manager dialog box, go to the ‘Check for Updates’ tab and update GeneMANIA if necessary. Note: if the ‘Check for Updates’ tab is blank, it means that you have installed the most recent versions of all apps.

1. Create the network

• In Cytoscape, go to Control Panel and locate and select the Network Tab in the Control Panel
• Make sure the GeneMANIA search provider is selected in the Network Search Bar.
• Choose Homo sapiens from the list of supported organisms
• Copy and paste the gene list STAD_MutSig.txt in the field.
• Locate the “More Options…” button at the right side of the field and only select ‘Physical interactions’ as ‘Interaction Networks’ and set 0 to the ‘Max Resultant Genes’.

• Click the “Search Network” button

• Note: The network may look slightly different compared to below screenshot if the underlying database has been updated since the screenshot was taken.

Screenshot of the output:

1. Explore the functions in the GeneMANIA Results Panel.
   • Go to ‘Results Panel’ located at the right side and select the GeneMANIA tab. Choose the ‘Functions’ tab to visualize the list of enriched GO gene-sets. Question Can you see which genes are included in these gene-sets? Hint: you can click on a function of your choice to see corresponding nodes highlighted in yellow.
2. Improve the visual style:

• Color nodes by function.
  • In Control Panel, select the ‘Style’ tab and go to the ‘Node’ panel.
  • Expand the ‘Fill Color’ field using the down arrow and set ‘Column’ to ‘annotation name’ which is the top field (/!\ not ‘annotations’). Select one annotation of your choice by clicking on the white space and choose a color. Repeat for 2 annotation names. For the current example, we have selected ‘BMP signaling pathway’, ‘Notch signaling pathway’ and’Wnt signaling pathway’. Hint: the annotation names are displayed in alphabetical order.
• Edge width (optional). In Control Panel, go to the ‘Edge’ panel. Expand the ‘Width’ field using the down arrow. A grah is displayed. Double click on the graph to select it and move the left and right handles up. Look at the changes on the network (suggested values are approximately 7 for the left handle and approximately 18 for the right handle). Click on OK.

1. Create a subnetwork containing CTNNB1 and connected genes
   • Locate CTNNB1, use the “First neighbors of selected nodes” icon (has the shape of 2 houses) in the toolbar to highlight genes connected to CTNNB1
   • Create a subnetwork using the approriate icon.
   • How many nodes do contain this subnetwork? Hint: Go to Control Panel, Network and look at the number of nodes corresponding to your subnetwork.

Note: on the above the screenshot, the CCBN1 color of the edges were changed by bypassing the style (select the edges, right click, edit, ByPass Style, Set Bypass to Selected Edges , change color in Byp. column)

---

Optional part 1: Launch a GeneMANIA search using the “Local Search” option (for big networks)

• In Cytoscape , open the GeneMANIA app and select ‘GeneMANIA Local Search’. Copy and paste the MutSig genes in the ‘Genes of Interest’ field.
• In Advanced Options, select only ‘Physical interactions’ as ‘Interaction Networks’ and set 0 in the “Find the top” 0 “related genes”.
• Click on ‘Start’.

!Warnings:

• If you use it for the first time and you haven’t installed data as it was said in the installation instructions, only install “CORE” data as the full data may take 1 hour to download.

Note: There are 2 ways to perform a GeneMANIA search. The first option using the Network search bar from the Control Panel is doing a seach by calling and connecting the GeneMANIA server (same as the website:https://genemania.org/). The other option as just showed here is to select GeneMANIA from the Apps menu and click on ‘Local Search…’. This option will use a database that is installed locally on your computer when you first use GeneMANIA. As it does not imply any connection to the server, this option is the best choice for large query, e.g input gene list size greater than 100 or resulting network containing more than 200 nodes.

The network and predicted functions should be the same as the ones obtained in part 2. Feel free to explore the network or follow the same steps as part 2.

---

Optional part 2: Use STRING from the Network Search Bar

STRING (Search Tool for the Retrieval of Interacting Genes/Proteins) is a biological database and web resource of known and predicted protein–protein interactions.

• In Cytoscape, go to Control Panel and locate and select the Network Tab in the Control Panel
• Make sure ‘STRING protein query’ is selected in the Network Search Bar.
• Type CTNNB1 in the search field.
• Click the “Search Network” button
• Explore the network!

–

Congratulations! You have reached the end of the integrated assignment.