#########################################################################
#
# 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/>.
############################################################################
#' isotopeAnnotation
#'
#' Finds and evaluate isotope candidates for each ion mass.
#'
#'
#' @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 = peaks).
#' \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 isoNumber Integer. Number of isotopes to be found.
#' @param tolerance Integer. Mass tolerance for the candidates in scans or ppms.
#' @param charge Integer. Charge of the patterns to be found.
#' @param scoreThreshold Numeric. Score value to consider a ion mass a good isotope candidate. Only the ions that have this number or greater will undergo the following isotope searching stages.
#' @param toleranceUnits String. Must be 'ppm' or 'scan'. If ToleranceUnits is 'scan' then ImageVector must be the mass channels vector of the rMSI image (rMSIObj$mass).
#' @param imageVector Numeric Vector. The mass channels vector of the imaging dataset containing all the scans.
#'
#' @return A list containing the results of the test and other kind of information.
#'
isotopeAnnotation <- function(PeakMtx, isoNumber = 2, tolerance = 30, charge = 1, scoreThreshold = 0.8, toleranceUnits = "ppm", imageVector = NULL)
{
result <- list()
##### Type check #####
{
if(!(is.matrix(PeakMtx$intensity) & (length(PeakMtx$mass)>0)))
{
stop("Error: PeakMtx has not the rMSIprocPeakMatrix format. Look at help for more information.")
}
if(isoNumber < 1)
{
stop("Error: IsoNumber must be 1 or a greater integer")
}
if ((scoreThreshold > 1) | (scoreThreshold < 0))
{
stop("Error: ScoreThreshold must be between 0 and 1")
}
if(toleranceUnits == "ppm")
{
InScans <- FALSE
imageVector <- c(1,2)
} else if(toleranceUnits == "scan")
{
InScans <- TRUE
} else
{
stop("Error: ToleranceUnits must be 'ppm' or 'scan'")
}
if((toleranceUnits == "scan") & (length(PeakMtx$mass) > length(imageVector)))
{
stop("Error: If ToleranceUnits is 'scan' then ImageVector must be the mass channels vector of the rMSI image. (rMSIObj$mass)")
}
if((toleranceUnits == "ppm") & (tolerance>500))
{
writeLines("Maximum allowed tolerance is 500. Replacing the value")
tolerance <- 500
}
if(!(isoNumber%%1 == 0))
{
stop("Error: IsoNumber must be an integer")
}
if((toleranceUnits == "scan") & !(tolerance%%1 == 0))
{
stop("Error: If ToleranceUnits is 'scan' then Tolerance must be integer")
}
if((toleranceUnits == "scan") & (tolerance>length(imageVector)))
{
stop("Error: Not enough scans in your mass axis, reduce the tolerance")
}
}
##### Calling the C++ method #####
result <- C_isotopeAnnotator(length(PeakMtx$mass), #number of mass peaks
length(imageVector), #number of massChannels
sum(PeakMtx$numPixels), #number of pixels
isoNumber, #number of isotopes
PeakMtx$intensity, #peak matrix
PeakMtx$mass, #matrix mass axis
imageVector, #image mass axis
tolerance, #tolerance in ppm or scans
scoreThreshold, #score threshold
InScans, #tolerance units
charge)
##### Output format #####
result <- DeisotopingOutputFormat(result, scoreThreshold)
return(result)
}
#' DeisotopingOutputFormat
#'
#' Gives format to the output.
DeisotopingOutputFormat <- function(r, ScoreThreshold)
{
##### Names the results list #####
NullVec <- matrix(nrow = length(r)-1, ncol = length(r[[1]]))
Nullcnt <- 0
for(i in 2:length(r))
{
names(r[[i]]) <- signif(r[[1]],digits = 8)
Nullcnt <- 0
for(j in 1:length(r[[i]]))
{
if(!is.null(r[[i]][[j]]))
{
row.names(r[[i]][[j]]) <- c("M+N mass","Abs ppm error","M+0 index", "M+N index","Final score","Morphology score","Intensity score", "Mass error score", "Isotopic ratio", "Number of C atoms")
colnames(r[[i]][[j]]) <- (1:(dim(r[[i]][[j]])[2]))
}else
{
Nullcnt <- Nullcnt + 1
NullVec[i-1,Nullcnt] <- j
}
}
r[[i]] <- r[[i]][-(NullVec[i-1,1:Nullcnt])]
}
##### Removes the mass vector and the empty isotope lists #####
r <- r[-1]
flagRemList <- FALSE
RemList <- c()
for(i in 1:length(r))
{
if(length(r[[i]])==0)
{
flagRemList <- TRUE
RemList <- c(RemList,i)
}
}
if(flagRemList)
{
r <- r[-RemList]
}
if(length(r) == 0)
{
cat(paste("Number of monoisotopic ions found: 0\n"))
return()
}
##### Creats the isotopic & monoisotopic ions index vectors #####
MonoVec <- c()
CompVec <- c()
for(i in (1:(length(r))))
{
for(j in (1:length(r[[i]])))
{
if((!is.null(r[[i]][[j]])) & (as.numeric(r[[i]][[j]][5,which.max(r[[i]][[j]][5,])]) >= ScoreThreshold))
{
if(i == 1)
{
MonoVec <- c(MonoVec, r[[i]][[j]][3,which.max(r[[i]][[j]][5,])])
}
}
for(k in (1:ncol(r[[i]][[j]])))
{
if(!is.nan(r[[i]][[j]][5, k]))
{
if(as.numeric(r[[i]][[j]][5, k]) > ScoreThreshold )
{
CompVec <- c(CompVec, r[[i]][[j]][4, k])
}
}
}
}
}
names(r) <- paste("M",1:length(r),sep = "")
r$isotopicPeaks <- (unique(CompVec))
r$monoisotopicPeaks <- (unique(MonoVec))
r$monoisotopicPeaks <- setdiff(r$monoisotopicPeaks, r$isotopicPeaks)
MonoVec <- r$monoisotopicPeaks
cat(paste("Number of monoisotopic ions found:", length(MonoVec),"\n"))
return(r)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.