R/mzpick.R

In blosloos/enviPick: Peak Picking for High Resolution Mass Spectrometry Data

mzpick <-
function(           MSlist,
                    minpeak=4,
                    drtsmall=20,
                    drtfill=10,
                    drttotal=200,
                    recurs=4,
                    weight=2,
                    SB=3,
                    SN=2,
                    minint=1E4,
                    maxint=1E7,
                    ended=2,
					get_mass = "wmean",
                    progbar=FALSE,
                    from=FALSE,
                    to=FALSE
                    ){

    ############################################################################
    # check inputs #############################################################
	if(all(is.na(match(get_mass, c("mean", "wmean"))))) stop("Wrong wmean argument")
    if(!length(MSlist)==8) stop("This is not an MSlist object")
    if(minpeak <= 0) stop("minpeak must be > 0!")
    if(drtsmall <= 0 || drtsmall<=0) stop("drt must be > 0!")
    if(!is.loaded("pickpeak")) stop(".Call to pickpeak failed!")
    if(!is.loaded("gapfill")) stop(".Call to gapfill failed!")
    if(!is.loaded("picklist")) stop(".Call to picklist failed!")
    if(minint >= maxint){stop("Revise your minint & maxint settings!")}
    if(ended < 1){stop("Wrong ended argument!")}
    if(!is.integer(recurs)){recurs <- ceiling(recurs); recurs <- as.integer(recurs)}
    if(!is.integer(ended)){ended <- ceiling(ended); ended <- as.integer(ended)}
    if(!is.integer(minpeak)){minpeak <- ceiling(minpeak); minpeak <- as.integer(minpeak)}
    if(!is.list(MSlist)) stop("Invalid MSlist argument!")
    if(!MSlist[[1]][[3]]) stop("MSlist invalid. Use mzMLtoMSlist, mzpart and mzclust first.")
    n <- if(!to) nrow(MSlist[[6]]) else to
    m <- if(!from)1 else from   
    ############################################################################
    # reset if rerun ###########################################################
    MSlist[[7]] <- MSlist[[8]] <- 0
    MSlist[[4]][[2]][, 7] <- rep(0,length(MSlist[[4]][[2]][, 4]))
    ############################################################################
    if(progbar){
      prog <- winProgressBar("EIC peak detection",min=m,max=n);
      setWinProgressBar(prog, 0, title = "EIC peak detection", label = NULL);
    }
    startat <- peaknumb <- c(0)
    MSlist[[4]][[2]][, 4] <- seq(1, length(MSlist[[4]][[2]][, 4]), 1); 
    for(i in m:n){    
		if(MSlist[[6]][i, 3] >= minpeak){
		  if(progbar == TRUE){setWinProgressBar(prog, i, title = "EIC  peak detection", label = NULL);}
		  # same RT, interpolate ###################################################
		  out1 <- .Call("gapfill",
			as.numeric(MSlist[[4]][[2]][MSlist[[6]][i, 1]:MSlist[[6]][i, 2], 3]),
			as.numeric(MSlist[[4]][[2]][MSlist[[6]][i, 1]:MSlist[[6]][i, 2], 2]),
			as.integer(order(MSlist[[4]][[2]][MSlist[[6]][i,1]:MSlist[[6]][i, 2], 3], decreasing = FALSE)),
			as.numeric(MSlist[[4]][[2]][MSlist[[6]][i, 1]:MSlist[[6]][i, 2], 1]),
			as.numeric(MSlist[[4]][[2]][MSlist[[6]][i, 1]:MSlist[[6]][i, 2], 4]),
			as.numeric(MSlist[[4]][[1]]),
			as.numeric(drtfill),
			PACKAGE = "enviPick"
			)
		  out1 <- matrix(out1, ncol = 10)
		  colnames(out1) <- c("m/z", "intens", "RT", "index", "intens_filt", "1pick", "pickcrit", "baseline", "intens_corr", "2pick")               
		  # Filter step ############################################################
		  out1[, 5] <- out1[, 2]   # yet to be implemented!
		  # 1st peak detection & baseline subtraction & 2nd peak detection ##########  
		  out2 <- .Call("pickpeak",          
			as.numeric(out1),
			as.numeric(drtsmall),
			as.numeric(drttotal),
			as.integer(minpeak),
			as.integer(recurs),
			as.numeric(weight),   # weight
			as.numeric(SB),       # SB 
			as.numeric(SN),       # SN                                                                           
			as.numeric(minint),   # minimum intensity
			as.numeric(maxint),   # maximum intensity threshold
			as.integer(ended),
			as.integer(2),
			as.integer(1),
			as.integer(1),
			as.integer(0),
			PACKAGE = "enviPick"
		  )     
		  out2 <- matrix(out2, ncol = 10);
		  colnames(out2) <- c("m/z","intens","RT","index","intens_filt","1pick","pickcrit","baseline","intens_corr","2pick");
		  # assign final entries ####################################################
		  if(!all(out2[, 10] == 0)){      
			peaknumb <- peaknumb + max(out2[,10]);
			out2[,10] <- out2[, 10] + startat;         
			for(k in 1:length(out2[, 10])){ # very unelegant approach to remove gap-filling -> revise!
				if(out2[k, 10] != startat){
					MSlist[[4]][[2]][out2[k, 4], 7] <- out2[k, 10]
				}
			}        
			startat <- c(max(out2[,10]));
			MSlist[[4]][[2]][MSlist[[6]][i, 1]:MSlist[[6]][i, 2],] <-
				MSlist[[4]][[2]][MSlist[[6]][i, 1]:MSlist[[6]][i, 2],][order(
					MSlist[[4]][[2]][MSlist[[6]][i, 1]:MSlist[[6]][i, 2], 7], decreasing = FALSE),];   
		  }
		}
    }              
    ############################################################################
	maxit <- max(MSlist[[4]][[2]][,7]);
    # generate peak ID table ###################################################
    if(maxit){
		if(progbar == TRUE){setWinProgressBar(prog, i, title = "Generate index", label = NULL);}
		index <- .Call("indexed",
			as.integer(MSlist[[4]][[2]][, 7]),
			as.numeric(MSlist[[4]][[2]][, 2]),
			as.integer(minpeak),
			as.numeric(maxint),
			as.integer(maxit),
			PACKAGE="enviPick"
		)
		if(any(index[,2]!=0)){
			index <- index[index[,2] != 0,, drop = FALSE];
			partID <- .Call("partID",
				as.integer(index),
				as.integer(length(MSlist[[4]][[2]][, 7])),
				PACKAGE = "enviPick"  
			)
			MSlist[[4]][[2]][, 7] <- partID
			colnames(index) <- c("start_ID", "end_ID", "number_peaks")
			MSlist[[7]]<-index
		}
	}
    ############################################################################
    maxit <- max(MSlist[[4]][[2]][, 7])
	# generate peaklist ########################################################
    if(maxit){
 	  peaklist <- matrix(0, ncol = 11, nrow = maxit)
      colnames(peaklist) <- c("m/z", "var_m/z", "max_int", "sum_int", "RT", "minRT", "maxRT", "part_ID", "EIC_ID", "peak_ID", "Score")
      if(progbar == TRUE){setWinProgressBar(prog, i, title = "Generate peak table", label = NULL);}
      for(i in 1:length(MSlist[[7]][, 1])){
		if(get_mass == "mean"){
			peaklist[i, 1] <- mean(MSlist[[4]][[2]][MSlist[[7]][i, 1]:MSlist[[7]][i, 2], 1])
        }
		if(get_mass == "wmean"){
			peaklist[i,1] <- weighted.mean(
				x = (MSlist[[4]][[2]][MSlist[[7]][i, 1]:MSlist[[7]][i, 2], 1]),
				w = (MSlist[[4]][[2]][MSlist[[7]][i, 1]:MSlist[[7]][i, 2], 2])
			)
        }		
		peaklist[i, 2] <- var(MSlist[[4]][[2]][MSlist[[7]][i, 1]:MSlist[[7]][i, 2], 1])
        peaklist[i, 3] <- max(MSlist[[4]][[2]][MSlist[[7]][i, 1]:MSlist[[7]][i, 2], 2])
        peaklist[i, 4] <- sum(MSlist[[4]][[2]][MSlist[[7]][i, 1]:MSlist[[7]][i, 2], 2])
        peaklist[i, 5] <- (MSlist[[4]][[2]][MSlist[[7]][i, 1]:MSlist[[7]][i,2], 3][MSlist[[4]][[2]][MSlist[[7]][i, 1]:MSlist[[7]][i, 2], 2]==max(MSlist[[4]][[2]][MSlist[[7]][i, 1]:MSlist[[7]][i, 2], 2])[1]])[1]
		peaklist[i, 6] <- min(MSlist[[4]][[2]][MSlist[[7]][i, 1]:MSlist[[7]][i, 2], 3])
        peaklist[i, 7] <- max(MSlist[[4]][[2]][MSlist[[7]][i, 1]:MSlist[[7]][i, 2], 3])        
        peaklist[i, 8:10] <- MSlist[[4]][[2]][MSlist[[7]][i, 1], 5:7]
        peaklist[i, 11] <- 0
      }
      peaklist <- peaklist[order(peaklist[, 3],decreasing = TRUE),, drop = FALSE];
      MSlist[[8]] <- peaklist;
	  MSlist[[3]][[6]] <- length(peaklist[, 1])
    }else{
	  MSlist[[3]][[6]] <- "No peaks picked"
	}
    if(progbar == TRUE){close(prog)}
    ############################################################################
    MSlist[[1]][[5]] <- TRUE;
    MSlist[[3]][[7]] <- length(MSlist[[4]][[2]][MSlist[[4]][[2]][,7]!=0,2])  
    ############################################################################
	MSlist[[2]][[2]][22] <- as.character(minpeak)
	MSlist[[2]][[2]][23] <- as.character(drtsmall) 	
	MSlist[[2]][[2]][24] <- as.character(drtfill) 	
	MSlist[[2]][[2]][25] <- as.character(drttotal) 	
	MSlist[[2]][[2]][26] <- as.character(recurs) 	
	MSlist[[2]][[2]][27] <- as.character(weight) 	
	MSlist[[2]][[2]][28] <- as.character(SB) 	
	MSlist[[2]][[2]][29] <- as.character(SN) 	
	MSlist[[2]][[2]][30] <- as.character(minint) 	
	MSlist[[2]][[2]][31] <- as.character(maxint) 	
	MSlist[[2]][[2]][32] <- as.character(ended) 	
	MSlist[[2]][[2]][33] <- as.character(from) 	
	MSlist[[2]][[2]][34] <- as.character(to) 	
    ############################################################################
    return(MSlist);

}