R/NMR_metabolite_identification.R

In specmine: Metabolomics and Spectral Data Analysis and Mining

nmr_identification = function(dataset, ppm.tol, frequency_scores, solvent_scores, organism_scores,
                              method='Match_uniq', per.sample=FALSE, tresh_zero=0, alpha=10e-4,
                              verbose = TRUE){
  
  if(!method %in% c('Match_uniq', 'Hyper', 'Hyper_uniq')) {
    stop('Invalid method. Valid methods: Match_uniq, Hyper or Hyper_uniq')
  }
  
  if(is.null(frequency_scores) | is.null(solvent_scores) | is.null(organism_scores)) {
    stop('Arguments missing: frequency_scores, solvent_score or organism_scores. Please check help page of the function to know how to set them')
  }
  
  res = list()
  
  if(method %in% c('Match_uniq', 'Hyper_uniq')) {
    if (verbose) message("Calculating uniqueness scores...")
    uniq_scores = uniqueness_scores(0)
    if (verbose) message("Done")
  } else {
    uniq_scores = NULL
  }
  
  if(method %in% c('Hyper', 'Hyper_uniq')) {
    n_peaks_total = length(unique(unlist(spectra_list)))
  } else {
    n_peaks_total = NULL
  }
  
  if (verbose) message("Getting compounds for organism scores...")
  all_orgs = get_OrganismsCodes()
  cpds_groups = list()
  for(group in names(organism_scores)) {
    if (!group %in% c("other", "not_in_kegg")) {
      if (verbose) message("-- Getting compounds from ", group)
    }
    if(is_organism(group)) {
      cpds_groups[[group]] = compounds_in_organism(group)
    } else if (is_group(group, all_orgs)) {
      cpds_groups[[group]] = compounds_in_group(group, all_orgs)
    }
  }
  if (verbose) message("Done")
  
  if(!per.sample) {
    samples_peaks = as.numeric(rownames(dataset$data))
    return(
      identification_nmr_peaks(
        sample_peaks = samples_peaks,
        method = method,
        ppm.tol = ppm.tol,
        frequency_scores = frequency_scores,
        solvent_scores = solvent_scores,
        organism_scores = organism_scores,
        cpds_groups = cpds_groups,
        uniq_scores = uniq_scores,
        n_peaks_total = n_peaks_total,
        alpha = alpha,
        verbose = verbose
      )
    )
  }
  
  res = list()
  for(samp in colnames(dataset$data)) {
    if (verbose) message("-SAMPLE ", samp)
    sample_peaks = as.numeric(names(dataset$data[dataset$data[, samp] > tresh_zero, 1]))
    res[[samp]] = identification_nmr_peaks(
      sample_peaks = sample_peaks,
      method = method,
      ppm.tol = ppm.tol,
      frequency_scores = frequency_scores,
      solvent_scores = solvent_scores,
      organism_scores = organism_scores,
      cpds_groups = cpds_groups,
      uniq_scores = uniq_scores,
      n_peaks_total = n_peaks_total,
      alpha = alpha,
      verbose = verbose
    )
  }
  return(res)
}

