R/calc_dbe.R

In ume: Ultrahigh-Resolution Mass Spectrometry Data Evaluation for Complex Organic Matter

#' @title Calculate Double Bond Equivalent (DBE)
#'
#' @description Calculates the Double Bond Equivalent (DBE) for a given neutral molecular formula.
#' DBE is a measure of unsaturation, representing the total number of rings and pi bonds
#' in a molecule. This function uses the `masses` data table to determine valence information
#' for each element in the input molecular formula.
#'
#' @name calc_dbe
#' @family calculations
#' @inheritParams main_docu
# @param element_names Specify whether element symbols in mfd are in 'lower_case' (default) or 'upper_case'
#' @import data.table
#'
#' @return A numeric vector of the same length as the number of rows in `mfd`,
#' where each entry represents the calculated DBE for the corresponding molecular formula.
#' The result vector is named 'dbe'.
#'
#' @examples
#' # Example with user-defined data
#' calc_dbe("C6H10O6")
#' calc_dbe("C6H10Br2")
#' calc_dbe(c("C3[13C1]H10O4", "C6H10O6"))
#'
#' # Example with demo data from UME package
#' calc_dbe(mfd = mf_data_demo)
#'
#' @details
#' This function computes DBE based on the molecular formula specified in `mfd`.
#' `mfd` can be a data.table or a character string or character vector of molecular formula strings.
#'
#' For each isotope in the formula, DBE is calculated as the sum of (valence - 2)
#' multiplied by the count of that isotope, divided by 2, and then adding 1.
#' Elements with a valence of 2 are excluded from the DBE calculation.
#'
#' The function will stop and print an error if any elements in `mfd` have missing valence information
#' in `masses`.
#' @export

calc_dbe <- function(mfd,
                     masses = ume::masses,
                     verbose = FALSE,
                     ...) {

  full_name <- new_name <- isotope_name <- NULL

# Check if input is a character or character vector of molecular formulas
  if(is.character(mfd)){
    mfd <- convert_molecular_formula_to_data_table(mfd)
  }

  # Check for empty input
  if (nrow(mfd) == 0) {
    stop("Input data.table 'mfd' is empty. Please provide a non-empty data.table.")
  }

  isotope_info <- get_isotope_info(mfd)
  setnames(mfd, isotope_info$orig_name, isotope_info$label, skip_absent = TRUE)

# Create a list of element assignments
  # tmp <- melt(data = isotope_info, id.vars = c("symbol", "nm", "valence", "hill_order"),
  #             measure.vars = c("element", "symbol", "isotope", "label", "orig_name"),
  #             variable.name = "isotope_id", value.name = "isotope_name")
  # tmp$isotope_id <- NULL
  # tmp[, full_name:=paste0(nm, tolower(symbol))]
  # tmp <- unique(tmp)
  #
  # cols_to_update <- grep("^[0-9]", tolower(names(mfd)), value = TRUE, invert = TRUE)
  # cols_to_update <- cols_to_update[cols_to_update %in% tmp$isotope_name]
  #
  # if(length(cols_to_update)>0){
  # new_names <- tmp[isotope_name %in% cols_to_update, .(nm=min(nm)), isotope_name]
  # new_names[, new_name:=paste0(nm, isotope_name)]
  #
  # setnames(mfd, cols_to_update, new_names$new_name)
  # mfd
  # }

  # Test if all element names and valences exist in data.table "masses"
  # ume::masses supplies the (most common) valence of each isotope at normal conditions
  # Retrieve information on isotopes that are contained in mfd and
  # remove all elements with valence = 2 (doesn't affect dbe calculation)

  # tmp$isotope_name <- NULL
  # tmp <- unique(tmp)

  elements_included <- isotope_info[!is.na(valence) & valence != 2,]

    if(length(elements_included)==0 | nrow(isotope_info[is.na(valence)]) > 0)  {
      stop("Name or valence of elements in mfd not found in masses.Rdata (consider revising argument 'element_names').")
    } else {
      .msg("Included elements for DBE calculation:\n '%s'", paste(elements_included$label, collapse = "', '"))
    }

  # DBE calculation
    x_sum <- 0
    for (i in unname(unlist(elements_included[, label]))){
      x <- 0
      x = (mfd[, get(i)] * (elements_included[label == i, valence] - 2)) # first part of DBE calculation
      x_sum <- x_sum + x
    }

    dbe <- as.vector(x_sum/2 + 1) # second part of DBE calculation

    return(dbe)
  }