R/helperFunctions.R

In protViz/bibliospec: Reading and Analyzing Mass Spectrometric Search Results

#' Get Calibration models for spectrum
#' 
#' @param spectraMeta - get it by calling spectraMeta <- bibliodb$getSpectraMeta()
#' @param iRTpeptides - 
#' 
#' @export
#' @import MASS
#' @import ggplot2
#' 
#' @examples 
#' data(iRTpeptides)
#' 
getCalibrationModelsRT <- function(spectraMeta, iRTpeptides){
  iRTSpectraMeta <- merge(spectraMeta, iRTpeptides, by.x="peptideModSeq", by.y="peptide")
  plot = ggplot2::ggplot(iRTSpectraMeta, aes( x = retentionTime, y = iRT)) + geom_point() + facet_wrap( ~ fileName) + stat_smooth(method = "lm", col = "red")
  myrlm <- function(x){ 
    x <- unique(subset(x, select = c("iRT", "retentionTime")))
    model <- MASS::rlm(iRT ~ retentionTime, x) 
  }
  iRTmodelList <- plyr::dlply(iRTSpectraMeta, .(fileName), myrlm)
  return(list(iRTSpectraMeta=iRTSpectraMeta, iRTmodelList=iRTmodelList, plot = plot))
}
#' Apply Calibration models
#'
#' @param spectraMeta - get it by calling spectraMeta <- bibliodb$getSpectraMeta()
#' @param modelList - get by calling getCalibrationModelsRT
#' @import MASS
#' @export
#' 
applyRTCalibrationModels <- function(spectraMeta, modelList){
  psmall <- dlply(spectraMeta, ~fileName)
  stopifnot(names(modelList)==names(psmall))
  
  rtCalibrate <- function(mylm, data){
    preddata <- predict(mylm,newdata = data.frame("retentionTime" = data$retentionTime))
    data$iRT <- preddata
    data
  }
  iRTCalibrated <- mapply(rtCalibrate, modelList, psmall,SIMPLIFY = FALSE)
  spectraMetaWithIRT <- plyr::rbind.fill(iRTCalibrated)
  spectraMetaWithIRT
}
#'  For each precursor select top PSM
#'  
#' @param x extract best scoring PSM
#' @param score_larger_better indicate if large score is better
#'  
#' @export
#'  
getTOPhit <- function(x, score_larger_better=FALSE){
  
  id <- if(!score_larger_better){
    x$RefSpectraId[which.min(x$score)]
  }else{
    x$RefSpectraId[which.max(x$score)]
  }
  
  
  c(meanIRT = mean(x$iRT),
    medianIRT = stats::median(x$iRT),
    sdIRT = stats::sd(x$iRT),
    RefSpectraId = id,
    topOfNrPrecursors = nrow(x))
}