## ----mms-case-study-1, include = FALSE----------------------------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    fig.align = "center"
)

.DT <- function(x) {
    if (knitr::is_html_output()) {
        DT::datatable(
            x,
            options = list(scrollX = TRUE, pageLength = 10, dom = "ftip"),
            rownames = FALSE
        )
    } else {
        knitr::kable(x)
    }
}

## ----setup, message = FALSE, warning=FALSE------------------------------------
library(struct)
library(MetMashR)
library(ggplot2)

## ----step0--------------------------------------------------------------------
# the studies to import
study_ids <- c("ST001039", "ST001048", "ST001317")

# read each source
sources <- lapply(
    study_ids,
    function(id) read_source(mwb_study_source(source = id, tag = id))
)

# use MetMashR to merge all sources
combine_step <- combine_sources(
    source_list = list(),
    keep_cols = "metabolite_name",
    source_col = "source"
)

# apply the merge
combine_step <- model_apply(combine_step, sources)

# extract the merged table
combined <- predicted(combine_step)

# report number of records
nrow(combined$data)

## ----cache-setup--------------------------------------------------------------
# Function that returns an rds_cache by name. Caches for steps shared by both
# variants are keyed by name/InChIKey, not by the service that produced the
# InChIKey. The chebi_name and cts_lite caches are specific to the ChEBI variant.
cached <- function(name, variant = FALSE) {
    prefix <- if (variant) "mms_case_study_chebi_variant_" else "mms_case_study_"
    rds_cache(source = file.path(
        system.file("cached", package = "MetMashR"),
        paste0(prefix, name, ".rds")
    ))
}

## ----step0.51-----------------------------------------------------------------
tidy_dictionary <- list(
    # trailing replicate/fragment label, e.g. "_R1", "_R2"
    list(pattern = "_[A-Za-z][0-9]+$", replace = ""),
    # remove short (<=3 character) batch/replicate codes
    list(pattern = "\\s\\([A-Za-z0-9]{1,3}\\)\\s*$", replace = "")
)

## ----step0.515----------------------------------------------------------------
# flag names reporting two "/"-joined candidate identifications
ambiguous_name_step <- compute_column(
    input_columns = "name_tidied",
    output_column = "ambiguous_name",
    # compute_column passes fcn a 1-column data.frame, not a bare vector
    fcn = function(x) grepl("\\)/", x[[1]])
)

# split into one record per candidate name, at a "/" immediately following a
# closing parenthesis; names without this pattern are left as a single record
split_ambiguous_name_step <- split_records(
    column_name = "name_tidied",
    separator = "(?<=\\))/"
)

## ----mms-case-study-2, step0.52-----------------------------------------------
# matches any trailing "(...)", not just abbrev(expansion). 
# e.g. "LysoPC (16:0)" -> "16:0". OK for this dataset; check before reusing 
# elsewhere
abbrev_dictionary <- list(
    list(
        pattern = "\\s\\([^()]+\\)\\s*$",
        replace = function(inp) {
            m <- regmatches(
                inp$x,
                regexec("\\s\\(([^()]+)\\)\\s*$", inp$x)
            )
            if (length(m[[1]]) > 1) trimws(m[[1]][2]) else inp$x
        }
    )
)

## ----mms-case-study-3, step0.53-----------------------------------------------
# dictionary to convert lyso-lipids to LIPID MAPS shorthand
lyso_dictionary <- list(
    list(pattern = "^Lyso-?PC", replace = "LPC", perl = TRUE, ignore.case = TRUE),
    list(pattern = "^Lyso-?PE", replace = "LPE", perl = TRUE, ignore.case = TRUE),
    list(pattern = "^Lyso-?PI", replace = "LPI", perl = TRUE, ignore.case = TRUE),
    list(pattern = "^Lyso-?PG", replace = "LPG", perl = TRUE, ignore.case = TRUE),
    list(pattern = "^Lyso-?PS", replace = "LPS", perl = TRUE, ignore.case = TRUE),
    list(pattern = "^Lyso-?PA", replace = "LPA", perl = TRUE, ignore.case = TRUE)
)

## ----mms-case-study-4, step0.54-----------------------------------------------
# workflow steps to apply the cleaning steps to the metabolite names
normalise_step <- list(
    normalise_strings(
        search_column = "metabolite_name",
        output_column = "name_tidied",
        dictionary = tidy_dictionary
    ),
    ambiguous_name_step,
    split_ambiguous_name_step,
    normalise_strings(
        search_column = "name_tidied",
        output_column = "name_normalised",
        dictionary = c(
            abbrev_dictionary,
            racemic_dictionary, # included in MetMashR
            lyso_dictionary
        )
    )
)

## ----step1--------------------------------------------------------------------
# query pubchem by synonym and return inchikey (shared)
pubchem_translate <- pubchem_id_exchange(
    query_column = "name_normalised",
    input_type = "synonyms",
    output_type = "inchikey",
    cache = cached("pubchem_id_exchange"),
    cache_mode = "offline"
)

# flag disagreement, prioritise the second service, drop records with no
# inchikey and merge duplicates. inchikey_col is the column returned by the
# second service and tag its name in InChIKey_source
translate_common <- function(inchikey_col, tag) {
    list(
        # flag where the first 14 characters of the inchikey match but the
        # first 26 do not (the final protonation character is not compared)
        compute_column(
            input_columns = c(inchikey_col, "inchikey_pubchem_id_exchange"),
            output_column = "ambiguous_inchikey",
            fcn = function(x) {
                a <- x[[1]]
                b <- x[[2]]
                ambiguous <- substr(a, 1, 14) == substr(b, 1, 14) &
                    substr(a, 1, 26) != substr(b, 1, 26)
                ambiguous[is.na(ambiguous)] <- FALSE
                ambiguous
            }
        ),
        # select the second service over pubchem where both exist
        prioritise_columns(
            column_names = c(inchikey_col, "inchikey_pubchem_id_exchange"),
            output_name = "InChIKey",
            source_name = "InChIKey_source",
            source_tags = c(tag, "pubchem")
        ),
        # exclude any results without an inchikey
        filter_na(column_name = "InChIKey", mode = "exclude"),
        # merge duplicate records
        combine_records(
            group_by = "InChIKey",
            default_fcn = fuse_unique(separator = " || "),
            fcns = list(
                metabolite_name = fuse_unique(separator = " || "),
                source = fuse_unique(separator = " || "),
                ambiguous_inchikey = function(x) any(x, na.rm = TRUE),
                ambiguous_name = function(x) any(x, na.rm = TRUE)
            )
        )
    )
}

# RefMet variant
translate_refmet <- c(
    list(
        pubchem_translate,
        # query refmet by synonym
        mwb_refmet_lookup(
            query_column = "name_normalised",
            cache = cached("refmet"),
            cache_mode = "offline",
            delay = 0.5
        )
    ),
    translate_common("inchi_key_refmet", "refmet")
)

# ChEBI variant
translate_chebi <- c(
    list(
        pubchem_translate,
        # query chebi by name and return a chebi accession, name and inchikey
        chebi_lookup(
            query_column = "name_normalised",
            search_by = "name",
            cache = cached("chebi_name", variant = TRUE),
            cache_mode = "offline",
            delay = 1
        )
    ),
    translate_common("inchikey_chebi", "chebi")
)

## ----step2--------------------------------------------------------------------
# PubChem properties and ClassyFire classes, shared by both variants
attribute_core <- list(
    # query the inchikey with pubchem and return various properties
    pubchem_property_lookup(
        query_column = "InChIKey",
        search_by = "inchikey",
        property = c(
            "IUPACName", "MolecularWeight", "MolecularFormula",
            "InChI", "CanonicalSMILES", "Charge"
        ),
        cache = cached("pubchem_property"),
        cache_mode = "offline",
        delay = 0.4
    ),
    # query the SMILES (from pubchem) with classyfire and return class
    # information.
    classyfire_batch_lookup(
        query_column = "ConnectivitySMILES_pubchem",
        output_items = c("kingdom", "superclass", "class"),
        output_fields = "name",
        cache = cached("classyfire"),
        cache_mode = "offline",
        delay = 3
    )
)

# RefMet variant
attribute_step <- c(
    attribute_core,
    list(
        # query the inchikey with MWB and return other cross reference
        # identifiers
        mwb_compound_lookup(
            input_item = "inchi_key",
            query_column = "InChIKey",
            output_item = c("hmdb_id", "chebi_id", "lm_id", "pubchem_cid"),
            cache = cached("mwb_compound"),
            cache_mode = "offline",
            delay = 0.5
        ),
        # query kegg with the ChEBI id and retrieve the kegg compound
        # identifier.
        kegg_lookup(
            get = "compound",
            from = "chebi",
            query_column = "chebi_id_mwb",
            cache = cached("kegg"),
            cache_mode = "offline"
        )
    )
)

# ChEBI variant: no MWB lookup. The ChEBI id comes from chebi_lookup() and
# CTS-Lite is added.
attribute_step_chebi <- c(
    attribute_core,
    list(
        # kegg_lookup() expects the bare ChEBI number, not "CHEBI:nnnn"
        compute_column(
            input_columns = "chebi_id_chebi",
            output_column = "chebi_number",
            fcn = function(x) gsub("CHEBI:", "", x[[1]], fixed = TRUE)
        ),
        kegg_lookup(
            get = "compound",
            from = "chebi",
            query_column = "chebi_number",
            cache = cached("kegg"),
            cache_mode = "offline"
        ),
        cts_lite_lookup(
            query_column = "InChIKey",
            cache = cached("cts_lite", variant = TRUE),
            cache_mode = "offline"
        )
    )
)

## ----step3--------------------------------------------------------------------
# create model sequence
curation_step <- list(
    # extract inchikey connectivity/skeleton only (first 14 characters) -
    # stereo-blind grouping, see prose above
    compute_column(
        input_columns = "InChIKey",
        output_column = "inchikey_skeleton",
        # compute_column passes fcn a 1-column data.frame, not a bare vector
        fcn = function(x) substr(x[[1]], 1, 14)
    ),
    # combine records with the same inchikey skeleton. ambiguous_inchikey
    # must also flag skeleton-only merges of records that still have
    # different full InChIKeys
    combine_records(
        group_by = "inchikey_skeleton",
        default_fcn = fuse_unique(separator = " || "),
        fcns = list(
            ambiguous_inchikey = function(x) {
                any(x, na.rm = TRUE) ||
                    dplyr::n_distinct(dplyr::pick(InChIKey)[[1]], na.rm = TRUE) > 1
            },
            ambiguous_name = function(x) any(x, na.rm = TRUE)
        )
    ),
    # remove any duplicate records
    unique_records()
)

## ----full-workflow,warning=FALSE----------------------------------------------
# create model sequences
workflow_refmet <- Reduce(
    `+`,
    c(normalise_step,
      translate_refmet,
      attribute_step,
      curation_step)
)

workflow_chebi <- Reduce(
    `+`,
    c(normalise_step,
      translate_chebi,
      attribute_step_chebi,
      curation_step)
)

# apply model sequences to merged sources and extract the results
result_refmet <- predicted(model_apply(workflow_refmet, combined))
result_chebi <- predicted(model_apply(workflow_chebi, combined))

# report number of records
c(
    refmet = nrow(result_refmet$data),
    chebi = nrow(result_chebi$data)
)

## ----result-table-------------------------------------------------------------
preview_cols <- c(
    "metabolite_name", "InChIKey", "IUPACName_pubchem",
    "MolecularFormula_pubchem", "MolecularWeight_pubchem",
    "hmdb_id_mwb", "chebi_id_mwb", "pubchem_cid_mwb", "compound_kegg", "lm_id_mwb",
    "kingdom.name_cfb", "superclass.name_cfb", "class.name_cfb",
    "ambiguous_inchikey", "ambiguous_name", "source"
)

## ----result-table-refmet------------------------------------------------------
.DT(head(result_refmet$data[, preview_cols], 10))

## ----result-table-chebi-------------------------------------------------------
chebi_cols <- c(
    "metabolite_name", "InChIKey", "IUPACName_pubchem",
    "MolecularFormula_pubchem", "MolecularWeight_pubchem",
    "chebi_id_chebi", "compound_kegg",
    "kingdom.name_cfb", "superclass.name_cfb", "class.name_cfb",
    "compound_name_cts", "annotation_type_count_cts",
    "ambiguous_inchikey", "ambiguous_name", "source"
)
.DT(head(result_chebi$data[, chebi_cols], 10))

## ----comparison, out.width = "70%", fig.width = 8, fig.height = 6-------------
# url to the content
zenodo_url <- paste0(
    "https://zenodo.org/api/records/8226097/files/",
    "MDM_Suppl_vSubmitted.xlsx/content"
)

# prepare BiocFileCache object
gold_db <- BiocFileCache_database(
    source = zenodo_url,
    resource_name = "MMS_zenodo_suppl",
    bfc_fun = cache_as_is,
    import_fun = function(path) {
        # the Zenodo API URL has no .xlsx extension; openxlsx needs one
        tmp <- tempfile(fileext = ".xlsx")
        file.copy(path, tmp, overwrite = TRUE)
        on.exit(unlink(tmp))
        openxlsx::read.xlsx(tmp, sheet = "Table S3", colNames = TRUE)
    }
)

# import as an annotation_source
gold <- read_source(gold_db)

# restrict to the three MWB studies, then derive the same
# inchikey skeleton used in Step 3
gold_workflow <- 
    filter_labels(
        column_name = "Reference", 
        labels = study_ids, mode = "include") + 
    compute_column(
        input_columns = "InChIKey", 
        output_column = "inchikey_skeleton",
        fcn = function(x) substr(x[[1]], 1, 14))

# apply the model to the imported source
gold <- predicted(model_apply(gold_workflow, gold))

# set labels that get used by the Venn chart
gold$name <- "MMS (Villalba et al.)"
gold$tag <- "MMS (Villalba et al.)"

# Both workflow results are already annotation_sources with an
# inchikey_skeleton column from Step 3 so they can be used directly
result_refmet$name <- "MetMashR (RefMet)"
result_refmet$tag <- "MetMashR (RefMet)"
result_chebi$name <- "MetMashR (ChEBI + CTS-Lite)"
result_chebi$tag <- "MetMashR (ChEBI + CTS-Lite)"

# compare the inchikey skeleton columns using a venn diagram
C <-
    annotation_venn_chart(
        factor_name = "inchikey_skeleton",
        labels = TRUE,
        legend = TRUE,
        fill_colour = ".group"
)
chart_plot(C, result_refmet, result_chebi, gold)

## ----mismatch-check, echo = FALSE---------------------------------------------
gold_keys <- unique(gold$data$inchikey_skeleton)
variant_keys <- list(
    "MetMashR (RefMet)" = unique(result_refmet$data$inchikey_skeleton),
    "MetMashR (ChEBI + CTS-Lite)" = unique(result_chebi$data$inchikey_skeleton)
)

overlap <- do.call(rbind, lapply(names(variant_keys), function(nm) {
    k <- variant_keys[[nm]]
    n_shared <- length(intersect(k, gold_keys))
    data.frame(
        variant = nm,
        n_skeletons = length(k),
        shared_with_published = n_shared,
        pct_of_variant = round(100 * n_shared / length(k), 1),
        pct_of_published = round(100 * n_shared / length(gold_keys), 1),
        only_in_variant = length(setdiff(k, gold_keys)),
        only_in_published = length(setdiff(gold_keys, k))
    )
}))

both <- Reduce(intersect, variant_keys)
either <- Reduce(union, variant_keys)

## ----mismatch-table, echo = FALSE---------------------------------------------
.DT(overlap)

## ----mms-case-study-5, warning=FALSE------------------------------------------
sessionInfo()