R/library_visualizer.R

In daniellyz/MergeION: Merge Scans from Multiple LC-MS/MS Raw Data Files

#' Visualize selected mass spectra in the spectral library
#'
#' The function plots mass spectra of selected compound in a library
#'
#' @param library A list generated by the function library_generator() or library_similarity()$SELECTED
#' @param id Compound id of the spectra to be plotted (both most recent MS1 and MS2 spectra will be plotted!)
#' @param png.out Boolean. True if plotted spectra are exported as png images!
#'
#' @examples
#'
#' data(DRUG_THERMO_LIBRARY)
#'
#' # Search library using query command lines:
#' query = library_manager(library2,query=c("IONMODE=Positive","RT=1.2"), logical="AND", rt_search=6)
#'
#' # Visualize the first compound found:
#' library_visualizer(query$SELECTED, id = query$ID_SELECTED[1])
#'
#' @export
#'
#' @importFrom MSnbase fData readMgfData
#' @importFrom graphics plot title legend abline text
#' @importFrom grDevices png dev.off

library_visualizer<-function(library, id = library$metadata$ID[1], png.out=F, show.legend = T){

  options(stringsAsFactors = FALSE)
  options(warn=-1)
  max_display = 5 # Display text of 5 most abundant mass peaks and higher than 5%

  #################
  ### Check inputs:
  #################

  if (missing(library)){
    stop("Please provide the output of library_generator() or a .mgf file as input library!")
  }

  if (is.character(library)){
    if (file_ext(library)!="mgf"){
      stop("The file extension of your input library must be mgf!")
    }}

  if (is.list(library)){
    if (length(library)==2 & "complete" %in% names(library)){
      library = library$complete
    }
    if (length(library)!=2 || (!is.list(library$sp)) || !is.data.frame(library$metadata)){
      stop("Please make sure your input library is a valid output of library_generator()!")
    }}

  ##################################
  ### Reading from spectral library:
  ##################################

  if (is.character(library)){ # If input is a mgf file name
    library=readMGF2(library)}

  if (nrow(library$metadata)==0){stop("The library for visualization is empty!")}

  library1 = library_manager(library, query = paste0("ID =", id))$SELECTED

  metadata = library1$metadata
  spectrum_list = library1$sp

  if (nrow(metadata)==0){stop("The library does not contain the selected ID!")}

  ##############
  ### Visualize:
  ##############

  if (png.out){png(paste0("QUERY_ID_",id,".png"), width = 700, height = 480)}

  metadata1 = metadata[as.numeric(metadata$MSLEVEL)==1,]
  metadata2 = metadata[as.numeric(metadata$MSLEVEL)==2,]

  if ((nrow(metadata1)>0) & (nrow(metadata2)>0)){
    par(mfrow = c(2,1), mar=c(2,2,2,2))  # Joint plot MS1 MS2
  } else {par(mfrow = c(1,1), mar=c(2,2,2,2))}

  if (nrow(metadata1)>0){

    scan1 = max(as.numeric(metadata1$SCANS)) # Scan selected
    ind1 = which(as.numeric(metadata$SCANS)==scan1)

    metadata1 = metadata[ind1,]
    spectrum1 = spectrum_list[[ind1]]
    prec_mz = round(as.numeric(metadata[ind1,"PEPMASS"]),3)
    adduct1 = metadata$ADDUCT[ind1]
    cpd1 = "Unknown"
    if ("COMPOUND" %in% colnames(metadata)){cpd1 = metadata$COMPOUND[ind1]}

    xrange = c(prec_mz-2,prec_mz+8)
    ranges = which((spectrum1[,1]>=prec_mz-2) & (spectrum1[,1]<=prec_mz+8))
    yrange = c(10, max(spectrum1[ranges,2])*1.2)

    plot(spectrum1[,1],spectrum1[,2], type = "h", xlim = xrange, ylim = yrange,
         xlab = "m/z", ylab = "Intensity", font.lab=2)

    kkk = min(nrow(spectrum1), max_display)
    max_pics = order(spectrum1[,2], decreasing = T)[1:kkk]
    max_pics = intersect(max_pics,which(spectrum1[,2]>=max(spectrum1[,2])*0.05))
    text_pics = spectrum1[max_pics,1]
    int_pics = spectrum1[max_pics,2]*1.1
    text(text_pics, int_pics, as.character(round(text_pics,3)))

    title(main = paste0("ID: ",id, " ; SCAN: ", scan1, " (", cpd1, ")"),
          font = 3, cex.main = 1.5)
    if (show.legend){
    legend("topleft", bty = "n", cex = 1.2, text.font = 1.5,
           legend = paste0(
             "Precursor: ", prec_mz,  "\n",
             "Adduct: ", adduct1, "\n",
             "MS Level: 1"))}
  }

  if (nrow(metadata2)>0){

    scan2 = max(as.numeric(metadata2$SCANS)) # Scan selected
    ind2 = which(as.numeric(metadata$SCANS)==scan2)

    metadata2 = metadata[ind2,]
    spectrum2 = spectrum_list[[ind2]]
    prec_mz = round(as.numeric(metadata[ind2,"PEPMASS"]),3)
    adduct2 = metadata$ADDUCT[ind2]
    cpd2 = "Unknown"
    if ("COMPOUND" %in% colnames(metadata)){cpd2 = metadata$COMPOUND[ind2]}

    xrange = c(50,prec_mz+8)
    ranges = which((spectrum2[,1]>=50) & (spectrum2[,1]<=prec_mz+8))
    yrange = c(10, max(spectrum2[ranges,2])*1.2)

    spectrum2 = spectrum2[ranges,]
    if (length(ranges)==1){spectrum = matrix(spectrum, ncol=2)}

    plot(spectrum2[,1],spectrum2[,2], type = "h", xlim = xrange, ylim = yrange,
         xlab = "m/z", ylab = "Intensity", font.lab=2)

    kkk = min(nrow(spectrum2), max_display)
    max_pics = order(spectrum2[,2], decreasing = T)[1:kkk]
    max_pics = intersect(max_pics,which(spectrum2[,2]>=max(spectrum2[,2])*0.05))
    text_pics = spectrum2[max_pics,1]
    int_pics = spectrum2[max_pics,2]*1.1
    text(text_pics, int_pics, as.character(round(text_pics,3)))

    title(main = paste0("ID: ", id, " ; SCAN: ", scan2, " (", cpd2, ")"),
          font = 3, cex.main = 1.5)

    if (show.legend){
      legend("topleft", bty = "n", cex = 1.2, text.font = 1.5,
            legend = paste0(
            "Precursor: ", prec_mz,  "\n",
            "Adduct: ", adduct2, "\n",
            "MS Level: 2"))}

    legend("topleft", bty = "n", cex = 1.2, text.font = 1.5,
            legend = paste0(
            "Precursor: ", prec_mz,  "\n",
            "Adduct: ", adduct2, "\n",
            "MS Level: 2"))
    abline(v = prec_mz, col = "red")
   }

    if (png.out){dev.off()}
}