inst/ShinyApp/DeltaMS/data/DeltaMSfunctions/plotMassSpectrum_doubleLabel.R

In Pohnert-Lab/DeltaMS: Isotopologue finder

# Jan 2017
# Plot MassSepectrum  
# plots isotopologue peaks from a label report to PDF named outputfile
# Highligts the plots of the selected isotopic ratio (?)
# Distance between Peaks in isotopologue group (28.01.2016)
plotMassSpectrum_doubleLabel <- function(isoLabelReport, rawFiles, iRatio, errRatio, disPeak, maxNumDelta, RTwin, maxLT, infoList, classes, labeledSamples, labIso, outputfile) {
  # isoLabelReport = Output of getIsoLabelReport_doubleLabel().
  # rawFiles = Names of the raw files
  # iRatio = Ratio of the added isotopes
  # errRatio = Relative error of the isotope ratio
  # #  If the found peaks correspond to the predetermined isotope ratio, they are marked pink.
  # disPeak = distance between peaks # 28.01.2016 # Peaks with distance greater "disPeak*maxNumDelta" are not plotted
  # maxNumDelta = number of peaks 
  # RTwin = timeframe of RT 
  # #maxLT = Maximum number of peaks of a group to be plotted. # 14-02-2017 
  # labIso = Labels of Isotopes; used as legend
  # classes = names of classes (intern)
  # labelSamples = class of the labeled sample (intern)
  if (is.null(maxNumDelta))  {
    maxNumDelta <- 10000;
  } 
  if (is.null(disPeak))  {
    disPeak <- 1;
  } 
  #  Maximum number of peaks to plot:  maxLT (14-02-17) 
  if (is.null(maxLT))  {
    maxLT <- 2;
  } 
  
  # isotopic ratio # 18.01.2016
  if (is.null(iRatio))  {
    iRatio <- c(0.5, 0.5);
  } 
  if (is.null(errRatio))  {
    errRatio <-0.05;
  } 
  if (is.null(RTwin))  {
    RTwin <-10;
  } 
  # legende # 10.05.2017
  if (is.null(labIso))  { 
    labIso <- c("I1", "I2");
  } 
  if (is.null(infoList))  {
    return("Error, Parameter 'infoList' is need! Perform first the function 'InfoOsoList'")
  } 
  className <-  as.matrix(classes)  
# fixed configuration paramter
# Color of the "correct" peaks   
  col1 <- "darkblue"
#  
  #  Classification of Isotopic pattern
  colI1 <- "darkorchid4"
  colI2 <- "darkorchid1"
  
  # Color of the "not-correct" peaks   
  colN1 <- "gray40"
  colN2 <- "gray80" # "lightblue" #
  
  
  # Color of the "correct" peaks   
  col1 <- "darkblue"
  col2 <- "blue" # "lightblue" #
  #  Classification of Isotopic pattern
  colI1 <- "darkorchid4"
  colI2 <-  "orchid1" #"darkorchid1"
  
#  colfuncN <- colorRampPalette(c(colN1, colN2))
  colfuncI <- colorRampPalette(c(colI1, colI2))
 # colfunc <- colorRampPalette(c(col1, col2))
  
## end: fixed ...  
  

  # Number of Isotopic pattern
  numPat <- length(iRatio)
  
  # Number of the calculated Isotopic pattern
  calcPat <- dim(infoList$iPat)[2]
  
  # Number useed
  curPat <- min(c(numPat, calcPat, maxLT))
  
  
  # Sort the data for the output
  outIDX <- order(infoList$infoMat[,1], infoList$infoMat[,2],infoList$infoMat[,3], infoList$infoMat[,4],infoList$infoMat[,5], decreasing = c(FALSE, FALSE, TRUE, TRUE,TRUE), method="radix")

  # Count raw-files
  numRaw <- length(rawFiles)
  listRaw <- list()
  for (i in 1:numRaw){
    listRaw[[i]] <- xcmsRaw(rawFiles[i])
  }
  # cCol <- max(c(3, numRaw))
  cCol <- 3
  numPlot <- cCol * 5
  #
  pdf(outputfile, width = 8.5, height = 11);
 # x11()
  numGroups = length(isoLabelReport[[1]]);
  par(mfrow=c(5,cCol));
  iz <- 1
  
  # for all isotopic patterns
  for (ic in 1 : numPat) {
    if (ic > curPat) { #Currently possible number
      break
    }
    if (ic == 1) {
 #     par(mfrow=c(4,3));
      par(mfrow=c(5,cCol));
    }
    toPlot <- (as.matrix(iRatio[[ic]]))
    rWord <-  c("m/z")
    rNames <- c("m/z")
    
    legTxt <- paste(ic,"x",labIso[1])
    
    for (i in 1:(ic)) {
      rNames <- cbind(rNames, paste(rWord, "+", format(disPeak*i, digits = 4)) )
      if (i == ic) {
        legTxt <- cbind(legTxt, paste(ic,"x",labIso[2]))
      }else{
        legTxt <- cbind(legTxt, paste(ic-i,"x",labIso[1], ",", i,"x",labIso[2]))
      }
    }
    
    # Output of the pattern to be identified
    YLIM <- c(0, 1)
    mp = barplot(t(toPlot), beside = TRUE, legend.text = legTxt, 
                 main = paste("Level ", ic, "\n", "Isotopic pattern"), 
                 axisnames = FALSE, ylim = YLIM , col  =  colfuncI(ic+1)); #= c(colI1, colI2)); #, col = c("dark red", "grey")
    text(seq(1.5, (ic)*2+1.5, by = 2), pos=2, par("usr")[3]-0.15, xpd = TRUE, labels = rNames, srt = 45);
    iz <- iz + 1
 #  Format of the output-file new line or not ?     
 #   for (k in 2:cCol) {
#      plot(c(0, 1), c(0, 1), type="l")
#      lines(c(0, 1), c(1, 0), type="l")
#      iz <- iz +1
#    }

    for (i in 1:numGroups) {
       # Each group starts with a new line
       if (!(iz %% cCol == 0)) {
          iz <- iz+ cCol - (iz %% cCol)
       } 
       if (iz %% numPlot == 1) {
          par(mfrow=c(5,cCol));
       }
       # current index
       cIdx <- outIDX[i]
#       if (i == 6 ){
#          print(i)
#       }
       #
       # Determination of the data for the output
       toPlot = cbind(isoLabelReport$meanAbsU[[cIdx]]);
       sds = cbind(isoLabelReport$cvTotalU[[cIdx]]*isoLabelReport$meanAbsU[[cIdx]]);
       YLIM = c(0, max(toPlot));
       rat12 <- toPlot/sum(toPlot)
       
       #
       numPeaks = dim(toPlot)[1];

    mzPeaks <- isoLabelReport$isotopologue[[cIdx]]
    rtPeaks <- isoLabelReport$rt[[cIdx]]
    rownames(toPlot) = as.character(round(isoLabelReport$isotopologue[[cIdx]], 3));
    
    # RT-window as parameter ?    
    for (k in 1:numRaw) {
      iz <- iz + 1
      idRT <- which(listRaw[[k]]@scantime>rtPeaks[1]-0.05*RTwin & listRaw[[k]]@scantime<rtPeaks[numPeaks]+0.05*RTwin)
      if (length(idRT) < 1) {
        idRT <- which(listRaw[[k]]@scantime>rtPeaks[1]-0.2*RTwin & listRaw[[k]]@scantime<rtPeaks[numPeaks]+0.2*RTwin)
        if (length(idRT) < 1) {
          idRT <- which(listRaw[[k]]@scantime>rtPeaks[1]-0.5*RTwin & listRaw[[k]]@scantime<rtPeaks[numPeaks]+0.5*RTwin)
        }
        if (length(idRT) < 1) {
          idRT <- which(listRaw[[k]]@scantime>rtPeaks[1]-0.75*RTwin & listRaw[[k]]@scantime<rtPeaks[numPeaks]+0.75*RTwin)
        }
      }
      if (length(idRT) < 1) {
        plot(c(0, 1), c(0, 1), type="l")
        lines(c(0, 1), c(1, 0), type="l")
        next
      }
      df <- disPeak

      ## handle last spectrum
      if (idRT[1] == length(listRaw[[k]]@scanindex)) {
        followingScanIndex <- length(listRaw[[k]]@env$mz)
      } else {
        followingScanIndex <- listRaw[[k]]@scanindex[idRT[1]+1]
      }
      
      idx <- (listRaw[[k]]@scanindex[idRT[1]]+1):min(followingScanIndex,
                                            length(listRaw[[k]]@env$mz), na.rm=TRUE)
      mzrange <- c(mzPeaks[1]-df, mzPeaks[numPeaks]+df)
      if (length(mzrange) >= 2) {
        mzrange <- range(mzrange)
        idx <- idx[listRaw[[k]]@env$mz[idx] >= mzrange[1] & listRaw[[k]]@env$mz[idx] <= mzrange[2]]
      }
      if ((((infoList$infoMat[cIdx,1] == 1) && infoList$infoMat[cIdx,2] == (ic+1) && infoList$infoMat[cIdx,4] == 1))
          || (((infoList$infoMat[cIdx,1] == 1) && infoList$infoMat[cIdx,5] == (ic+1) ))) {
        
          points <- cbind(listRaw[[k]]@env$mz[idx], listRaw[[k]]@env$intensity[idx])
#          title = paste("#",k, ": ", round(listRaw[[k]]@scantime[idRT[1]]/60, 2),
#                    " min (scan ",idRT[1], ")", sep = "")
          title = paste("#", i, ": ", className[k],"\n", "#",k, ": ", round(listRaw[[k]]@scantime[idRT[1]]/60, 2),
                        " min (scan ",idRT[1], ")", sep = "")
          
          if (length(idx) > 0) {
              plot(points, type="h", main = title, xlab="m/z", ylab="Intensity")
          } else {
            plot(c(0, 1), c(0, 1), type="l", main = title, xlab="m/z", ylab="Intensity")
            lines(c(0, 1), c(1, 0), type="l")
          }
      }    
      
#     
      
      scanVal <- getScan(listRaw[[k]], idRT[1], mzrange = c(mzPeaks[1]-df, mzPeaks[numPeaks]+df))
      if (length(scanVal) == 0) {
        next
      }
      
#      hv1 <- scanVal[,1]- mzPeaks[1]
#      i1 <- scanVal[which.min(abs(hv1)),2]
#      hv2 <- scanVal[,1]- mzPeaks[2]
#      i2 <- scanVal[which.min(abs(hv2)),2]
      iA <- matrix(data=0, nrow=1, ncol=numPeaks)
      for (kk in 1:numPeaks){
        hvA <- scanVal[,1]- mzPeaks[kk]
        iA[kk] <- scanVal[which.min(abs(hvA)),2]
        
      }
      if ((((infoList$infoMat[cIdx,1] == 1) && infoList$infoMat[cIdx,2] == (ic+1) && infoList$infoMat[cIdx,4] == 1))
          || (((infoList$infoMat[cIdx,1] == 1) && infoList$infoMat[cIdx,5] == (ic+1) ))) {
        iz <- iz + 1
        if ((infoList$infoMat[cIdx,3] == 1) ) { # || (infoList$iPat[cIdx,ic] == 1)) { # No search of current pattern in other groups
            for (k in 1: numPeaks) {
             text(mzPeaks[k], iA[k], rownames(toPlot)[k], col = c(col1))
             lines(c(mzPeaks[k], mzPeaks[k]), c(0,iA[k]),  col = c(col1), type="l", lwd = 2)
           }
        } else {
           for (k in 1: numPeaks) {
             text(mzPeaks[k], iA[k], rownames(toPlot)[k])
             lines(c(mzPeaks[k], mzPeaks[k]), c(0,iA[k]),  type="l", lwd = 1)
           }
        }
      }
    }     
    
    
  }
  }
  ################################################
  # the remaining groups of peaks
  notUse <- 1 # makes no sense
  if (!notUse) {
  iz <- 1
  maxDist <- max(infoList$infoMat[,1]) # Maximum distance of the peaks (factor)
  for (ic in 2 : maxDist) {
    if (ic == 2) {
      #     par(mfrow=c(4,3));
      par(mfrow=c(5,cCol));
    }
    toPlot <- (as.matrix(iRatio[[1]]))
    rWord <-  c("m/z")
    rNames <- c("m/z")
    
    for (i in 1) {
      rNames <- cbind(rNames, paste(rWord, "+", format(disPeak*ic, digits = 4)) )
    }
    
    YLIM <- c(0, 1)
    mp = barplot(t(toPlot), beside = TRUE, legend.text = c("I1", "I2"), 
                 main = paste("Distance ", ic, "\n", rNames[2]), 
                 axisnames = FALSE, ylim = YLIM , col = c(colI1, colI2)); #, col = c("dark red", "grey")
    text(seq(1.5, (ic)*2+1.5, by = 2), pos=2, par("usr")[3]-0.15, xpd = TRUE, labels = rNames, srt = 45);
    iz <- iz + 1
    
    #  Format of the output-file new line or not ?     
#    for (k in 2:cCol) {
#      plot(c(0, 1), c(0, 1), type="l")
#      lines(c(0, 1), c(1, 0), type="l")
#      iz <- iz +1
#    }
    
    for (i in 1:numGroups) {
      # Each group starts with a new line
      if (!(iz %% cCol == 0)) {
        iz <- iz+ cCol - (iz %% cCol)
      } 
      if (iz %% numPlot == 1) {
        par(mfrow=c(5,cCol));
      }
      # current index
      cIdx <- outIDX[i]
      if (infoList$infoMat[cIdx,1] != ic) {
        next
      }
      # Determination of the data for the output
      toPlot = cbind(isoLabelReport$meanAbsU[[cIdx]]);
      sds = cbind(isoLabelReport$cvTotalU[[cIdx]]*isoLabelReport$meanAbsU[[cIdx]]);
      YLIM = c(0, max(toPlot));
      rat12 <- toPlot/sum(toPlot)
      
      #
      numPeaks = dim(toPlot)[1];

      mzPeaks <- isoLabelReport$isotopologue[[cIdx]]
      rtPeaks <- isoLabelReport$rt[[cIdx]]
      rownames(toPlot) = as.character(round(isoLabelReport$isotopologue[[cIdx]], 3));
      
      # RT-window as parameter ?    
      for (k in 1:numRaw) {
        iz <- iz + 1
        idRT <- which(listRaw[[k]]@scantime>rtPeaks[1]-0.05*RTwin & listRaw[[k]]@scantime<rtPeaks[numPeaks]+0.05*RTwin)
        if (length(idRT) < 1) {
          idRT <- which(listRaw[[k]]@scantime>rtPeaks[1]-0.2*RTwin & listRaw[[k]]@scantime<rtPeaks[numPeaks]+0.2*RTwin)
          if (length(idRT) < 1) {
            idRT <- which(listRaw[[k]]@scantime>rtPeaks[1]-0.5*RTwin & listRaw[[k]]@scantime<rtPeaks[numPeaks]+0.5*RTwin)
          }
          if (length(idRT) < 1) {
            idRT <- which(listRaw[[k]]@scantime>rtPeaks[1]-0.75*RTwin & listRaw[[k]]@scantime<rtPeaks[numPeaks]+0.75*RTwin)
          }
        }
        if (length(idRT) < 1) {
          plot(c(0, 1), c(0, 1), type="l")
          lines(c(0, 1), c(1, 0), type="l")
          next
        }
        df <- disPeak

        ## handle last spectrum
        if (idRT[1] == length(listRaw[[k]]@scanindex)) {
          followingScanIndex <- length(listRaw[[k]]@env$mz)
        } else {
          followingScanIndex <- listRaw[[k]]@scanindex[idRT[1]+1]
        }
        
        idx <- (listRaw[[k]]@scanindex[idRT[1]]+1):min(followingScanIndex,
                                                       length(listRaw[[k]]@env$mz), na.rm=TRUE)
        mzrange <- c(mzPeaks[1]-df, mzPeaks[numPeaks]+df)
        if (length(mzrange) >= 2) {
          mzrange <- range(mzrange)
          idx <- idx[listRaw[[k]]@env$mz[idx] >= mzrange[1] & listRaw[[k]]@env$mz[idx] <= mzrange[2]]
        }

          points <- cbind(listRaw[[k]]@env$mz[idx], listRaw[[k]]@env$intensity[idx])
#          title = paste("#",k, ": ", round(listRaw[[k]]@scantime[idRT[1]]/60, 2),
#                        " min (scan ",idRT[1], ")", sep = "")
          title = paste("#", i, ": ", className[k],"\n", "#",k, ": ", round(listRaw[[k]]@scantime[idRT[1]]/60, 2),
                        " min (scan ",idRT[1], ")", sep = "")
          
  #        plot(points, type="h", main = title, xlab="m/z", ylab="Intensity")
          
          if (length(idx) > 0) {
            plot(points, type="h", main = title, xlab="m/z", ylab="Intensity")
          } else {
            plot(c(0, 1), c(0, 1), type="l", main = title, xlab="m/z", ylab="Intensity")
            lines(c(0, 1), c(1, 0), type="l")
          }
          
 
        scanVal <- getScan(listRaw[[k]], idRT[1], mzrange = c(mzPeaks[1]-df, mzPeaks[numPeaks]+df))
        if (length(scanVal) == 0) {
          next
        }
        
        
#        hv1 <- scanVal[,1]- mzPeaks[1]
#        i1 <- scanVal[which.min(abs(hv1)),2]
#        hv2 <- scanVal[,1]- mzPeaks[2]
#        i2 <- scanVal[which.min(abs(hv2)),2]
        iA <- matrix(data=0, nrow=1, ncol=numPeaks)
        for (k in 1:numPeaks){
          hvA <- scanVal[,1]- mzPeaks[k]
          iA[k] <- scanVal[which.min(abs(hvA)),2]
          
        }
        for (k in 1: numPeaks) {
          text(mzPeaks[k], iA[k], rownames(toPlot)[k])
         lines(c(mzPeaks[k], mzPeaks[k]), c(0,iA[k]),  type="l", lwd = 1)
         }
      }     
    }
  }
  } # bDel
  
  
  
  ################################################
  
  dev.off();
}