inst/doc/MetCirc.R

In MetCirc: Navigating mass spectral similarity in high-resolution MS/MS metabolomics data

## ----knitr, include=FALSE, cache=FALSE----------------------------------------
library("knitr")

## ----eval=FALSE---------------------------------------------------------------
#  if (!requireNamespace("BiocManager", quietly=TRUE))
#      install.packages("BiocManager")
#  BiocManager::install("MetCirc")

## ----eval=TRUE----------------------------------------------------------------
library(MetCirc)

## ----eval=TRUE,echo=c(1:2)----------------------------------------------------
## load data from Li et al., 2015
data("sd02_deconvoluted", package = "MetCirc")

## load binnedMSP
data("sd01_outputXCMS", package = "MetCirc")
data("binnedMSP", package = "MetCirc")
## load similarityMat
data("similarityMat", package = "MetCirc")

## ----eval=TRUE----------------------------------------------------------------
## load idMSMStissueproject
data("idMSMStissueproject", package = "MetCirc")

## ----eval=TRUE,echo=c(1:2)----------------------------------------------------
## load data from GNPS
data("convertMsp2Spectra", package = "MetCirc")

## ----eval=TRUE----------------------------------------------------------------
## get all MS/MS spectra
id_uniq <- unique(sd02_deconvoluted[, "id"])

## obtain precursor m/z from id_uniq
prec_mz_l <- lapply(strsplit(as.character(id_uniq), split=" _ "), "[", 2)
prec_mz_l <- lapply(prec_mz_l, as.numeric)

## obtain m/z from fragments per precursor m/z
mz_l <- lapply(id_uniq, function(x) sd02_deconvoluted[sd02_deconvoluted[, "id"] == x, "mz"])
## obtain corresponding intensity values 
int_l <- lapply(id_uniq, function(x) sd02_deconvoluted[sd02_deconvoluted[, "id"] == x, "intensity"])

## obtain retention time by averaging all retention time values
rt_l <- lapply(id_uniq, function(x) sd02_deconvoluted[sd02_deconvoluted[, "id"] == x, "rt"])
rt_l <- lapply(rt_l, mean)

## create list of spectrum2 objects
spN_l <- lapply(1:length(mz_l), function(x) new("Spectrum2", rt=rt_l[[x]], precursorMz=prec_mz_l[[x]], mz=mz_l[[x]], intensity=int_l[[x]]))

## combine list of spectrum2 objects to MSpectra object
spectra_li <- MSnbase::MSpectra(spN_l, elementMetadata=S4Vectors::DataFrame(show=rep(TRUE, length(spN_l)))) 

## ----eval=TRUE----------------------------------------------------------------
## get all MS/MS spectra
tissue <- tissue[tissue[, "id"] %in% compartmentTissue[, "mz_rt_pcgroup"], ]
id_uniq <- unique(tissue[, "id"])

## obtain precursor m/z from id_uniq
prec_mz_l <- lapply(strsplit(as.character(id_uniq), split="_"), "[", 1)
prec_mz_l <- lapply(prec_mz_l, as.numeric)

## obtain m/z from fragments per precursor m/z
mz_l <- lapply(id_uniq, function(x) tissue[tissue[, "id"] == x, "mz"])
## obtain corresponding intensity values

## obtain retention time by averaging all retention time values
int_l <- lapply(id_uniq, function(x) tissue[tissue[, "id"] == x, "intensity"])
rt_l <- lapply(id_uniq, function(x) tissue[tissue[, "id"] == x, "rt"])
rt_l <- lapply(rt_l, mean)

## create list of spectrum2 objects
spN_l <- lapply(1:length(mz_l), function(x) new("Spectrum2", rt=rt_l[[x]], precursorMz=prec_mz_l[[x]], mz=mz_l[[x]], intensity=int_l[[x]]))

