.GRcompAlignOneEIC<-function(ieic,lrefs,tabeic,eicfile=NULL,whichrt="rtcor",
eicParams){
require(mgcv)
if(is.null(eicfile)) eicfile=paste(eicParams$dirEic,tabeic$GrpEic[which(tabeic$Id==ieic)],
eicParams$addEic,".rda",sep="")
load(eicfile)
tmpeic=dfeic[dfeic$eic==ieic,]
rtref=codaref=rep(NA,length(lrefs));names(rtref)=names(codaref)=lrefs
cmax=drt=matrix(NA,nrow=length(lrefs),ncol=ncol(m),dimnames=list(lrefs,colnames(m)))
if(any(!lrefs%in%dfeic$samp)) return(list(cmax=cmax,drt=drt,rtref=rtref,codaref=codaref))
meic = .GRconvEIC(tmpeic, whichrt = whichrt, bw = eicParams$nsmo1 * eicParams$bw, delrt = eicParams$bw/2)
m=meic[[1]]
xrt=meic[[2]]
lenout=2^ceiling(log2(nrow(m)*2))
for(iref in lrefs[lrefs%in%colnames(m)]){
rec=apply(m,2,.GRcompVals,m[,iref],lenout,21)
cmax[iref,]=apply(rec,2,max)
drt[iref,]=(which.max(rec[,iref])-apply(rec,2,which.max))*median(diff(rt))
rtref[iref]=weighted.median(xrt,m[,iref])
codaref[iref]= .GRcodadw2(tmpeic$y[tmpeic$samp==iref])
}
return(list(cmax=cmax,drt=drt,rtref=rtref,codaref=codaref))
}
################
.GRcompAlignGrpEIC<-function(grpeic,lrefs,drtmax=0.3,tabeic=NULL,eicfile=NULL,whichrt="rtcor",eicParams,verbose=FALSE){
require(limma)
if(is.null(eicfile)) eicfile=paste(eicParams$dirEic,grpeic,
eicParams$addEic,".rda",sep="")
if(verbose) cat("* ",eicfile,": ")
load(eicfile)
if(!is.null(tabeic)) leics=tabeic$Id[tabeic==grpeic] else leics=unique(dfeic$eic)
leics=leics[leics%in%dfeic$eic]
ares=list()
for(ieic in leics){
if(verbose) cat(ieic," ",sep=" ")
tmpeic=dfeic[dfeic$eic==ieic,]
if("ineic" %in%names(tmpeic)) tmpeic=tmpeic[tmpeic[,"ineic"]==1,]
meic = .GRconvEIC(tmpeic, whichrt = whichrt, bw = eicParams$nsmo1 * eicParams$bw, delrt = eicParams$bw/2)
m=meic[[1]]
xrt=meic[[2]]
rtref=codaref=rep(NA,length(lrefs));names(rtref)=names(codaref)=lrefs
cmax=drt=matrix(NA,nrow=length(lrefs),ncol=ncol(m),dimnames=list(lrefs,colnames(m)))
nmax=ceiling(drtmax/median(diff(xrt))+3)
if(nmax%%2==0) nmax=nmax+1
lenout=2^ceiling(log2(nrow(m)*2))
for(iref in lrefs[lrefs%in%colnames(m)]){
rec=apply(m,2,.GRcompVals,m[,iref],lenout,nmax)
cmax[iref,]=apply(rec,2,max)
drt[iref,]= (apply(rec,2,which.max)-which.max(rec[,iref]))*median(diff(xrt))
rtref[iref]=weighted.median(xrt,m[,iref])
codaref[iref]= .GRcodadw2(tmpeic$y[tmpeic$samp==iref])
}
ares[[ieic]]=list(cmax=cmax,drt=drt,rtref=rtref,codaref=codaref)
}
if(verbose) cat("\n")
return(ares)
}
################
.GRcorGrpEIC<-function(lgrpeic,shmat,tabeic=NULL,byGrp=FALSE,dosave=FALSE,
whichrt="rtcor",newrt="rtcor2",eicParams,verbose=FALSE,retres=TRUE,ncl=1){
.inGRcorGrpEIC<-function(grpeic,shmat,tabeic=NULL,byGrp=FALSE,dosave=FALSE,
whichrt="rtcor",newrt="rtcor2",eicParams,verbose=TRUE,retres=TRUE){
eicfile=paste(eicParams$dirEic,grpeic,eicParams$addEic,".rda",sep="")
if(verbose) cat("* ",eicfile)
load(eicfile)
if(!is.null(tabeic)) leics=tabeic$Id[tabeic$GrpEic==grpeic] else leics=unique(dfeic$eic)
leics=leics[leics%in%dfeic$eic]
eicinf=attr(dfeic,"eic")
eicinf=eicinf[eicinf$Id%in%leics,]
rts=as.numeric(rownames(shmat))
addrt=apply(shmat,2,function(x) approx(rts,x,eicinf$rtap)$y)
rownames(addrt)=eicinf$Id
addrtv=as.vector(addrt)
names(addrtv)=paste(rep(rownames(addrt),ncol(addrt)),rep(colnames(addrt),each=nrow(addrt)),sep=";;;")
addnewrt=addrtv[paste(dfeic$eic,dfeic$samp,sep=";;;")]
if(any(is.na(addnewrt))){
if(verbose) cat(" some missing!")
addnewrt[is.na(addnewrt)]=0
}
addnewrt=dfeic[,whichrt]+addnewrt
dfeic[,newrt]=addnewrt
if(dosave) save(file=eicfile,dfeic,eicinfos,sampinfos)
if(verbose) cat("\n")
if(retres) invisible(dfeic) else return(nrow(eicinf))
}
###################
if(ncl!=1){
require("snowfall")
ncl=max(1,min(ncl,parallel:::detectCores()))
}
if(!is.null(tabeic) & is.null(lgrpeic)) lgrpeic=unique(tabeic$GrpEic)
lfiles=paste(ifelse(is.null(eicParams$dirEic),"./",eicParams$dirEic),
lgrpeic,
ifelse(is.null(eicParams$addEic),"./",eicParams$addEic),".rda",sep="")
lgrpeic=unique(lgrpeic[file.exists(lfiles)])
if(!is.null(tabeic)) lgrpeic=unique(tabeic$GrpEic[tabeic$GrpEic%in%lgrpeic])
d0=proc.time()[3]
cat("Started at ",date()," on ",ncl," processors\n",sep="")
if(ncl==1){
allr=lapply(lgrpeic,.inGRcorGrpEIC,shmat=shmat,tabeic=tabeic,byGrp=byGrp,dosave=dosave,whichrt=whichrt,newrt=newrt,eicParams,
verbose=verbose,retres=retres)
if(retres) totop=sum(sapply(allr,function(x) length(unique(dfeic$eic)))) else totop=sum(unlist(allr))
}
if(ncl>1){
sfInit(parallel=TRUE, cpus=ncl, type='SOCK')
sfExport( ".inGRcorGrpEIC", local=TRUE )
sfLibrary(GRMeta)
sfLibrary(limma)
allr=sfClusterApplyLB(lgrpeic,.inGRcorGrpEIC,shmat=shmat,tabeic=tabeic,byGrp=byGrp,dosave=TRUE,
whichrt=whichrt,newrt=newrt,eicParams,verbose=FALSE,retres=FALSE)
totop=sum(unlist(allr))
sfStop()
}
d1=proc.time()[3]
cat("\nCompleted at ",date()," -> took ",round(d1-d0,1)," secs to process ",length(lgrpeic)," EIC groups\n",sep="")
cat(" * ",totop," EICs in all, ",round(totop/length(lgrpeic),3)," on av.\n")
return(invisible(allr))
}
################
.GRcorGrpEICo<-function(grpeic,shmat,tabeic=NULL,eicfile=NULL,byGrp=FALSE,dosave=FALSE,
whichrt="rtcor",newrt="rtcor2",eicParams,verbose=FALSE,ncl=1){
if(is.null(eicfile)) eicfile=paste(eicParams$dirEic,grpeic,
eicParams$addEic,".rda",sep="")
if(verbose) cat("* ",eicfile)
load(eicfile)
if(!is.null(tabeic)) leics=tabeic$Id[tabeic==grpeic] else leics=unique(dfeic$eic)
leics=leics[leics%in%dfeic$eic]
eicinf=attr(dfeic,"eic")
eicinf=eicinf[eicinf$Id%in%leics,]
rts=as.numeric(rownames(shmat))
addrt=apply(shmat,2,function(x) approx(rts,x,eicinf$rtap)$y)
rownames(addrt)=eicinf$Id
addrtv=as.vector(addrt)
names(addrtv)=paste(rep(rownames(addrt),ncol(addrt)),rep(colnames(addrt),each=nrow(addrt)),sep=";;;")
addnewrt=addrtv[paste(dfeic$eic,dfeic$samp,sep=";;;")]
if(any(is.na(addnewrt))){
if(verbose) cat(" some missing!")
addnewrt[is.na(addnewrt)]=0
}
addnewrt=dfeic[,whichrt]+addnewrt
dfeic[,newrt]=addnewrt
if(dosave) save(file=eicfile,dfeic,eicinfos,sampinfos)
if(verbose) cat("\n")
return(invisible(dfeic))
}
.GRcompAlignMGrpEIC<-function(tabeic,lrefs,drtmax=0.3,whichrt="rtcor",eicParams,ncl=1){
#
# if(is.null(names(lfiles))) names(lfiles)=1:length(lfiles)
# lfiles=lfiles[file.exists(lfiles)]
# if(length(lfiles)==0) return(NULL)
#
require(limma)
leicgrp=unique(tabeic$GrpEic)
leicgrp=leicgrp[file.exists(paste(eicParams$dirEic,leicgrp,eicParams$addEic,".rda",sep=""))]
print(leicgrp)
if(ncl!=1){
require("snowfall")
ncl=max(1,min(ncl,parallel:::detectCores()))
}
d0=proc.time()[3]
cat("Started at ",date(),sep="")
acef=list()
if(ncl==1){
cat(" on 1 processor\n",sep="")
acef=lapply(leicgrp,.GRcompAlignGrpEIC,
lrefs,tabeic,drtmax=drtmax,eicfile=NULL,whichrt=whichrt,eicParams,verbose=TRUE)
}
if(ncl>1){
cat(" on ",ncl," processors\n",sep="")
sfInit(parallel=TRUE, cpus=ncl, type='SOCK')
sfLibrary(GRMeta)
sfLibrary(limma)
#if(local)
# sfExport( ".GRcompAlignGrpEIC", local=TRUE )
acef=sfClusterApplyLB(leicgrp,.GRcompAlignGrpEIC,
lrefs,tabeic,drtmax=drtmax,eicfile=NULL,whichrt=whichrt,eicParams)
sfStop()
}
acef=unlist(acef,recursive=FALSE)
d1=proc.time()[3]
cat("\nCompleted at ",date()," -> took ",round(d1-d0,1)," secs - ",round((d1-d0)/length(leicgrp),1)," secs per file\n",sep="")
return(invisible(acef))
}
.GRcompAlignEICs<-function(ares,newrts=NULL,lsamps=NULL,lrefs=NULL,rtlim=NULL,drtmax=0.5){
if(is.null(lsamps)) lsamps=unique(unlist(lapply(ares,function(x) colnames(x[[1]]))))
if(is.null(lrefs)) lrefs=unique(unlist(lapply(ares,function(x) rownames(x[[1]]))))
adrt=lapply(ares,function(x){
x=x$drt[match(lrefs,rownames(x$drt)),match(lsamps,colnames(x$drt))]
dimnames(x)=list(lrefs,lsamps);x
})
acmax=lapply(ares,function(x){
x=x$cmax[match(lrefs,rownames(x$cmax)),match(lsamps,colnames(x$cmax))]
dimnames(x)=list(lrefs,lsamps);x
})
artref=sapply(ares,function(x) x$rtref[lrefs])
acodaref=sapply(ares,function(x) x$codaref[lrefs])
if(is.null(newrts)){
toadd=max(2,ceiling(drtmax/eicParams$bw))
newrts=(floor(min(artref,na.rm=T)/(2*eicParams$bw))-toadd):(ceiling(max(artref,na.rm=T)/(2*eicParams$bw))+toadd)
newrts=newrts*eicParams$bw*2
}
doborderu=doborderb=FALSE
if(!is.null(rtlim) & any(newrts<rtlim[1])){
lbe=which(newrts<rtlim[1])
lbe2=which(newrts>=rtlim[1])[1:3]
doborderb=TRUE
}
if(!is.null(rtlim) & any(newrts>rtlim[2])){
lup=which(newrts>rtlim[2])
lup2=rev(which(newrts<=rtlim[2]))[1:3]
doborderu=TRUE
}
###################
# for each sample
allmpr=list()
for(isamp in lsamps){
cat(".")
allmpr[[isamp]]=matrix(NA,nrow=length(newrts),ncol=length(lrefs),dimnames=list(newrts,lrefs))
mrt=sapply(adrt,function(x) x[,isamp])
cmax=sapply(acmax,function(x) x[,isamp])
cmax[which(cmax<0.2)]=0
for(iref in lrefs){
df=na.omit(data.frame(y=mrt[iref,],x=artref[iref,],w=cmax[iref,],cw=acodaref[iref,]))
df$w[df$w<0.2]=0
if(sum(df$w>0)>11){
m<-try(mgcv:::gam(y~s(x),data=df,weights=w),TRUE)
if("try-error"%in%class(m)) m<-try(mgcv:::gam(y~(x),data=df,weights=w),TRUE)
pr=mgcv:::predict.gam(m,data.frame(x=newrts))
if(any(abs(pr)>drtmax)) pr[which(abs(pr)>drtmax)]=drtmax*sign(pr[which(abs(pr)>drtmax)])
if(doborderb) pr[lbe]=median(pr[lbe2])
if(doborderu) pr[lup]=median(pr[lup2])
allmpr[[isamp]][,iref]=pr
}
}
}
cat("\n")
shmat=sapply(allmpr,function(x){
pr=apply(x,1,median)
if(any(abs(pr)>drtmax)) pr[which(abs(pr)>drtmax)]=drtmax*sign(pr[which(abs(pr)>drtmax)])
if(doborderb) pr[lbe]=median(pr[lbe2])
if(doborderu) pr[lup]=median(pr[lup2])
pr=mgcv:::gam(pr~s(newrts))$fit
if(doborderb) pr[lbe]=pr[lbe2[which.min(lbe2)]]
if(doborderu) pr[lup]=pr[lup2[which.max(lup2)]]
pr
})
rownames(shmat)=newrts
return(list(ShMat=shmat,RT=newrts,AllPr=allmpr,Ref=lrefs,RTref=artref))
}
#
# iref
#
#
# mpr=matrix(0,nrow=length(rts1),ncol=length(lsamp),dimnames=list(rts1,lsamp))
# for(isamp in lsamp[lsamp!=iref]){
# df=na.omit(data.frame(y=msh[,isamp],x=refrt,w=mcor[,isamp]))
# df$w[df$w<0.2]=0
# if(sum(df$w>0)>11){
# m<-try(mgcv:::gam(y~s(x),data=df,weights=w),TRUE)
# if("try-error"%in%class(m)) m<-try(mgcv:::gam(y~(x),data=df,weights=w),TRUE)
# pr=mgcv:::predict.gam(m,data.frame(x=rts1))*delrt
# # pr[rts1<min(df$x,na.rm=T) | rts1>max(df$x,na.rm=T)]=NA
# mpr[,isamp]=pr
# }
# }
#
#
#
#
# leics=rownames(tousemat)
# lsamp=colnames(tousemat)
# mcor=msh=matrix(NA,ncol=length(lsamp),nrow=length(leics),dimnames=list(leics,lsamp))
# refrt=rep(NA,length(leics));names(refrt)=leics
# leics2=names(which(tousemat[,as.numeric(iref)]))
# d0=proc.time()[3]
#
# if(serial) cat("Reference ",iref," (",length(leics2),")",sep="")
# for(i in 1:length(leicf)){
# load(leicf[i])
# if(serial) cat(".")
# leic=as.character(unique(dfeic$Eic))
# leic=leic[leic%in%leics2]
# for(ix in leic){
# tt=dfeic[dfeic$Eic==ix,]
# vmint=filP$Mint[ix]
# if(is.na(vmint)) vmint=median(filP$Mint)
# re<-filtreScan(csc=as.vector(tt[,"y"]>vmint),idx=tt$Samp,alpha=filP$LSc,perc=filP$Perc)
# tt$isin=re
#
# lrt=rts1[which(rts1>=(min(tt[,whichrt])-delrt) & rts1<=(max(tt[,whichrt])+delrt))]
# if(max(lrt)<max(tt[,whichrt])) lrt=c(lrt,max(lrt)+(1:ceiling((max(tt[,whichrt])-max(lrt))/delrt))*delrt)
# if(min(lrt)>min(tt[,whichrt])) lrt=c(min(lrt)-(ceiling((min(lrt)-min(tt[,whichrt]))/delrt):1)*delrt,lrt)
# m=sapply(unique(tt$Samp), function(ik) approx(tt[,whichrt][tt$Samp==ik],tt$y[tt$Samp==ik],lrt)$y)
# dimnames(m)=list(lrt,unique(tt$Samp))
# m[is.na(m)]=0
#
# lenout=2^ceiling(log2(nrow(m)*2))
# rec=apply(m,2,compVals,m[,iref],lenout,maxsh)
# # matplot(rec,typ="l",main=paste(iref,ix),col=mfiles$cols[as.numeric(colnames(rec))],lwd=2,lty=1)
# mcor[ix,]=apply(rec,2,max)[as.character(lsamp)]
# msh[ix,]=apply(rec,2,which.max)[as.character(lsamp)]-which.max(rec[,iref])
# if(sum(tt$Samp==iref & tt$isin)==0) print(c(i,ix))
# refrt[ix]=weighted.median(tt[,whichrt][tt$Samp==iref & tt$isin],tt$y[tt$Samp==iref & tt$isin])
# }
# }
# d1=proc.time()[3]
# if(serial) cat(":",round(d1-d0,1),"secs.\n Compute GAMs")
# top=paste("Reference ",iref," (",length(leics2),"): ",round(d1-d0,1),"secs.",sep="")
# # save(file=paste("totot",iref,".rda",sep=""),refrt,msh,mcor,lsamp,mfiles,leicf,tousemat,maxsh,whichrt,filP,border,rts1,delrt)
#
# mpr=matrix(0,nrow=length(rts1),ncol=length(lsamp),dimnames=list(rts1,lsamp))
# for(isamp in lsamp[lsamp!=iref]){
# df=na.omit(data.frame(y=msh[,isamp],x=refrt,w=mcor[,isamp]))
# df$w[df$w<0.2]=0
# if(sum(df$w>0)>11){
# m<-try(mgcv:::gam(y~s(x),data=df,weights=w),TRUE)
# if("try-error"%in%class(m)) m<-try(mgcv:::gam(y~(x),data=df,weights=w),TRUE)
# pr=mgcv:::predict.gam(m,data.frame(x=rts1))*delrt
# # pr[rts1<min(df$x,na.rm=T) | rts1>max(df$x,na.rm=T)]=NA
# mpr[,isamp]=pr
# }
# }
# ### smooth border
# # save(file=paste("totot",iref,".rda",sep=""),refrt,mpr)
#
# lb1=which(rts1>min(refrt,na.rm=T) & rts1<=(min(refrt,na.rm=T)+border+2*delrt))
# l2rep=1:lb1[(floor(length(lb1)/2))]
# mpr[l2rep,]=matrix(rep(apply(mpr[lb1,],2,median),each=length(l2rep)),ncol=length(lsamp))
# lb1=which(rts1<max(refrt,na.rm=T) & rts1>=(max(refrt,na.rm=T)-border-2*delrt))
# l2rep=nrow(mpr):rev(lb1)[(floor(length(lb1)/2))]
# mpr[l2rep,]=matrix(rep(apply(mpr[lb1,],2,median),each=length(l2rep)),ncol=length(lsamp))
# if(serial) cat("\n")
#
# return(list(sht=msh,pr=mpr,cor=mcor,rts=rts1,ref=iref,top=top))
# }
#
#
# #######
# ## comp alignment based on EIC
# compAlignMEIC<-function(tabeic,lrefs,whichrt="rtcor",maxsh=50,params,border=0.2,ncl=1){
#
# require(mgcv)
#
# d0=proc.time()[3]
#
# if(missing(lrefs)) stop('Sid of the ref samples')
#
#
# if(ncl>1){
# if(ncl>parallel:::detectCores()) ncl=min(ncl,parallel:::detectCores()-1,length(lrefs))
# if(sfIsRunning()) sfStop()
# sfInit(parallel=TRUE, cpus=ncl, type="SOCK")
# sfExport("filtreScan","compVals", "compAlignEIC","compAlignEICCl", local=TRUE )
# allr=sfClusterApplyLB(1:length(lrefs),compAlignEICCl,lrefs,mfiles,leicf,tousemat,maxsh,whichrt,filP,border)
# names(allr)=lrefs
# sfStop()
# allr=allr[sapply(allr,length)>0]
# cat(sapply(allr,function(x) x$top),sep="\n")
# }
#
# if(ncl==1){
# allr=list()
# for(iref in lrefs){
# resp<-try(compAlignEIC(iref=iref,mfiles,leicf,tousemat,maxsh=maxsh,whichrt=whichrt,filP=filP,border=border,serial=TRUE),TRUE)
# if(class(resp)!="try-error") allr[[iref]]<-resp
# }
# }
# d1=proc.time()[3]
#
# cat("\n Took",round(d1-d0,1),"secs. to process ",length(allr)," refs\n")
#
# return(list(allres=allr,whichrt=whichrt,border=border))
# }
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.