MetaPhOR was developed to enable users to assess metabolic dysregulation using transcriptomic-level data (RNA-sequencing and Microarray data) and produce publication-quality figures. A list of differentially expressed genes (DEGs), which includes fold change and p value, from DESeq2 (Love, Huber, and Anders (2014)) or limma (Ritchie et al. (2015)), can be used as input, with sample size for MetaPhOR, and will produce a data frame of scores for each KEGG pathway. These scores represent the magnitude and direction of transcriptional change within the pathway, along with estimated p-values (Rosario et al. (2018)). MetaPhOR then uses these scores to visualize metabolic profiles within and between samples through a variety of mechanisms, including: bubble plots, heatmaps, and pathway models.

This command line can be used to install and load MetaPhOR.

if (!require(“BiocManager”, quietly = TRUE))

install.packages(“BiocManager”)
BiocManager::install(“MetaPhOR”)

```
library(MetaPhOR)
## Warning: replacing previous import 'utils::findMatches' by
## 'S4Vectors::findMatches' when loading 'AnnotationDbi'
```

Minimal data preparation is required to run MetaPhOR. DEGs may be loaded into R in the form of .csv or .tsv for use with this package. The DEG file must contain columns for log fold change, adjusted p-value, and HUGO gene names. By default, MetaPhOR assumes DESeq2 header stylings: “log2FoldChange” and “padj”. In any function that assumes these headers, however, the user can define column names for these values. Below is a sample DEG file resulting from limma:

```
exdegs <- read.csv(system.file("extdata/exampledegs.csv",
package = "MetaPhOR"),
header = TRUE)
```

X | logFC | AveExpr | t | P.Value | adj.P.Val | B |
---|---|---|---|---|---|---|

TATDN1 | 0.0996779 | 6.112478 | 6.074387 | 0e+00 | 0.0000523 | 10.736686 |

ZNF706 | 0.1055866 | 6.446976 | 5.845307 | 0e+00 | 0.0000958 | 9.237467 |

DCAF13 | 0.0882855 | 6.393943 | 5.513711 | 1e-07 | 0.0002918 | 7.531541 |

TRMT12 | 0.1093866 | 6.071099 | 5.434470 | 1e-07 | 0.0003550 | 7.377948 |

IGF2BP1 | 0.9256589 | 4.454204 | 5.631783 | 0e+00 | 0.0002065 | 7.273836 |

MAP6D1 | 0.1619178 | 5.719428 | 5.342911 | 1e-07 | 0.0004612 | 7.090578 |

“pathwayAnalysis” first assigns scores and their absolute values using log fold change and p value to each gene (Rosario et al. (2018)). These transcript-level scores, along with sample size, are then utilized to calculate both scores (directional change) and absolute value scores (magnitude of change) (Rosario et al. (2018)) for each KEGG Pathway (Kanehisa and Goto (2000)). We then utilize a bootstrapping method, to randomly calculate 100,000 scores per pathway, based on the number of genes in that pathway and model the distribution. This distribution can then be used to evaluate where the actual score for that pathway sits in relation to the distribution, and can assign a p-value to the achieved score.

For example, if the polyamine biosynthetic pathway contains 13 genes, we can get a score for the sum of the 13 genes that exist within that pathway. Using bootstrapping, we can then randomly sample, with replacement, 13 genes to create scores, 100,000 times. We use these random samples (100,000) to generate a distribution, and we can calculate a p-value dependent on where the score that consists of the 13 genes that actually exist within the pathway falls within the distribution.

Taken together, the scores and p-values resulting from “pathwayAnalysis” provide a measure for both the biological and statistical significance of metabolic dysregulation.

**Note: A seed MUST be set before utilizing this function to ensure reproducible
results by bootstrapping. It is NECESSARY that the seed remain the same
throughout an analysis.**

**Note: Bootrapping is performed, by default, with 100,000 iterations of
resampling for optimal power. The number of iterations can be decreased for
improved speed. We use 50,000 iterations for improved speed of examples.***

pathwayAnalysis() requires:

- The file path to the DEG list of interest
- The name of the column containing HUGO gene names
- The sample size of the DEG analysis
- The number of iterations of resampling to be performed during bootstrapping
- Correct headers for fold change and p value columns (as indicated above)

A partial output of the pathway analysis function can be seen as follows:

```
set.seed(1234)
brca <- pathwayAnalysis(DEGpath = system.file("extdata/BRCA_DEGS.csv",
package = "MetaPhOR"),
genename = "X",
sampsize = 1095,
iters = 50000,
headers = c("logFC", "adj.P.Val"))
```

Scores | ABSScores | ScorePvals | ABSScorePvals | |
---|---|---|---|---|

Cardiolipin.Metabolism | 0.0293332 | 0.0299031 | 0.22996 | 0.60408 |

Cardiolipin.Biosynthesis | 0.0031671 | 0.0031671 | 0.43074 | 0.91878 |

Cholesterol.Biosynthesis | 0.2400430 | 0.2614076 | 0.06334 | 0.35952 |

Citric.Acid.Cycle | 0.0314352 | 0.1469312 | 0.48340 | 0.92848 |

Cyclooxygenase.Arachidonic.Acid.Metabolism | 0.0004634 | 0.0080444 | 0.52228 | 0.92286 |

Prostaglandin.Biosynthesis | -0.0392858 | 0.1046875 | 0.79502 | 0.35816 |

The metabolic profile determined by pathway analysis can be easily visualized using “bubblePlot.” Scores are plotted on the x-axis, while absolute value scores are plotted on the y-axis. Each point represents a KEGG pathway, where point size represents p-value (the smaller the p value, the larger the point) and point color is dictated by scores. Negative scores, which indicate transcriptional downregulation, are blue, and positive scores, which indicate transcriptional upregulation, are red. The top ten points, either by smallest p value or greatest dysregulation by score, are labeled with their pathway names. The plot demonstrates which pathways are the most statistically and biologically relevant.

bubblePlot() requires:

- The output of pathwayAnalysis(), as a data frame
- An indication which values to use, in order to label points: either “Pval” or “LogFC”
- Optional: Numeric value for point label text size (default = .25)

**Bubble Plot Labeled by P Value**

```
pval <- bubblePlot(scorelist = brca,
labeltext = "Pval",
labelsize = .85)
plot(pval)
```