CoRe Metabolomics

Christina Schmidt

Heidelberg Universiy

Dimitrios Prymidis

Cologne Universiy

A Consumption-Release (CoRe) metabolomics experiment usually refers to a cell culture experiment where metabolomics is performed on the cell culture media.

In this tutorial we showcase how to use MetaProViz:

• To process raw peak data and identify outliers.
  
• To perform differential metabolite analysis (DMA) to generate Log2Distance and statistics and perform pathway analysis using Over Representation Analysis (ORA) on the results.
  
• To do metabolite clustering analysis (MCA) to find clusters of metabolites with similar behaviors and perform pathway analysis using ORA on each cluster.
  
• To use specific visualizations to aid biological interpretation of the results.
  
  
  First if you have not done yet, install the required dependencies and load the libraries:

# 1. Install Rtools if you haven’t done this yet, using the appropriate version (e.g.windows or macOS). 
# 2. Install the latest development version from GitHub using devtools
#devtools::install_github("https://github.com/saezlab/MetaProViz")

library(MetaProViz)

#dependencies
library(tidyverse)

#Please install the Biocmanager Dependencies:
#BiocManager::install("clusterProfiler")
#BiocManager::install("EnhancedVolcano")

1. Loading the example data

Here we choose an example datasets, which is publicly available on metabolomics workbench project PR001418 including metabolic profiles of human renal epithelial cells HK2 and cell renal cell carcinoma (ccRCC) cell lines cultured in Plasmax cell culture media. Here we use the integrated raw peak data as example data using the trivial metabolite name in combination with the KEGG ID as the metabolite identifiers.

As part of the MetaProViz package you can load the example data into your global environment using the function toy_data():

1. CoRe experiment (CoRe)
The raw data are available via metabolomics workbench study ST002226 were exometabolomics of HK2 and ccRCC cell lines 786-O, 786-M1A, 786-M2A, OS-RC-2, OS-LM1 and RFX-631 were performed.

MetaProViz::ToyData(data="CoRe")

Preview of the DF CoRe including columns with sample information and metabolite ids with their measured values.

	Conditions	Biological_Replicates	GrowthFactor	valine-d8	hipppuric acid-d5	2-hydroxyglutarate	2-ketoglutarate	3-Dehydro-L-threonate
MS51-06	HK2	1	249.2817	780552871	3127630257	84950547	169244159	1807489245
MS51-07	HK2	2	249.2817	802602348	3256031922	60753859	151064767	695228424
MS51-08	HK2	3	249.2817	831984796	3308009345	73718363	171281531	791442407
MS51-09	HK2	4	249.2817	822744518	3209731571	65933166	112033043	315209589
MS51-10	HK2	5	249.2817	805565867	3297793480	68183576	170902744	615035217
MS51-11	786-O	1	297.3423	841873509	3418515398	75941661	215553005	501977089
MS51-12	786-O	2	297.3423	825462965	3218049751	62100210	195308040	484681099

2. Additional information mapping the trivial metabolite names to KEGG IDs and selected pathways (MappingInfo)

MetaProViz::ToyData(data="MappingInfo")

Preview of the DF Pathways including the trivial metabolite identifiers used in the experiment as well as KEGG IDs and pathway information.

	HMDB	KEGG.ID	KEGGCompound	Pathway
N-acetylaspartate	HMDB0000812	C01042	N-Acetyl-L-aspartate	Alanine, aspartate and glutamate metabolism
argininosuccinate	HMDB0000052	C03406	N-(L-Arginino)succinate	Alanine, aspartate and glutamate metabolism
N-acetylaspartylglutamate	HMDB0001067	C12270	N-Acetylaspartylglutamate	Alanine, aspartate and glutamate metabolism
tyrosine	HMDB0000158	C00082	L-Tyrosine	Amino acid metabolism
asparagine	HMDB0000168	C00152	L-Asparagine	Amino acid metabolism

3. KEGG pathways that are loaded via KEGG API using the package KEGGREST and can be used to perform pathway analysis. (KEGG_Pathways)

#This will use KEGGREST to query the KEGG API to load the pathways:
MetaProViz::LoadKEGG()
#> Cached file loaded from: C:\Users\chris\AppData\Local/Cache/KEGG_Metabolite.rds

Preview of the DF KEGG_Pathways.

term	Metabolite	MetaboliteID	Description
Glycolysis / Gluconeogenesis - Homo sapiens (human)	Pyruvate	C00022	Glycolysis / Gluconeogenesis - Homo sapiens (human)
Glycolysis / Gluconeogenesis - Homo sapiens (human)	Acetyl-CoA	C00024	Glycolysis / Gluconeogenesis - Homo sapiens (human)
Glycolysis / Gluconeogenesis - Homo sapiens (human)	D-Glucose	C00031	Glycolysis / Gluconeogenesis - Homo sapiens (human)
Glycolysis / Gluconeogenesis - Homo sapiens (human)	Acetate	C00033	Glycolysis / Gluconeogenesis - Homo sapiens (human)
Glycolysis / Gluconeogenesis - Homo sapiens (human)	Oxaloacetate	C00036	Glycolysis / Gluconeogenesis - Homo sapiens (human)

2. Run MetaProViz Analysis

Currently, MetaProViz contains four different modules, which include different methods and can be used independently from each other or in combination (see introduction for more details). Here we will go trough each of those modules and apply them to the example data.

Pre-processing

MetaProViz includes a pre-processing module with the function Preprocessing() that has multiple parameters to perform customize data processing.
Feature_Filtering applies the 80%-filtering rule on the metabolite features either on the whole dataset (=“Standard”) (Bijlsma et al. 2006) or per condition (=“Modified”) (Wei et al. 2018). This means that metabolites are removed were more than 20% of the samples (all or per condition) have no detection. In case of the CoRe experiment, the blank samples are ignored during feature filtering, since often metabolites are released from a cell and not naturally present in the culture media leading to no detection in the blank. With the parameter Feature_Filt_Value we enable the adaptation of the stringency of the filtering based on the experimental context. For instance, patient tumour samples can contain many unknown subgroups due to gender, age, stage etc., which leads to a metabolite being detected in only 50% (or even less) of the tumour samples, hence in this context it could be considered to change the Feature_Filt_Value from the default (=0.8). If Feature_Filtering = "None", no feature filtering is performed. In the context of Feature_Filtering it is also noteworthy that the function Pool_Estimation() can be used to estimate the quality of the metabolite detection and will return a list of metabolites that are variable across the different pool measurements (pool = mixture of all experimental samples measured several times during the LC-MS run) . Variable metabolite in the pool sample should be removed from the data.
The parameter TIC_Normalization refers to Total Ion Count (TIC) normalisation, which is often used with LC-MS derived metabolomics data. If TIC_Normalization = TRUE, each feature (=metabolite) in a sample is divided by the sum of all intensity value (= total number of ions) for the sample and finally multiplied by a constant ( = the mean of all samples total number of ions). Noteworthy, TIC normalisation should not be used with small number of features (= metabolites), since TIC assumes that on “average” the ion count of each sample is equal if there were no instrument batch effects (Wulff and Mitchell 2018).
The parameter MVI refers to Missing Value Imputation (MVI) and if MVI = TRUE half minimum (HM) missing value imputation is performed per feature (= per metabolite). Here it is important to mention that HM has been shown to perform well for missing vales that are missing not at random (MNAR) (Wei et al. 2018).
Lastly, the function Preprocessing() performs outlier detection and adds a column “Outliers” into the DF, which can be used to remove outliers. The parameter HotellinsConfidence can be used to choose the confidence interval that should be used for the Hotellins T2 outlier test (Hotelling 1931).

Since our example data contains pool samples, we will do Pool_Estimation() before applying the Preprocessing() function. This is important, since one should remove the features (=metabolites) that are too variable prior to performing any data transformations such as TIC as part of the Preprocessing() function.
It is worth mentioning that the Coefficient of variation (CV) is calculated by dividing the standard deviation (SD) by the mean. Hence CV depends on the SD, which in turn works for normally distributed data.

Pool_Estimation_result<- MetaProViz::PoolEstimation(InputData = Media[,-c(1:3)],
                                                    SettingsFile_Sample = Media[,1:3],
                                                    SettingsInfo = c(PoolSamples = "Pool", Conditions="Conditions"),
                                                    CutoffCV = 100)

Pool_Estimation_result_DF_CV <-Pool_Estimation_result[["DF"]][["CV"]]

#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Bin width defaults to 1/30 of the range of the data. Pick better value with
#> `binwidth`.

Preview of the Pool_Estimation result.

Metabolite	CV	HighVar	MissingValuePercentage
valine-d8	1.744198	FALSE	0
hipppuric acid-d5	1.253983	FALSE	0
2-hydroxyglutarate	7.199141	FALSE	0
2-ketoglutarate	2.766989	FALSE	0
3-Dehydro-L-threonate	9.222780	FALSE	0

The results from the Pool_Estimation() is a table that has the CVs. If there is a high variability, one should consider to remove those features from the data. For the example data nothing needs to be removed. If you have used internal standard in your experiment you should specifically check their CV as this would indicate technical issues (here valine-d8 and hippuric acid-d5).

Now we will apply the Preprocessing() function to the example data and have a look at the output produced. You will notice that all the chosen parameters and results are documented in messages. All the results data tables, the Quality Control (QC) plots and outlier detection plots are returned and can be easily viewed. Importantly, here we are able to specify that we have a CoRe experiment setting the parameter CoRe=TRUE, in which case a few additional data processing steps are applied:
1. Blank sample: This refers to media samples where no cells have been cultured in, which will be used as blank. In detail, the mean of the blank sample of a feature (= metabolite) will be substracted from the values measured in each sample for the same feature. In the column “Condition” of the Experimental_design DF, you will need to label your blank samples with “blank”.
2. Growth factor or growth rate: This refers to the different conditions and is either based on cell count or protein quantification at the start of the experiment (t0) and at the end of the experiment (t1) resulting in the growth factor (t1/t0). Otherwise, one can experimentally estimate the growth rate of each condition. Ultimately, this measure is used to normalize the data, since the amount of growth will impact the consumption and release of metabolites from the media and hence we need to account for this. If you do not have this information, this will be set to 1, yet be aware that this may affect the results.

You can pass these additional information via the parameter Input_SettingsInfo, by passing the column name for the CoRe_norm_factor in the Input_SettingsFile and the condition name for the CoRe_media in the Input_data file.

#Prepare the input:
Media_input <- Media%>%
  subset(!Conditions=="Pool", select = -c(1:3))#remove pool samples and remove the information columns

Media_Metadata <- Media%>%
  subset(!Conditions=="Pool", select = c(1:3))#remove pool samples and keep the information columns only

PreProcessing_res <-  MetaProViz::PreProcessing(InputData=Media_input,
                                                SettingsFile_Sample =Media_Metadata,
                                                SettingsInfo = c(Conditions = "Conditions",
                                                                       Biological_Replicates = "Biological_Replicates",
                                                                       CoRe_norm_factor = "GrowthFactor",
                                                                       CoRe_media = "blank"),
                                                FeatureFilt = "Modified",
                                                FeatureFilt_Value = 0.8,
                                                TIC = TRUE,# As we have raw data we will perform total ion count norm
                                                MVI=TRUE, #We assume the values are not missing at random and perform half minimum MVI
                                                MVI_Percentage=50,
                                                HotellinsConfidence = 0.99,# We perform outlier testing using 0.99 confidence interval
                                                CoRe = TRUE) 

# Now we can have a look at the results table:
Media_Preprocessed <-  PreProcessing_res[["DF"]][["Preprocessing_output"]]

