R/annotation.R

In prafols/rMSIproc: MSI Processing in C++

#########################################################################
#     
#     Copyright (C) 2018 Lluc Semente Fernandez
# 
#     This program is free software: you can redistribute it and/or modify
#     it under the terms of the GNU General Public License as published by
#     the Free Software Foundation, either version 3 of the License, or
#     (at your option) any later version.
# 
#     This program is distributed in the hope that it will be useful,
#     but WITHOUT ANY WARRANTY; without even the implied warranty of
#     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#     GNU General Public License for more details.
# 
#     You should have received a copy of the GNU General Public License
#     along with this program.  If not, see <http://www.gnu.org/licenses/>.
############################################################################

#' rMSIannotation
#' 
#' Searches for isotopic ions in the peak matrix and evaluates them using morphology, intensity and mass error criteria.
#' This algorithm only searches carbon-based isotopic ions.
#' Using the infromation generated by the isotopes test, the algorithm then proceeds to search adduct pairs between the monoisotopic ions.
#' 
#' @param PeakMtx List. An rMSIprocPeakMatrix. Must contain at least the following categories: \itemize{
#' \item PeakMtx$intensity. A matrix containg the intensities of the peaks for each pixel (rows = pixels, cols = ions).  
#' \item PeakMtx$mass. A vector containg the masses of each peak. Must be in the same order with the columns of the intensity marix.   
#' \item PeakMtx$numPixels. Number of pixels (rows in your matrix).   
#' }
#' @param iso.number Integer. Number of isotopes for each monoisotopic ion to be found.
#' @param iso.tolerance Integer. Mass tolerance for the distance between isotope candidates in scans or ppms. (1.003355 Da +/- iso.tolerance)
#' @param iso.charge Integer. Charge of the patterns to be found.
#' @param iso.scoreThreshold Numeric between 0 and 1. ILS value to consider a pair of ions isotopes.
#' @param iso.toleranceUnits String. The units of iso.tolerance. Must be 'ppm' or 'scan'. If ToleranceUnits is 'scan' then ImageVector must be the mass channels vector of the rMSI image (rMSIObj$mass).
#' @param iso.imageVector Numeric Vector. The mass channels vector of the imaging dataset containing all the scans.
#' 
#' @param addu.adductDataFrame Data frame containing the adducts to be searched. By defalut (+H, +Na and +K). It must have two columns.\itemize{
#'  \item $name. Column containing the names as strings of the elements or molecules that form adducts (i.e. "Na", "K", "H") 
#'  \item $mass. Masses of the adduct forming elements.
#'  }
#' @param addu.tolerance Integer. Mass error tolerance in ppm.
#' @param csv.results Boolean. If TRUE, the annotations will be writen in CSV files under the name A.csv, B.csv and C.csv. By default FALSE.
#' @param csv.path Path. Path to the folder where the CSV files shoul be stored. By default the working directory.
#' @return List. 
#' \itemize{
#'   \item $A: Data frame of the pairs of M+0 ions cataloged as adducts.
#'   \item $B: Data frame of the pairs of M+0 ions and non-isotopic ions cataloged as adducts.
#'   \item $C: Data frame of the all the M+0 ions.
#'   \item $isotopicTestData: All the test results computed to determine the monoisotopic ions.
#'   \item $monoisotopicPeaks: Peak matrix column indeces of monoisotopic peaks. 
#'   \item $isotopicPeaks: Peak matrix column indeces of isotopic peaks.
#' }
#' @export
#' 
rMSIannotation <- function(PeakMtx, 
                           iso.number = 2, 
                           iso.tolerance = 30, 
                           iso.charge = 1,
                           iso.scoreThreshold = 0.75,
                           iso.toleranceUnits = "ppm",
                           iso.imageVector = NULL,
                           addu.tolerance = 50,
                           addu.adducts = data.frame(mass = c(38.963706,22.98976,1.007825),
                                                    name = c("K","Na","H")),
                           csv.results = FALSE,
                           csv.path = getwd())
{
  isotopeObj <- list(NULL)
  adductObj <- list(NULL)
    
  isotopeObj <- isotopeAnnotation(PeakMtx, isoNumber = iso.number, tolerance = iso.tolerance, 
                                  charge = iso.charge, scoreThreshold = iso.scoreThreshold,
                                  toleranceUnits = iso.toleranceUnits, imageVector = iso.imageVector)
  
  if(length(isotopeObj$monoisotopicPeaks) == 0)
  {
    cat(paste("Procedure cancelled \n"))
    return() 
  }
  
  if(length(isotopeObj$monoisotopicPeaks) >= 2)
  {
    adductObj <- adductAnnotation(isotopeObj = isotopeObj, PeakMtx = PeakMtx,
                                  adductDataFrame = addu.adducts, tolerance = addu.tolerance)
  }

  results <- list()
  results <- annotationOutpuFormat(isotopeObj, adductObj, PeakMtx$mass)
  results$parameters <- list(isotope_number = iso.number,
                            isotope_charge = iso.charge,
                            isotope_mass_tolerance = iso.tolerance,
                            isotope_charge = iso.charge,
                            isotope_score_threshold = iso.scoreThreshold,
                            isotope_tolerance_units = iso.toleranceUnits,
                            adduct_mass_tolerance = addu.tolerance,
                            adduct_ions = addu.adducts)
  
  if(csv.results)
  {
    if(!is.null(adductObj[[1]]))
    {
      if(!is.null(results$A)) {write.csv(x = results$A, file = paste(csv.path,"/A.csv",sep=""),row.names = FALSE)}
      if(!is.null(results$B)) {write.csv(x = results$B, file = paste(csv.path,"/B.csv",sep=""),row.names = FALSE)}
    }
    write.csv(x = results$C, file = paste(csv.path,"/C.csv",sep=""),row.names = FALSE)
  }
  return(results)
}


