In MetClassNet/MetClassNetR: Create and analize networks built from metabolomics Data

knitr::opts_chunk$set(echo = TRUE)

Reading mzTab-M

library(MSnbase)

fl <- "/vol/bioinvindex/Submissions/MTBLS1582/MTBLS1582/Height_NRGneg_2020471011.mzTab"
mtbls1582mztab <- MzTab(fl)

## Extract objects from mzTab-M
mt <- unlist(metadata(mtbls1582mztab))
sm <- smallMolecules(mtbls1582mztab)
mf <- moleculeFeatures(mtbls1582mztab)

## Extract assay and file names
assaynames <- mt[grep("^assay\\[[0-9]*\\]$", names(mt))]
filenames <- mt[grep("^ms_run\\[[0-9]*\\]-location$", names(mt))]

## Extract feature definitions
featid <- mf$SMF_ID
rt <- mf$retention_time_in_seconds
mz <- mf$exp_mass_to_charge

## Extract some chemistry
smiles <- sm$smiles
inchi <- sm$inchi

## Extract data matrix
## and replace actual assay names

ecols <- grep("abundance_assay", names(mf))
e <- as.matrix(mf[, ecols])
colnames(e) <- assaynames[sub("abundance_", "", colnames(e))]

library(MSnbase) 
library(QFeatures) 

fl <- "/vol/bioinvindex/Submissions/MTBLS1582/MTBLS1582/Height_NRGneg_2020471011.mzTab"
mtbls1582mztab <- MzTab(fl)

## Extract objects from mzTab-M
mt <- unlist(metadata(mtbls1582mztab))
mf <- moleculeFeatures(mtbls1582mztab)

## Extract assay names
assaynames <- mt[grep("^assay\\[[0-9]*\\]$", names(mt))]

## Extract data matrix
## and replace actual assay names

ecols <- grep("abundance_assay", names(mf))

## Fix assay names with information from Metadata MTD section
colnames(mf)[ecols] <- assaynames[sub("abundance_", "", colnames(mf[, ecols]))]

## Convert to QFeature
data_qF <- readQFeatures(mf, ecol = ecols, fname = "SMF_ID", name = "pos")

#head(assay(data_qF[["pos"]]))
#head(rowData(data_qF[["pos"]]))

## MS/MS parts
library(MsBackendMsp)
be <- MsBackendMsp()

## Import a single file.
fl2 <- "/vol/bioinvindex/Submissions/MTBLS1582/MTBLS1582/20203121416_spectra_NRGneg.msp"
msms <- backendInitialize(be, fl2)

## Paranoid check
if (length(msms) != nrow(e)) {
  stop("Different number of features in MS1 and MS2")
}

write.csv(cbind("m/z"=mz, "dummy"=999, "RT"=rt)[1:111,],
          file="MTBLS1582-envihomolog.csv",
          row.names=FALSE)

# 
# https://www.envihomolog.eawag.ch/
# F=results-peaks-MTBLS1582-envihomolog.csv cat $F | grep -v '"0"'| egrep 'mz|["/]611["/]' | cut -d "," -f 2 --complement  | column -t -s, | tr -d \" | expand

# F=results-peaks-MTBLS1582-envihomolog.csv cat $F | grep -v '"0"'| egrep 'mz|["/]611["/]' | cut -d "," -f 2 --complement  | column -t -s, | tr -d \"

library(nontarget)

# (0.2) list of adducts/isotopes - package enviPat ############
data(adducts);
data(isotopes);

## Use MTBLS1582 peaklist
peaklist <- data.frame(mass=mz, intensity=e[,1], rt=rt)

######################################################
# (2.1) run grouping of peaks for different adducts ##
# of the same candidate molecule #####################
adduct<-adduct.search(peaklist,
  adducts,
  rttol=2,
  mztol=5, ppm=TRUE,
  use_adducts=c("M-H", "M+FA-H", "M+Hac-H"), 
                ion_mode="negative");
# (2.2) plot results #################################
plotadduct(adduct);

# (4.1) Screen for homologue series ##################
homol <- homol.search(peaklist,
                      isotopes,     elements=c("C","H","O"), use_C=TRUE,
                      minmz=5,  maxmz=50,
                      minrt=5,  maxrt=60,
                      ppm=TRUE,
                      mztol=5,  rttol=5,
                      minlength=4,
                      mzfilter=FALSE,
                      spar=.45,     R2=.98,
                      plotit=FALSE)
# (4.2) Plot results #################################
plothomol(homol,xlim=FALSE,ylim=FALSE,plotlegend=TRUE);

library(tibble)
library(dplyr)
library(igraph)

library(visNetwork) # for visNetwork() and friends
library(networkD3)  # for saveNetwork()

nl <- as_tibble(peaklist)

el <- as_tibble(homol[["Peaks in homologue series"]]) %>% 
    filter(`to ID` != "0") %>% select (c("peak ID", "to ID", 
                                         "m/z increment", "RT increment"))

g <- el %>% select (c("peak ID", "to ID")) %>% as.matrix %>% graph_from_edgelist

## Write to File
write_graph(g, "/tmp/g.gml", "gml")

## Some Plotting
data <- toVisNetworkData(g)
vn <- visNetwork(nodes = data$nodes, 
                 edges = data$edges)
vn
saveNetwork(vn, "/tmp/vn.html")

sessionInfo()

MetClassNet/MetClassNetR documentation built on June 30, 2023, 2:12 p.m.