### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ###
#'
#' Weighting around the precursor window -> to be improved with larger windo etc...
#' TO DO: align properly the deconv/peak picking with the EIC integration!!!
#'
#' @param eic vector of mz-mzPrec
#' @param eicbl vector of mz-mzPrec for the blank samples
#' @param parDeco weight
#' @param imain for the plot tile
#' @param doPlot weight
#' @keywords internal
#'
#' @export
.MGrefFct0<-function(eic,eicbl=NULL,parDeco,imain=NULL,doPlot=T,lExcl=NULL){
reroi=.MGrefineROIs(eic[which(!is.na(eic[,"y"])),],parDeco,minRTwin = parDeco$minRTwin,drt=parDeco$psdrt)
if(is.null(reroi)) return(NULL)
l2k=which(reroi[,"intensity"]>=parDeco$minHeightMS1 & (reroi[,"rtmax"]-reroi[,"rtmin"])>parDeco$minRTwin)
if(length(l2k)==0) return(NULL)
reroi=reroi[l2k,,drop=F]
nspan=floor(parDeco$span/2)*2+1
minNoise=parDeco$minNoiseMS1
llmatpks=list()
for(iiroi in 1:nrow(reroi)){
#########
iroi=reroi[iiroi,]
mzr=range(iroi[c("mzmin","mzmax")])
rtr0=rtr=range(iroi[c("rtmin","rtmax")])
rtr=rtr+c(-3,3)*parDeco$psdrt*(parDeco$span+1) ## extra padding
xeic=eic[which(eic[,"rt"]>=rtr[1] & eic[,"rt"]<=rtr[2]),]
l2excl=which(xeic[,"mz"]<mzr[1] | xeic[,"mz"]>mzr[2])
if(length(l2excl)) xeic[l2excl,c("mz","y")]=NA
xeic=xeic[order(xeic[,"scan"],!is.na(xeic[,"mz"])),]
xeic=xeic[which(!(duplicated(xeic[,"scan"]) & is.na(xeic[,"mz"]))),]
irmz=range(xeic[,"mz"],na.rm = T)
##### dual smooth/peak detec
# drt = parDeco$psdrt ; span=parDeco$span;bw=parDeco$bw;minNoise=parDeco$minNoiseMS1;minHeight=parDeco$minHeightMS1;sbslr=2
newx=.MGfill(xeic,parDeco$psdrt,nspan,parDeco$bw,minNoise)
# y=newx[,"y2"]
# n2pad=floor(parDeco$drt/parDeco$psdrt/2)
# y=c(rep(y[1],n2pad),y,rep(rev(y)[1],n2pad))
# bsl=GRMeta:::.GRbslrf(1:length(y),y,NoXP = NULL)
# bsl$fit[bsl$fit<minNoise]=minNoise
# bslscore <- (y - bsl$fit)/max(bsl$sigma, 10^-3)
# bslscore[which(abs(bslscore) > 10)] = sign(bslscore[which(abs(bslscore) > 10)]) * 10
# #newx$bslc=bslscore[nspan*5+1+(1:length(newx$y))]
# newx=data.frame(cbind(newx,bsl=bsl$fit[n2pad+(1:nrow(newx))],bslc=bslscore[n2pad+(1:nrow(newx))]))
# noise.local=(newx$y2/newx$bsl >=2 & newx$bslc>=1)*1;noise.local[newx$bslc< -1]=-1
#
# bsl=newx$bsl;bslscore=newx$bslc;x=newx[,"rt"];y=newx[,"y2"];span=nspan;minNoise = minNoise*1.01;v2alim = 0.8;snr.thresh = parDeco$sbslr
#### quick integration
pks=.MGsimpleIntegr2(x=newx[,"rt"],y=newx[,"y2"],bsl = newx$bsl,bslscore = newx$bslc,snr.thresh = parDeco$sbslr,span=nspan,minNoise = minNoise*1.01,v2alim = 0.8)
# pks=.MGsimpleIntegr(x=newx[,"rt"],y=newx[,"y2"],noise.local =newx$bslc,snr.thresh = 1,span=nspan,minNoise = minNoise*1.01,v2alim = 0.8)
if(nrow(pks)==0) next
##### check blank from and align to newx$x
newx$ybl=rep(parDeco$minNoiseMS1,nrow(newx))
if(!is.null(eicbl)){
tmpbl=eicbl
l2excl=which(tmpbl[,"mz"]<mzr[1] | tmpbl[,"mz"]>mzr[2])
if(length(l2excl)) tmpbl[l2excl,c("mz","y")]=NA
tmpbl=tmpbl[order(tmpbl[,"scan"],!is.na(tmpbl[,"mz"])),]
tmpbl=tmpbl[which(!(duplicated(tmpbl[,"scan"]) & is.na(tmpbl[,"mz"]))),]
tmpbl=tmpbl[tmpbl[,"blid"]%in%names(which(tapply(!is.na(tmpbl[,"y"]),tmpbl[,"blid"],sum)>10)),,drop=F]
if(nrow(tmpbl)>1){
tmpbl[which(is.na(tmpbl[,"y"]) | tmpbl[,"y"]<minNoise),"y"]=minNoise
tmpbl<-try(do.call("cbind",tapply(1:nrow(tmpbl),tmpbl[,"blid"],function(xbl){
ybl = .GRasysm(tmpbl[xbl,"y"], p = 0.5, lambda = 10^5)
ybl[ybl <= minNoise | is.na(ybl)] = minNoise
ybl=approx(tmpbl[xbl,"rt"],ybl,newx$rt,yleft = ybl[1],yright = rev(ybl)[1])$y})),TRUE)
if(!"try-error"%in%class(tmpbl)) newx$ybl=apply(tmpbl,1,median)
}
}
## compute the SBR around the apex
pks$pk.sbr=round(sapply(pks$pk.loc,function(i) max((newx$y2/newx$ybl)[which(abs(newx$rt-i)<=(2*parDeco$psdrt))])),3)
pks$pk.snr=round(sapply(pks$pk.loc,function(i) max((newx$y2/newx$bsl)[which.min(abs(newx$rt-i))])),3)
sbrlim=ifelse(any(lExcl>=irmz[1] & lExcl<=irmz[2]),0,parDeco$sbr)
l2k=as.numeric(names(which(tapply(pks$pk.sbr,pks$pk.cl,max)>=sbrlim &
tapply(pks$pk.snr,pks$pk.cl,max)>=parDeco$sbslr &
tapply(pks$pk.int*2,pks$pk.cl,max)>parDeco$minHeightMS1)))
pks=pks[pks$pk.cl%in%l2k,,drop=F]
if(nrow(pks)==0) next
pks$pk.cl=as.numeric(factor(pks$pk.cl))
## compute coda/roi stats
pks$pk.roi=rep(1,nrow(pks))
if(nrow(pks)>1)
for(i in 2:(nrow(pks))) pks$pk.roi[i]=pks$pk.roi[i-1]+ifelse((pks$pk.left[i]-pks$pk.right[i-1])>=parDeco$drt,1,0)
pks=cbind(pks,do.call("rbind",tapply(1:nrow(pks),pks$pk.roi,function(x){
nrtr0=range(c(pks$pk.left[x],pks$pk.right[x]))+c(-.49,.49)*parDeco$drt
nrtr0[1]=max(nrtr0[1],min(rtr0))
nrtr0[2]=min(nrtr0[2],max(rtr0))
l1=which(newx$rt>=nrtr0[1] & newx$rt<=nrtr0[2])
nrtr=range(c(pks$pk.left[x],pks$pk.right[x]))
l2=which(newx$rt>=nrtr[1] & newx$rt<=nrtr[2])
data.frame(roi.id=sprintf("R%.4f@%.2f-%.2f",median(newx[l2,"mz"],na.rm=T),min(newx[l2,"rt"],na.rm=T),max(newx[l2,"rt"],na.rm=T)),
roi.coda=.GRcoda(newx$y2[l1]),roi.mz50=unname(median(newx[l2,"mz"],na.rm=T)),
roi.mz10=unname(quantile(newx[l2,"mz"],.1,na.rm=T)),roi.mz90=unname(quantile(newx[l2,"mz"],.9,na.rm=T)))
}))[as.character(pks$pk.roi),,drop=FALSE])
## integrate peaks
newx$y2v=newx$y2-newx$bsl
newx$y2v[newx$y2v<0]=0
idx=1
matpks=list()
for(idx in 1:nrow(pks)){
tmp=newx[which(newx$rt>=pks$pk.left[idx] & newx$rt<=pks$pk.right[idx]),]
lnna=which(!is.na(tmp$y))
if(length(lnna)==0) next
lnna=lnna[order(abs(lnna-which.max(tmp$y2)),-tmp$y[lnna])]
if(length(lnna)>nspan) lnna=lnna[1:nspan]
iap=lnna[which.max(tmp$y[lnna])]
matpks[[idx]]=data.frame(id=idx,
pk.id=sprintf("P%.5f@%.3f",tmp$mz[iap],tmp$rt[iap]),
pk.rtmin=min(tmp$rt),
pk.rtmax=max(tmp$rt),
pk.rt=tmp$rt[which.max(tmp$y2)],
pk.rtap=tmp$rt[iap],
pk.int.ap=tmp$y[iap],
pk.int.apsm=max(tmp$y2),
pk.intmin=tmp$y2[1],
pk.bslmin=tmp$bsl[1],
pk.intmax=rev(tmp$y2)[1],
pk.bslmax=rev(tmp$bsl)[1],
pk.npts=sum(!is.na(tmp$y)),
pk.mz=round(tmp$mz[iap],6),
pk.mzmin=round(min(tmp$mz[iap],na.rm = T),6),
pk.mzmax=round(max(tmp$mz[iap],na.rm = T),6),
pk.mzmed=suppressWarnings(round(matrixStats:::weightedMedian(tmp$mz,tmp$y,na.rm=T),6)),
pk.area.sm=round(unname(.GRgetArea(tmp$rt,tmp$y2)[1]),3),
pk.area=round(.GRgetArea(tmp$rt[!is.na(tmp$y)],tmp$y[!is.na(tmp$y)])[1],3),
pk.snr=tmp$y[iap]/tmp$bsl[iap])
}
if(length(matpks)==0) next
matpks=do.call("rbind",matpks)
matpks=cbind(matpks,pk.sbr=pks[matpks$id,c("pk.sbr")],pks[matpks$id,grep("^roi.",names(pks))])
matpks$id=paste0("R",pks$pk.roi,"C",pks$pk.cl)[matpks$id]
matpks$cl.id=tapply(1:nrow(matpks),matpks$id,function(x) sprintf("C%.5f@%.3f",median(matpks$pk.mz[x]),median(matpks$pk.rtap[x])))[matpks$id]
llmatpks[[iiroi]]=matpks
###
if(doPlot){
# pks=ire[[1]]
# newx=ire[[2]]
yl=pretty(range(sqrt(c(0,xeic[,"y"],newx$y)),na.rm=T))
# yl=pretty(sqrt(range(pretty(yl^2))))
xl=pretty(rtr0,n=11)
par(mfrow=c(2,1))
par(mar=c(1,5,2,0))
plot(xeic[,"rt"],sqrt(xeic[,"y"]),pch=16,main=paste(imain,iiroi,round(.GRcoda(newx$y2),2)),xlab="",ylab="Height*10^6",
xlim=range(xl),ylim=range(yl),axes=F)
axis(2,at=yl,yl^2/1000000,las=2,pos=min(xl))
# axis(1,at=xl,pos=min(yl))
# abline(v=c(unique(matpks$pk.rtmax),unique(matpks$pk.rtmin)),col="grey",lwd=2)
lines(newx$rt,sqrt(newx$y2),col=2)
lines(newx$rt,sqrt(newx$bsl),col=4,lwd=2,lty=1)
lines(newx$rt,sqrt(newx$bsl*parDeco$sbslr),col=4,lwd=1,lty=4)
lines(newx$rt,sqrt(newx$ybl*parDeco$sbr),col=3,lwd=1,lty=2)
lines(newx$rt,sqrt(newx$ybl),col=3,lwd=2,lty=1)
abline(h=sqrt(parDeco$minHeightMS1))
points(matpks$pk.rtap,sqrt(matpks$pk.int.ap),pch=8,col=factor(matpks$cl.id),cex=2)
abline(v=matpks$pk.rtmin+parDeco$psdrt/2,lty=3)
abline(v=matpks$pk.rtmax-parDeco$psdrt/2,lty=4)
################
par(mar=c(4,5,0,0))
yl=pretty(range(c(range(matpks$pk.mzmed,matpks$pk.mz)*(1+c(-1,1)*2.1*10^-6),quantile(xeic[,"mz"],c(.1,.9),na.rm=T))))
xl=pretty(rtr0,n=11)
plot(xeic[,"rt"],xeic[,"mz"],pch=16,main=NULL,xlab="Retention time",ylab="",xlim=range(xl),ylim=range(yl),axes=F)
axis(2,at=yl,las=2,pos=min(xl))
nv=ceiling(diff(range(yl))/(median(yl)*10^-6))
abline(h=(-nv:nv)*median(yl)*10^-6+median(yl),lty=3,col="grey")
axis(1,at=xl,pos=min(yl))
# abline(v=c(unique(matpks$pk.rtmax),unique(matpks$pk.rtmin)),col="grey",lwd=2)
points(matpks$pk.rtap,matpks$pk.mz,pch=19,col=factor(matpks$cl.id),cex=2)
abline(v=matpks$pk.rtmin+parDeco$psdrt/2,lty=3)
abline(v=matpks$pk.rtmax-parDeco$psdrt/2,lty=4)
par(mfrow=c(1,1),mar=c(5.1 ,4.1 ,4.1, 2.1))
}
# reroi2=c(reroi2,list(pks))
}
allmatpks=do.call("rbind",llmatpks)
return(allmatpks)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.