## combine list of spectrum2 objects to MSpectra object, 
## use SPL, LIM, ANT, STY for further analysis
spectra_tissue <- MSnbase::MSpectra(spN_l, 
    elementMetadata=S4Vectors::DataFrame(compartmentTissue[, c("SPL", "LIM", "ANT", "STY")])) 

## ----eval=TRUE----------------------------------------------------------------
spectra_msp <- convertMsp2Spectra(msp2spectra)

## ----eval=TRUE----------------------------------------------------------------
## similarity Matrix 
similarityMat <- compare_Spectra(spectra_tissue[1:100,], fun=normalizeddotproduct, binSize=0.01)

## ----cluster,eval=TRUE,fig.show='hide'----------------------------------------
## load package amap
hClustMSP <- hcluster(similarityMat, method="spearman")
## visualise clusters
plot(hClustMSP, labels = FALSE, xlab="", sub="")

## ----truncateGroup1,eval=TRUE,echo=TRUE---------------------------------------
## define three clusters
cutTreeMSP <- cutree(hClustMSP, k=3)
## extract feature identifiers that belong to module 1
module1 <- names(cutTreeMSP)[as.vector(cutTreeMSP) == "1"] 
## create a new similarity matrix that contains only highly related features
similarityMat_module1 <- similarityMat[module1, module1]

## ----eval=TRUE----------------------------------------------------------------
linkDf <- createLinkDf(similarityMatrix=similarityMat, spectra=spectra_tissue,
    condition=c("SPL", "LIM", "ANT", "STY"), lower=0.5, upper=1)

## ----eval=FALSE---------------------------------------------------------------
#  selectedFeatures <- shinyCircos(similarityMat, spectra=spectra_tissue, condition=c("LIM", "ANT", "STY"))

## ----circos,eval=TRUE,results='hide',message=FALSE,fig.show='hide'------------
## order similarity matrix according to retention time
condition <- c("SPL", "LIM", "ANT", "STY")
simM <- orderSimilarityMatrix(similarityMatrix=similarityMat, 
    spectra=spectra_tissue, type="retentionTime", group=FALSE)
groupname <- rownames(simM)
inds <- MetCirc:::spectraCond(spectra_tissue, 
    condition=condition)
inds_match <- lapply(inds, function(x) {inds_match <- match(groupname, x)
    inds_match <- inds_match[!is.na(inds_match)]; x[inds_match]})
inds_cond <- lapply(seq_along(inds_match), 
    function(x) {
        if (length(inds_match[[x]]) > 0) {
            paste(condition[x], inds_match[[x]], sep="_")
        } else character()
})
inds_cond <- unique(unlist(inds_cond))

## create link matrix
linkDf <- createLinkDf(similarityMatrix=simM, spectra=spectra_tissue, 
    condition=c("SPL", "LIM", "ANT", "STY"), lower=0.9, upper=1)
## cut link matrix (here: only display links between groups)
linkDf_cut <- cutLinkDf(linkDf=linkDf, type="inter")

## set circlize paramters
circos.par(gap.degree=0, cell.padding=c(0, 0, 0, 0), 
            track.margin=c(0, 0))

## here set indSelected arbitrarily
indSelected <- 14
selectedFeatures <- inds_cond[indSelected]

## actual plotting
plotCircos(inds_cond, linkDf_cut, initialize=TRUE, featureNames=TRUE, 
    cexFeatureNames=0.2, groupSector=TRUE, groupName=FALSE, 
    links=FALSE, highlight=TRUE)

highlight(groupname=inds_cond, ind=indSelected, linkDf=linkDf_cut)

## plot without highlighting
plotCircos(inds_cond, linkDf_cut, initialize=TRUE, featureNames=TRUE, 
    cexFeatureNames=0.2, groupSector=TRUE, groupName=FALSE, links=TRUE, 
    highlight=FALSE)

## ----session,eval=TRUE,echo=FALSE---------------------------------------------
sessionInfo()