<jupyter_start><jupyter_text>Introduction to pathway analysis in metabolomicsThis pathway analysis tutorial forms part of the Hands-on Data Analysis for Metabolic Profiling course at Imperial College London.The topics that are covered in this tutorial are:- Compound identifier mapping- Over-representation analysis in metabolomics- Gene set enrichment analysis in metabolomics- Pathway database choice- Pathway network visualisationYou should not focus on understanding every line of code in this tutorial, but rather it should help provide an overview of the key concepts in metabolomics pathway analysis.Please make sure you run each cell in order to ensure the code works correctly.<jupyter_code># Import the required packages
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import gseapy
from sklearn.decomposition import PCA

# Import the pathway analysis functions
# This file contains custom functions to run t-tests and ORA - feel free to check it out and re-use the code!
import pathwayfunctions<jupyter_output><empty_output><jupyter_text>Metabolomics data and pathway importAll pathway analysis methods require at least a dataset and a set of pathways. First we are going to import the KEGG human pathways and the processed metabolomics dataset we will be working with. It is important to ensure the data has been processed prior to pathway analysis. This includes:1. Sample normalisation (e.g. Probabilistic Quotient Normalisation)2. Missing value imputation3. Log transformation4. Centering and scaling the features to mean = 0 and standard deviation = 15. Outlier detection (e.g. using PCA)For more details about pre-processing metabolomics data, see here: https://www.intechopen.com/books/metabolomics-fundamentals-and-applications/processing-and-visualization-of-metabolomics-data-using-r<jupyter_code># Import pathway database
KEGG_human_pathways = pd.read_csv("KEGG_human_pathways_compounds_R98.csv", dtype=str, index_col=0)
Reactome_human_pathways = pd.read_csv("Reactome_human_pathways_compounds_R76.csv", dtype=str, index_col=0)

# Import procsesed datasets
covid_data = pd.read_csv("Su_COVID_metabolomics_processed.csv", index_col=0)
covid_data_KEGG = pd.read_csv("Su_COVID_metabolomics_processed_KEGG.csv", index_col=0)
covid_data_CHEBI = pd.read_csv("Su_COVID_metabolomics_processed_CHEBI.csv", index_col=0)<jupyter_output><empty_output><jupyter_text>The dataset we will be working with is from **Su et al. "Multi-Omics Resolves a Sharp Disease-State Shift between Mild and Moderate COVID-19".** It is an untargeted mass spectrometry dataset generated by Metabolon containing 133 COVID-19 plasma samples of varying WHO ordinal scale severity and 133 control plasma samples. A total of 1050 compounds were profiled.<NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, et al. Multi-Omics Resolves a Sharp Disease-State Shift between Mild and Moderate COVID-19. Cell. 2020;183: 1479-1495.e20. [doi:10.1016/j.cell.2020.10.037](doi:10.1016/j.cell.2020.10.037)<jupyter_code># Inspect the processed data abundance matrix
# Note the Group and WHO_status columns at the far right of the matrix
covid_data.head(6)<jupyter_output><empty_output><jupyter_text>Let's perform PCA to visualise the differences between the COVID and control sample groups.<jupyter_code>sns.set_style("ticks")

# Perform a 5-component PCA
pca = PCA(n_components=5).fit_transform(covid_data.iloc[:, :-2])
pca_df = pd.DataFrame(pca)
pca_df.columns = ["PC " + str(i) for i in range(1, 6)]
pca_df["Group"] = covid_data["Group"].values

# Use the pairplot function to visualise the transformed data
sns.pairplot(pca_df.iloc[:, np.r_[0:4, 5]],
             hue="Group",
             diag_kind="auto")
