annotate moulinetteJF/retraitement_MSpurity_V4.R @ 0:67733206be53 draft

" master branch Updating"

## read and process mspurity W4M files
## create a summary of fragment for each precursor and a graphics of peseudo spectra + correlation on which checking of fragment is based on
## V3 try to identify and process multiple files for 1 precursor which may occur if different collision energy are used
## V4 elimination of correlation = NA. Correlation is done with precursor, if precursor is not present correlation with most intense peak

##require(msPurity)

###########################################################################################

plot_pseudoSpectra <- function(x, fid, sumInt, vmz, corAbsInt, refcol, cNAME) {
  ## function to compute sum of intensities among scans for all m/z kept (cor > r_threshold & minimum number of scans) 
  ## and plot pseudo spectra
  ## x dataframe scan X fragments with scans number in the 1st column and ions in next with intensities
  ## fid file id when several a precursor has been detected in several files
  # du fait de la difference de nombre de colonne entre la dataframe qui inclue les scans en 1ere col, mzRef se decale de 1
  refcol <- refcol-1
  ## compute relative intensities max=100%
  relInt <- sumInt[-1]
  relInt <- relInt/max(relInt)
  
  ## define max value on vertical axis (need to increase in order to plot the label of fragments)
  ymax <- max(relInt)+0.2*max(relInt)
  
  par(mfrow=c(2,1))
  plot(vmz, relInt, type="h", ylim=c(0,ymax), main = cNAME)
  ## low correl coef. will be display in grey
  corLow <- which(round(corAbsInt,2) < r_threshold)
  
  lbmzcor <- paste(as.character(vmz),"(r=",round(corAbsInt,2),")", sep="")
  
  if (length(corLow) > 0) {
    text(vmz[corLow], relInt[corLow], lbmzcor[corLow], cex=0.5, col = "grey", srt = 90 , adj=0)
    if (length(vmz) - length(corLow) > 1) text(vmz[-c(refcol,corLow)], relInt[-c(refcol,corLow)], lbmzcor[-c(refcol,corLow)], cex=0.6, col = 1, srt = 90 , adj=0)
  } else {
    if (length(vmz) > 1) text(vmz[-c(refcol)], relInt[-c(refcol)], lbmzcor[-c(refcol)], cex=0.6, col = 1, srt = 90 , adj=0)
  }
  
  text(vmz[refcol], relInt[refcol], lbmzcor[refcol], cex=0.8, col = 2, srt = 90, adj=0 )
  
  ## prepare result file
  corValid <- (round(corAbsInt,2) >= r_threshold)
  cpRes <- data.frame(rep(cNAME, length(vmz)),
                      rep(fid, length(vmz)),
                      vmz,corAbsInt,sumInt[-1],relInt,corValid)
  
  colnames(cpRes) <- c("compoundName","fileid","fragments_mz","CorWithPrecursor","AbsoluteIntensity","relativeIntensity","corValid")
  return(cpRes)
  
}

