stat.getMSEA_Metabolon: Metabolite set enrichment analysis (MSEA) using pathway...

View source: R/stat.getMSEA_Metabolon.r

stat.getMSEA_MetabolonR Documentation

Metabolite set enrichment analysis (MSEA) using pathway knowledge curated by Metabolon

Description

A function that returns the pathway enrichment score for all perturbed metabolites in a patient's full metabolomic profile.

Usage

stat.getMSEA_Metabolon(data, class.labels, pathway_knowledgebase = "Metabolon")

Arguments

data

- A data matrix, where rows are metabolites and columns are sample IDs.

class.labels

- A vector of binary indicators, where 1 is phenotype1 and 0 is phenotype2.

pathway.knowledgebase

- The filename of the .gmt file associated with the pathway knowledge desired. Currently only "Metabolon" is offered, though "KEGG", "WikiPathways", "SMPDB" and/or "Reactome" can be added in future versions.

Examples

data(Miller2015)
# Get class labels for diagnostic class "diagClass"
diagClass = "Citrullinemia"
cohorts = list()
cohorts[[diagClass]] = colnames(Miller2015[,which(Miller2015[1,]==diagClass)])
cohorts[["ref"]] = colnames(Miller2015[,which(Miller2015[1,]=="No biochemical genetic diagnosis")])
class.labels = unlist(cohorts)
class.labels[-which(class.labels %in% cohorts[[diagClass]])] = 0
class.labels[which(class.labels %in% cohorts[[diagClass]])] = 1
class.labels = as.numeric(class.labels)

# Create the data matrix 
data_mx = Miller2015[-1,which(colnames(Miller2015) %in% unlist(cohorts)))]

# Generate a .gmt file.
population = rownames(data_mx)
paths.hsa = list.dirs(path=system.file("extdata", package="CTDext"), full.names = FALSE)
paths.hsa = paths.hsa[-which(paths.hsa %in% c("", "RData", "allPathways", "Pathway_GMTs"))]
sink(system.file("extdata/Pathway_GMTs/Metabolon.gmt", package="CTDext"))
for (p in 1:length(paths.hsa)) {
  load(system.file(sprintf("extdata/RData/%s.RData", paths.hsa[p]), package="CTDext"))
  pathway.compounds = V(ig)$label[which(V(ig)$shape=="circle")]
  pathCompIDs = unique(tolower(pathway.compounds[which(pathway.compounds %in% population)]))
  cat(sprintf("%s         %s\n", paths.hsa[p], paste(pathCompIDs, collapse="    ")))
}
sink()
print("test")
pathway.data = stat.getMSEA_Metabolon(data_mx, class.labels, pathway_knowledgebase = "Metabolon")

BRL-BCM/CTDext documentation built on May 7, 2022, 5:31 a.m.