plt.show()<jupyter_output><empty_output><jupyter_text>The PCA pairplot shows the transformed data from different principal components (PCs) plotted against each other. In this dataset, the strongest separation between the groups is captured by PC1 and PC2.  Where can I download pathways?For the purposes of this tutorial, we will be using KEGG and Reactome pathway databases provided as files in this folder. It is important to use as up to date resources as possible when running any type of pathway analysis. If you wish to download these pathway collections yourself you can do so by:**KEGG**:- Non-programmatic access: Purchasing a KEGG license and downloading pathways via the KEGG FTP site (https://www.kegg.jp/kegg/download/)- Programmatic access: Downloading KEGG pathways via the KEGG API (https://www.kegg.jp/kegg/rest/keggapi.html)**Reactome**:- Non-programmatic access: Downloading the Reactome pathways via the Reactome downloads site (https://reactome.org/download-data). You will need the https://reactome.org/download/current/ChEBI2Reactome_All_Levels.txt file.You may need to write additional code to properly parse the files into a useable format. The files we will use in this tutorial have already been parsed in the desired format.<jupyter_code># Let's inspect the KEGG pathways file:
KEGG_human_pathways.head()<jupyter_output><empty_output><jupyter_text>This CSV file contains:- **Index column**: contians the KEGG pathway identifiers- **Pathway name column**: contains the KEGG pathway names- **Compound columns**: Each row of the file contains columns each containing the KEGG compound identifiers for each pathway  Compound identifier conversionThe first step in the pathway analysis workflow is to convert the compound names or identifiers to the identifier type of the pathway database used. In this example, we are using the KEGG database and will be converting componund names to KEGG compound identifiers.The MetaboAnalyst name conversion tool allows you to convert between several metabolite identifier types and names.https://www.metaboanalyst.ca/MetaboAnalyst/upload/ConvertView.xhtml. If you wish to do this in a programmatic manner, you can use the MetaboAnalyst R package. Have a go at converting the compound names of the Su_compound_names.txt file to KEGG ID's. Inspect the results. As there are 900+ compounds profiled in this dataset, just look at a few examples. The full converted dataset is provided below.What percentage of the compounds in the dataset are you able to map? Which class(es) of compounds are particularly difficult to map?For the purposes of this tutorial, the compound names have already been converted to KEGG and ChEBI identifiers in the covid_data_KEGG and covid_data_CHEBI matrices respectively:<jupyter_code>covid_data_KEGG.head(5)<jupyter_output><empty_output><jupyter_text>The column identifiers beginning with "C0" represent KEGG compounds  Over-representation analysis (ORA)Over-representation analysis is one of the commonest pathway analysis methods in metabolomics. It tests the null hypothesis that a subset of compounds of interest (e.g. those differentially abundant between study groups) are no more present in a certain biological pathway than background compounds. We can build a 2-by-2 contingency table such as the one below, which differentiates between the number of compounds in the dataset which are in a pathway, those that are differentially abundant, and vice versa.       A test statistic and P-value can be determined for each pathway using this contingency table by using a **right-sided Fisher's exact test**. \begin{equation}p = 1 - \sum_{i=0}^{k-1} \frac{\binom{M}{i}\binom{N-M}{n-i}}{\binom{N}{n}}\end{equation}The essential inputs of ORA are:- **Background set** (N)- **List of metabolites of interest** (n)- **Pathway set** (M represents one pathway)The Venn diagram depicts ORA, in which we are testing the significance of the overlap between the pathway and the compounds of interest, **k**.   ORA step 1: Determining the list of compounds of interestAs input to ORA, we require a list of compounds of interest, which are usually metabolites that are differentially abundant between study groups.We will obtain this list by using a series of t-tests. The function below computes a two-sided t-test for each metabolite in the dataset, testing for a significance in metabolite abundance levels between the COVID and healthy groups. The P-values are then adjusted for multiple-correction testing using the Benjamini-Hochberg FDR method.The pathwayfunctions.t_tests() function takes 3 parameters:- **Abundance matrix**. This is either covid_data_KEGG or covid_data_CHEBI. - **Sample group/classes column**. In the COVID dataframe, this column is called "Group"- **Multiple testing correction method** (string) e.g. "fdr_bh"<jupyter_code># Use the t-test function to compute t-tests 
# Note that for the abundance matrix parameter, we use covid_data_KEGG.iloc[:, :-2] to provide the abundance matrix
# with the final two columns removed (WHO_status and Group) as these are string columns and not metabolites
t_test_res = pathwayfunctions.t_tests(covid_data_KEGG.iloc[:, :-2], covid_data_KEGG["Group"], "fdr_bh")

# Select differentially abundant metabolites by filtering on the adjusted P-value at ≤ 0.05
DA_metabolites = t_test_res[t_test_res["P-adjust"] <= 0.05]["Metabolite"].tolist()
print(t_test_res)
print("There are", len(DA_metabolites), "differentially abundant metabolites in this dataset")<jupyter_output><empty_output><jupyter_text>Background setThe other essential component of ORA is the background set, or metabolite "universe". Our assay-specific background set consists of all the compounds profiled in our dataset:<jupyter_code># Select all the columns in the data matrix apart from the last two (metadata) columns (Group and WHO_status)
assay_specific_background_set = covid_data_KEGG.iloc[:, :-2].columns.tolist()
print("The size of the background set is:", len(assay_specific_background_set), "compounds")<jupyter_output><empty_output><jupyter_text>ORA step 2: Running ORANow we are ready to perform ORA with KEGG pathways. The function (code in the pathwayfunctions.py file) loops through each pathway in the pathway database and performs a one-sided Fisher's exact test to compute a test statistic and P-value for the enrichment of each pathway.The pathwayfunctions.over_representation_analysis() function takes 3 parameters:- **List of differentially abundant metabolites** (the compound identifiers must match those in the pathways)- **Background set**: list of compound identifiers forming the background set- **Pathways dataframe**: pandas DataFrame containing pathway identifiers and their respective compounds<jupyter_code>ora_res = pathwayfunctions.over_representation_analysis(DA_metabolites, assay_specific_background_set, KEGG_human_pathways)<jupyter_output><empty_output><jupyter_text>Inspect the ORA results. The 'hits' column shows the fraction of DA metabolites present in the pathway (note the pathway is only composed of compounds which are also present in the background set). The 'coverage' column shows the fraction of pathway metabolites in the background list based on the total number of metabolites in a specific KEGG pathway.<jupyter_code># Look at the top 10 pathways ranked by P-value
ora_res.sort_values(by="P-value").head(10)<jupyter_output><empty_output><jupyter_text>What do you think about the level of pathway coverage in this dataset? Would you recommend pathway analysis for datasets with fewer than 100 metabolites profiled? Let's plot the top 10 pathways ordered by P-value:<jupyter_code>plt.style.use("seaborn")

# we can use a colourmap to colour the bars corresponding to the P-values
cmap = sns.color_palette("mako_r", as_cmap=True)
norm = plt.Normalize(vmin=0, vmax=0.4)

# plot the top 10 pathways
top_10_ora_res = ora_res.sort_values(by="P-value").iloc[0:10, :]
plt.barh(top_10_ora_res["Pathway_name"], top_10_ora_res["P-value"], color=cmap(norm(top_10_ora_res["P-value"])))
plt.xlabel("P-value")
plt.ylabel("Pathway name")

# plot a vertical line to represent P ≤ 0.05
plt.axvline(0.05, color="black", linestyle='--')
plt.title("Top 10 enriched pathways (assay-specific background set)")
plt.show()<jupyter_output><empty_output><jupyter_text>What happens if we use a non-assay specific background set instead?This is a common pitfall of metabolomics pathway analysis, and the default background set often consists of all compounds in a pathway database, for example.Here we will use all KEGG compounds as the non-assay specific background set:<jupyter_code># create the non-specific background set by selecting all the unique compounds mapped to KEGG pathways
non_specific_background_set = list(set([x for x in KEGG_human_pathways.iloc[:, 1:].values.flatten() if x is not np.nan]))<jupyter_output><empty_output><jupyter_text>Now we will perform ORA again with the non-specific background set<jupyter_code># perform ORA using non-specific background set
ora_res_non_specific_background = pathwayfunctions.over_representation_analysis(DA_metabolites, non_specific_background_set, KEGG_human_pathways)<jupyter_output><empty_output><jupyter_text>Plot the results of ORA using the assay-specific vs. non-assay-specific background set<jupyter_code># plotting code
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))
plt.style.use("seaborn")
top_10_ora_res_non_specific = ora_res_non_specific_background.sort_values(by="P-value").iloc[0:10, :]
ax1.barh(top_10_ora_res_non_specific["Pathway_name"], top_10_ora_res_non_specific["P-value"])
ax1.set_xlabel("P-value")
ax1.set_ylabel("Pathway name")
ax1.set_title("Top 10 enriched pathways \n (non-specific background set)")