#> For Consumption Release experiment we are using the method from Jain M.  REF: Jain et. al, (2012), Science 336(6084):1040-4, doi: 10.1126/science.1218595.
#> Here we apply the modified 80%-filtering rule that takes the class information (Column `Conditions`) into account, which additionally reduces the effect of missing values. REF: Yang et. al., (2015), doi: 10.3389/fmolb.2015.00004)
#> filtering value selected: 0.8
#> 3 metabolites where removed: N-acetylaspartylglutamate, hypotaurine, S-(2-succinyl)cysteine
#> NA values were found in Control_media samples for metabolites. For metabolites including NAs MVI is performed unless all samples of a metabolite are NA.
#> Metabolites with high NA load (>20%) in Control_media samples are: dihydroorotate.
#> Metabolites with only NAs (=100%) in Control_media samples are: hydroxyphenylpyruvate. Those NAs are set zero as we consider them true zeros
#> Missing value imputation is performed, as a complementary approach to address the missing value problem, where the missing values are imputing using the `half minimum value`. REF: Wei et. al., (2018), Reports, 8, 663, doi:https://doi.org/10.1038/s41598-017-19120-0
#> Total Ion Count (TIC) normalization is used to reduce the variation from non-biological sources, while maintaining the biological variation. REF: Wulff et. al., (2018), Advances in Bioscience and Biotechnology, 9, 339-351, doi:https://doi.org/10.4236/abb.2018.98022
#> CoRe data are normalised by substracting mean (blank) from each sample and multiplying with the CoRe_norm_factor
#> Identification of outlier samples is performed using Hotellin's T2 test to define sample outliers in a mathematical way (Confidence = 0.99 ~ p.val < 0.01) REF: Hotelling, H. (1931), Annals of Mathematical Statistics. 2 (3), 360–378, doi:https://doi.org/10.1214/aoms/1177732979.
#> HotellinsConfidence value selected: 0.99
#> There are possible outlier samples in the data
#> Filtering round 1 Outlier Samples: MS51-06
#> Filtering round 2 Outlier Samples: MS51-09

Preview of the pre-processing results, which has an additional column Outlier including the results of Hotellins T2.

	Conditions	Biological_Replicates	GrowthFactor	Outliers	valine-d8	hipppuric acid-d5	2-hydroxyglutarate	2-ketoglutarate	3-Dehydro-L-threonate
MS51-06	HK2	1	249.2817	Outlier_filtering_round_1	3531526217	-1673750986	5513004267	30614207219	368500681526
MS51-07	HK2	2	249.2817	no	1958125674	1886704843	-1209262732	24617648085	78317130942
MS51-08	HK2	3	249.2817	no	-4532562455	-39893563477	767684868	26771101645	88956838993
MS51-09	HK2	4	249.2817	Outlier_filtering_round_2	-527252000	-38568863508	-530294707	13989640920	-18139183358
MS51-10	HK2	5	249.2817	no	-9262935498	-36751484785	-390772056	26968551969	49440162385

In the output table you can now see the column “Outliers” and for the Condition HK2 CCM, we can see that based on Hotellin’s T2 test, samples were detected as outliers in the first and second round of filtering.
As part of the Preprocessing() function several plots are generated and saved. Additionally, the ggplots are returned into the list to enable further modifiaction using the ggplot syntax. These plots include plots showing the outliers for each filtering round and other QC plots.

As part of the MetaProViz visualization module one can easily further customize the PCA plot and adapt color and shape for the information of interest. You can see more below for the VizPCA() function.
Before we proceed, we will remove the outlier:

Media_Preprocessed <-Media_Preprocessed%>%
  subset(!Outliers=="Outlier_filtering_round_1")

In metabolomics, sometimes samples are injected (=measured) several times, which can be termed as analytical replicates. The MetaProViz pre-processing module includes the function ReplicateSum(), which will summarize those and save the results.

DMA

Differential Metabolite Analysis (DMA) between two conditions (e.g. Tumour versus Healthy) usually calculates the Log2FC, p-value, adjusted p-value and t-value. Yet, in a CoRe experiment the normalized metabolite values can be either a negative value, if the metabolite has been consumed from the media, or a positive value, if the metabolite has been released from the cell into the culture media. Since we can not calculate a Log2FC using negative values, we calculate the absolute difference between the mean of Condition 1 versus the mean of Condition 2. The absolute difference is log2 transformed in order to make the values comparable between the different metabolites, resulting in the Log2Dist. The result doesn’t consider whether one product is larger than the other; it only looks at the magnitude of their difference. To reflect the direction of change between the two conditions we multiply with -1 if C1 < C2. By setting the paramteter CoRe = TRUE, instead of calclulating the Log2FC, the Log2 Distance is calculated.
With the different parameters STAT_pval and STAT_padj one can choose the statistical tests such as t.test, wilcoxon test, limma, annova, kruskal walles, etc. (see function reference for more information).
As input one can use the pre-processed data we have generated using the Preprocessing module, but here one can of course use any DF including metabolite values, even though we recommend to normalize the data and remove outliers prior to DMA. Moreover, we require the Input_SettingsFile_Sample including the sample metadata with information which condition a sample corresponds to. Additionally, we enable the user to provide a Plot_SettingsFile_Metab containing the metadata for the features (metabolites), such as KEGG ID, pathway, retention time, etc.

By defining the numerator and denominator as part of the Input_SettingsInfo parameter, it is defined which comparisons are performed:
1. one_vs_one (single comparison): numerator=“Condition1”, denominator =“Condition2”
2. all_vs_one (multiple comparison): numerator=NULL, denominator =“Condition”
3. all_vs_all (multiple comparison): numerator=NULL, denominator =NULL (=default)

As input we will use the pre-processed data we have generated using the Preprocessing module, but here one can of course use any DF including metabolite values and information about the conditions that should be compared (even though we recommend to normalize the data and remove outliers prior to DMA).

In the example data we have seven different cell lines, healthy (HK2) and cancer (ccRCC: 786-M1A, 786-M2A, 786-O, OSRC2, OSLM1B and RFX631) and hence we can perform multiple different comparisons. The results can be automatically saved and all the results are returned in a list with the different data frames. If parameter Plot=TRUE, an overview Volcano plot is generated and saved.

# Perform multiple comparison All_vs_One using annova:
DMA_Annova <-  MetaProViz::DMA(InputData=Media_Preprocessed[,-c(1:6)],
                               SettingsFile_Sample=Media_Preprocessed[,c(1:4)],
                               SettingsInfo = c(Conditions="Conditions", Numerator=NULL, Denominator = "HK2"),
                               StatPval ="aov",
                               StatPadj="fdr",
                               SettingsFile_Metab = MappingInfo, 
                               CoRe=TRUE)

#Inspect the DMA results tables:
DMA_786M1A_vs_HK2 <- DMA_Annova[["DMA"]][["786-M1A_vs_HK2"]]
Shapiro <- DMA_Annova[["ShapiroTest"]][["DF"]][["Shapiro_result"]]

#> There are no NA/0 values
#> For the condition HK2 82.35 % of the metabolites follow a normal distribution and 17.65 % of the metabolites are not-normally distributed according to the shapiro test. You have chosen aov, which is for parametric Hypothesis testing. `shapiro.test` ignores missing values in the calculation.
#> For the condition 786-O 95.71 % of the metabolites follow a normal distribution and 4.29 % of the metabolites are not-normally distributed according to the shapiro test. You have chosen aov, which is for parametric Hypothesis testing. `shapiro.test` ignores missing values in the calculation.
#> For the condition 786-M1A 97.14 % of the metabolites follow a normal distribution and 2.86 % of the metabolites are not-normally distributed according to the shapiro test. You have chosen aov, which is for parametric Hypothesis testing. `shapiro.test` ignores missing values in the calculation.
#> For the condition 786-M2A 88.57 % of the metabolites follow a normal distribution and 11.43 % of the metabolites are not-normally distributed according to the shapiro test. You have chosen aov, which is for parametric Hypothesis testing. `shapiro.test` ignores missing values in the calculation.
#> For the condition OSRC2 92.86 % of the metabolites follow a normal distribution and 7.14 % of the metabolites are not-normally distributed according to the shapiro test. You have chosen aov, which is for parametric Hypothesis testing. `shapiro.test` ignores missing values in the calculation.
#> For the condition OSLM1B 85.71 % of the metabolites follow a normal distribution and 14.29 % of the metabolites are not-normally distributed according to the shapiro test. You have chosen aov, which is for parametric Hypothesis testing. `shapiro.test` ignores missing values in the calculation.
#> For the condition RFX631 97.14 % of the metabolites follow a normal distribution and 2.86 % of the metabolites are not-normally distributed according to the shapiro test. You have chosen aov, which is for parametric Hypothesis testing. `shapiro.test` ignores missing values in the calculation.
#> For 67.65% of metabolites the group variances are equal.
#> No condition was specified as numerator and HK2 was selected as a denominator. Performing multiple testing `all-vs-one` using aov.

1. Preview of the Shaprio results for the different conditions.

Code	Metabolites with normal distribution [%]	Metabolites with not-normal distribution [%]	Shapiro p.val(2-hydroxyglutarate)	Shapiro p.val(2-ketoglutarate)
HK2	82.35	17.65	0.6833007	0.0446492
786-O	95.71	4.29	0.4938675	0.3823712
786-M1A	97.14	2.86	0.9979050	0.1384610
786-M2A	88.57	11.43	0.5546558	0.5369470
OSRC2	92.86	7.14	0.8899005	0.2007242
OSLM1B	85.71	14.29	0.4643014	0.9022803
RFX631	97.14	2.86	0.9292099	0.0247568

2. Preview of the DMA results for the comparison of 786-M1A versus HK2 cells.

Metabolite	Log2(Distance)	p.adj	t.val	CoRe_specific	CoRe	Mean_786-M1A	Mean_HK2	MS51-16	MS51-17	MS51-18	MS51-19	MS51-20	MS51-07	MS51-08	MS51-09	MS51-10	HMDB	KEGG.ID	KEGGCompound	Pathway
aconitate	-32.22764	0.0000000	5029068800	Released	Released	1240145450.47062	6269214250.16407	1801837769	1235254549	1163829715	818770348	1181034871	7455333818	5658710336	5862998752	6099814094	HMDB0000072	C00417	cis-Aconitate	Citrate cycle (TCA cycle)
arginine	34.44492	0.9999860	-23385974128	Consumed	Consumed	-349377950011.591	-372763924139.179	-239724147819	-332286298627	-640453325380	-317923893286	-216502084946	-247877938391	-198621970371	-527002953935	-517552833859	HMDB0000517	C00062	L-Arginine	Amino acid metabolism
aspartate	-34.07452	0.0000000	18090607147	Consumed in 786-M1A and Released HK2	Released/Consumed	-730678742.798673	17359928404.2693	1151442869	-1794575260	-2464906916	-1239561110	694206704	18576102375	19746440364	14563233762	16553937116	HMDB0000191	C00049	L-Aspartate	Amino acid metabolism
betaine	39.11101	0.0015506	-593725992183	Released	Released	625658719313.162	31932727130.3549	727410421212	529916930520	696363835606	763335853575	411266555653	254572823562	7711585466	-435055677655	300502177148	HMDB0000043	C00719	Betaine	Not assigned
carbamoyl phosphate	29.60718	0.5445135	-817802132	Released in 786-M1A and Consumed HK2	Released/Consumed	747379583.468335	-70422548.1409482	458618834	546774968	1903694954	730697132	97112031	130196291	-425231689	479500454	-466155248	HMDB0001096	C00169	Carbamoyl phosphate	Purine metabolism

Using the DMA results, we can now use the MetaProViz visualization module and generate further customized Volcano plots VizVolcano(). You can see some examples below.

