R/mzclust.R

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

mzclust <-
function(           MSlist,
                    dmzdens=10,
                    ppm=TRUE,
                    drtdens=60,
                    minpeak=4,
                    maxint=1E6,
                    progbar=FALSE,
                    merged=TRUE,
                    from=FALSE,
                    to=FALSE
                     ){

    ############################################################################
    if(!is.logical(ppm)){stop("invalid ppm argument, TRUE or FALSE")}
    if(!is.logical(merged)){stop("invalid merged argument, TRUE or FALSE")}
    if(!is.logical(progbar)){stop("invalid progbar argument, TRUE or FALSE")}
    if(!length(MSlist)==8){stop("This is not an MSlist object")}
    if(minpeak<1){stop("invalid minpeak argument")}
    if(dmzdens<=0 || drtdens<=0){stop("invalid: some argument <= 0")}
    if(!is.loaded("getEIC")){stop(".Call to getEIC failed")}
    if(!MSlist[[1]][[2]]){stop("MSlist invalid. Use mzMLtoMSlist and mzpart to upload and partition .mzML data first.")}
    if(ppm){ppm2=1}else{ppm2=2};
    if(merged){merged2=1}else{merged2=2};
    if(!to){n=length(MSlist[[5]][,1])}else{n=to}
    if(!from){m=1}else{m=from}
    ############################################################################
    # reset if rerun ###########################################################
    MSlist[[6]]<-0;
    MSlist[[7]]<-0;
    MSlist[[8]]<-0;
    MSlist[[4]][[2]][,6]<-rep(0,length(MSlist[[4]][[2]][,4]));
    MSlist[[4]][[2]][,7]<-rep(0,length(MSlist[[4]][[2]][,4]));
    ############################################################################
    # cluster ################################################################
    if(progbar==TRUE){  prog <- winProgressBar("EIC  clustering", min = m, max = n);
                        setWinProgressBar(prog, 0, title = "EIC  clustering", label = NULL);}    
    startat<-c(0);    
    for(k in m:n){ 
      if(progbar==TRUE){setWinProgressBar(prog, k, title = paste0("EIC  clustering: ID", k, " @", (MSlist[[5]][k,2] - MSlist[[5]][k,1] + 1)," measurements"), label = NULL)}
      if((MSlist[[5]][k,2] - MSlist[[5]][k,1] + 1) > 1){      
        clusters <- .Call("getEIC",
                          as.numeric(MSlist[[4]][[2]][MSlist[[5]][k,1]:MSlist[[5]][k,2],1]),       # mz
                          as.numeric(MSlist[[4]][[2]][MSlist[[5]][k,1]:MSlist[[5]][k,2],3]),       # RT
                          as.numeric(MSlist[[4]][[2]][MSlist[[5]][k,1]:MSlist[[5]][k,2],2]),       # intens
                          as.integer(order(MSlist[[4]][[2]][MSlist[[5]][k,1]:MSlist[[5]][k,2],2], decreasing = TRUE)),  # intensity order
                          as.integer(order(MSlist[[4]][[2]][MSlist[[5]][k,1]:MSlist[[5]][k,2],3], decreasing = FALSE)),  # RT order                          
                          as.numeric(dmzdens),
                          as.integer(ppm2),
                          as.numeric(drtdens),                                                   
                          as.integer(merged2),
                          PACKAGE = "enviPick"
                        )
        clusters[,10] <- clusters[,10] + startat;             
        MSlist[[4]][[2]][(MSlist[[5]][k,1]):(MSlist[[5]][k,2]),6] <- clusters[,10]                       
        MSlist[[4]][[2]][(MSlist[[5]][k,1]):(MSlist[[5]][k,2]),] <- (MSlist[[4]][[2]][(MSlist[[5]][k,1]):(MSlist[[5]][k,2]),][order(clusters[,10],decreasing=FALSE),])
        startat <- c(max(clusters[,10]))
      }else{
        MSlist[[4]][[2]][(MSlist[[5]][k,1]):(MSlist[[5]][k,2]),6] <- startat;   
        startat <- startat + 1;   
      }
    }    
    if(progbar == TRUE){close(prog)};
    maxit <- max(MSlist[[4]][[2]][,6]);
    if(maxit > 0){
		index <- .Call("indexed",
			as.integer(MSlist[[4]][[2]][,6]),
			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];
			colnames(index) <- c("start_ID", "end_ID", "number_peaks");
			MSlist[[6]] <- index;
			partID <- .Call("partID",
				as.integer(index),
				as.integer(length(MSlist[[4]][[2]][,6])),
				PACKAGE = "enviPick"  
			)
			MSlist[[4]][[2]][,6] <- partID;
		}
	}
    ############################################################################
    MSlist[[1]][[3]]<-TRUE;
    MSlist[[1]][[4]]<-FALSE;
    MSlist[[1]][[5]]<-FALSE;
    MSlist[[3]][[4]]<-length(index[,3]);
    MSlist[[3]][[5]]<-length(MSlist[[4]][[2]][MSlist[[4]][[2]][,6]!=0,6]);
    MSlist[[4]][[2]][,4]<-seq(1,length(MSlist[[4]][[2]][,4]),1);
    ############################################################################
    MSlist[[2]][[2]][14]<-as.character(dmzdens)
    MSlist[[2]][[2]][15]<-as.character(ppm)
    MSlist[[2]][[2]][16]<-as.character(drtdens)
	MSlist[[2]][[2]][17]<-as.character(minpeak)
    MSlist[[2]][[2]][18]<-as.character(maxint)	
    MSlist[[2]][[2]][19]<-as.character(merged)
    MSlist[[2]][[2]][20]<-as.character(from)
    MSlist[[2]][[2]][21]<-as.character(to)	
	############################################################################
    return(MSlist)
}