R/library_visualizer.R

In daniellyz/MergeION2: Merging and Networking 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 input_library A list generated by the function library_generator() or the name of mgf/msp/RData spectral library file
#' @param id Compound id of the spectra to be plotted (both most recent MS1 and MS2 spectra will be plotted!)
#' @param type Character. "Complete" if raw specra are visualized. "Consensus" if consensus spectra are visualized
#' @param query_spectrum  Two-column data matrix. Optional for mirror plot. Two columns represent m/z and intensity of query tandem spectrum.
#' @param add.legend Boolean. If TRUE, figure legend is added.
#' 
#' @examples
#' data(DRUG_THERMO_LIBRARY)
#'
#' # Search library using query command lines:
#'  query = process_query(library2,query=c("IONMODE=Positive"),  rt_search=6, sepString = "@xxx@")
#' # Visualize the first compound found:
#'  library_visualizer(query$SELECTED, id = "1")
#' 
#' # Plot against query spectrum:
#' query_spectrum = cbind(c(71.084, 89.6, 455.178), c(90,40,39))
#' library_visualizer(library2, id = "1", query_spectrum,  type = "complete") 
#'
#' @export
#'
#' @importFrom MSnbase fData readMgfData
#' @importFrom graphics plot title legend abline text
#' @importFrom grDevices png dev.off
#' @importFrom OrgMassSpecR SpectrumSimilarity

library_visualizer<-function(input_library, id = input_library$metadata$ID[1], type = "complete", query_spectrum=NULL, add.legend = T){

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

  #####################################
  ### Reading from spectral library:###
  #####################################
  
  input_library = library_reader(input_library)
  type = type[1]
  
  if (!type %in% c("complete", "consensus")){
    stop("Please choose the visualized library type between complete and consensus!")
  }
  
  if (is.null(input_library$consensus)){
    type = "complete"
    message("No consensus library is available! The entire library is used for visualization!")
  }
  
  if (type == "complete"){input_library = input_library$complete}
  if (type == "consensus"){input_library = input_library$consensus}
  
  library1 = process_query(input_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 without query spectrum:####
  #########################################

  if (is.null(query_spectrum)){
    
    metadata1 = metadata[as.numeric(metadata$MSLEVEL)==1,,drop = FALSE]
    metadata2 = metadata[as.numeric(metadata$MSLEVEL)==2,,drop = FALSE]

    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))}

    # Plot MS1 spectrum:
    
    if (nrow(metadata1)>0){

      scan1 = max(as.numeric(metadata1$SCANS)) # Scan selected
      ind1 = which(as.numeric(metadata$SCANS)==scan1)
      if (is.na(scan1)){
        scan1 = max(metadata1$SCANS) # Scan selected
        ind1 = which(metadata$SCANS==scan1)
      }
      
      metadata1 = metadata[ind1,]
      spectrum1 = spectrum_list[[ind1]]
      prec_mz = round(as.numeric(metadata[ind1,"PEPMASS"]),4)
      adduct1 = metadata$ADDUCT[ind1]
      cpd1 = "Unknown"
      if ("COMPOUND" %in% colnames(metadata)){cpd1 = metadata$COMPOUND[ind1]}

      xrange = c(prec_mz-30,prec_mz+60)
      ranges = which((spectrum1[,1]>=prec_mz-30) & (spectrum1[,1]<=prec_mz+60))
      yrange = c(-max(spectrum1[ranges,2])*0.01, 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)
      abline(h=0)
      
      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,4)))

      title(main = paste0("ID: ",id, " ; SCAN: ", scan1, " (", cpd1, ")"),
          font = 3, cex.main = 1.5)
      
      if (add.legend){
      legend("topleft", bty = "n", cex = 1.2, text.font = 1.5,
           legend = paste0(
             "Precursor: ", prec_mz,  "\n",
             "Adduct: ", adduct1, "\n",
             "MS Level: 1"))}
    }
  
  # Plot MS2 spectrum:
    
  if (nrow(metadata2)>0){

    scan2 = max(as.numeric(metadata2$SCANS)) # Scan selected
    ind2 = which(as.numeric(metadata$SCANS)==scan2)
    
    if (is.na(scan2)){
      scan2 = max(metadata2$SCANS) # Scan selected
      ind2 = which(metadata$SCANS==scan2)
    }
    
    metadata2 = metadata[ind2,]
    spectrum2 = spectrum_list[[ind2]]
    prec_mz = round(as.numeric(metadata[ind2,"PEPMASS"]),4)
    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(max(spectrum2[ranges,2])*0.01, max(spectrum2[ranges,2])*1.2)

    spectrum2 = spectrum2[ranges,1:2,drop=FALSE]
    #spectrum2 = apply(spectrum2[,1:2,drop=FALSE], 2, as.numeric)

    plot(spectrum2[,1],spectrum2[,2], type = "h", xlim = xrange, ylim = yrange,
         xlab = "m/z", ylab = "Intensity", font.lab=2)
    abline(h=0)
    
    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,4)))

    if (add.legend){
    legend("topleft", bty = "n", cex = 1.2, text.font = 1.5,
          legend = paste0(
          "Precursor: ", prec_mz,  "\n",
          "Adduct: ", adduct2, "\n",
          "MS Level: 2"))}
   }
  }
  
  #########################################
  ### Visualize against query spectrum:####
  #########################################
  
  if (!is.null(query_spectrum)){
    
    if (ncol(query_spectrum)<2){
      stop("Spectrum must have 2 columns m/z and intensity!")
    }
    
    dat = query_spectrum[,1:2]
    dat[,2]=dat[,2]/max(dat[,2])*100 # Force normalization
    
    scan_max = max(as.numeric(metadata$SCANS)) # Most recent Scan selected
    ind = which(as.numeric(metadata$SCANS)==scan_max)
    if (length(ind)==0){ind = sample(1:nrow(metadata), 1)}
    
    metadata = metadata[ind,,drop=FALSE]
    spectrum = spectrum_list[[ind]]
    
    bottom.label = paste0("Library spectrum of ID = ", metadata$ID)
    xlim = c(min(dat[,1])*0.8, max(dat[,1])*1.2)
    SpectrumSimilarity(dat, spectrum, top.label="Query spectrum",bottom.label= bottom.label,xlim=xlim, print.graphic = T)
  }
}