R/isotope_distributions.R

In jasenfinch/mzAnnotation: Tools for Putative Annotation of Electrospray Ionisation High Resolution Mass Spectrometry Data

#' Isotopic distribution calculator
#' @param MF a molecular formula
#' @param charge total charge
#' @param limit the minimum relative abundance threshold
#' @param element_table a table containing the elemental information. Defaults to `elements()`.
#' @return A tibble containing the isotope distribution for the molecular formula.
#' @examples 
#' isotopeDistribution(
#'   'C4H5O5',
#'   charge = -1)
#' @importFrom stringr str_split str_extract str_replace_all
#' @importFrom dplyr bind_cols group_by summarise right_join rename bind_rows desc
#' @export

isotopeDistribution <- function(MF,charge, limit = 0.00009 , element_table = elements()) {
  element_table <- element_table %>%
    mutate(Name = str_c(round(AtomicMass),Element))
  atomFrequencies <- tibble(Element = names(count.elements(MF)),Frequency = count.elements(MF))
  
  Isotopes <- filter(element_table,Element %in% atomFrequencies$Element, RelativeAbundance != 1) %>%
    mutate(Isotope = str_c(round(AtomicMass),Element))
  
  Isotopes <- map(unique(Isotopes$Element),~{
    element <- .
    numberEle <- atomFrequencies$Frequency[atomFrequencies$Element == element]
    availableIsotopes <- Isotopes %>%
      filter(Element == element) 
    t <- map(availableIsotopes$Isotope, ~{
      abund <- availableIsotopes$RelativeAbundance[availableIsotopes$Isotope == .]
      relativeAbundance <- 1
      for(i in 2:(numberEle + 1)) {
        relativeAbundance[i] <- relativeAbundance[i - 1] * abund * (numberEle - (i - 1) + 1) / (i - 1)
      }
      relativeAbundance <- relativeAbundance[-1]
      iso <- tibble(Element = rep(element,length(relativeAbundance)),`Relative Abundance` = relativeAbundance,Isotope = rep(.,length(relativeAbundance)), Frequency = 1:length(relativeAbundance)) %>%
        mutate(Name = str_c(Isotope,Frequency,sep = ' '))
    }) %>%
      bind_rows()
    return(t) 
  }) %>%
    bind_rows() %>%
    filter(`Relative Abundance` > limit)
  
  combinations <- map(2:length(unique(Isotopes$Isotope)),~{
    numberCombinations <- .
    com <- Isotopes$Name %>% 
      combn(m = numberCombinations) %>%
      t() %>%
      {
        suppressMessages(as_tibble(.,.name_repair = 'unique'))
      }
    differentCombinations <- map(com,~{
      unlist(map(str_split(.,' '),~{
        .[[1]][1]
      }))
    }) %>%
      bind_cols() %>%
      apply(1,function(x){
        T %in% (table(x) > 1)
      })
    com <- com %>%
      filter(differentCombinations == FALSE) %>%
      apply(1,str_c,collapse = '; ')
    return(com)
  }) %>%
    unlist() %>%
    map(~{
      isos <- unlist(str_split(.,'; '))
      isos <- Isotopes %>%
        filter(Name %in% isos)
      tibble(Name = ., `Relative Abundance` = prod(isos$`Relative Abundance`))
    }) %>%
    bind_rows() %>%
    filter(`Relative Abundance` > limit)
  
  Isotopes <- Isotopes %>%
    select(Name, `Relative Abundance`) %>%
    bind_rows(combinations) %>%
    rename(Isotope = Name) %>% 
    bind_rows(tibble(Isotope = NA, `Relative Abundance` = 1)) %>%
    arrange(desc(`Relative Abundance`)) %>%
    mutate(Probability = `Relative Abundance`/sum(`Relative Abundance`))
  
  mzs <- map(Isotopes$Isotope,~{
    if(is.na(.)){
      calcAccurateMass(MF,charge = charge)
    } else {
      Name <- unlist(map(str_split(., '; '),~{unlist(map(str_split(.,' '),~{.[1]}))}))
      Frequency <- as.integer(unlist(map(str_split(., '; '),~{unlist(map(str_split(.,' '),~{.[2]}))})))
      isoDat <- tibble(Name = Name,Frequency = Frequency) %>%
        mutate(Element = str_replace_all(Name,'[:digit:]','')) %>%
        left_join(select(element_table,Name,Element,AtomicMass),by = c('Name' = 'Name', 'Element' = 'Element')) %>%
        mutate(AtomicMass = AtomicMass * Frequency)
      isoMass <- sum(isoDat$AtomicMass)
      newMF <- isoDat %>%
        group_by(Element) %>%
        summarise(Frequency = sum(Frequency)) %>%
        right_join(atomFrequencies,by = c('Element' = 'Element'),suffix = c('1','2'))
      newMF$Frequency1[is.na(newMF$Frequency1)] <- 0
      newMF <- newMF %>%
        mutate(Frequency = Frequency2 - Frequency1,MF = str_c(Element,Frequency)) %>%
        select(MF) %>% 
        unlist() %>%
        str_c(collapse = '') %>%
        calcAccurateMass(charge = charge)
      mz <- isoMass + newMF
      
      return(mz)
    }
  }) %>%
    unlist() %>%
    round(5)
  
  Isotopes <- Isotopes %>%
    mutate(`m/z` = mzs) %>%
    select(Isotope,`m/z`,`Relative Abundance`,Probability)
  return(Isotopes)
}