## ----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)

## ----flow-legend, echo = FALSE, fig.wide = TRUE, out.width = "100%", fig.alt = "Legend: blue boxes are external lookups, peach boxes compute a column, red boxes filter records, yellow boxes merge records and green boxes split records."----
knitr::include_graphics("mms_figures/mms_legend.png")

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

# read the stored copy of each study table
sources <- lapply(study_ids, function(id) {
    read_source(rds_database(
        source = system.file(
            paste0("extdata/vignette_snapshots/MWB_", id, ".rds"),
            package = "MetMashR"
        ),
        tag = id,
        .writable = FALSE
    ))
})

# 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 the rds_cache for a lookup in one variant
# ("refmet" or "chebi"). Every lookup has a separate cache per variant so that
# the two workflows do not share any retrieved responses.
cached <- function(variant, name) {
    rds_cache(source = file.path(
        system.file("cached", package = "MetMashR"),
        paste0("mms_case_study_", variant, "_", 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 = "(?<=\\))/"
)

## ----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
        }
    )
)

## ----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)
)

## ----flow-normalise, echo = FALSE, fig.wide = TRUE, out.width = "100%", fig.cap = "Normalisation workflow, the same in both variants.", fig.alt = "Normalisation flow diagram: normalise_strings producing name_tidied, compute_column for ambiguous_name, split_records splitting name_tidied, and normalise_strings producing name_normalised."----
knitr::include_graphics("mms_figures/mms_normalise.png")

## ----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
        )
    )
)

## ----flow-step1, echo = FALSE, fig.wide = TRUE, out.width = "100%", fig.cap = "MMS Step 1 workflow. The variants differ only in the second name-translation lookup.", fig.alt = "Step 1 flow diagram: pubchem_id_exchange, then a split into mwb_refmet_lookup for the RefMet variant and chebi_lookup for the ChEBI variant, rejoining at compute_column for ambiguous_inchikey, then prioritise_columns, filter_na and combine_records."----
knitr::include_graphics("mms_figures/mms_step1_translation.png")

## ----step1--------------------------------------------------------------------
# query pubchem by synonym and return inchikey (same step in both variants,
# each with its own cache)
pubchem_translate <- function(variant) {
    pubchem_id_exchange(
        query_column = "name_normalised",
        input_type = "synonyms",
        output_type = "inchikey",
        cache = cached(variant, "pubchem_id_exchange"),
        cache_mode = "offline"
    )
}

# flag disagreement, prioritise refmet/chebi over pubchem, drop records with
# no inchikey and merge duplicates. inchikey_col is the column returned by
# refmet/chebi 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 refmet/chebi over pubchem where both exist (first column
        # listed takes priority)
        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("refmet"),
        # query refmet by synonym
        mwb_refmet_lookup(
            query_column = "name_normalised",
            cache = cached("refmet", "mwb_refmet"),
            cache_mode = "offline",
            delay = 0.5
        )
    ),
    translate_common("inchi_key_refmet", "refmet")
)

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

## ----flow-step2, echo = FALSE, fig.wide = TRUE, out.width = "100%", fig.cap = "MMS Step 2 workflow. The variants differ in how the ChEBI identifier for the KEGG lookup is obtained, and only the ChEBI variant queries CTS-Lite.", fig.alt = "Step 2 flow diagram: pubchem_property_lookup and classyfire_batch_lookup, then a split into mwb_compound_lookup for the RefMet variant and compute_column for chebi_number for the ChEBI variant, rejoining at kegg_lookup, followed by cts_lite_lookup for the ChEBI variant only."----
knitr::include_graphics("mms_figures/mms_step2_attributes.png")

## ----step2--------------------------------------------------------------------
# PubChem properties and ClassyFire classes (same steps in both variants,
# each with its own caches)
attribute_core <- function(variant) 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", "ConnectivitySMILES", "Charge"
        ),
        cache = cached(variant, "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(variant, "classyfire"),
        cache_mode = "offline",
        delay = 3
    )
)

# RefMet variant
attribute_step <- c(
    attribute_core("refmet"),
    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("refmet", "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("refmet", "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("chebi"),
    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("chebi", "kegg"),
            cache_mode = "offline"
        ),
        cts_lite_lookup(
            query_column = "InChIKey",
            cache = cached("chebi", "cts_lite"),
            cache_mode = "offline"
        )
    )
)

## ----flow-step3, echo = FALSE, fig.wide = TRUE, out.width = "100%", fig.cap = "MMS Step 3 workflow, the same in both variants.", fig.alt = "Step 3 flow diagram: compute_column for inchikey_skeleton, combine_records by inchikey_skeleton and unique_records."----
knitr::include_graphics("mms_figures/mms_step3_curation.png")

## ----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
applied_refmet <- model_apply(workflow_refmet, combined)
applied_chebi <- model_apply(workflow_chebi, combined)
result_refmet <- predicted(applied_refmet)
result_chebi <- predicted(applied_chebi)

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

## ----record-counts------------------------------------------------------------
record_counts <- function(W) {
    # output of the step before the InChIKey filter
    steps <- seq_len(length(W))
    i <- which(vapply(steps, function(k) is(W[k], "filter_na"), logical(1)))[1]
    looked_up <- predicted(W[i - 1])$data
    final <- predicted(W)$data
    keys <- unique(looked_up$InChIKey[!is.na(looked_up$InChIKey)])
    c(
        input_records = nrow(combined$data),
        names_looked_up = nrow(looked_up),
        names_without_inchikey = sum(is.na(looked_up$InChIKey)),
        distinct_full_inchikeys = length(keys),
        distinct_14_character_groups = nrow(final)
    )
}
data.frame(
    refmet = record_counts(applied_refmet),
    chebi = record_counts(applied_chebi)
)

## ----worked-example-----------------------------------------------------------
worked_example <- function(W, inchikey_col, tag) {
    steps <- seq_len(length(W))
    i <- which(vapply(steps, function(k) {
        is(W[k], "prioritise_columns")
    }, logical(1)))[1]
    # lookups before prioritising, and the selected InChIKey after
    before <- predicted(W[i - 1])$data
    after <- predicted(W[i])$data
    pc <- before$inchikey_pubchem_id_exchange
    other <- before[[inchikey_col]]
    # first example of each outcome
    pick <- c(
        both_agree = which(!is.na(pc) & !is.na(other) & pc == other &
            before$metabolite_name != before$name_normalised)[1],
        both_differ = which(!is.na(pc) & !is.na(other) & pc != other)[1],
        pubchem_only = which(!is.na(pc) & is.na(other))[1],
        other_only = which(is.na(pc) & !is.na(other))[1],
        unresolved = which(is.na(pc) & is.na(other))[1]
    )
    pick <- pick[!is.na(pick)]
    out <- data.frame(
        outcome = names(pick),
        original_name = before$metabolite_name[pick],
        name_sent = before$name_normalised[pick],
        inchikey_pubchem = pc[pick],
        inchikey_other = other[pick],
        selected = after$InChIKey[pick],
        selected_from = after$InChIKey_source[pick]
    )
    names(out)[5] <- paste0("inchikey_", tag)
    out
}
.DT(worked_example(applied_refmet, "inchi_key_refmet", "refmet"))

## ----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",
    "inchikey_cts", "match_level_cts",
    "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-------------
# stored copy of Table S3 from the supplementary file of Zenodo record 8226097
gold_db <- rds_database(
    source = system.file(
        "extdata/vignette_snapshots/MMS_zenodo_8226097_TableS3.rds",
        package = "MetMashR"
    ),
    .writable = FALSE
)

# 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()