Nothing
peaksDataset <- function(fns=dir(, "[Cc][Dd][Ff]"), verbose=TRUE,
mz=seq(50, 550), rtDivide=60, rtrange=NULL){
rawdata <- vector("list", length(fns))
rawrt <- vector("list", length(fns))
for(i in 1:length(fns)){
if(verbose)
cat(" Reading ", fns[i], "\n")
a <- xcmsRaw(fns[i])
D1 <- cbind.data.frame(mz=a@mzrange[1]:a@mzrange[2],
a@env$profile)
D2 <- data.frame(mz=mz)
if(a@mzrange[1]-range(mz)[1] < 3 & a@mzrange[2]-range(mz)[2] <
3){
mrg <- merge(D1, D2, by='mz', all=TRUE)
}else{
mrg <- merge(D1, D2, by='mz', all=FALSE)
}# not a very clean approach. The problem is in the raw m/z;
# the possible cause is the centroid approx. diring the
# acquisition of the cgromatogram se all=TRUE funziona con
# mzrange= 50-550; se all=FALSE funziona con mzrange= 100-400
mrg[is.na(mrg)] <- 0
rownames(mrg) <- 1:nrow(mrg)
colnames(mrg) <- 1:ncol(mrg)
a@env$profile <- data.matrix(mrg)
a@mzrange <- range(mz) ## allows mzrange param without errors
if(is.null(mz))
mz <- seq(a@mzrange[1], a@mzrange[2])
this.mz <- seq(a@mzrange[1], a@mzrange[2])
rawrt[[i]] <- a@scantime/rtDivide
nc <- length(rawrt[[i]])
if(length(rtrange) == 2){
w <- which(rawrt[[i]] >= rtrange[1] & rawrt[[i]] <= rtrange[2])
rawrt[[i]] <- rawrt[[i]][w]
}else{
w <- seq(1, nc)
}
rawdata[[i]] <- a@env$profile[, w]
rm(a)
if(length(this.mz) != length(mz)){
d <- matrix(0, nrow=length(mz), ncol=length(w))
d[match(this.mz, mz), ] <- rawdata[[i]]
rawdata[[i]] <- d
}
}
gc()
nm <- lapply(fns, function(u){
sp <- strsplit(u, split="/")[[1]]
sp[length(sp)]
})
nm <- sub(".CDF", "", nm)
names(rawdata) <- names(rawrt) <- nm
new("peaksDataset", rawdata=rawdata,
rawrt=rawrt, mz=mz, files=fns)
}
.plotpD<-function(object,runs=1:length(object@rawdata),mzind=1:nrow(object@rawdata[[1]]),mind=NULL,plotSampleLabels=TRUE,calcGlobalMax=FALSE,peakCex=0.8,plotPeaks=TRUE,plotPeakBoundaries=FALSE,plotPeakLabels=FALSE,plotMergedPeakLabels=TRUE,mlwd=3,usePeaks=TRUE,plotAcrossRuns=FALSE,overlap=F,rtrange=NULL,cols=NULL,thin=1,max.near=median(object@rawrt[[1]]),how.near=50,scale.up=1,...) {
# only specific one of ma or runs
if(is.null(rtrange))
rtrange<-c(min(object@rawrt[[1]]),max(object@rawrt[[1]]))
if(is.factor(cols))
cols<-c("black","green","blue","red","grey","orange")[cols]
if(is.null(cols))
cols<-rep(c("black","green","blue","red"),each=length(object@rawdata)/4)
if(length(cols)!=length(runs))
cols<-rep("black",length(runs))
if (usePeaks) {
datalist<-object@peaksdata
rtlist<-object@peaksrt
} else {
datalist<-object@rawdata
rtlist<-object@rawrt
}
if(calcGlobalMax) {
mx<-0
for(i in 1:length(runs)) {
ind<-runs[i]
mxind<-which(abs(object@rawrt[[ind]]-max.near)<how.near)
b<-colSums(matrix(object@rawdata[[ind]][mzind,],nrow=length(mzind)))
if (mx < max(b[mxind]))
mx<-max(b[mxind])
}
}
if( length(object@peaksdata)!=length(object@rawdata) )
plotPeaks<-FALSE
for(i in 1:length(runs)) {
ind<-runs[i]
b<-colSums(matrix(object@rawdata[[ind]][mzind,],nrow=length(mzind)))
if (!calcGlobalMax) {
mxind<-which(abs(object@rawrt[[ind]]-max.near)<how.near)
mx<-max(b[mxind])
}
b<-b/mx*scale.up
if (!overlap) {
if(i==1) {
plot(object@rawrt[[ind]],b,type="l",xlab="Retention Time",ylab="",xlim=rtrange,ylim=c(-.5,length(runs)),col=cols[i],xaxs="i",axes=FALSE,...)
axis(1)
} else {
lines(object@rawrt[[ind]],b+i-1,col=cols[i],...)
}
if (plotSampleLabels)
text(rtrange[1]+.1*diff(rtrange),i-.5,names(object@rawdata)[ind],col="black")
} else {
if(i==1) {
plot(object@rawrt[[ind]],b,type="l",xlab="Retention Time",ylab="",xlim=rtrange,ylim=c(0,1),col=cols[i],xaxs="i",axes=FALSE)
axis(1)
} else {
lines(object@rawrt[[ind]],b,col=cols[i])
}
}
if(plotPeaks & !overlap) {
rt<-rtlist[[ind]]; d<-min(diff(rt))*.4
for(j in 1:length(rt)) {
if(rt[j] > rtrange[1] & rt[j] < rtrange[2]) {
if (plotPeakBoundaries) {
st<-object@rawrt[[ind]][ object@peaksind.start[[ind]][j] ]
en<-object@rawrt[[ind]][ object@peaksind.end[[ind]][j] ]
#cat(ind,st,en,rt[j],"\n")
if(j%%2==0)
lines(c(st,st,rt[j],rt[j],rt[j],en,en),c(i-1+.025,i-1,i-1,i-1+.3,i-1,i-1,i-1+.025),col="blue")
else
lines(c(st,st,rt[j],rt[j],rt[j],en,en),c(i-1+.025,i-1,i-1,i-1+.3,i-1,i-1,i-1+.025)-.025,col="red")
#lines(c(rt[j],rt[j]),c(i-1,i-1+.3))
} else {
lines(c(rt[j]-d,rt[j]+d),c(i-1,i-1))
lines(c(rt[j],rt[j]),c(i-1,i-1+.3))
}
}
}
if(plotPeakLabels)
text(rt,i-1+.15,labels=1:length(rt),adj=0,cex=peakCex)
}
}
if(!is.null(mind) & !overlap) {
rws<-as.integer(seq(1,nrow(mind),by=thin))
for(i in rws) {
if( !is.na(mind[i,1]) & plotMergedPeakLabels )
text(rtlist[[runs[1]]][mind[i,1]],0,srt=90,i,col="blue",pos=1,cex=peakCex*.8,adj=1)
for(j in 1:(ncol(mind)-1)) {
ind1<-mind[i,j]
ind2<-mind[i,j+1]
if(!is.na(ind1) & !is.na(ind2)) {
rt1<-rtlist[[runs[j]]][ind1]
rt2<-rtlist[[runs[j+1]]][ind2]
#cat(ind1,ind2,runs[j],runs[j+1],rt1,rtrange[1],rt2,rtrange[2],"\n")
if( (rt1 > rtrange[1] | rt2 > rtrange[1]) & (rt1 < rtrange[2] | rt2 < rtrange[2])) {
if(j > 1)
lines(c(rt1,rt2),c(j+.4-1,j+.95-1),lwd=mlwd,lty=1,col="darkgrey")
else
lines(c(rt1,rt2),c(j+.4-1,j+.95-1),lwd=mlwd,lty=1,col="darkgrey")
}
}
}
if (plotAcrossRuns) {
if (ncol(mind) <= 2) next
for(j in 1:(ncol(mind)-2)) {
ind1<-mind[i,j]
ind2<-mind[i,j+1]
ind3<-mind[i,j+2]
if(!is.na(ind1) & is.na(ind2) & !is.na(ind3)) {
rt1<-rtlist[[runs[j]]][ind1]
rt3<-rtlist[[runs[j+2]]][ind3]
#cat(rt1,rtrange[1],rt3,rtrange[2],"\n")
if( (rt1 > rtrange[1] | rt3 > rtrange[1]) & (rt1 < rtrange[2] | rt3 < rtrange[2]))
lines(c(rt1,rt3),c(j+.4-1,j+.95),lwd=mlwd,lty=1,col="green")
}
}
if (ncol(mind) <= 3) next
for(j in 1:(ncol(mind)-3)) {
ind1<-mind[i,j]
ind2<-mind[i,j+1]
ind3<-mind[i,j+2]
ind4<-mind[i,j+3]
if(!is.na(ind1) & is.na(ind2) & is.na(ind3) & !is.na(ind4)) {
rt1<-rtlist[[runs[j]]][ind1]
rt4<-rtlist[[runs[j+3]]][ind4]
#cat(rt1,rtrange[1],rt4,rtrange[2],"\n")
if( (rt1 > rtrange[1] | rt4 > rtrange[1]) & (rt1 < rtrange[2] | rt4 < rtrange[2]))
lines(c(rt1,rt4),c(j+.4-1,j+1.95),lwd=mlwd,lty=1,col="red")
}
}
}
}
}
}
#setGeneric("plot", function(object, ...) standardGeneric("plot"))
#setGeneric("plot",function(x,y,...){standardGeneric("plot")})
#setMethod("plot",signature("peaksDataset","missing"),function(x,y,...) plotpD(x,...))
setMethod("plot","peaksDataset",function(x,y,...) .plotpD(x,...))
setMethod("show","peaksDataset",
function(object) {
cat("An object of class \"",class(object),"\"\n",sep="")
cat(length(object@rawdata), "samples:",names(object@rawdata),"\n",sep=" ")
cat(length(object@mz), "m/z bins - range: (",range(object@mz),")\n",sep=" ")
cat("scans:",sapply(object@rawdata,ncol),"\n",sep=" ")
cat("peaks:",sapply(object@peaksdata,ncol),"\n",sep=" ")
})
addAMDISPeaks<-function(object,fns=dir(,"[Eu][Ll][Uu]"),verbose=TRUE,...) {
if(length(fns)!=length(object@rawdata))
stop("Number of files must be the same as the number of runs (and must match).")
mn<-1e6; mx<--1
for(i in 1:length(object@rawrt)) {
rng<-range(object@rawrt[[i]])
if(mn > rng[1]) mn<-rng[1]
if(mx < rng[2]) mx<-rng[2]
}
if (verbose)
cat("Reading retention time range:",c(mn,mx),"\n")
for(i in 1:length(object@files)) {
if (verbose)
cat("Reading",fns[i],"...")
csin<-parseELU(fns[i],mz=object@mz,rtrange=c(mn,mx),...)
thisind<-rep(NA,nrow(csin$tab))
thisind.st<-rep(NA,nrow(csin$tab))
thisind.en<-rep(NA,nrow(csin$tab))
thisrt<-rep(NA,nrow(csin$tab))
thisdata<-matrix(0,nrow=length(object@mz),ncol=nrow(csin$tab))
for(j in 1:length(thisind)) {
thisind[j]<-which.min(abs(object@rawrt[[i]]-csin$tab[j,2]))
thisind.st[j]<-csin$tab[j,3]
thisind.en[j]<-csin$tab[j,4]
thisrt[j]<-object@rawrt[[i]][thisind[j]]
take<-which(csin$peaks[,j]>0)
#thisdata[take,j]<-sqrt(object@rawdata[[i]][take,thisind[j]])
thisdata[take,j]<-object@rawdata[[i]][take,thisind[j]]
}
if (verbose)
cat(" Done.\n")
object@peaksind[[i]]<-thisind
thisind.st[thisind.st==0]<-1
object@peaksind.start[[i]]<-thisind.st
object@peaksind.end[[i]]<-thisind.en
object@peaksrt[[i]]<-thisrt
object@peaksdata[[i]]<-thisdata
}
names(object@peaksdata)<-names(object@peaksrt)<-names(object@peaksind)<-names(object@rawdata)
new("peaksDataset",object)
}
addChromaTOFPeaks<-function(object,fns=dir(,"[Tt][Xx][Tx]"),rtDivide=60,verbose=TRUE,...) {
if(length(fns)!=length(object@rawdata))
stop("Number of files must be the same as the number of runs (and must match).")
mn<-1e6; mx<--1
for(i in 1:length(object@rawrt)) {
rng<-range(object@rawrt[[i]])
if(mn > rng[1]) mn<-rng[1]
if(mx < rng[2]) mx<-rng[2]
}
if (verbose)
cat("Reading retention time range:",c(mn,mx),"\n")
for(i in 1:length(object@files)) {
if (verbose)
cat("Reading",fns[i],"...")
csin<-parseChromaTOF(fns[i],mz=object@mz,rtrange=c(mn,mx),rtDivide=rtDivide,...)
thisrt<-csin$tab$rt
thisind<-rep(NA,length(thisrt))
thisdata<-matrix(0,nrow=length(object@mz),ncol=nrow(csin$tab))
for(j in 1:length(thisind)) {
thisind[j]<-which.min(abs(object@rawrt[[i]]-csin$tab$rt[j]/rtDivide))
thisrt[j]<-object@rawrt[[i]][thisind[j]]
take<-which(csin$peaks[,j]>0)
#thisdata[take,j]<-sqrt(object@rawdata[[i]][take,thisind[j]])
thisdata[take,j]<-object@rawdata[[i]][take,thisind[j]]
}
if (verbose)
cat(" Done.\n")
object@peaksind[[i]]<-thisind
object@peaksrt[[i]]<-thisrt
object@peaksdata[[i]]<-thisdata
}
names(object@peaksdata)<-names(object@peaksrt)<-names(object@peaksind)<-names(object@rawdata)
new("peaksDataset",object)
}
setMethod("plotImage","peaksDataset",
function(object,run=1,rtrange=c(11,13),main=NULL,mzrange=c(50,200),SCALE=log2,...) {
mz<-object@mz
for(i in run) {
rt<-object@rawrt[[i]]
d<-object@rawdata[[i]]
cols<-colorpanel(50,"green","blue","red")
rind<-which(rt>rtrange[1]&rt<rtrange[2])
mind<-which(mz>mzrange[1]&mz<mzrange[2])
ff<-t(d[mind,rind])
ff[ff==0]<-min(ff[ff>0])
if( is.null(main) )
tit<-object@files[i]
else
tit<-main
#cat("min=",min(SCALE(ff),na.rm=TRUE)," max=",max(SCALE(ff),na.rm=T),"\n")
image(rt[rind],mz[mind],SCALE(ff),col=cols,xlab="retention time",ylab="m/z",xlim=rtrange,main=tit,...)
#abline(h=c(72.5,73.5,146.5,147.5),col="white",lty=3)
#if (!is.null(pts)) {
# points(pts[,1],pts[,2],pch=19)
# matplot(t(pts[,8:9]),t(pts[,c(2,2)]),type="l",col="black",lty=1,add=TRUE)
#}
}
return(invisible(ff))
})
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.