ORA using the DMA results

Over Representation Analysis (ORA) is a pathway enrichment analysis (PEA) method that determines if a set of features (=metabolic pathways) are over-represented in the selection of features (=metabolites) from the data in comparison to all measured features (metabolites) using the Fishers exact test. The selection of metabolites are usually the most altered metabolites in the data, which can be selected by the top and bottom t-values. Given that for CoRe data it is important to consider weather a metabolite was consumed or released, it is sensible to perform ORA on each metabolite cluster.
Of course, there are many other PEA methods such as the well known GSEA. Here we do not aim to provide an extensive tool for different methods to perform pathway enrichment analysis and only focus on ORA since we can apply this to perform standard pathway enrichment as well as pathway enrichment on clusters of metabolites. If you are interested in using different pathway enrichment methods please check out specialized tools such as decopupleR (Badia-I-Mompel et al. 2022).

Here we will use the KEGG pathways (Kanehisa and Goto 2000). Before we can perform ORA on the DMA results, we have to ensure that the metabolite names match with the KEGG IDs or KEGG trivial names. In general, the PathwayFile requirements are column “term”, “Metabolite” and “Description”, and the Input_data requirements are column “t.val” and column “Metabolite”.

#Since we have performed multiple comparisons (all_vs_HK2), we will run ORA for each of this comparison
DM_ORA_res<- list()

comparisons <- names(DMA_Annova[["DMA"]])
for(comparison in comparisons){
  #Ensure that the Metabolite names match with KEGG IDs or KEGG trivial names. 
  DMA <- DMA_Annova[["DMA"]][[comparison]]
  DMA <- DMA[complete.cases(DMA),-1]%>%#we remove metabolites that do not have a KEGG ID/KEGG pathway
    remove_rownames()%>%
    column_to_rownames("KEGGCompound")#We use the KEGG trivial names to match with the KEGG pathways
  
  #Perform ORA: Here we use 
  DM_ORA_res[[comparison]] <- MetaProViz::ClusterORA(InputData=DMA,
                                                     SettingsInfo=c(ClusterColumn="CoRe_specific", PathwayTerm= "term", PathwayFeature= "Metabolite"),
                                                     RemoveBackground=FALSE,#we do not have any background
                                                     PathwayFile=KEGG_Pathways,
                                                     PathwayName="KEGG",
                                                     minGSSize=3,
                                                     maxGSSize=1000)
}
#> Number of metabolites in cluster `Released in 786-M1A and Consumed HK2`: 3
#> 
#> Number of metabolites in cluster `Released`: 19
#> Number of metabolites in cluster `Consumed in 786-M1A  and Released HK2`: 9
#> Number of metabolites in cluster `Consumed`: 28
#> Number of metabolites in cluster `Released in 786-M2A and Consumed HK2`: 4
#> Number of metabolites in cluster `Released`: 21
#> Number of metabolites in cluster `Consumed in 786-M2A  and Released HK2`: 7
#> Number of metabolites in cluster `Consumed`: 27
#> Number of metabolites in cluster `Released in 786-O and Consumed HK2`: 4
#> Number of metabolites in cluster `Released`: 23
#> Number of metabolites in cluster `Consumed in 786-O  and Released HK2`: 5
#> Number of metabolites in cluster `Consumed`: 27
#> Number of metabolites in cluster `Consumed`: 29
#> Number of metabolites in cluster `Released`: 21
#> Number of metabolites in cluster `Consumed in OSLM1B  and Released HK2`: 7
#> Number of metabolites in cluster `Released in OSLM1B and Consumed HK2`: 2
#> Number of metabolites in cluster `Released in OSRC2 and Consumed HK2`: 4
#> Number of metabolites in cluster `Released`: 21
#> Number of metabolites in cluster `Consumed in OSRC2  and Released HK2`: 7
#> Number of metabolites in cluster `Consumed`: 27
#> Number of metabolites in cluster `Released in RFX631 and Consumed HK2`: 1
#> Number of metabolites in cluster `Released`: 22
#> Number of metabolites in cluster `Consumed`: 30
#> Number of metabolites in cluster `Consumed in RFX631  and Released HK2`: 6

#Lets check how the results look like:
MC_ORA_786M1A_vs_HK2_Consumed <- DM_ORA_res[["786-M1A_vs_HK2"]][["DF"]][["Consumed"]]

Preview of the ORA results for the comparison of 786-M1A versus HK2 cells focusing on pathways enriched in consumed metabolites.

GeneRatio	BgRatio	pvalue	p.adjust	qvalue	Metabolites_in_pathway	Count	Metabolites_in_Pathway	Percentage_of_Pathway_detected
12/26	17/51	0.0453148	0.6780045	0.6717073	L-Arginine/L-Asparagine/L-Glutamine/L-Histidine/L-Lysine/L-Methionine/L-Phenylalanine/L-Serine/L-Threonine/L-Tryptophan/L-Tyrosine/L-Valine	12	52	23.08
4/26	4/51	0.0598239	0.6780045	0.6717073	Linoleate/(9Z)-Octadecenoic acid/Hexadecanoic acid/Octadecanoic acid	4	74	5.41
13/26	18/51	0.0248808	0.6780045	0.6717073	L-Arginine/L-Asparagine/L-Cystine/L-Glutamine/L-Histidine/L-Lysine/L-Methionine/L-Phenylalanine/L-Serine/L-Threonine/L-Tryptophan/L-Tyrosine/L-Valine	13	40	32.50
3/26	3/51	0.1248499	0.7074830	0.7009119	(9Z)-Octadecenoic acid/Hexadecanoic acid/Octadecanoic acid	3	58	5.17
8/26	11/51	0.0980791	0.7074830	0.7009119	L-Asparagine/L-Glutamine/L-Methionine/L-Phenylalanine/L-Serine/L-Threonine/L-Tryptophan/L-Valine	8	18	44.44

MCA

Metabolite Clustering Analysis (MCA) is a module, which includes different functions to enable clustering of metabolites into groups based on logical regulatory rules. This can be particularly useful if one has multiple conditions and aims to find patterns in the data.

MCA_CoRe

This metabolite clustering method is based on logical regulatory rules to sort metabolites into metabolite clusters. Here you additionally need intracellular samples corresponding to the CoRe samples.
Here we will define if a feature (= metabolite) is assigned into:
1. “UP”, which means a metabolite is significantly up-regulated in the underlying comparison.
2. “DOWN”, which means a metabolite is significantly down-regulated in the underlying comparison.
3. “No Change”, which means a metabolite does not change significantly in the underlying comparison and/or is not defined as up-regulated/down-regulated based on the Log2FC threshold chosen.

Therebye “No Change” is further subdivided into four states:
1. “Not Detected”, which means a metabolite is not detected in the underlying comparison.
2. “Not Significant”, which means a metabolite is not significant in the underlying comparison.
3. “Significant positive”, which means a metabolite is significant in the underlying comparison and the differential metabolite abundance is positive, yet does not meet the threshold set for “UP” (e.g. Log2FC >1 = “UP” and we have a significant Log2FC=0.8).
4. “Significant negative”, which means a metabolite is significant in the underlying comparison and the differential metabolite abundance is negative, yet does not meet the threshold set for “DOWN”.

Lastly, we also take into account the CoRe direction, meaning if a metabolite was:
1. “Released”, which means is released into the media in both conditions of the underlying comparison.
2. “Consumed”, which means is consumed from the media in both conditions of the underlying comparison.
3. “Released/Consumed”, which means is consumed/released in one condition, whilst the opposite occurs in the second condition of the underlying comparison.
4. “Not Detected”, which means a metabolite is not detected in the underlying comparison.
This definition is done individually for each comparison and will impact in which metabolite cluster a metabolite is sorted into.
Since we have two comparisons (Intracellular and CoRe), we can choose between different Background settings, which defines which features will be considered for the clusters (e.g. you could include only features (= metabolites) that are detected in both comparisons, removing the rest of the features).The background methods backgroundMethod are the following from 1.1. - 1.4. from most restrictive to least restrictive:
1.1. Intra&CoRe: Most stringend background setting and will lead to a small number of metabolites.
1.2. CoRe: Focus is on the metabolite abundance of the CoRe.
1.3. Intra: Focus is on the metabolite abundance of intracellular.
1.4. Intra|CoRe: Least stringent background method, since a metabolite will be included in the input if it has been detected on one of the two conditions.

Lastly, we will get clusters of metabolites that are defined by the metabolite change in the two conditions. For example, if Alanine is “UP” based on the thresholds in both comparisons it will be sorted into the cluster “Core_UP”. As there are three 6-state6-state4 transitions between the comparisons, the flows are summarised into smaller amount of metabolite clusters using different Regulation Groupings (RG): 1. RG1_All
2. RG2_Significant taking into account genes that are significant (UP, DOWN, significant positive, significant negative)
3. RG3_SignificantChange only takes into account genes that have significant changes (UP, DOWN).

In order to define which group a metabolite is assigned to, we set two different thresholds. For intracellular those are based on the differential metabolite abundance (Log2FC) and the significance (e.g. p.adj). For the CoRe data this is based on the Log2 Distance and the significance (e.g. p.adj). For Log2FC we recommend a threshold of 0.5 or 1, whilst for the Log2 Distance one should check the distance ranges and base the threhold on this.

Regulatory rules:

#Example of all possible flows:
MCA_CORE <- MetaProViz::MCA_rules(Method="CoRe")

Metabolite Clustering Analysis: CoRe.