# adding labels for the p-values 
for index, data in enumerate(top_10_ora_res_non_specific["P-value"]):
    ax1.text(data+0.000001, index, s=round(top_10_ora_res_non_specific["P-value"].tolist()[index], 6))

# Ensures pathways are the same in both results (non specific background set can have additional pathways)
intersect = (set(ora_res["Pathway_ID"].tolist()) & set(ora_res_non_specific_background["Pathway_ID"].tolist()))
ora_res_all = ora_res_non_specific_background.loc[ora_res_non_specific_background["Pathway_ID"].isin(intersect)]

# convert p-values to -log10 scale
ora_res_pvals = np.negative(np.log10(ora_res["P-value"].tolist()))
ora_res_all_pvals = np.negative(np.log10(ora_res_all["P-value"].tolist()))

sns.regplot(x=ora_res_pvals, y=ora_res_all_pvals, ax=ax2, line_kws={"color": "red"})
ax2.set_ylabel("-log10 P-value (Non-assay specific background set)")
ax2.set_xlabel("-log10 P-value (Assay-specific background set)")
ax2.set_title("Comparison of P-values")
ax2.axhline(y=-np.log10(0.05), linewidth=1, color='black', linestyle='--')
ax2.axvline(x=-np.log10(0.05), linewidth=1, color='black', linestyle='--')
ax2.set(ylim=(0, 10), xlim=(0, 5))