identification_nmr_peaks = function(sample_peaks, method, ppm.tol, frequency_scores, solvent_scores,
                                    organism_scores, cpds_groups,
                                    uniq_scores = NULL, n_peaks_total = NULL, alpha = 10e-4,
                                    verbose = TRUE){
  
  if (verbose) message("Matching samples to reference peaks...")
  if(method == "Hyper") {
    res = Hyper_method(sample_peaks, ppm.tol, n_peaks_total, alpha)
  } else if(method == "Hyper_uniq") {
    res = Hyper_uniq_method(sample_peaks, ppm.tol, uniq_scores, n_peaks_total, alpha)
  } else {
    res = Match_uniq_method(sample_peaks, ppm.tol, uniq_scores)
  }
  if (verbose) message("Done")
  
  if (verbose) message("Calculating Final scores...")
  score = c()
  score_frequency = c()
  score_solvents = c()
  score_organisms = c()
  for(i in 1:dim(res$results_table)[1]) {
    reference_id = res$results_table$SPCMNS[i]
    
    score_fr = score_freq(reference_id, frequency_scores)
    score_solv = score_solvent(reference_id, solvent_scores)
    score_org = score_organism(res$results_table$SPCMNM[i], cpds_groups, organism_scores)
    score_final = (res$results_table$match_score[i] + score_fr + score_solv + score_org) / 4
    
    score = c(score, score_final)
    score_frequency = c(score_frequency, score_fr)
    score_solvents = c(score_solvents, score_solv)
    score_organisms = c(score_organisms, score_org)
  }
  
  if(method == 'Hyper') {
    res$results_table = cbind(res$results_table[, c('SPCMNM', 'Name', 'SPCMNS')],
                              score,
                              res$results_table[, c('match_score')],
                              score_frequency, score_solvents, score_organisms,
                              res$results_table[, c('n.peaks.matched', 'detailed_results_id')])
    colnames(res$results_table)[5] = 'match_score'
  }
  else if(method == 'Hyper_uniq') {
    res$results_table = cbind(res$results_table[, c('SPCMNM', 'Name', 'SPCMNS')],
                              score,
                              res$results_table[, c('match_score', 'hypergeometric_score', 'uniqueness_score')],
                              score_frequency, score_solvents, score_organisms,
                              res$results_table[, c('n.peaks.matched', 'detailed_results_id')])
  }
  else {
    res$results_table = cbind(res$results_table[, c('SPCMNM', 'Name', 'SPCMNS')],
                              score,
                              res$results_table[, c('match_score', 'ratio', 'uniqueness_score')],
                              score_frequency, score_solvents, score_organisms,
                              res$results_table[, c('n.peaks.matched', 'detailed_results_id')])
  }
  colnames(res$results_table)[4] = 'Final_Score'
  res$results_table = res$results_table[order(res$results_table$Final_Score, decreasing = TRUE), ]
  if (verbose) message("Done.")
  return(res)
}

Hyper_method = function(peaks, ppm.tol, n_peaks_total, alpha=10e-4){
  res = list()
  Name=c()
  spcmnm=c()
  ratio=c()
  uniqueness=c()
  match_score=c()
  p.fdr=c()
  n.peaks.matched=c()
  
  spcmns=c()
  j=0
  for (reference_id in names(spectra_list)){
    reference=spectra_list[[reference_id]]
    
    res.ref=match_Hyper(peaks, unique(reference$chemical.shift), n_peaks_total, ppm.tol)
    
    if(res.ref$matched_peaks!=0){
      spcmnm_ids=get_spcmnm_from_spcmns(reference_id)
      for (i in 1:length(spcmnm_ids)){
        spcmns=c(spcmns, reference_id)
        j=j+1
        id=paste(reference_id, j, sep="_")
        res$more_results[[id]]=res.ref[c("matched_peaks_ref", "matched_peaks_samp", "reference_peaks")]
        spcmnm=c(spcmnm, spcmnm_ids[i])
        Name=c(Name, codes[codes$SPCMNM==spcmnm_ids[i], "NAME"])
        match_score=c(match_score, res.ref$score)
        n.peaks.matched=c(n.peaks.matched, res.ref$matched_peaks)
      }
    }
  }
  
  p.fdr=p.adjust(match_score, "fdr")
  hypergeometric_score = rep(0, length(p.fdr))
  hypergeometric_score[p.fdr>alpha] = 0
  hypergeometric_score[p.fdr<alpha] = (1-(p.fdr[p.fdr<alpha]/alpha))
  
  res$results_table=data.frame(SPCMNM=spcmnm, Name, SPCMNS=spcmns, match_score=hypergeometric_score, n.peaks.matched, detailed_results_id=names(res$more_results), stringsAsFactors=FALSE)
  maintain_non_zeros = res$results_table$match_score!=0
  maintain_non_zeros_detailed = res$results_table$detailed_results_id[maintain_non_zeros]
  res$results_table = res$results_table[maintain_non_zeros,]
  res$more_results = res$more_results[maintain_non_zeros_detailed]
  
  return(res)
}