Intra	CoRe	Core_Direction	RG1_All	R2_Significant	RG3_Change
DOWN	DOWN	Released	Intra DOWN+ CoRe DOWN_Released	Both_DOWN (Released)	Both_DOWN (Released)
DOWN	Not Detected	Not Detected	Intra DOWN+ CoRe Not Detected	None	None
DOWN	Not Significant	Released	Intra DOWN+ CoRe Not Significant_Released	None	None
DOWN	Significant Negative	Released	Intra DOWN+ CoRe Significant Negative_Released	Both_DOWN (Released)	None
DOWN	Significant Positive	Released	Intra DOWN+ CoRe Significant Positive_Released	Opposite (Released UP)	None
DOWN	UP	Released	Intra DOWN+ CoRe UP_Released	Opposite (Released UP)	Opposite (Released UP)
UP	DOWN	Released	Intra UP+ CoRe DOWN_Released	Opposite (Released DOWN)	Opposite (Released DOWN)
UP	Not Detected	Not Detected	Intra UP+ CoRe Not Detected	None	None
UP	Not Significant	Released	Intra UP+ CoRe Not Significant_Released	None	None
UP	Significant Negative	Released	Intra UP+ CoRe Significant Negative_Released	Opposite (Released UP)	None
UP	Significant Positive	Released	Intra UP+ CoRe Significant Positive_Released	Both_UP (Released)	None
UP	UP	Released	Intra UP+ CoRe UP_Released	Both_UP (Released)	Both_UP (Released)
Not Detected	DOWN	Released	Intra Not Detected+ CoRe DOWN_Released	CoRe_DOWN (Released)	CoRe_DOWN (Released)
Not Detected	Not Detected	Not Detected	Intra Not Detected+ CoRe Not Detected	None	None
Not Detected	Not Significant	Released	Intra Not Detected+ CoRe Not Significant_Released	None	None
Not Detected	Significant Negative	Released	Intra Not Detected+ CoRe Significant Negative_Released	None	None
Not Detected	Significant Positive	Released	Intra Not Detected+ CoRe Significant Positive_Released	None	None
Not Detected	UP	Released	Intra Not Detected+ CoRe UP_Released	CoRe_UP (Released)	CoRe_UP (Released)
Significant Negative	DOWN	Released	Intra Significant Negative+ CoRe DOWN_Released	Both_DOWN (Released)	CoRe_DOWN (Released)
Significant Negative	Not Detected	Not Detected	Intra Significant Negative+ CoRe Not Detected	None	None
Significant Negative	Not Significant	Released	Intra Significant Negative+ CoRe Not Significant_Released	None	None
Significant Negative	Significant Negative	Released	Intra Significant Negative+ CoRe Significant Negative_Released	None	None
Significant Negative	Significant Positive	Released	Intra Significant Negative+ CoRe Significant Positive_Released	None	None
Significant Negative	UP	Released	Intra Significant Negative+ CoRe UP_Released	Opposite (Released UP)	CoRe_UP (Released)
Significant Positive	DOWN	Released	Intra Significant Positive+ CoRe DOWN_Released	Opposite (Released DOWN)	CoRe_DOWN (Released)
Significant Positive	Not Detected	Not Detected	Intra Significant Positive+ CoRe Not Detected	None	None
Significant Positive	Not Significant	Released	Intra Significant Positive+ CoRe Not Significant_Released	None	None
Significant Positive	Significant Negative	Released	Intra Significant Positive+ CoRe Significant Negative_Released	None	None
Significant Positive	Significant Positive	Released	Intra Significant Positive+ CoRe Significant Positive_Released	None	None
Significant Positive	UP	Released	Intra Significant Positive+ CoRe UP_Released	Both_UP (Released)	CoRe_UP (Released)
Not Significant	DOWN	Released	Intra Not Significant+ CoRe DOWN_Released	CoRe_DOWN (Released)	CoRe_DOWN (Released)
Not Significant	Not Detected	Not Detected	Intra Not Significant+ CoRe Not Detected	None	None
Not Significant	Not Significant	Released	Intra Not Significant+ CoRe Not Significant_Released	None	None
Not Significant	Significant Negative	Released	Intra Not Significant+ CoRe Significant Negative_Released	None	None
Not Significant	Significant Positive	Released	Intra Not Significant+ CoRe Significant Positive_Released	None	None
Not Significant	UP	Released	Intra Not Significant+ CoRe UP_Released	CoRe_UP (Released)	CoRe_UP (Released)
DOWN	DOWN	Consumed	Intra DOWN+ CoRe DOWN_Consumed	Both_DOWN (Consumed)	Both_DOWN (Consumed)
DOWN	Not Detected	Not Detected	Intra DOWN+ CoRe Not Detected	None	None
DOWN	Not Significant	Consumed	Intra DOWN+ CoRe Not Significant_Consumed	None	None
DOWN	Significant Negative	Consumed	Intra DOWN+ CoRe Significant Negative_Consumed	Both_DOWN (Consumed)	None
DOWN	Significant Positive	Consumed	Intra DOWN+ CoRe Significant Positive_Consumed	Opposite (Consumed UP)	None
DOWN	UP	Consumed	Intra DOWN+ CoRe UP_Consumed	Opposite (Consumed UP)	Opposite (Consumed UP)
UP	DOWN	Consumed	Intra UP+ CoRe DOWN_Consumed	Opposite (Consumed DOWN)	Opposite (Consumed DOWN)
UP	Not Detected	Not Detected	Intra UP+ CoRe Not Detected	None	None
UP	Not Significant	Consumed	Intra UP+ CoRe Not Significant_Consumed	None	None
UP	Significant Negative	Consumed	Intra UP+ CoRe Significant Negative_Consumed	Opposite (Consumed UP)	None
UP	Significant Positive	Consumed	Intra UP+ CoRe Significant Positive_Consumed	Both_UP (Consumed)	None
UP	UP	Consumed	Intra UP+ CoRe UP_Consumed	Both_UP (Consumed)	Both_UP (Consumed)
Not Detected	DOWN	Consumed	Intra Not Detected+ CoRe DOWN_Consumed	CoRe_DOWN (Consumed)	CoRe_DOWN (Consumed)
Not Detected	Not Detected	Not Detected	Intra Not Detected+ CoRe Not Detected	None	None
Not Detected	Not Significant	Consumed	Intra Not Detected+ CoRe Not Significant_Consumed	None	None
Not Detected	Significant Negative	Consumed	Intra Not Detected+ CoRe Significant Negative_Consumed	None	None
Not Detected	Significant Positive	Consumed	Intra Not Detected+ CoRe Significant Positive_Consumed	None	None
Not Detected	UP	Consumed	Intra Not Detected+ CoRe UP_Consumed	CoRe_UP (Consumed)	CoRe_UP (Consumed)
Significant Negative	DOWN	Consumed	Intra Significant Negative+ CoRe DOWN_Consumed	Both_DOWN (Consumed)	CoRe_DOWN (Consumed)
Significant Negative	Not Detected	Not Detected	Intra Significant Negative+ CoRe Not Detected	None	None
Significant Negative	Not Significant	Consumed	Intra Significant Negative+ CoRe Not Significant_Consumed	None	None
Significant Negative	Significant Negative	Consumed	Intra Significant Negative+ CoRe Significant Negative_Consumed	None	None
Significant Negative	Significant Positive	Consumed	Intra Significant Negative+ CoRe Significant Positive_Consumed	None	None
Significant Negative	UP	Consumed	Intra Significant Negative+ CoRe UP_Consumed	Opposite (Consumed UP)	CoRe_UP (Consumed)
Significant Positive	DOWN	Consumed	Intra Significant Positive+ CoRe DOWN_Consumed	Opposite (Consumed DOWN)	CoRe_DOWN (Consumed)
Significant Positive	Not Detected	Not Detected	Intra Significant Positive+ CoRe Not Detected	None	None
Significant Positive	Not Significant	Consumed	Intra Significant Positive+ CoRe Not Significant_Consumed	None	None
Significant Positive	Significant Negative	Consumed	Intra Significant Positive+ CoRe Significant Negative_Consumed	None	None
Significant Positive	Significant Positive	Consumed	Intra Significant Positive+ CoRe Significant Positive_Consumed	None	None
Significant Positive	UP	Consumed	Intra Significant Positive+ CoRe UP_Consumed	Both_UP (Consumed)	CoRe_UP (Consumed)
Not Significant	DOWN	Consumed	Intra Not Significant+ CoRe DOWN_Consumed	CoRe_DOWN (Consumed)	CoRe_DOWN (Consumed)
Not Significant	Not Detected	Not Detected	Intra Not Significant+ CoRe Not Detected	None	None
Not Significant	Not Significant	Consumed	Intra Not Significant+ CoRe Not Significant_Consumed	None	None
Not Significant	Significant Negative	Consumed	Intra Not Significant+ CoRe Significant Negative_Consumed	None	None
Not Significant	Significant Positive	Consumed	Intra Not Significant+ CoRe Significant Positive_Consumed	None	None
Not Significant	UP	Consumed	Intra Not Significant+ CoRe UP_Consumed	CoRe_UP (Consumed)	CoRe_UP (Consumed)
DOWN	DOWN	Released/Consumed	Intra DOWN + CoRe DOWN_Released/Consumed	Both_DOWN (Released/Consumed)	Both_DOWN (Released/Consumed)
DOWN	Not Detected	Not Detected	Intra DOWN + CoRe Not Detected	None	None
DOWN	Not Significant	Released/Consumed	Intra DOWN + CoRe Not Significant_Released/Consumed	None	None
DOWN	Significant Negative	Released/Consumed	Intra DOWN + CoRe Significant Negative_Released/Consumed	Both_DOWN (Released/Consumed)	None
DOWN	Significant Positive	Released/Consumed	Intra DOWN + CoRe Significant Positive_Released/Consumed	Opposite (Released/Consumed UP)	None
DOWN	UP	Released/Consumed	Intra DOWN + CoRe UP_Released/Consumed	Opposite (Released/Consumed UP)	Opposite (Released/Consumed UP)
UP	DOWN	Released/Consumed	Intra UP + CoRe DOWN_Released/Consumed	Opposite (Released/Consumed DOWN)	Opposite (Released/Consumed DOWN)
UP	Not Detected	Not Detected	Intra UP + CoRe Not Detected	None	None
UP	Not Significant	Released/Consumed	Intra UP + CoRe Not Significant_Released/Consumed	None	None
UP	Significant Negative	Released/Consumed	Intra UP + CoRe Significant Negative_Released/Consumed	Opposite (Released/Consumed UP)	None
UP	Significant Positive	Released/Consumed	Intra UP + CoRe Significant Positive_Released/Consumed	Both_UP (Released/Consumed)	None
UP	UP	Released/Consumed	Intra UP + CoRe UP_Released/Consumed	Both_UP (Released/Consumed)	Both_UP (Released/Consumed)
Not Detected	DOWN	Released/Consumed	Intra Not Detected + CoRe DOWN_Released/Consumed	CoRe_DOWN (Released/Consumed)	CoRe_DOWN (Released/Consumed)
Not Detected	Not Detected	Not Detected	Intra Not Detected + CoRe Not Detected	None	None
Not Detected	Not Significant	Released/Consumed	Intra Not Detected + CoRe Not Significant_Released/Consumed	None	None
Not Detected	Significant Negative	Released/Consumed	Intra Not Detected + CoRe Significant Negative_Released/Consumed	None	None
Not Detected	Significant Positive	Released/Consumed	Intra Not Detected + CoRe Significant Positive_Released/Consumed	None	None
Not Detected	UP	Released/Consumed	Intra Not Detected + CoRe UP_Released/Consumed	CoRe_UP (Released/Consumed)	CoRe_UP (Released/Consumed)
Significant Negative	DOWN	Released/Consumed	Intra Significant Negative + CoRe DOWN_Released/Consumed	Both_DOWN (Released/Consumed)	CoRe_DOWN (Released/Consumed)
Significant Negative	Not Detected	Not Detected	Intra Significant Negative + CoRe Not Detected	None	None
Significant Negative	Not Significant	Released/Consumed	Intra Significant Negative + CoRe Not Significant_Released/Consumed	None	None
Significant Negative	Significant Negative	Released/Consumed	Intra Significant Negative + CoRe Significant Negative_Released/Consumed	None	None
Significant Negative	Significant Positive	Released/Consumed	Intra Significant Negative + CoRe Significant Positive_Released/Consumed	None	None
Significant Negative	UP	Released/Consumed	Intra Significant Negative + CoRe UP_Released/Consumed	Opposite (Released/Consumed UP)	CoRe_UP (Released/Consumed)
Significant Positive	DOWN	Released/Consumed	Intra Significant Positive + CoRe DOWN_Released/Consumed	Opposite (Released/Consumed DOWN)	CoRe_DOWN (Released/Consumed)
Significant Positive	Not Detected	Not Detected	Intra Significant Positive + CoRe Not Detected	None	None
Significant Positive	Not Significant	Released/Consumed	Intra Significant Positive + CoRe Not Significant_Released/Consumed	None	None
Significant Positive	Significant Negative	Released/Consumed	Intra Significant Positive + CoRe Significant Negative_Released/Consumed	None	None
Significant Positive	Significant Positive	Released/Consumed	Intra Significant Positive + CoRe Significant Positive_Released/Consumed	None	None
Significant Positive	UP	Released/Consumed	Intra Significant Positive + CoRe UP_Released/Consumed	Both_UP (Released/Consumed)	CoRe_UP (Released/Consumed)
Not Significant	DOWN	Released/Consumed	Intra Not Significant + CoRe DOWN_Released/Consumed	CoRe_DOWN (Released/Consumed)	CoRe_DOWN (Released/Consumed)
Not Significant	Not Detected	Not Detected	Intra Not Significant + CoRe Not Detected	None	None
Not Significant	Not Significant	Released/Consumed	Intra Not Significant + CoRe Not Significant_Released/Consumed	None	None
Not Significant	Significant Negative	Released/Consumed	Intra Not Significant + CoRe Significant Negative_Released/Consumed	None	None
Not Significant	Significant Positive	Released/Consumed	Intra Not Significant + CoRe Significant Positive_Released/Consumed	None	None
Not Significant	UP	Released/Consumed	Intra Not Significant + CoRe UP_Released/Consumed	CoRe_UP (Released/Consumed)	CoRe_UP (Released/Consumed)