plt.show()
# Look at the results of the top 10 enriched pathways using the non-assay-specific background set. 
top_10_ora_res_non_specific<jupyter_output><empty_output><jupyter_text>What do you notice about the P-values obtained using the assay-specific background set as opposed to the non-assay specific background set containing all KEGG compounds?  Comparing pathway databasesHow do pathway analysis results differ when using different pathway databases? Let's compare ORA results obtained using KEGG and Reactome human pathway databases:<jupyter_code># Re-run and save the ORA KEGG results for comparison 
ora_KEGG = pathwayfunctions.over_representation_analysis(DA_metabolites,
                                                         assay_specific_background_set,
                                                         KEGG_human_pathways)<jupyter_output><empty_output><jupyter_text>First we need to determine differential metabolites for the dataset with Reactome identifiers (ChEBI). We perform a series of t-tests on the metabolites and keep only those with an adjusted P-value ≤ 0.05.<jupyter_code>t_test_res_reactome = pathwayfunctions.t_tests(covid_data_CHEBI.iloc[:, :-2], covid_data_CHEBI["Group"], "fdr_bh")
DA_metabolites_reactome = t_test_res_reactome[t_test_res_reactome["P-adjust"] <= 0.05]["Metabolite"].tolist()
print(len(DA_metabolites_reactome), "differential metabolites using Reactome with a total background set size of:", len(covid_data_CHEBI.iloc[:, :-2].columns.tolist()))<jupyter_output><empty_output><jupyter_text>Now compute ORA with the Reactome pathways. We first take the columns from the "covid_data_CHEBI" file to use as the Reactome assay-specific background set.<jupyter_code>reactome_assay_specific_bgset = covid_data_CHEBI.iloc[:, :-2].columns.tolist()
ora_Reactome = pathwayfunctions.over_representation_analysis(DA_metabolites_reactome, reactome_assay_specific_bgset, Reactome_human_pathways)<jupyter_output><empty_output><jupyter_text>Plot the ORA results from using KEGG vs. Reactome pathways:<jupyter_code>fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6), sharex=True)
plt.style.use("seaborn")
top_10_ora_res_kegg = ora_KEGG.sort_values(by="P-value").iloc[0:10, :]
ax1.barh(top_10_ora_res_kegg["Pathway_name"], top_10_ora_res_kegg["P-value"], color="tab:cyan")
ax1.set_xlabel("P-value")
ax1.set_ylabel("Pathway name")
ax1.set_title("Top 10 enriched pathways (KEGG)")
ax1.axvline(0.05, color="black", linestyle='--')

top_10_ora_res_reactome = ora_Reactome.sort_values(by="P-value").iloc[0:10, :]
ax2.barh(top_10_ora_res_reactome["Pathway_name"], top_10_ora_res_reactome["P-value"], color="tab:red")
ax2.set_xlabel("P-value")
ax2.set_ylabel("Pathway name")
ax2.set_title("Top 10 enriched pathways (Reactome)")
plt.tight_layout()
ax2.axvline(0.05, color="black", linestyle='--')