Hyper_uniq_method = function(peaks, ppm.tol, uniq_scores, n_peaks_total, alpha=10e-4){
  res = list()
  Name=c()
  spcmnm=c()
  ratio=c()
  uniqueness=c()
  match_score=c()
  p.fdr=c()
  n.peaks.matched=c()
  
  spcmns=c()
  j=0
  for (reference_id in names(spectra_list)){
    reference=spectra_list[[reference_id]]
    
    res.ref=match_Hyper(peaks, unique(reference$chemical.shift), n_peaks_total, ppm.tol)
    
    if(res.ref$matched_peaks!=0){
      spcmnm_ids=get_spcmnm_from_spcmns(reference_id)
      for (i in 1:length(spcmnm_ids)){
        spcmns=c(spcmns, reference_id)
        j=j+1
        id=paste(reference_id, j, sep="_")
        res$more_results[[id]]=res.ref[c("matched_peaks_ref", "matched_peaks_samp", "reference_peaks")]
        spcmnm=c(spcmnm, spcmnm_ids[i])
        Name=c(Name, codes[codes$SPCMNM==spcmnm_ids[i], "NAME"])
        match_score=c(match_score, res.ref$score)
        n.peaks.matched=c(n.peaks.matched, res.ref$matched_peaks)
      }
    }
  }
  
  p.fdr=p.adjust(match_score, "fdr")
  hypergeometric_score = rep(0, length(p.fdr))
  hypergeometric_score[p.fdr>alpha] = 0
  hypergeometric_score[p.fdr<alpha] = (1-(p.fdr[p.fdr<alpha]/alpha))
  
  spcmns_uniq_scores = c(NA, length(hypergeometric_score))
  Hyper_uniq_score = rep(0, length(hypergeometric_score))
  for(i in 1:length(hypergeometric_score)){
    if(hypergeometric_score[i]!=0){
      Hyper_uniq_score[i] = (hypergeometric_score[i] + uniq_scores[spcmns[i]]) / 2
      spcmns_uniq_scores[i] = uniq_scores[spcmns[i]]
    }
  }
  
  res$results_table=data.frame(SPCMNM=spcmnm, Name, SPCMNS=spcmns, match_score=Hyper_uniq_score, hypergeometric_score, uniqueness_score=spcmns_uniq_scores, n.peaks.matched, detailed_results_id=names(res$more_results), stringsAsFactors=FALSE)
  maintain_non_zeros = res$results_table$match_score!=0
  maintain_non_zeros_detailed = res$results_table$detailed_results_id[maintain_non_zeros]
  res$results_table = res$results_table[maintain_non_zeros,]
  res$more_results = res$more_results[maintain_non_zeros_detailed]
  
  return(res)
}

Match_uniq_method = function(peaks, ppm.tol, uniq_scores){
  res = list()
  Name=c()
  spcmnm=c()
  ratio=c()
  uniqueness=c()
  match_score=c()
  p.fdr=c()
  n.peaks.matched=c()
  
  spcmns=c()
  j=0
  for (reference_id in names(spectra_list)){
    reference=spectra_list[[reference_id]]
    
    res.ref=match_Match_uniq(peaks, unique(reference$chemical.shift), uniq_scores[[reference_id]], ppm.tol)
    
    if(res.ref$matched_peaks!=0){
      spcmnm_ids=get_spcmnm_from_spcmns(reference_id)
      for (i in 1:length(spcmnm_ids)){
        spcmns=c(spcmns, reference_id)
        j=j+1
        id=paste(reference_id, j, sep="_")
        res$more_results[[id]]=res.ref[c("matched_peaks_ref", "matched_peaks_samp", "reference_peaks")]
        spcmnm=c(spcmnm, spcmnm_ids[i])
        Name=c(Name, codes[codes$SPCMNM==spcmnm_ids[i], "NAME"])
        ratio=c(ratio, res.ref$ratio)
        uniqueness=c(uniqueness, res.ref$uniqueness)
        match_score=c(match_score, res.ref$score)
        n.peaks.matched=c(n.peaks.matched, res.ref$matched_peaks)
      }
    }
  }
  
  res$results_table=data.frame(SPCMNM=spcmnm, Name, SPCMNS=spcmns, match_score, ratio, uniqueness_score=uniqueness, n.peaks.matched, detailed_results_id=names(res$more_results), stringsAsFactors=FALSE)
  
  return(res)
}