Now we can load the corresponding pre-processed intracellular example data for the comparison of 786M-1A versus HK2 (For the detailed pre-processing please see the vignette “Standard Metabolomics”).

#Load the Pre-processed intracellular data:
MetaProViz::ToyData(data="Standard_DMA")

#Perform metabolite clustering:
MCA_CoRe_res <- MetaProViz::MCA_CoRe(InputData_Intra =Intra_DMA_786M1A_vs_HK2 ,
                                     InputData_CoRe = DMA_786M1A_vs_HK2,
                                     SettingsInfo_Intra=c(ValueCol="Log2FC",StatCol="p.adj", StatCutoff= 0.05, ValueCutoff=1),
                                     SettingsInfo_CoRe=c(DirectionCol="CoRe", ValueCol="Log2(Distance)",StatCol="p.adj", StatCutoff= 0.05, ValueCutoff=28),
                                     FeatureID= "Metabolite",
                                     BackgroundMethod="Intra&CoRe",
                                     FolderPath=NULL)

#Lets check how the results look like:
MCA_res <- MCA_CoRe_res[["MCA_CoRe_Results"]]
ClusterSummary <- MCA_CoRe_res[["MCA_CoRe_Summary"]]

MetaProViz::MCA_CoRe for the comparison of 786-M1A versus HK2 cells in intracellular and CoRe samples.

Metabolite	Intra_DF_Cutoff	Intra_DF_Cutoff_Specific.x	CoRe_DF_Detected	CoRe_DF_Cutoff	CoRe_DF_Cutoff_Specific	BG_Method	RG1_All	RG2_Significant	RG3_Change
adenosine	No Change	Not Significant	FALSE	No Change	Not Detected	FALSE	Background = FALSE	Background = FALSE	Background = FALSE
ADP	No Change	Not Significant	FALSE	No Change	Not Detected	FALSE	Background = FALSE	Background = FALSE	Background = FALSE
betaine	No Change	Significant Negative	TRUE	UP	UP	TRUE	Intra Significant Negative + CoRe UP_Released	Opposite (Released UP)	CoRe_UP (Released)
creatine	No Change	Significant Positive	TRUE	UP	UP	TRUE	Intra Significant Positive + CoRe UP_Released	Both_UP (Released)	CoRe_UP (Released)

MetaProViz::MCA_CoRe Summary of number of metabolites per cluster.

Regulation Grouping	SiRCle cluster Name	Number of Features
RG2_Significant	Background = FALSE	151
RG2_Significant	None	18
RG2_Significant	Both_UP (Released)	3
RG2_Significant	Both_DOWN (Released/Consumed)	3
RG2_Significant	Opposite (Released DOWN)	1
RG2_Significant	Both_DOWN (Released)	4
RG2_Significant	Opposite (Released/Consumed DOWN)	1
RG2_Significant	Opposite (Released UP)	1
RG2_Significant	CoRe_DOWN (Consumed)	4
RG2_Significant	CoRe_UP (Consumed)	1
RG2_Significant	Both_DOWN (Consumed)	6
RG2_Significant	CoRe_UP (Released)	3
RG2_Significant	Opposite (Consumed DOWN)	2
RG2_Significant	CoRe_DOWN (Released)	1
RG3_Change	Background = FALSE	151
RG3_Change	None	18
RG3_Change	Both_UP (Released)	1
RG3_Change	CoRe_DOWN (Released/Consumed)	4
RG3_Change	CoRe_DOWN (Released)	2
RG3_Change	Both_DOWN (Released)	4
RG3_Change	CoRe_UP (Released)	6
RG3_Change	CoRe_DOWN (Consumed)	8
RG3_Change	CoRe_UP (Consumed)	1
RG3_Change	Opposite (Consumed DOWN)	2
RG3_Change	Both_DOWN (Consumed)	2

ORA on each metabolite cluster

As explained in detail above, Over Representation Analysis (ORA) is a pathway enrichment analysis (PEA) method. As ORA is based on the Fishers exact test it is perfect to test if a set of features (=metabolic pathways) are over-represented in the selection of features (= clusters of metabolites) from the data in comparison to all measured features (all metabolites). In detail, MC_ORA() will perform ORA on each of the metabolite clusters using all metabolites as the background.

Pathway Input for MetaProViz::MC_ORA.

	HMDB	KEGG.ID	KEGGCompound	Pathway
N-acetylaspartate	HMDB0000812	C01042	N-Acetyl-L-aspartate	Alanine, aspartate and glutamate metabolism
argininosuccinate	HMDB0000052	C03406	N-(L-Arginino)succinate	Alanine, aspartate and glutamate metabolism
N-acetylaspartylglutamate	HMDB0001067	C12270	N-Acetylaspartylglutamate	Alanine, aspartate and glutamate metabolism
tyrosine	HMDB0000158	C00082	L-Tyrosine	Amino acid metabolism
asparagine	HMDB0000168	C00152	L-Asparagine	Amino acid metabolism
glutamate	HMDB0000148	C00025	L-Glutamate	Amino acid metabolism

MC_ORA_result<- MetaProViz::ClusterORA(InputData=MCA_CoRe_res[["MCA_CoRe_Results"]]%>%column_to_rownames("Metabolite"),
                                       SettingsInfo=c(ClusterColumn="RG2_Significant", 
                                                        BackgroundColumn="BG_Method", 
                                                        PathwayTerm= "Pathway", #This is the column name including the pathways names
                                                        PathwayFeature= "Metabolite"),
                                       RemoveBackground=TRUE,
                                       PathwayFile=MappingInfo%>%rownames_to_column("Metabolite"),
                                       PathwayName="KEGG",
                                       minGSSize=3,
                                       maxGSSize=1000 ,
                                       SaveAs_Table= "csv")
#> Number of metabolites in cluster `None`: 18
#> Number of metabolites in cluster `Both_UP (Released)`: 3
#> Number of metabolites in cluster `Both_DOWN (Released/Consumed)`: 3
#> Number of metabolites in cluster `Opposite (Released DOWN)`: 1
#> Number of metabolites in cluster `Both_DOWN (Released)`: 4
#> Number of metabolites in cluster `Opposite (Released/Consumed DOWN)`: 1
#> Number of metabolites in cluster `Opposite (Released UP)`: 1
#> Number of metabolites in cluster `CoRe_DOWN (Consumed)`: 4
#> Number of metabolites in cluster `CoRe_UP (Consumed)`: 1
#> Number of metabolites in cluster `Both_DOWN (Consumed)`: 6
#> Number of metabolites in cluster `CoRe_UP (Released)`: 3
#> Number of metabolites in cluster `Opposite (Consumed DOWN)`: 2
#> Number of metabolites in cluster `CoRe_DOWN (Released)`: 1

#Lets check how the results look like:
Both_UP_Released <- MC_ORA_result[["DF"]][["Both_UP (Released)"]]

MetaProViz::MC_ORA results for the RG2_Significant cluster Both_UP (Released).

ID	GeneRatio	BgRatio	pvalue	p.adjust	qvalue	Metabolites_in_pathway	Count	Metabolites_in_Pathway	Percentage_of_Pathway_detected
Arginine and proline metabolism	1/3	3/48	0.1795791	0.5043941	0.5043941	creatine	1	7	14.29
Citrate cycle (TCA cycle)	1/3	6/48	0.3362627	0.5043941	0.5043941	2-ketoglutarate	1	10	10.00
Amino acid metabolism	1/3	18/48	0.7652636	0.7652636	0.7652636	glutamate	1	35	2.86
Fatty acyl carnitines	NA	NA	NA	NA	NA	NA	0	11	0.00

Here we see that the pathways have a low amount of genes included that were also part of the cluster and the pathways are not significant. This is due to multiple factors, first we only start with a small number of metabolites with KEGG IDs and secondly we only included metabolites if they where detected in both, intracellular and CoRe samples (parameter backgroundMethod="Intra&CoRe"). Hence, by for example setting parameter backgroundMethod="Intra|CoRe", we will obtain larger metabolite clusters.

3. Run MetaProViz Visualisation

The big advantages of the MetaProViz visualization module is its flexible and easy usage, which we will showcase below and that the figures are saved in a publication ready style and format. For instance, the x- and y-axis size will always be adjusted for the amount of samples or features (=metabolites) plotted, or in the case of Volcano plot and PCA plot the axis size is fixed and not affected by figure legends or title. In this way, there is no need for many adjustments and the figures can just be dropped into the presentation or paper and are all in the same style.

All the VizPlotName() functions are constructed in the same way. Indeed, with the parameter Plot_SettingsInfo the user can pass a named vector with information about the metadata column that should be used to customize the plot by colour, shape or creating individual plots, which will all be showcased for the different plot types. Via the parameter Plot_SettingsFile the user can pass the metadata DF, which can be dependent on the plot type for the samples and/or the features (=metabolites). In case of both the parameter is named Plot_SettingsFile_Sample and Plot_SettingsFile_Metab.

In each of those Plot_Settings, the user can label color and/or shape based on additional information (e.g. Pathway information, Cluster information or other other demographics like gender). Moreover, we also enable to plot individual plots where applicable based on those MetaData (e.g. one plot for each metabolic pathway).
For this we need a metadata table including information about our samples that could be relevant to e.g. color code:

MetaData_Sample <- Media_Preprocessed[,c(1:2)]%>%
   mutate(Status = case_when(Conditions=="HK2" ~ 'Healthy',
                               TRUE ~ 'Cancer'))

Metadata table including additional information about our Samples.

	Conditions	Biological_Replicates	Status
MS51-07	HK2	2	Healthy
MS51-08	HK2	3	Healthy
MS51-09	HK2	4	Healthy
MS51-10	HK2	5	Healthy
MS51-11	786-O	1	Cancer
MS51-12	786-O	2	Cancer
MS51-13	786-O	3	Cancer
MS51-14	786-O	4	Cancer
MS51-15	786-O	5	Cancer
MS51-16	786-M1A	1	Cancer
MS51-17	786-M1A	2	Cancer
MS51-18	786-M1A	3	Cancer
MS51-19	786-M1A	4	Cancer
MS51-20	786-M1A	5	Cancer
MS51-21	786-M2A	1	Cancer
MS51-23	786-M2A	2	Cancer
MS51-24	786-M2A	3	Cancer
MS51-25	786-M2A	4	Cancer
MS51-26	OSRC2	1	Cancer
MS51-27	OSRC2	2	Cancer
MS51-28	OSRC2	3	Cancer
MS51-29	OSRC2	4	Cancer
MS51-30	OSRC2	5	Cancer
MS51-31	OSLM1B	1	Cancer
MS51-32	OSLM1B	2	Cancer
MS51-33	OSLM1B	3	Cancer
MS51-34	OSLM1B	4	Cancer
MS51-35	OSLM1B	5	Cancer
MS51-36	RFX631	1	Cancer
MS51-37	RFX631	2	Cancer
MS51-38	RFX631	3	Cancer
MS51-39	RFX631	4	Cancer
MS51-40	RFX631	5	Cancer

Moreover, we can use MetaData for our features (=Metabolites), which we loaded with the MappingInfo and we can also add the information on which cluster a metabolite was assigned to in the MetaProViz::MCA() analysis above:

MetaData_Metab <-MappingInfo

Metadata table including additional information about the Metabolites.

HMDB	KEGG.ID	KEGGCompound	Pathway
HMDB0001067	C12270	N-Acetylaspartylglutamate	Alanine, aspartate and glutamate metabolism
HMDB0000158	C00082	L-Tyrosine	Amino acid metabolism
HMDB0000148	C00025	L-Glutamate	Amino acid metabolism
HMDB0250980	NA	NA	Amino acid metabolism
NA	NA	NA	Amino acid metabolism
HMDB0004207	NA	NA	Amino acid metabolism

Noteworthy, here we can also use the KEGG pathways we used for the pathway analysis.

PCA plots

Principal component analysis (PCA) is a dimensionality reduction method that reduces all the measured features (=metabolites) of one sample into a few features in the different principal components, whereby each principal component can explain a certain percentage of the variance between the different samples. Hence, this enables interpretation of sample clustering based on the measured features (=metabolites).
As mentioned above, PCA plots can be quite useful for quality control, but of course it offers us many more opportunities, which will be showcased here.

As input, we need a DF that contains the samples as rownames and the features (=metabolites) as column names:

Input_PCA <- Media_Preprocessed[,-c(1:4)] #remove columns that include Metadata such as cell type,...

Input_data for MetaProViz::VizPCA(), with samples as rownames and metabolites as column names.

	valine-d8	hipppuric acid-d5	2-hydroxyglutarate	2-ketoglutarate	3-Dehydro-L-threonate	acetylcarnitine	acetylcholine
MS51-07	1958125674	1886704843	-1209262732	24617648085	78317130942	54840869029	281972444
MS51-08	-4532562455	-39893563477	767684868	26771101645	88956838993	-112289337769	654251353
MS51-09	-527252000	-38568863508	-530294707	13989640920	-18139183358	340901904708	1304844375
MS51-10	-9262935498	-36751484785	-390772056	26968551969	49440162385	206100411537	873125674
MS51-11	23691497374	89855056020	3905415435	50824311101	42493117379	-288960790882	745540459
MS51-12	28311034766	66171788552	402018301	46919285080	42865593269	-437681160359	-69825257

Now lets check out the standard plot:

MetaProViz::VizPCA(InputData=Input_PCA)

Figure: Standard Settings.

Next, we can interactively choose shape and color using the additional information of interest from our Metadata. Especially for complex data, such as patient data, it can be valuable to use different demographics (e.g. age, gender, medication,…) for this. First lets check if we have any batch effect by colour coding for the biological replicates, which would be the case if the replicates cluster together.

MetaProViz::VizPCA(SettingsInfo= c(color="Biological_Replicates"),
                   SettingsFile_Sample = MetaData_Sample ,
                   InputData=Input_PCA,
                   PlotName = "Batch Effect")

Figure: Do we have a batch effect?

Given the biological replicates are numeric, we can also set color_scale to continuous:

MetaProViz::VizPCA(SettingsInfo= c(color="Biological_Replicates"),
                   SettingsFile_Sample = MetaData_Sample ,
                   InputData=Input_PCA,
                   ColorScale = "continuous",
                   PlotName = "Batch Effect (continuous color scale)")

Figure: Do we have a batch effect?

Next, we can colour code for condition and use the biological replicates in the shape parameter:

MetaProViz::VizPCA(SettingsInfo= c(color="Conditions", shape="Biological_Replicates"),
                   SettingsFile_Sample = MetaData_Sample ,
                   InputData=Input_PCA,
                   PlotName = "Sample Conditions")

Figure: Do the samples cluster for the conditions?

The different cell lines we have are either control or cancerous, so we can display this too.

MetaProViz::VizPCA(SettingsInfo=  c(color="Status"),
                   SettingsFile_Sample = MetaData_Sample ,
                   InputData=Input_PCA,
                   PlotName = "Sample Status")

Figure: Do the samples cluster for the Cell status?

Heatmaps

Clustered heatmaps can be useful to understand the patterns in the data, which will be showcased on different examples.
As input, we need a DF that contains the samples as rownames and the features (=metabolites) as column names:

Input_Heatmap <-   Media_Preprocessed[,-c(1:4)] #remove columns that include Metadata such as cell type,...

Input for MetaProViz::VizHeatmap(), with samples as rownames and metabolites as column names.

	valine-d8	hipppuric acid-d5	2-hydroxyglutarate	2-ketoglutarate	3-Dehydro-L-threonate
MS51-07	1958125674	1886704843	-1209262732	24617648085	78317130942
MS51-08	-4532562455	-39893563477	767684868	26771101645	88956838993
MS51-09	-527252000	-38568863508	-530294707	13989640920	-18139183358
MS51-10	-9262935498	-36751484785	-390772056	26968551969	49440162385
MS51-11	23691497374	89855056020	3905415435	50824311101	42493117379
MS51-12	28311034766	66171788552	402018301	46919285080	42865593269

Now we can generate an overview heatmap. Since we plot all metabolites the metabolite names are not plotted since this would get too crowded (You can enforce this by changing the parameter enforce_FeatureNames = TRUE).

MetaProViz::VizHeatmap(InputData = Input_Heatmap, 
                       PlotName = "Overview")

Overview heatmap.

Overview heatmap.

Here we can add as many sample metadata information as needed at the same time:

MetaProViz::VizHeatmap(InputData = Input_Heatmap,
                       SettingsFile_Sample = MetaData_Sample,
                       SettingsInfo = c(color_Sample = list("Conditions","Biological_Replicates", "Status")),
                       PlotName = "Colour Samples")

Colour for sample metadata.

Colour for sample metadata.

Moreover, we can also add metabolite metadata information:

MetaProViz::VizHeatmap(InputData = Input_Heatmap,
                       SettingsFile_Sample = MetaData_Sample,
                       SettingsInfo = c(color_Metab = list("Pathway")),
                       SettingsFile_Metab =  MappingInfo, 
                       PlotName = "Colour Metabolites")

Colour for metabolite metadata.

Colour for metabolite metadata.

Lastly, by generate individual plot for e.g. each pathway or the metabolite clusters by adding individual (individual_Sample or individual_Metab) to Plot_SettingsInfo. At the same time we can still maintain the metadata information for both, the samples and the metabolites. Together this can help us to draw biological conclusions about the different pathways: Indeed, we can observe for the D-Amino acid metabolism many metabolites fall into the MCA-Cluster Core_DOWN, meaning in comparison to HK2 cells we have a negative Log2FC for 786-O and 786-M1A.

# individual: One individual plot for each pathway, col annotation: Colour for samples
MetaProViz::VizHeatmap(InputData = Input_Heatmap, 
                       SettingsFile_Sample = MetaData_Sample,
                       SettingsInfo = c(individual_Metab = "Pathway",
                                        color_Sample = list("Conditions","Biological_Replicates"),
                                        color_Metab = list("Pathway")),
                       SettingsFile_Metab =  MetaData_Metab,
                       PlotName = "Pathway")

You can also choose to make individual plots for any Sample Metadata using individual_Sample (e.g. in patients you may want to plot male and female separately). Moreover, you can also use both at the same time.

Volcano plot

In general,we have three different Plot_Settings, which will also be used for other plot types such as lollipop graphs.
1. "Standard" is the standard version of the plot, with one dataset being plotted.
2. "Conditions" here two or more datasets will be plotted together. How datasets can be plotted together depends on the plot type.
3. "PEA" stands for Pathway Enrichment Analysis, and is used if the results of an GSE analysis should be plotted as here the figure legends will be adapted. You can find an example for this in the vignette Standard Metabolomics

Here we will look at all the different options we have to display our results from the different analysis, which will help us to interpret our results as this can be sometimes difficult to do from the many data tables.
Just a quick reminder, how the input data look like:
1. Results of Differential Metabolite Analysis (DMA): Log2(Distance) and stats:

Input_data for MetaProViz::VizVolcano() are for example differential analysis results from MetaProViz::DMA().

Metabolite	Log2(Distance)	p.adj	t.val	CoRe_specific	CoRe
2-hydroxyglutarate	31.49006	0.3729158	-3016157146	Released in 786-M1A and Consumed HK2	Released/Consumed
2-ketoglutarate	34.69337	0.0000000	-27780823966	Released	Released
3-Dehydro-L-threonate	33.18913	0.9902860	-9793175378	Released	Released
acetylcarnitine	-38.27316	0.1061020	332176666710	Consumed in 786-M1A and Released HK2	Released/Consumed
acetylcholine	-29.84850	0.0035677	966707924	Consumed in 786-M1A and Released HK2	Released/Consumed
acetylornithine	32.26766	0.0196175	-5170506138	Consumed	Consumed
aconitate	-32.22764	0.0000000	5029068800	Released	Released

Standard

Here we will first look into the results from the differential analysis (see section DMA above) for the comparison of 786-M1A_vs_HK2:

# Run with default parameter --> only need to provide Input_data and the title we like
MetaProViz::VizVolcano(InputData=DMA_786M1A_vs_HK2%>%column_to_rownames("Metabolite"),
                       x= "Log2(Distance)")

Figure: Standard figure displaying DMA results.

If you seek to plot the metabolite names you can change the paramter SelectLab from its default (SelectLab="") to NULL and the metabolite names will be plotted randomly.

# Run with default parameter --> only need to provide Input_data and the title we like
MetaProViz::VizVolcano(InputData=DMA_786M1A_vs_HK2%>%column_to_rownames("Metabolite"),
                       x= "Log2(Distance)",
                       SelectLab = NULL)

Figure: Standard figure displaying DMA results.

With the parameter SelectLab you can also pass a vector with Metabolite names that should be labeled:

# Run with default parameter --> only need to provide Input_data and the title we like
MetaProViz::VizVolcano(InputData=DMA_786M1A_vs_HK2%>%column_to_rownames("Metabolite"),
                       x= "Log2(Distance)",
                       SelectLab = c("histidine", "phenylalanine", "lactate"))

Figure: Standard figure displaying DMA results.

As explained above, when analyzing CoRe data it is important to take into account if a metabolite is consumed or released. we can use this information to colour code and or shape the metabolites on the plot.
For this we need to add this information into the Metadata_Metabolite file:

# colour for consumption and release: For this we need to add this information into the Metadata_Metabolite file
MetaData_Metab <- merge(MappingInfo%>%rownames_to_column("Metabolite"), DMA_786M1A_vs_HK2[,c(1,5:6)], by="Metabolite", all.y=TRUE)%>%
  column_to_rownames("Metabolite")

Metadata table including additional information about the Metabolites.

	HMDB	KEGG.ID	KEGGCompound	Pathway	CoRe_specific	CoRe
3-Dehydro-L-threonate	NA	NA	NA	Amino acid metabolism	Released	Released
acetylcarnitine	HMDB0000201	C02571	O-Acetylcarnitine	Fatty acyl carnitines	Consumed in 786-M1A and Released HK2	Released/Consumed
acetylornithine	HMDB0003357	C00437	N-Acetylornithine	Arginine and proline metabolism	Consumed	Consumed
glutamate	HMDB0000148	C00025	L-Glutamate	Amino acid metabolism	Released	Released
glutamine	HMDB0000641	C00064	L-Glutamine	Amino acid metabolism	Consumed	Consumed
glycerylphosphorylcholine	HMDB0252858	NA	NA	Not assigned	Released	Released

Now we can make the different plots:

#Now we need to add our Plot_SettingsFile and the Plot_SettingsInfo:
MetaProViz::VizVolcano(PlotSettings="Standard",
                       SettingsInfo= c(color="CoRe_specific"),
                       SettingsFile_Metab= MetaData_Metab,
                       InputData=DMA_786M1A_vs_HK2%>%column_to_rownames("Metabolite"),
                       x= "Log2(Distance)",
                       PlotName= "786M1A versus HK2",
                       Subtitle= "Results of DMA. Colour coded for consumption/release" )