plt.show()<jupyter_output><empty_output><jupyter_text>Pathway databases are always evolving, so an important consideration is to use the lastest version of a database when running your analysis.   Your turn! 📚 Use ORA to find pathways enriched in severe vs. mild COVIDWe previously looked at pathway analysis to find out which pathways differ between COVID and non-COVID patients. Now it's your turn to use the code presented in this notebook to run an ORA analysis.  First, the metabolite abundance matrix is filtered only for those patients with WHO status 1-2 (mild) or 5-7 (severe):<jupyter_code>covid_KEGG_mild_severe = covid_data_KEGG[covid_data_KEGG["WHO_status"].isin(["1-2", "5-7"])]
# There are 45 mild covid patients and 28 severe covid patients
# This bit of code is dropping the old group column (COVID vs. non-COVID) and adding a new column (mild vs. severe)
covid_KEGG_mild_severe = covid_KEGG_mild_severe.drop("Group", axis=1)
covid_KEGG_mild_severe["Group"] = ["Severe" if i == "5-7" else "Mild" for i in covid_KEGG_mild_severe["WHO_status"]]
# This is doing the same thing as above - filtering the samples - but using the Reactome compound identifiers instead
covid_Reactome_mild_severe = covid_data_CHEBI[covid_data_CHEBI["WHO_status"].isin(["1-2", "5-7"])]
covid_Reactome_mild_severe = covid_Reactome_mild_severe.drop("Group", axis=1)
covid_Reactome_mild_severe["Group"] = ["Severe" if i == "5-7" else "Mild" for i in covid_Reactome_mild_severe["WHO_status"]]
covid_KEGG_mild_severe.head()<jupyter_output><empty_output><jupyter_text>Complete the code in the following cells. You will need to use the functions:pathwayfunctions.t_tests(abundance_matrix, group_column, multiple_correction_method)pathwayfunctions.over_representation_analysis(DA_metabolite_list, background_list, pathway_set)Your code will need to:1. Use the t-test function to run t-tests on all the metabolites2. Filter the t-test results for only those metabolites with adjusted P-value ≤ a threshold3. Obtain the background set for the dataset, taking the column names (=metabolites) present in the dataset matrix4. Run ORA<jupyter_code># STEP ONE
# Use the t-test function to compute t-tests 
# Use "fdr_bh" (Benjamini-Hochberg method) as the parameter for multiple testing correction for the t-test function
t_testres_covid_severity = 

# Select differentially abundant metabolites by filtering on the adjusted P-value at ≤ 0.05
DA_metabolites_covid_severity = 

print(len(DA_metabolites_covid_severity))

# Obtain the background set 
assay_specific_background_set_covid_severity = 
# STEP TWO
# Run ORA using the 
# - list of DA metabolites
# - background list
# - pathways: KEGG (KEGG_human_pathways) or Reactome (Reactome_human_pathways)
ora_res_covid_severity =<jupyter_output><empty_output><jupyter_text>If you are stuck, try to re-use the code snippets from the previous ORA analysis above, changing the variable names to those relating to the filtered abundance matrix (containing only severe and mild COVID patients): **covid_KEGG_mild_severe**<jupyter_code># Use the pandas .sort_values() method to sort the ORA results dataframe by pathway p-values
ora_res_covid_severity.sort_values(by="P-value")<jupyter_output><empty_output><jupyter_text>Run the following cell to plot your results<jupyter_code># Plot the top 10 enriched pathways
plt.style.use("seaborn")

# we can use a colourmap to colour the bars corresponding to the P-values
cmap = sns.color_palette("YlOrBr_r", as_cmap=True)
norm = plt.Normalize(vmin=0, vmax=0.4)

# plot the top 10 pathways
top_10_ora_res = ora_res_covid_severity.sort_values(by="P-value").iloc[0:10, :]
plt.barh(top_10_ora_res["Pathway_name"], top_10_ora_res["P-value"], color=cmap(norm(top_10_ora_res["P-value"])))
plt.xlabel("P-value")
plt.ylabel("Pathway name")

