#' get the xcmsRaw lighter
#'
#' @name eicXcmsRaw
#' @rdname eicXcmsRaw-methods
#' @import xcms
#' @export
setGeneric(name="eicXcmsRaw",def=function(object,eicmat,inthr) standardGeneric("eicXcmsRaw"))
#' @rdname eicXcmsRaw-methods
#' @aliases eicXcmsRaw
setMethod("eicXcmsRaw","xcmsRaw",
eicXcmsRaw<-function(object,eicmat,inthr=0){
ascans=do.call("rbind",lapply(1:nrow(eicmat),function(i)
.GRrawMat(object,mzrange = eicmat[i,c("mzmin","mzmax")],rtrange = eicmat[i,c("rtmin","rtmax")]*60)[,1:4,drop=F]))
ascans=ascans[order(ascans[,1],ascans[,2]),,drop=FALSE]
ascans=ascans[!duplicated(ascans[,4]),,drop=FALSE]
ascans=ascans[which(ascans[,"y"]>inthr),,drop=FALSE]
scrange=min(ascans[,1]):max(ascans[,1])
sc2time=object@scantime[scrange]
tabsc=table(factor(ascans[,1],levels = scrange))
scanidx=c(0,cumsum(tabsc))
scanidx=as.integer(scanidx[-length(scanidx)])
newtic=round(as.numeric(tapply(ascans[,"y"],ascans[,"scan"],sum)[as.character(scrange)]),1)
newtic[is.na(newtic)]=0
ob<-new("xcmsRaw")
ob@env <- new.env(parent = .GlobalEnv)
ob@env$mz<-as.numeric(round(ascans[,"mz"],6))
ob@env$intensity<-as.numeric(round(ascans[,"y"],1))
ob@scanindex<-scanidx
ob@scantime<-sc2time
ob@acquisitionNum<-1:length(ob@scanindex)
ob@filepath<-object@filepath
ob@mzrange<-range(ob@env$mz)
ob@profmethod<-object@profmethod
ob@tic<-newtic
ob@profparam<-list()
ob<-xcms:::remakeTIC(ob)
return(ob)
}
)
######################################################################################
#' @title Remove artefact in xcmsRaw object
#' get the xcmsRaw lighter
#'
#' @name satremXcmsRaw
#' @rdname satremXcmsRaw-methods
#' @param object xcmsRaw object
#' @param reso resolution
#' @param typ 1/2 2 correrspond to the extended model
#' @param inthr intensity threshold
#' @param rtlim min/max retention time in min.
#' @import xcms
#' @export
setGeneric(name="satremXcmsRaw",def=function(object,reso,typ,inthr,rtlim) standardGeneric("satremXcmsRaw"))
#' @rdname satremXcmsRaw-methods
#' @aliases satremXcmsRaw
setMethod("satremXcmsRaw","xcmsRaw",
satremXcmsRaw<-function(object,reso=70000,typ=2,inthr=0,rtlim=c(-Inf,Inf)){
scsten=cbind(object@scanindex+1,c(object@scanindex[-1],length(object@env$intensity)))
rownames(scsten)=1:length(object@scanindex)
newdat=list()
nrm1=nrm2=0
lsc=which((object@scantime/60)>=rtlim[1] & (object@scantime/60)<=rtlim[2])
cat(" -- processing",length(lsc),"scans\n")
for(ix in lsc){
re<-xcms:::getScan(object,as.numeric(ix))
l2excl=which(re[,2]<inthr)
nrm1=nrm1+length(l2excl)
if(length(l2excl)>0) re=re[-l2excl,,drop=F]
l2excl=satremoval(re,reso,typ)
nrm2=nrm2+length(l2excl)
if(length(l2excl)>0) re=re[-l2excl,,drop=F]
if(nrow(re)==0) next
newdat[[ix]]=cbind(re,scan=ix)
}
ntot=diff(range(scsten[as.character(lsc),]))+1
cat(" -- removed ",nrm1,"/",nrm2," out of",ntot," duples - avg.=",round(100*nrm1/ntot,2),"/",round(100*nrm2/ntot,2)," perc. per scan\n",sep="")
newdat=do.call("rbind",newdat)
newdat=newdat[order(newdat[,"scan"],newdat[,1]),,drop=F]
scrange=min(newdat[,"scan"]):max(newdat[,"scan"])
sc2time=object@scantime[scrange]
tabsc=table(factor(newdat[,"scan"],levels = scrange))
scanidx=c(0,cumsum(tabsc))
scanidx=as.integer(scanidx[-length(scanidx)])
newtic=round(as.numeric(tapply(newdat[,"intensity"],newdat[,"scan"],sum)[as.character(scrange)]),1)
newtic[is.na(newtic)]=0
ob<-new("xcmsRaw")
ob@env <- new.env(parent = .GlobalEnv)
ob@env$mz<-as.numeric(round(newdat[,"mz"],6))
ob@env$intensity<-as.numeric(round(newdat[,"intensity"],1))
ob@scanindex<-scanidx
ob@scantime<-sc2time
ob@acquisitionNum<-1:length(ob@scanindex)
ob@filepath<-object@filepath
ob@mzrange<-range(ob@env$mz)
ob@profmethod<-object@profmethod
ob@tic<-newtic
ob@profparam<-list()
ob<-xcms:::remakeTIC(ob)
return(ob)
}
)
######################################################################################
#' Artfact removal in duples
#' @param re matrix of mz/int
#' @param reso mass resolution
#' @param typ 1/2 extended
#' @return ids of entries to be remove
#' @export
satremoval<-function(re,reso,typ=2){
# print(reso)
if(is.na(reso)) return(c())
m=re[,1]
v=re[,2]
if(typ==1) unique(unlist(lapply(1:length(m),function(i){
m0=m[i]
gamma=m0/reso
v1=v[i]*(gamma^2)/(gamma^2+(m-m0)^2)
which(v<(0.99*v1))})))
if(typ==2) unique(unlist(lapply(1:length(m),function(i){
m0=m[i]
gamma=m0/reso
v1=v[i]*(gamma^2)/(gamma^2+(m-m0)^2)
gamma2=m0/(reso*0.05)
v2=0.05*v[i]*(gamma2^2)/(gamma2^2+(m-m0)^2)
# l=which(v1<(0.05*v[i]))
# v1[l]=v2[l]
which(v<(0.99*v1) | v<(0.99*v2))})))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.