## function for extraction of fragments corresponding to precursors detected by MSPurity
Xtract_fragments <- function(mzref, rtref, cNAME, r_threshold = 0.85, seuil_ra = 0.1, tolmz = 0.01, tolrt = 60) {
  
  ## filter precursor in the precursors file based on mz and rt in the compound list
  cat("processing ",cNAME,"\n")
  selprec <- which((abs(prec$precurMtchMZ - mzref) <= tolmz) & (abs(prec$precurMtchRT - rtref) <= tolrt))
  
  ## check if there is the precursor in the file
  if (length(selprec) > 0) {
    
    sprecini <- prec[selprec,]
    
    ## check if fragments corresponding to precursor are found in several files (collision energy)
    ## this lead to a processing for each fileid 
    mf <- levels(as.factor(sprecini$fileid))
    nbf <- length(mf)
    if (nbf > 1) cat(" several files detected for this compounds :\n")
    
    for (f in 1:nbf) {
      
      sprec <- sprecini[sprecini$fileid == mf[f],]
      
      ## selection of fragment in the fragments file with the grpid common in both fragments and precursors
      selfrgt <- levels(as.factor(sprec$grpid))
      sfrgt <- frgt[frgt$grpid %in% selfrgt & frgt$fileid == mf[f],]
      
      ## filter fragments on relative intensity seuil_ra = user defined parameter (MSpurity flags could be used here)
      sfrgtfil <- sfrgt[sfrgt$ra > seuil_ra,]
      
      mznominal <- round(x = sfrgtfil$mz, mzdecimal)
      sfrgtfil <- data.frame(sfrgtfil, mznominal)
      
      ## creation of cross table row=scan col=mz X=ra
      vmz <- levels(as.factor(sfrgtfil$mznominal))
      #fscan <- levels(as.factor(sfrgtfil$acquisitionNum))
      
      cat(" fragments :",vmz)
      # dsIntra <- matrix(NA,nrow = length(vscan), ncol = length(vmz))
      # rownames(dsIntra) <- fscan
      # dsIntra <- data.frame(fscan,dsIntra)
      # colnames(dsIntra) <- c("fscan",vmz)
      
      ## mz of precursor in data precursor to check correlation with
      mzPrec <- paste("mz",round(mean(sprec$mz),mzdecimal),sep="")
      
      for (m in 1:length(vmz)) {
        
        ## absolute intensity
        cln <- c(which(colnames(sfrgtfil) == "acquisitionNum"), which(colnames(sfrgtfil) == "i"))
        int_mz <- sfrgtfil[sfrgtfil$mznominal == vmz[m], cln]
        colnames(int_mz)[2] <- paste("mz", vmz[m], sep="")
        
        ## average intensities of mass in duplicate scans
        compScans <- aggregate(x = int_mz, by = list(int_mz[[1]]), FUN = mean)
        int_mz <- compScans[,-1]
        
        if (m == 1) {
          dsAbsInt <- int_mz
          #dsRelInt <- ra_mz
        } else {
          dsAbsInt <- merge(x = dsAbsInt, y = int_mz, by.x = 1, by.y = 1, all.x=TRUE, all.y=TRUE)
          #dsRelInt <- merge(x = dsRelInt, y = ra_mz, by.x = 1, by.y = 1, all.x=TRUE, all.y=TRUE)
          
        }
      }
      ## for debug
      ## write.table(x = dsAbsInt,file=paste(cNAME,"dsAbsInt.txt",sep=""), row.names = FALSE, sep="\t")
      
      ## elimination of mz with less than minNumberScan scans (user defined parameter)
      xmz <- rep(NA,ncol(dsAbsInt)-1)
      sumInt <- rep(NA,ncol(dsAbsInt))
      nbxmz <- 0
      NbScanCheck <- min(nrow(dsAbsInt),minNumberScan)
      
      for (j in 2:ncol(dsAbsInt)) {
        sumInt[j] <- sum(dsAbsInt[j],na.rm = TRUE)
        if (sum(!is.na(dsAbsInt[[j]])) < NbScanCheck) {
          nbxmz <- nbxmz + 1
          xmz[nbxmz] <- j
        }
      }
      
      xmz <- xmz[-which(is.na(xmz))]
      if (length(xmz) > 0) {
        dsAbsInt <- dsAbsInt[,-c(xmz)]
        sumInt <- sumInt[-c(xmz)]
        ## liste des mz keeped decale de 1 avec dsAbsInt
        vmz <- as.numeric(vmz[-c(xmz-1)])
      }
 
      ## reference ion for correlation computing = precursor OR maximum intensity ion in precursor is not present
      refcol <- which(colnames(dsAbsInt) == mzPrec)
      if (length(refcol) == 0) refcol <- which(sumInt == max(sumInt, na.rm = TRUE))
      
      pdf(file=paste(cNAME,"_processing_file",mf[f],".pdf",sep=""), width = 8, height = 11 ); par(mfrow=c(3,2))
      
      ## pearson correlations between absolute intensities computing
      corAbsInt <- rep(NA,length(vmz))
      
      if (length(refcol) > 0) {
        for (i in 2:length(dsAbsInt)) {
          corAbsInt[i-1] <- cor(x = dsAbsInt[[refcol]], y=dsAbsInt[[i]], use = "pairwise.complete.obs", method = "pearson")
          plot(dsAbsInt[[refcol]],dsAbsInt[[i]], 
               xlab = colnames(dsAbsInt)[refcol], ylab=colnames(dsAbsInt)[i], 
               main=paste(cNAME," corr coeff r=",round(corAbsInt[i-1],2),sep=""))
        }
        ## plot pseudo spectra
        resCompByfile <- plot_pseudoSpectra(x= dsAbsInt, fid = mf[f], sumInt = sumInt, vmz = vmz, corAbsInt = corAbsInt, refcol = refcol, cNAME = cNAME )
        if (f == 1) resComp <- resCompByfile
      } else {
        resCompByfile <- NULL
        cat(" non detected in fragments file \n")
      }
      if (!is.null(resCompByfile)) resComp <- rbind(resComp,resCompByfile)
      cat("\n")
      dev.off()
    }
  } else {
    resComp <- NULL
    cat(" non detected in precursor file \n")
  }
  
  return(resComp)
}


#################################### start  #######################################################

## FOLDER AND FILES
setwd(dir = "C:/Users/jfmartin/Documents/PROJETS/PROG/tools_R/MSPurity/positif2/V2")

#load("Galaxy541-[msPurity.filterFragSpectra_on_data_540__RData].rdata")

## MSpurity precursors file
prec <- read.table(file = "msPurity.filterFragSpectra_on_data_60__peaklist_precursors.tsv",
                   header = TRUE, sep="\t")

## MSpurity fragments file
frgt <- read.table(file = "msPurity.filterFragSpectra_on_data_60__peaklist_fragments.tsv",
                   header = TRUE, sep="\t")

## list of compounds : col1=Name of molecule, col2=m/z, col3=retention time
compounds <- read.table(file = "compounds_pos.txt", sep="\t", header = TRUE)

## PARAMETERS
## tolerance for mz(dalton) rt(seconds) to match the standard in the compounds list with the precursor MSpurity file
tolmz <- 0.01
tolrt <- 20

##  relative intensity threshold
seuil_ra = 0.05
## nb decimal for mz
mzdecimal <- 0
## r pearson correlation threshold between precursor and fragment absolute intensity
r_threshold <- 0.85
## fragments are kept if there are found in a minimum number of scans
minNumberScan <- 8

for (i in 1:nrow(compounds)) {
  ## loop execution for all compounds in the compounds file
  resCor <- NULL
  resCor <- Xtract_fragments(mzref = compounds[[2]][i], 
                             rtref = compounds[[3]][i] * 60,
                             cNAME = compounds[[1]][i],
                             r_threshold, seuil_ra, tolmz, tolrt)
  if (i == 1 & !is.null(resCor)) resAll <- resCor else if (!is.null(resCor)) resAll <- rbind(resAll,resCor)

  }

write.table(x = resAll, file = "compound_fragments_result.txt", sep="\t", row.names = FALSE)

########################################  end call ################################################