# plot a vertical line to represent P ≤ 0.05
plt.axvline(0.05, color="black", linestyle='--')
plt.title("Top 10 enriched pathways in severe vs. mild COVID patients")
plt.show()<jupyter_output><empty_output><jupyter_text>- What are the top 5 enriched pathways? Do they differ from the top 5 enriched in COVID vs. non-COVID patients?- Try the analysis with Reactome vs. KEGG pathways. How do the results differ? Which would you prefer to use? The answers for this exercise can be found in the "ORA_exercise_answers.txt" file.  m(GSEA)Over-representation analysis is a fast and straightforward approach to pathway analysis, but has a few **limitations**. These include:- Having to pre-select metabolites of interest (based on a hard threshold such as P-value cutoff)- Not taking into account the level of differences in these metabolites between study groups (e.g. the fold-change)Functional class scoring approaches such as gene set enrichment analysis (GSEA)* have been developed to overcome these challenges. They take all measured metabolites as input, which allows small but coordinated changes in the metabolites present in a pathway to be detected. Here we will use gseapy to perform this. gseapy can be adapted for metabolomics data simply by using metabolite data and pathways as input. The statistical test and assumptions remain the same. First we will create a ranked list of metabolites (based on differences in abundance between COVID and control groups) as input using the t-test function above that we also used to determine the list of metabolites of interest for ORA. We will rank the metabolites based on the t-test statistic. Other common ranking metrics include fold change, signal-to-noise ratio, or P-value and fold change multiplied together (signed fold-change). *<NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102: 15545–15550. doi:10.1073/pnas.0506580102  m(GSEA) step 1: Creating the ranked list of metabolitesWe use the same t-test function as in ORA, but instead of using the results to filter the adjusted P-values on, we use the *t-statistic* for all the metabolites in the dataset to rank them.<jupyter_code>t_test_res = pathwayfunctions.t_tests(covid_data_KEGG.iloc[:, :-2], covid_data_KEGG["Group"], "fdr_bh", return_stat=True)
ranked_list = t_test_res[["Metabolite", "t-statistic"]]
ranked_list.sort_values(by="t-statistic")<jupyter_output><empty_output><jupyter_text>m(GSEA) step 2: Running GSEA Next we will run m(GSEA) using the gseapy package. The inputs for this are:- The **ranked list** we just created (in the form of a pandas series)- A **gmt format file containing the pathways** (see format details here: https://software.broadinstitute.org/cancer/software/gsea/wiki/index.php/Data_formats)- The minimum size for a pathway to be included in the analysis (optional parameter)The outputs are stored in the GSEA_Prerank folder The GMT file is very similar to the pathway files we used for ORA. The index column containing the pathway identifiers simply becomes the first column:<jupyter_code>pd.read_csv("KEGG_human_pathways_compounds_R98.gmt", sep="\t", header=None).head(3)
# Run GSEA prerank using our ranked list of metabolites
mgsea_res = gseapy.prerank(rnk=ranked_list,
               gene_sets='KEGG_human_pathways_compounds_R98.gmt',
               min_size=3, # change this as you wish
               verbose=True,
               outdir=None) # specify an output directory path if you want to download the full results<jupyter_output><empty_output><jupyter_text>Let's explore the m(GSEA) results:<jupyter_code># sort pathways by lowest FDR
mgsea_res_df = mgsea_res.res2d.sort_values(by="fdr").head(10)

# add a column with the full pathway names to the dataframe
mgsea_res_df["pathway_name"] = mgsea_res_df.index.map(dict(zip(KEGG_human_pathways.index, KEGG_human_pathways["Pathway_name"])))
col = mgsea_res_df.pop("pathway_name")
mgsea_res_df.insert(0, "pathway_name", col)
mgsea_res_df.head(10)<jupyter_output><empty_output><jupyter_text>What do the extra columns in this dataframe mean? The following terms are specific to GSEA-based analyses:- **"es", Enrichment score**: GSEA enrichment score calculated using a weighted Kolmogorov-Smirnov statistic- **"nes", Normalised enrichment score**: Enrichment score normalised by the pathway size- **"ledge_genes", Leading edge genes**: These are the genes (or metabolites) that contributed to the enrichment of the pathway Now we will create some plots to examine the mGSEA results. We can plot the standard GSEA enrichment plot using the gseaplot function:<jupyter_code>from gseapy.plot import gseaplot

# to save your figure, make sure that ofname corresponds to a valid file path on your computer

# Here we are plotting the pathway with the lowest FDR, "Sphingolipid signaling pathway - Homo sapiens (human)"
gseaplot(rank_metric=mgsea_res.ranking,
         term="Sphingolipid signaling pathway - Homo sapiens (human)",
         **mgsea_res.results["hsa04071"],
         ofname=None)<jupyter_output><empty_output><jupyter_text>What does the m(GSEA) plot show?The m(GSEA) algorithm walks down the rank-ordered list of metabolites, each time adding a value if a metabolite in the pathway of interest is encountered (the magnitude of this value is based on the correlation of the metabolite with the phenotype (i.e. our t-statistic)), or subtracting a value if a metabolite encountered is not in the pathway of interest. - The enrichment score is recorded as the maximum deivation from 0 (i.e. the lowest point on the green curve). In this example the enrichment score is negative. As our ranking metric (t-test statistic) compares the means of the COVID samples relative to the controls, this means that a negative enrichment score suggests that the sphingolipid pathway is enriched in COVID patients. - The blue bars at the far right in the middle of the plot indicate where in the ranked metabolite list the metabolites in the pathway appear. If the bars cluster towards the right of the plot, this indicates they appear at the bottom of the ranked metabolite list and the pathway is negatively enriched. Alternatively, if the bars cluster at the left side of the plot, this means the metabolites appear at the top of the ranked list and the pathway is enriched. - The bottom part of the plot (ranked list metric) shows how the ranking metric changes as the algorithm walks down the list of genes. A positive ranking metric indicates correlation with the first phenotype (controls), whereas a negative ranking metric indicates a correlation with the second phenotype (COVID).- The leading edge metabolites (or genes) are those metabolites in the pathway being tested that appear before the maximum deviation from 0 in the case of a positive enrichment score, or after in the case of a negative enrichment score. These are usually the most interesting metabolites in a pathway, as they contribute strongly to the enrichment of the pathway. Let's compare the top 10 mGSEA results to ORA:<jupyter_code># plotting code
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6), sharex=True)
plt.style.use("seaborn")
top_10_ora_res = ora_res.sort_values(by="P-value").iloc[0:10, :]
ax1.barh(top_10_ora_res["Pathway_name"], top_10_ora_res["P-value"], color="tab:purple")
ax1.set_xlabel("P-value")
ax1.set_ylabel("Pathway name")
ax1.set_title("Top 10 enriched pathways (ORA)")
ax1.axvline(0.05, color="black", linestyle='--')