annotationOutpuFormat <- function(isotopeObj, adductObj, massAxis)
{
  if(!is.null(adductObj[[1]]))
  {
    if(!is.null(adductObj$A))
    {
      A <- adductObj$A
    }
    
    if(!is.null(adductObj$B))
    {
      B <- adductObj$B
    }
  }
  
  C <- data.frame(MonoisotopicMass = rep(0, times = length(isotopeObj$monoisotopicPeaks)),
                  ILS = rep(0, times = length(isotopeObj$monoisotopicPeaks)),
                  IsotopicIntensityRatio = rep(0, times = length(isotopeObj$monoisotopicPeaks)),
                  EstimatedCarbonAtoms = rep(0, times = length(isotopeObj$monoisotopicPeaks)),
                  MonoisotopicIndex = rep(0, times = length(isotopeObj$monoisotopicPeaks)))
  
  ord <- c()
  for (i in 1:length(isotopeObj$monoisotopicPeaks)) 
  {
    ord <- c(ord,which.min(abs(massAxis[sort(isotopeObj$monoisotopicPeaks)[i]]-as.numeric(names(isotopeObj$M1)))))
  }
  
  for(i in 1:length(ord))
  {
    C$MonoisotopicMass[i]       <- massAxis[isotopeObj$M1[[ord[i]]][3]]
    C$MonoisotopicIndex[i]      <- isotopeObj$M1[[ord[i]]][3,which.max(isotopeObj$M1[[ord[i]]][5,])]
    C$ILS[i]                    <- isotopeObj$M1[[ord[i]]][5,which.max(isotopeObj$M1[[ord[i]]][5,])]
    C$IsotopicIntensityRatio[i] <- isotopeObj$M1[[ord[i]]][9,which.max(isotopeObj$M1[[ord[i]]][5,])]
    C$EstimatedCarbonAtoms[i]   <- isotopeObj$M1[[ord[i]]][10,which.max(isotopeObj$M1[[ord[i]]][5,])]
  }
  
  C$MonoisotopicMass       <- trunc(C$MonoisotopicMass) + signif(C$MonoisotopicMass-trunc(C$MonoisotopicMass), digits = 4) 
  C$MonoisotopicIndex      <- trunc(C$MonoisotopicIndex) + signif(C$MonoisotopicIndex-trunc(C$MonoisotopicIndex), digits = 4)
  C$ILS                    <- trunc(C$ILS) + signif(C$ILS-trunc(C$ILS), digits = 4)
  C$IsotopicIntensityRatio <- trunc(C$IsotopicIntensityRatio) + signif(C$IsotopicIntensityRatio-trunc(C$IsotopicIntensityRatio), digits = 4)
  C$EstimatedCarbonAtoms   <- round(C$EstimatedCarbonAtoms)
  
  IsotopicTestData <- isotopeObj[-((length(isotopeObj)-1):length(isotopeObj))]
  
  if(!is.null(adductObj[[1]]))
  {
    if(is.null(adductObj$A))
    {
      results <- list(B = B,C = C, isotopicTestData = IsotopicTestData, monoisotopicPeaks = isotopeObj$monoisotopicPeaks, isotopicPeaks = isotopeObj$isotopicPeaks)
    }
    
    if(is.null(adductObj$B))
    {
      results <- list(A = A,C = C, isotopicTestData = IsotopicTestData, monoisotopicPeaks = isotopeObj$monoisotopicPeaks, isotopicPeaks = isotopeObj$isotopicPeaks)
    }
    
    if(!is.null(adductObj$B) & !is.null(adductObj$A))
    {
      results <- list(A = A,B = B,C = C, isotopicTestData = IsotopicTestData, monoisotopicPeaks = isotopeObj$monoisotopicPeaks, isotopicPeaks = isotopeObj$isotopicPeaks)
    }
  }
  else
  {
    results <- list(C = C, isotopicTestData = IsotopicTestData, monoisotopicPeaks = isotopeObj$monoisotopicPeaks, isotopicPeaks = isotopeObj$isotopicPeaks)
  }
  
return(results)
}