match_Hyper = function(sample_peaks, refMetab_peaks, n_peaks_total, PPMTOL=0.03){
  l.samp=length(sample_peaks)
  l.ref=length(refMetab_peaks)
  
  res=list()
  
  matched_peaks_samp=c()
  matched_peaks_ref=c()
  
  sample_peaks=sort(sample_peaks)
  refMetab_peaks=sort(refMetab_peaks)
  
  matched_peaks=0
  n_ref_peaks=l.ref
  n_notRef_peaks=n_peaks_total-n_ref_peaks
  n_samp_peaks=l.samp
  
  if (l.samp<l.ref){
    last_pR=1
    for(pS in 1:length(sample_peaks)){
      if (last_pR>l.ref) break
      for(pR in last_pR:l.ref){
        inf_limit=refMetab_peaks[pR]-PPMTOL
        if(round(sample_peaks[pS],3) < round(inf_limit,3)) break
        sup_limit=refMetab_peaks[pR]+PPMTOL
        if(round(sample_peaks[pS],3) <= round(sup_limit,3)){
          matched_peaks = matched_peaks + 1
          last_pR=pR+1
          
          matched_peaks_samp=c(matched_peaks_samp, round(sample_peaks[pS],3))
          matched_peaks_ref=c(matched_peaks_ref, round(refMetab_peaks[pR],3))
          break
        }
      }
    }
  }
  else{
    last_pS=1
    for(pR in 1:length(refMetab_peaks)){
      if (last_pS>l.samp) break
      for(pS in last_pS:l.samp){
        inf_limit=sample_peaks[pS]-PPMTOL
        if(round(refMetab_peaks[pR],3) < round(inf_limit,3)) break
        sup_limit=sample_peaks[pS]+PPMTOL
        if(round(refMetab_peaks[pR],3) <= round(sup_limit,3)){
          matched_peaks = matched_peaks + 1
          last_pS=pS+1
          
          matched_peaks_samp=c(matched_peaks_samp, round(sample_peaks[pS],3))
          matched_peaks_ref=c(matched_peaks_ref, round(refMetab_peaks[pR],3))
          break
        }
      }
    }
  }
  
  ht=dhyper(matched_peaks, n_ref_peaks, n_notRef_peaks, n_samp_peaks)
  
  res$matched_peaks_ref=matched_peaks_ref
  res$matched_peaks_samp=matched_peaks_samp
  res$reference_peaks=refMetab_peaks
  res$score=ht
  res$matched_peaks=matched_peaks
  
  return(res)
}

match_Match_uniq = function(sample_peaks, refMetab_peaks, uniq_score, PPMTOL=0.03){
  l.samp=length(sample_peaks)
  l.ref=length(refMetab_peaks)
  
  res=list()
  
  matched_peaks_samp=c()
  matched_peaks_ref=c()
  
  sample_peaks=sort(sample_peaks)
  refMetab_peaks=sort(refMetab_peaks)
  
  matched_peaks=0
  n_ref_peaks=l.ref
  n_samp_peaks=l.samp
  if (l.samp<l.ref){
    last_pR=1
    for(pS in 1:length(sample_peaks)){
      if (last_pR>l.ref) break
      for(pR in last_pR:l.ref){
        inf_limit=refMetab_peaks[pR]-PPMTOL
        if(round(sample_peaks[pS],3) < round(inf_limit,3)) break
        sup_limit=refMetab_peaks[pR]+PPMTOL
        if(round(sample_peaks[pS],3) <= round(sup_limit,3)){
          matched_peaks = matched_peaks + 1
          last_pR=pR+1
          
          matched_peaks_samp=c(matched_peaks_samp, round(sample_peaks[pS],3))
          matched_peaks_ref=c(matched_peaks_ref, round(refMetab_peaks[pR],3))
          break
        }
      }
    }
  }
  else{
    last_pS=1
    for(pR in 1:length(refMetab_peaks)){
      if (last_pS>l.samp) break
      for(pS in last_pS:l.samp){
        inf_limit=sample_peaks[pS]-PPMTOL
        if(round(refMetab_peaks[pR],3) < round(inf_limit,3)) break
        sup_limit=sample_peaks[pS]+PPMTOL
        if(round(refMetab_peaks[pR],3) <= round(sup_limit,3)){
          matched_peaks = matched_peaks + 1
          last_pS=pS+1
          
          matched_peaks_samp=c(matched_peaks_samp, round(sample_peaks[pS],3))
          matched_peaks_ref=c(matched_peaks_ref, round(refMetab_peaks[pR],3))
          break
        }
      }
    }
  }
  
  if(matched_peaks==0) score = 0
  else score = ((matched_peaks/n_ref_peaks) + uniq_score) / 2
  
  res$matched_peaks_ref=matched_peaks_ref
  res$matched_peaks_samp=matched_peaks_samp
  res$reference_peaks=refMetab_peaks
  res$score=score
  res$ratio=matched_peaks/n_ref_peaks
  if(matched_peaks!=0) res$uniqueness=uniq_score
  else res$uniqueness=NA
  res$matched_peaks=matched_peaks
  
  return(res)
}

uniqueness_scores = function(ppm_tolerance){
  
  all_peaks = c()
  for(spec in names(spectra_list)) all_peaks = c(all_peaks, unique(round(spectra_list[[spec]]$chemical.shift, digits=3)))
  peaks_table = table(all_peaks)
  
  uniq_values_peaks = c()
  for(peak in unique(all_peaks)){
    overlap = sum(peaks_table[as.character(all_peaks[all_peaks>=peak-ppm_tolerance & all_peaks<=peak+ppm_tolerance])])
    uniq_values_peaks = c(uniq_values_peaks, 1/overlap)
  }
  names(uniq_values_peaks) = as.character(unique(all_peaks))
  
  res = c()
  for(spec in names(spectra_list)){
    peaks = as.character(unique(round(spectra_list[[spec]]$chemical.shift, digits=3)))
    uniq_score_spec = mean(uniq_values_peaks[peaks])
    res = c(res, uniq_score_spec)
  }
  names(res) = names(spectra_list)
  
  return(res)
}

###DATA CONVERSION
##SPCMNS-->SPCMNM
get_spcmnm_from_spcmns=function(spcmns){
  metabolites=unlist(strsplit(spectra_options[spectra_options$SPCMNS==spcmns, "SPCMNM"], "; "))
  return(metabolites)
}

##SPCMNM-->KEGG
convert_spcmnm_to_kegg=function(spcmnm){
  keggs=unlist(strsplit(codes$KEGG[codes$SPCMNM==spcmnm], "; ")[[1]])
  return(keggs)
}

##CHEBI-->SPCMNM
convert_chebi_to_spcmnm=function(chebi){
  spcmnms=c()
  for(i in 1:length(codes$CHEBI)){
    chebi_ref=codes$CHEBI[i]
    if (chebi %in% unlist(strsplit(chebi_ref, "; ")[[1]])) spcmnms=c(spcmnms, codes$SPCMNM[i])
  }
  return(spcmnms)
}

###SCORES
##Frequency
score_freq=function(spcmns, scores=list('400'=0.5, '500'=1, '600'=0.5, '700'=0.5)){
  score=scores[[as.character(spectra_options[spectra_options$SPCMNS==spcmns, "FREQUENCY"])]]
  return(score)
}

##Solvent
score_solvent=function(spcmns, scores=list(CD3OD=1, D2O=0.8, Water=0.8, CDCl3=0.6, 'Acetone-d6'=0.6, 'Acetone'=0.6,
                                           "DMSO-d6"=0.4, "100%_DMSO"=0.4, "5%_DMSO"=0.4, C=0.2, C6D6=0.2,
                                           CD3CN=0.2, C2D2Cl4=0.2, CD2Cl2=0.2, CDC3OD=0.2, Ethanol=1)){
  score=scores[[spectra_options[spectra_options$SPCMNS==spcmns, "SOLVENT"]]]
  return(score)
}

##organism
score_organism=function(spcmnm, cpds_groups, scores=list('hsa'=1, 'Mammals'=0.8, 'Vertebrates'=0.7, 'Animals'=0.6, 'Eukaryotes'=0.3, 'other'=0.1, 'not_in_kegg'=0.6)){
  #To specify organism: Kegg code
  #Availale groups: "Eukaryotes", "Animals", "Vertebrates", "Mammals", "Birds", "Reptiles", "Amphibians", "Fishes", "Arthropods", "Insects",
  #"Nematodes", "Mollusks", "Cnidarians", "Plants", "Eudicots", "Monocots", "Green algae", "Red algae", "Fungi", "Protist", "Prokaryotes",
  #"Bacteria", "Archaea"
  
  cpds=convert_spcmnm_to_kegg(spcmnm)
  
  if(is.null(cpds)) return(scores[['not_in_kegg']])
  else{
    possible_scores=c(scores[['other']])
    for(cpd in cpds){
      for(group in names(scores)){
        if(!group%in%c('other', 'not_in_kegg')) if(cpd%in%cpds_groups[[group]]) possible_scores=c(possible_scores, scores[[group]])
      }
    }
    score=max(possible_scores)
  }
  return(score)
}