top_10_GSEA_res = mgsea_res_df.sort_values(by="pval").iloc[0:10, :]
ax2.barh(top_10_GSEA_res["pathway_name"], top_10_GSEA_res["pval"], color="tab:green")
ax2.set_xlabel("P-value")
ax2.set_ylabel("Pathway name")
ax2.set_title("Top 10 enriched pathways (GSEA)")
plt.tight_layout()
ax2.axvline(0.05, color="black", linestyle='--')

plt.show()<jupyter_output><empty_output><jupyter_text>Pathway visualisation  Annotating the KEGG metabolic network with iPath3 iPath3 is an open-source web application for pathway visualisation (https://pathways.embl.de/). Go to the iPath web server and select "Metabolism" from the main menu:  Paste the following information into the "Element selection" box: C00001 0000ff W20 0.5map00010 ff0000 W10R01148 800080 W10 After submitting the data you should see three separate entities highlighted including a compound, a pathway (glycolysis) and a reaction. - The first element corresponds to the KEGG compound/reaction/pathway ID- The second element is a hex colour code- The third element is the width of the node- The fourth element ranges from 0-1 and corresponds to the transparency of the node Let's try another submission with more compounds. This time, we will use the differentially abundant metabolites from the COVID dataset that we used as input for ORA, and visualise them on the KEGG metabolic network.<jupyter_code># Use the t-test function to compute t-tests 
t_test_res = pathwayfunctions.t_tests(covid_data_KEGG.iloc[:, :-2], covid_data_KEGG["Group"], "fdr_bh")

# Select differentially abundant metabolites by filtering on the adjusted P-value at ≤ 0.05
DA_metabolites = t_test_res[t_test_res["P-adjust"] <= 0.05]["Metabolite"].tolist()<jupyter_output><empty_output><jupyter_text>Now we will put the DA compunds into a **DataFrame**, along with specifications for the node width and colour. Putting these details into a DataFrame will make it easier to paste them all into iPath without having to write each individual element into the text box.<jupyter_code>ipath_params = pd.DataFrame(DA_metabolites)
ipath_params["Colour"] = "#ff0000"
ipath_params["Node width"] = "W20"
ipath_params["Transparency"] = 0.6
ipath_params.head()<jupyter_output><empty_output><jupyter_text>This is what our output dataframe looks like. We will now save it as a .tsv file and copy the data from here into iPath.<jupyter_code>ipath_params.to_csv("COVID_DA_metabolites_ipath_input.tsv", sep="\t")<jupyter_output><empty_output><jupyter_text>Try showing 'Advanced options' and experiment with some of the parameters:- Select 'Keep original colours'- Try changing the default node radius, opacity, and colours- Experiment with the species filter. E.g. type "HSA" for human pathways, or "SCE" for Saccaromyces cerevisiaeYou can export annotated pathway maps created using iPath to an .svg format, and then import these into InkScape or Illustrator to further edit the layers in detail and make publication-quality figures    Additional exercise: Using BioPython to visualise pathways Please note, you will need the BioPython and [ReportLab](https://pypi.org/project/reportlab/) packages to run this exercise.<jupyter_code>import Bio
from Bio.KEGG.REST import *
from Bio.KEGG.KGML import KGML_parser
from Bio.Graphics.KGML_vis import KGMLCanvas
from IPython.display import Image, HTML<jupyter_output><empty_output><jupyter_text>We can view KEGG pathway diagrams using the Image function:<jupyter_code>Image(kegg_get("hsa04071", "image").read())

# try replacing "hsa04071" with the name of another KEGG pathway of interest<jupyter_output><empty_output><jupyter_text>We can also view additional information about a pathway, such as the number of genes and compounds it contains:<jupyter_code>pathway = KGML_parser.read(kegg_get("hsa04071", "kgml"))
print(pathway)<jupyter_output><empty_output><jupyter_text>The BioPython KGML package allows you to parse KEGG pathway diagrams and add additional features to them. For example, we can highlight the differentially abundant metabolites (or genes) present in a pathway:<jupyter_code># If you are getting a 'URLError' please uncomment the following two lines and run the cell again!
# import ssl
# ssl._create_default_https_context = ssl._create_unverified_context

# hsa04071 is the sphingolipid signaling pathway
pathway = KGML_parser.read(kegg_get("hsa04071", "kgml"))
pathway_compounds = [cpd.name for cpd in pathway.compounds]

# We will use the list of differentially abundant compounds we used to perform ORA
DA_metabolites = t_test_res[t_test_res["P-adjust"] < 0.05]["Metabolite"].tolist()
DA_metabolites_kegg_format = ["cpd:"+i for i in DA_metabolites]

# and find which of these compounds are present in the sphingolipid signaling pathway
DA_metabolites_in_pathway = list(set(pathway_compounds) & set(DA_metabolites_kegg_format))
print("There are", len(DA_metabolites_in_pathway), "differentially abundant compounds in the pathway")

# we are highlighting all differentially abundant compounds in red circles
for size, element in zip(range(8, 8+len(pathway.compounds)), pathway.compounds):
    if element.name in DA_metabolites_in_pathway:
        for graphic in element.graphics:
            graphic.width = 20
            graphic.height = 20
            graphic.fgcolor = '#FF0000'
            graphic.bgcolor = '#FF0000'

canvas = KGMLCanvas(pathway,
                    import_imagemap=True,
                    show_genes=False,
                    label_compounds=False,
                    fontsize=14,
                    show_orthologs=False,
                    show_reaction_entries=False,
                    label_maps=False,
                    show_maps=False,
                    draw_relations=True)
# saves the pathway diagram as a pdf in the same folder as the notebook
canvas.draw("pathway_highlighted.pdf")<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup