R/annotate_spectra.R

In taxonomicallyinformedannotation/tima-r: Taxonomically Informed Metabolite Annotation

#' @title Annotate spectra
#'
#' @description This function annotates spectra
#'
#' @details It takes two files as input.
#'    A query file that will be matched against a library file.
#'
#' @include get_params.R
#' @include harmonize_adducts.R
#' @include import_spectra.R
#'
#' @param input Query file containing spectra. Currently an '.mgf' file
#' @param library Library containing spectra to match against.
#'    Can be '.mgf' or '.sqlite' (Spectra formatted)
#' @param polarity MS polarity. Must be 'pos' or 'neg'.
#' @param output Output file.
#' @param threshold Minimal similarity to report
#' @param ppm Relative ppm tolerance to be used
#' @param dalton Absolute Dalton tolerance to be used
#' @param qutoff Intensity under which ms2 fragments will be removed.
#' @param approx Perform matching without precursor match
#'
#' @return NULL
#'
#' @export
#'
#' @examples
#' \dontrun{
#' tima:::copy_backbone()
#' go_to_cache()
#' get_file(
#'   url = get_default_paths()$urls$examples$spectra_mini,
#'   export = get_params(step = "annotate_spectra")$files$spectral$raw
#' )
#' get_file(
#'   url = get_default_paths()$urls$examples$spectral_lib_mini$with_rt,
#'   export = get_default_paths()$data$source$libraries$spectra$exp$with_rt
#' )
#' annotate_spectra(
#'   library = get_default_paths()$data$source$libraries$spectra$exp$with_rt
#' )
#' unlink("data", recursive = TRUE)
#' }
annotate_spectra <- function(input = get_params(step = "annotate_spectra")$files$spectral$raw,
                             library = get_params(step = "annotate_spectra")$files$libraries$spectral,
                             polarity = get_params(step = "annotate_spectra")$ms$polarity,
                             output = get_params(step = "annotate_spectra")$files$annotations$raw$spectral$spectral,
                             threshold = get_params(step = "annotate_spectra")$annotations$thresholds$ms2$similarity$annotation,
                             ppm = get_params(step = "annotate_spectra")$ms$tolerances$mass$ppm$ms2,
                             dalton = get_params(step = "annotate_spectra")$ms$tolerances$mass$dalton$ms2,
                             qutoff = get_params(step = "annotate_spectra")$ms$thresholds$ms2$intensity,
                             approx = get_params(step = "annotate_spectra")$annotations$ms2approx) {
  stopifnot("Your input file does not exist." = file.exists(input))
  stopifnot("Polarity must be 'pos' or 'neg'." = polarity %in% c("pos", "neg"))
  ## Check if library file(s) exists
  stopifnot("Library file(s) do(es) not exist" = all(purrr::map(.x = library, .f = file.exists) |>
    unlist()))

  ## Not checking for ppm and Da limits, everyone is free.

  if (length(library) > 1) {
    library <- library[grepl(polarity, library, fixed = TRUE)]
  }

  log_debug("Loading spectra...")
  spectra <- input |>
    import_spectra(
      cutoff = qutoff,
      dalton = dalton,
      polarity = polarity,
      ppm = ppm
    )

  df_empty <- data.frame(
    feature_id = NA,
    candidate_spectrum_entropy = NA,
    candidate_adduct = NA,
    candidate_library = NA,
    candidate_structure_error_mz = NA,
    candidate_structure_name = NA,
    candidate_structure_inchikey_no_stereo = NA,
    candidate_structure_smiles_no_stereo = NA,
    candidate_structure_molecular_formula = NA,
    candidate_structure_exact_mass = NA,
    candidate_structure_xlogp = NA,
    candidate_score_similarity = NA,
    candidate_count_similarity_peaks_matched = NA
  )

  if (length(spectra) > 0) {
    log_debug("Loading spectral library")
    spectral_library <- unlist(library) |>
      purrr::map(
        .f = import_spectra,
        cutoff = qutoff,
        dalton = dalton,
        polarity = polarity,
        ppm = ppm,
        sanitize = FALSE
      ) |>
      purrr::map(
        .f = Spectra::applyProcessing,
        BPPARAM = BiocParallel::SerialParam()
      ) |>
      Spectra::concatenateSpectra()

    query_precursors <- spectra@backend@spectraData$precursorMz
    query_spectra <- spectra@backend@peaksData
    ## ISSUE see #148 find a way to have consistency in spectrum IDs
    query_ids <- spectra@backend@spectraData$acquisitionNum
    if (is.null(query_ids)) {
      query_ids <- spectra@backend@spectraData$spectrum_id
    }
    if (is.null(query_ids)) {
      query_ids <- spectra@backend@spectraData$SLAW_ID
    }
    if (is.null(query_ids)) {
      query_ids <- spectra@backend@spectraData$FEATURE_ID
    }
    rm(spectra)

    ## Fix needed
    lib_precursors <- spectral_library@backend@spectraData |>
      tidytable::transmute(precursor = tidytable::coalesce(tidytable::across(tidytable::any_of(
        c("precursorMz", "precursor_mz")
      )))) |>
      tidytable::pull()
    minimal <- pmin(lib_precursors - dalton, lib_precursors * (1 - (10^
      -6 * ppm)))
    maximal <- pmax(lib_precursors + dalton, lib_precursors * (1 + (10^
      -6 * ppm)))

    if (approx == FALSE) {
      log_debug("Reducing library size...")
      df_3 <- dplyr::inner_join(
        tidytable::tidytable(minimal, maximal, lib_precursors),
        tidytable::tidytable(val = unique(query_precursors)),
        by = dplyr::join_by(minimal <= val, maximal >= val)
      ) |>
        tidytable::distinct(minimal, .keep_all = TRUE)

      spectral_library <-
        spectral_library[lib_precursors %in% df_3$lib_precursors]

      ## Fix needed
      lib_precursors <- spectral_library@backend@spectraData |>
        tidytable::transmute(precursor = tidytable::coalesce(tidytable::across(tidytable::any_of(
          c("precursorMz", "precursor_mz")
        )))) |>
        tidytable::pull()
      minimal <- pmin(lib_precursors - dalton, lib_precursors * (1 - (10^
        -6 * ppm)))
      maximal <- pmax(lib_precursors + dalton, lib_precursors * (1 + (10^
        -6 * ppm)))
      rm(df_3)
    }

    lib_ids <- seq_along(spectral_library)
    spectral_library$spectrum_id <- lib_ids
    lib_spectra <- spectral_library@backend@peaksData
    safety <- lib_spectra[purrr::map(.x = lib_spectra, .f = length) != 0]
    if (length(safety) != 0) {
      log_debug("Annotating...")
      df_final <-
        calculate_entropy_score(
          lib_ids = lib_ids,
          lib_precursors = lib_precursors,
          lib_spectra = lib_spectra,
          query_ids = query_ids,
          query_precursors = query_precursors,
          query_spectra = query_spectra,
          dalton = dalton,
          ppm = ppm,
          threshold = threshold,
          approx = approx
        ) |>
        tidytable::as_tidytable()

      lib_adduct <- spectral_library@backend@spectraData$adduct
      if (is.null(lib_adduct)) {
        lib_adduct <- rep(NA_character_, length(spectral_library))
      }
      lib_inchikey <- spectral_library@backend@spectraData$inchikey
      if (is.null(lib_inchikey)) {
        lib_inchikey <- rep(NA_character_, length(spectral_library))
      }
      lib_inchikey2D <-
        spectral_library@backend@spectraData$inchikey_2D
      if (is.null(lib_inchikey2D)) {
        lib_inchikey2D <- rep(NA_character_, length(spectral_library))
      }
      lib_smiles <- spectral_library@backend@spectraData$smiles
      if (is.null(lib_smiles)) {
        lib_smiles <- rep(NA_character_, length(spectral_library))
      }
      lib_smiles2D <- spectral_library@backend@spectraData$smiles_2D
      if (is.null(lib_smiles2D)) {
        lib_smiles2D <- rep(NA_character_, length(spectral_library))
      }
      lib_library <- spectral_library@backend@spectraData$library
      if (is.null(lib_library)) {
        lib_library <- rep(NA_character_, length(spectral_library))
      }
      lib_mf <- spectral_library@backend@spectraData$formula
      if (is.null(lib_mf)) {
        lib_mf <- rep(NA_character_, length(spectral_library))
      }
      lib_mass <- spectral_library@backend@spectraData$exactmass
      if (is.null(lib_mass)) {
        lib_mass <- rep(NA_real_, length(spectral_library))
      }
      lib_name <- spectral_library@backend@spectraData$name
      if (is.null(lib_name)) {
        lib_name <- rep(NA_character_, length(spectral_library))
      }
      lib_xlogp <- spectral_library@backend@spectraData$xlogp
      if (is.null(lib_xlogp)) {
        lib_xlogp <- rep(NA_real_, length(spectral_library))
      }
      rm(spectral_library)

      df_meta <- tidytable::tidytable(
        "target_id" = lib_ids,
        "target_adduct" = lib_adduct,
        "target_inchikey" = lib_inchikey,
        "target_inchikey_no_stereo" = lib_inchikey2D,
        "target_smiles" = lib_smiles,
        "target_smiles_no_stereo" = lib_smiles2D,
        "target_library" = lib_library,
        "target_formula" = lib_mf,
        "target_exactmass" = lib_mass,
        "target_name" = lib_name,
        "target_xlogp" = lib_xlogp,
        "target_precursorMz" = lib_precursors
      )
      df_meta <- df_meta |>
        tima:::harmonize_adducts(adducts_colname = "target_adduct")
      rm(lib_precursors)

      df_final$candidate_spectrum_entropy <- as.numeric(df_final$candidate_spectrum_entropy)
      df_final$candidate_score_similarity <- as.numeric(df_final$candidate_score_similarity)
      df_final$candidate_count_similarity_peaks_matched <- as.integer(df_final$candidate_count_similarity_peaks_matched)

      df_final <- df_final |>
        tidytable::left_join(df_meta) |>
        tidytable::select(-target_id)

      df_final <- df_final |>
        tidytable::mutate(
          candidate_structure_error_mz = target_precursorMz - precursorMz,
          candidate_structure_inchikey_no_stereo = tidytable::if_else(
            condition = is.na(target_inchikey_no_stereo),
            true = target_inchikey |>
              gsub(
                pattern = "-.*",
                replacement = "",
                perl = TRUE
              ),
            false = target_inchikey_no_stereo
          ),
          candidate_structure_smiles_no_stereo = tidytable::coalesce(target_smiles_no_stereo, target_smiles)
        ) |>
        tidytable::select(tidyselect::any_of(
          c(
            "feature_id",
            "candidate_adduct" = "target_adduct",
            "candidate_library" = "target_library",
            "candidate_structure_error_mz",
            "candidate_structure_name" = "target_name",
            "candidate_structure_inchikey_no_stereo",
            "candidate_structure_smiles_no_stereo",
            "candidate_structure_molecular_formula" = "target_formula",
            "candidate_structure_exact_mass" = "target_exactmass",
            "candidate_structure_xlogp" = "target_xlogp",
            "candidate_spectrum_entropy",
            "candidate_score_similarity",
            "candidate_count_similarity_peaks_matched"
          )
        ))

      ## COMMENT AR: Not doing it because of thresholding
      ## df_final[is.na(df_final)] <- 0

      log_debug("Filtering results above threshold only...")
      df_final <- df_final |>
        tidytable::filter(candidate_score_similarity >= threshold) |>
        tidytable::arrange(tidytable::desc(candidate_score_similarity)) |>
        ## keep only the best result (per library for now)
        tidytable::distinct(
          feature_id,
          candidate_library,
          candidate_structure_inchikey_no_stereo,
          .keep_all = TRUE
        )

      log_debug(
        nrow(
          df_final |>
            ## else doesn't work if some are empty
            tidytable::distinct(
              candidate_structure_inchikey_no_stereo,
              candidate_structure_smiles_no_stereo
            )
        ),
        "Candidates were annotated on",
        nrow(df_final |>
          tidytable::distinct(feature_id)),
        "features, with at least",
        threshold,
        "similarity score."
      )
      if (nrow(df_final) == 0) {
        log_debug("No spectra were matched, returning an empty dataframe")
        df_final <- df_empty
      }
    } else {
      log_debug("No spectra left in the library,
              returning an empty dataframe")
      df_final <- df_empty
    }
    rm(
      query_precursors,
      query_spectra,
      query_ids,
      minimal,
      maximal
    )
  } else {
    log_debug("No spectra matched the given polarity,
              returning an empty dataframe")
    df_final <- df_empty
  }

  tima:::export_params(
    parameters = get_params(step = "annotate_spectra"),
    step = "annotate_spectra"
  )
  tima:::export_output(x = df_final, file = output[[1]])
  rm(df_final)

  return(output[[1]])
}