Figure: Standard figure displaying DMA results.

#If we want to use the shape instead of the colour for the cluster info, we can just change our Plot_SettingsInfo
MetaProViz::VizVolcano(PlotSettings="Standard",
                       SettingsInfo= c(shape="CoRe_specific"),
                       SettingsFile= MetaData_Metab,
                       InputData=DMA_786M1A_vs_HK2%>%column_to_rownames("Metabolite"),
                       x= "Log2(Distance)",
                       PlotName= "786M1A versus HK2",
                       Subtitle= "Results of DMA. Shape for consumption/release, color for significance." )

Figure: Standard figure displaying DMA results.

#Of course, we can also adapt both, color and shape for the same parameter:
MetaProViz::VizVolcano(PlotSettings="Standard",
                       SettingsInfo= c(shape="CoRe_specific", color="CoRe_specific"),
                       SettingsFile= MetaData_Metab,
                       InputData=DMA_786M1A_vs_HK2%>%column_to_rownames("Metabolite"),
                       x= "Log2(Distance)",
                       PlotName= "786M1A versus HK2",
                       Subtitle= "Results of DMA. Shape and color for consumption/release." )

Figure: Standard figure displaying DMA results.

Of course, here we may also want an individual plot for each of the consumption/release metabolites.

#individual plot for each metabolite behaviour:
MetaProViz::VizVolcano(PlotSettings="Standard",
                       SettingsInfo= c(individual="CoRe", shape="CoRe_specific"),
                       SettingsFile= MetaData_Metab,
                       InputData=DMA_786M1A_vs_HK2%>%column_to_rownames("Metabolite"),
                       x= "Log2(Distance)",
                       PlotName= "786M1A versus HK2",
                       Subtitle= "Results of DMA." )

Figure: Standard figure displaying DMA results.

Figure: Standard figure displaying DMA results.

Figure: Standard figure displaying DMA results.

Given that we also know, which metabolic pathway the metabolites correspond to, we can add this information into the plot. This is also a good example to showcase the flexibility of the visualisation function: Either you use the parameter Plot_SettingsFile= MetaData_Metab as above, but as we have the column “Pathway” also in our Input_data you can also pass Plot_SettingsFile= DMA_786-M1A_vs_HK2 or simply use the default Plot_SettingsFile=NULL, in which case the Plot_SettingsInfo information (here color) will be used from Input_data.

#Now we can use color for the pathways and shape for the metabolite clusters:
MetaProViz::VizVolcano(PlotSettings="Standard",
                       SettingsInfo= c(individual="CoRe", shape="CoRe_specific", color="Pathway"),
                       SettingsFile= MetaData_Metab,
                       InputData=DMA_786M1A_vs_HK2%>%column_to_rownames("Metabolite"),
                       x= "Log2(Distance)",
                       PlotName= "786M1A versus HK2",
                       Subtitle= "Results of DMA." )

Figure: Standard figure displaying DMA results colour coded for metabolic pathways and shaped for metabolic clusters.

Figure: Standard figure displaying DMA results colour coded for metabolic pathways and shaped for metabolic clusters.

Figure: Standard figure displaying DMA results colour coded for metabolic pathways and shaped for metabolic clusters.

Comparison

The parameter Plot_Settings="Compare" is helpful if you have performed multiple comparisons and seek to compare two of them in one plot:

#Make the plot
MetaProViz::VizVolcano(PlotSettings="Compare",
                       InputData=DMA_786M1A_vs_HK2%>%column_to_rownames("Metabolite"),
                       InputData2= DMA_Annova[["DMA"]][["786-O_vs_HK2"]]%>%column_to_rownames("Metabolite"),
                       ComparisonName= c(InputData="786M1A_vs_HK", InputData2= "786-O_vs_HK2"),
                       x= "Log2(Distance)",
                       PlotName= "786M1A vs HK2 compared to 7860 vs HK2",
                       Subtitle= "Results of DMA" )

Figure: Comparison.

Of course you have option to use shape or color to further customize your graph as well as make individual plots:

#Make the plot
MetaProViz::VizVolcano(PlotSettings="Compare",
                       SettingsInfo= c(color="Pathway"),
                       SettingsFile_Metab= MetaData_Metab,
                       InputData=DMA_786M1A_vs_HK2%>%column_to_rownames("Metabolite"),
                       InputData2= DMA_Annova[["DMA"]][["786-O_vs_HK2"]]%>%column_to_rownames("Metabolite"),
                       ComparisonName= c(InputData="786M1A_vs_HK", InputData2= "786-O_vs_HK2"),
                       x= "Log2(Distance)",
                       PlotName= "786M1A vs HK2 compared to 7860 vs HK2",
                       Subtitle= "Results of DMA" )

Figure: Comparison.

Now we do individual plots again:

MetaProViz::VizVolcano(PlotSettings="Compare",
                       SettingsInfo= c(individual="Pathway"),
                       SettingsFile_Metab= MetaData_Metab,
                       InputData=DMA_786M1A_vs_HK2%>%column_to_rownames("Metabolite"),
                       InputData2= DMA_Annova[["DMA"]][["786-O_vs_HK2"]]%>%column_to_rownames("Metabolite"),
                       ComparisonName= c(InputData="786M1A_vs_HK", InputData2= "786-O_vs_HK2"),
                       x= "Log2(Distance)",
                       PlotName= "786M1A vs HK2 compared to 7860 vs HK2",
                       Subtitle= "Results of DMA" )

PathwayEnrichmentAnalysis

If you have performed Pathway Enrichment Analysis (PEA) such as ORA or GSEA, we can also plot the results and add the information into the Figure legends.
Here we can for example use the results of the ORA we have performed on the differential expression results. Indeed for DMA_786M1A_vs_HK2 we have performed ORA on each cluster (consumed, released, consumed/released). Here, I will plot the ORA results of the metabolites that are released in both conditions, HK2 and 786-M1A.

#Prepare the Input:
#1. InputData=Pathway analysis input: Must have features as column names. Those feature names need to match features in the pathway analysis file SettingsFile_Metab. 
InputPEA <- DMA_786M1A_vs_HK2 %>%
  filter(!is.na(KEGGCompound)) %>%
  column_to_rownames("KEGGCompound")

#2. InputData2=Pathway analysis output: Must have same column names as SettingsFile_Metab for Pathway name
InputPEA2 <- MC_ORA_786M1A_vs_HK2_Consumed %>%
  dplyr::rename("term"="ID")

#3. SettingsFile_Metab= Pathways used for pathway analysis: Must have same column names as SettingsFile_Metab for Pathway name and feature names need to match features in the InputData. PEA_Feature passes this column name!

MetaProViz::VizVolcano(PlotSettings="PEA",
                       SettingsInfo= c(PEA_Pathway="term",# Needs to be the same in both, SettingsFile_Metab and InputData2.
                                       PEA_stat="p.adjust",#Column InputData2
                                       PEA_score="GeneRatio",#Column InputData2
                                       PEA_Feature="Metabolite"),# Column SettingsFile_Metab (needs to be the same as row names in InputData)
                       SettingsFile_Metab= KEGG_Pathways,#Must be the pathways used for pathway analysis
                       InputData= InputPEA, #Must be the data you have used as an input for the pathway analysis
                       InputData2= InputPEA2, #Must be the results of the pathway analysis
                       x= "Log2(Distance)",
                       PlotName= "KEGG",
                       Subtitle= "PEA" ,
                       SelectLab = NULL)

Session information

#> ─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────
#>  setting  value
#>  version  R version 4.3.2 (2023-10-31 ucrt)
#>  os       Windows 10 x64 (build 19045)
#>  system   x86_64, mingw32
#>  ui       RTerm
#>  language en
#>  collate  English_United Kingdom.utf8
#>  ctype    English_United Kingdom.utf8
#>  tz       Europe/Berlin
#>  date     2024-04-22
#>  pandoc   3.1.1 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────
#>  package          * version   date (UTC) lib source
#>  abind              1.4-5     2016-07-21 [1] CRAN (R 4.3.1)
#>  AnnotationDbi      1.64.1    2023-11-03 [1] Bioconductor
#>  ape                5.7-1     2023-03-13 [1] CRAN (R 4.3.2)
#>  aplot              0.2.2     2023-10-06 [1] CRAN (R 4.3.2)
#>  backports          1.4.1     2021-12-13 [1] CRAN (R 4.3.1)
#>  Biobase            2.62.0    2023-10-24 [1] Bioconductor
#>  BiocGenerics       0.48.1    2023-11-01 [1] Bioconductor
#>  BiocParallel       1.36.0    2023-10-24 [1] Bioconductor
#>  Biostrings         2.70.1    2023-10-25 [1] Bioconductor
#>  bit                4.0.5     2022-11-15 [1] CRAN (R 4.3.2)
#>  bit64              4.0.5     2020-08-30 [1] CRAN (R 4.3.2)
#>  bitops             1.0-7     2021-04-24 [1] CRAN (R 4.3.1)
#>  blob               1.2.4     2023-03-17 [1] CRAN (R 4.3.2)
#>  broom              1.0.5     2023-06-09 [1] CRAN (R 4.3.2)
#>  bslib              0.7.0     2024-03-29 [1] CRAN (R 4.3.3)
#>  cachem             1.0.8     2023-05-01 [1] CRAN (R 4.3.2)
#>  car                3.1-2     2023-03-30 [1] CRAN (R 4.3.2)
#>  carData            3.0-5     2022-01-06 [1] CRAN (R 4.3.2)
#>  cli                3.6.2     2023-12-11 [1] CRAN (R 4.3.2)
#>  clusterProfiler    4.10.1    2024-03-08 [1] Bioconductor 3.18 (R 4.3.3)
#>  codetools          0.2-19    2023-02-01 [2] CRAN (R 4.3.2)
#>  colorspace         2.1-0     2023-01-23 [1] CRAN (R 4.3.2)
#>  cowplot            1.1.3     2024-01-22 [1] CRAN (R 4.3.2)
#>  crayon             1.5.2     2022-09-29 [1] CRAN (R 4.3.2)
#>  data.table         1.14.10   2023-12-08 [1] CRAN (R 4.3.2)
#>  DBI                1.2.2     2024-02-16 [1] CRAN (R 4.3.3)
#>  desc               1.4.3     2023-12-10 [1] CRAN (R 4.3.2)
#>  digest             0.6.33    2023-07-07 [1] CRAN (R 4.3.2)
#>  DOSE               3.28.2    2023-12-10 [1] Bioconductor
#>  dplyr            * 1.1.4     2023-11-17 [1] CRAN (R 4.3.2)
#>  EnhancedVolcano    1.20.0    2023-10-24 [1] Bioconductor
#>  enrichplot         1.22.0    2023-10-24 [1] Bioconductor
#>  evaluate           0.23      2023-11-01 [1] CRAN (R 4.3.2)
#>  factoextra         1.0.7     2020-04-01 [1] CRAN (R 4.3.2)
#>  fansi              1.0.6     2023-12-08 [1] CRAN (R 4.3.2)
#>  farver             2.1.1     2022-07-06 [1] CRAN (R 4.3.2)
#>  fastmap            1.1.1     2023-02-24 [1] CRAN (R 4.3.2)
#>  fastmatch          1.1-4     2023-08-18 [1] CRAN (R 4.3.1)
#>  fgsea              1.28.0    2023-10-24 [1] Bioconductor
#>  forcats          * 1.0.0     2023-01-29 [1] CRAN (R 4.3.2)
#>  fs                 1.6.3     2023-07-20 [1] CRAN (R 4.3.2)
#>  generics           0.1.3     2022-07-05 [1] CRAN (R 4.3.2)
#>  GenomeInfoDb       1.38.5    2023-12-28 [1] Bioconductor 3.18 (R 4.3.2)
#>  GenomeInfoDbData   1.2.11    2024-01-08 [1] Bioconductor
#>  ggforce            0.4.1     2022-10-04 [1] CRAN (R 4.3.2)
#>  ggfortify        * 0.4.16    2023-03-20 [1] CRAN (R 4.3.2)
#>  ggfun              0.1.4     2024-01-19 [1] CRAN (R 4.3.2)
#>  ggplot2          * 3.4.4     2023-10-12 [1] CRAN (R 4.3.2)
#>  ggplotify          0.1.2     2023-08-09 [1] CRAN (R 4.3.2)
#>  ggpubr             0.6.0     2023-02-10 [1] CRAN (R 4.3.2)
#>  ggraph             2.1.0     2022-10-09 [1] CRAN (R 4.3.2)
#>  ggrepel            0.9.4     2023-10-13 [1] CRAN (R 4.3.2)
#>  ggsignif           0.6.4     2022-10-13 [1] CRAN (R 4.3.2)
#>  ggtree             3.10.0    2023-10-24 [1] Bioconductor
#>  glue               1.6.2     2022-02-24 [1] CRAN (R 4.3.2)
#>  GO.db              3.18.0    2024-01-08 [1] Bioconductor
#>  GOSemSim           2.28.0    2023-10-24 [1] Bioconductor
#>  graphlayouts       1.0.2     2023-11-03 [1] CRAN (R 4.3.2)
#>  gridExtra          2.3       2017-09-09 [1] CRAN (R 4.3.2)
#>  gridGraphics       0.5-1     2020-12-13 [1] CRAN (R 4.3.2)
#>  gson               0.1.0     2023-03-07 [1] CRAN (R 4.3.2)
#>  gtable             0.3.4     2023-08-21 [1] CRAN (R 4.3.2)
#>  hash               2.2.6.3   2023-08-19 [1] CRAN (R 4.3.2)
#>  HDO.db             0.99.1    2024-01-08 [1] Bioconductor
#>  highr              0.10      2022-12-22 [1] CRAN (R 4.3.2)
#>  hms                1.1.3     2023-03-21 [1] CRAN (R 4.3.2)
#>  htmltools          0.5.8.1   2024-04-04 [1] CRAN (R 4.3.3)
#>  httr               1.4.7     2023-08-15 [1] CRAN (R 4.3.2)
#>  igraph             1.6.0     2023-12-11 [1] CRAN (R 4.3.2)
#>  inflection         1.3.6     2022-06-15 [1] CRAN (R 4.3.2)
#>  IRanges            2.36.0    2023-10-24 [1] Bioconductor
#>  jquerylib          0.1.4     2021-04-26 [1] CRAN (R 4.3.2)
#>  jsonlite           1.8.8     2023-12-04 [1] CRAN (R 4.3.2)
#>  kableExtra         1.4.0     2024-01-24 [1] CRAN (R 4.3.2)
#>  KEGGREST           1.42.0    2023-10-24 [1] Bioconductor
#>  knitr              1.46      2024-04-06 [1] CRAN (R 4.3.3)
#>  labeling           0.4.3     2023-08-29 [1] CRAN (R 4.3.1)
#>  lattice            0.21-9    2023-10-01 [2] CRAN (R 4.3.2)
#>  lazyeval           0.2.2     2019-03-15 [1] CRAN (R 4.3.2)
#>  lifecycle          1.0.4     2023-11-07 [1] CRAN (R 4.3.2)
#>  lubridate        * 1.9.3     2023-09-27 [1] CRAN (R 4.3.2)
#>  magrittr           2.0.3     2022-03-30 [1] CRAN (R 4.3.2)
#>  MASS               7.3-60    2023-05-04 [2] CRAN (R 4.3.2)
#>  Matrix             1.6-1.1   2023-09-18 [2] CRAN (R 4.3.2)
#>  memoise            2.0.1     2021-11-26 [1] CRAN (R 4.3.2)
#>  MetaProViz       * 2.0.0     2024-04-18 [1] Github (saezlab/MetaProViz@5a365e2)
#>  munsell            0.5.1     2024-04-01 [1] CRAN (R 4.3.3)
#>  nlme               3.1-163   2023-08-09 [2] CRAN (R 4.3.2)
#>  patchwork          1.2.0     2024-01-08 [1] CRAN (R 4.3.2)
#>  pheatmap           1.0.12    2019-01-04 [1] CRAN (R 4.3.2)
#>  pillar             1.9.0     2023-03-22 [1] CRAN (R 4.3.2)
#>  pkgconfig          2.0.3     2019-09-22 [1] CRAN (R 4.3.2)
#>  pkgdown            2.0.7     2022-12-14 [1] CRAN (R 4.3.2)
#>  plyr               1.8.9     2023-10-02 [1] CRAN (R 4.3.2)
#>  png                0.1-8     2022-11-29 [1] CRAN (R 4.3.1)
#>  polyclip           1.10-6    2023-09-27 [1] CRAN (R 4.3.1)
#>  purrr            * 1.0.2     2023-08-10 [1] CRAN (R 4.3.2)
#>  qcc                2.7       2017-07-11 [1] CRAN (R 4.3.2)
#>  qvalue             2.34.0    2023-10-24 [1] Bioconductor
#>  R6                 2.5.1     2021-08-19 [1] CRAN (R 4.3.2)
#>  ragg               1.2.7     2023-12-11 [1] CRAN (R 4.3.2)
#>  rappdirs           0.3.3     2021-01-31 [1] CRAN (R 4.3.2)
#>  RColorBrewer       1.1-3     2022-04-03 [1] CRAN (R 4.3.1)
#>  Rcpp               1.0.11    2023-07-06 [1] CRAN (R 4.3.2)
#>  RCurl              1.98-1.13 2023-11-02 [1] CRAN (R 4.3.2)
#>  readr            * 2.1.4     2023-02-10 [1] CRAN (R 4.3.2)
#>  reshape2           1.4.4     2020-04-09 [1] CRAN (R 4.3.2)
#>  rlang              1.1.2     2023-11-04 [1] CRAN (R 4.3.2)
#>  rmarkdown          2.26      2024-03-05 [1] CRAN (R 4.3.3)
#>  RSQLite            2.3.4     2023-12-08 [1] CRAN (R 4.3.2)
#>  rstatix            0.7.2     2023-02-01 [1] CRAN (R 4.3.2)
#>  rstudioapi         0.16.0    2024-03-24 [1] CRAN (R 4.3.3)
#>  S4Vectors          0.40.2    2023-11-23 [1] Bioconductor
#>  sass               0.4.9     2024-03-15 [1] CRAN (R 4.3.3)
#>  scales             1.3.0     2023-11-28 [1] CRAN (R 4.3.2)
#>  scatterpie         0.2.1     2023-06-07 [1] CRAN (R 4.3.2)
#>  sessioninfo        1.2.2     2021-12-06 [1] CRAN (R 4.3.2)
#>  shadowtext         0.1.3     2024-01-19 [1] CRAN (R 4.3.2)
#>  stringi            1.8.3     2023-12-11 [1] CRAN (R 4.3.2)
#>  stringr          * 1.5.1     2023-11-14 [1] CRAN (R 4.3.2)
#>  svglite            2.1.3     2023-12-08 [1] CRAN (R 4.3.2)
#>  systemfonts        1.0.5     2023-10-09 [1] CRAN (R 4.3.2)
#>  textshaping        0.3.7     2023-10-09 [1] CRAN (R 4.3.2)
#>  tibble           * 3.2.1     2023-03-20 [1] CRAN (R 4.3.2)
#>  tidygraph          1.3.0     2023-12-18 [1] CRAN (R 4.3.2)
#>  tidyr            * 1.3.0     2023-01-24 [1] CRAN (R 4.3.2)
#>  tidyselect         1.2.1     2024-03-11 [1] CRAN (R 4.3.3)
#>  tidytree           0.4.6     2023-12-12 [1] CRAN (R 4.3.2)
#>  tidyverse        * 2.0.0     2023-02-22 [1] CRAN (R 4.3.2)
#>  timechange         0.2.0     2023-01-11 [1] CRAN (R 4.3.2)
#>  treeio             1.26.0    2023-10-24 [1] Bioconductor
#>  tweenr             2.0.2     2022-09-06 [1] CRAN (R 4.3.2)
#>  tzdb               0.4.0     2023-05-12 [1] CRAN (R 4.3.2)
#>  utf8               1.2.4     2023-10-22 [1] CRAN (R 4.3.2)
#>  vctrs              0.6.5     2023-12-01 [1] CRAN (R 4.3.2)
#>  viridis            0.6.5     2024-01-29 [1] CRAN (R 4.3.2)
#>  viridisLite        0.4.2     2023-05-02 [1] CRAN (R 4.3.2)
#>  withr              3.0.0     2024-01-16 [1] CRAN (R 4.3.2)
#>  xfun               0.43      2024-03-25 [1] CRAN (R 4.3.3)
#>  xml2               1.3.6     2023-12-04 [1] CRAN (R 4.3.2)
#>  XVector            0.42.0    2023-10-24 [1] Bioconductor
#>  yaml               2.3.8     2023-12-11 [1] CRAN (R 4.3.2)
#>  yulab.utils        0.1.4     2024-01-28 [1] CRAN (R 4.3.2)
#>  zlibbioc           1.48.0    2023-10-24 [1] Bioconductor
#> 
#>  [1] C:/Users/chris/AppData/Local/R/win-library/4.3
#>  [2] C:/Program Files/R/R-4.3.2/library
#> 
#> ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

Badia-I-Mompel, Pau, Jesús Vélez Santiago, Jana Braunger, Celina Geiss, Daniel Dimitrov, Sophia Müller-Dott, Petr Taus, et al. 2022. “decoupleR: Ensemble of Computational Methods to Infer Biological Activities from Omics Data.” Bioinformatics Advances 2 (1): vbac016. https://doi.org/10.1093/bioadv/vbac016.

Bijlsma, Sabina, Ivana Bobeldijk, Elwin R Verheij, Raymond Ramaker, Sunil Kochhar, Ian A Macdonald, Ben van Ommen, and Age K Smilde. 2006. “Large-Scale Human Metabolomics Studies: A Strategy for Data (Pre-) Processing and Validation.” Analytical Chemistry 78 (2): 567–74. https://doi.org/10.1021/ac051495j.

Hotelling, Harold. 1931. “The Generalization of Student’s Ratio.” The Annals of Mathematical Statistics 2 (3): 360–78. https://doi.org/10.1214/aoms/1177732979.

Kanehisa, M, and S Goto. 2000. “KEGG: Kyoto Encyclopedia of Genes and Genomes.” Nucleic Acids Research 28 (1): 27–30. https://doi.org/10.1093/nar/28.1.27.

Wei, Runmin, Jingye Wang, Mingming Su, Erik Jia, Shaoqiu Chen, Tianlu Chen, and Yan Ni. 2018. “Missing Value Imputation Approach for Mass Spectrometry-Based Metabolomics Data.” Scientific Reports 8 (1): 663. https://doi.org/10.1038/s41598-017-19120-0.

Wulff, Jacob E., and Matthew W. Mitchell. 2018. “A Comparison of Various Normalization Methods for LC/MS Metabolomics Data.” Advances in Bioscience and Biotechnology 09 (08): 339–51. https://doi.org/10.4236/abb.2018.98022.