#R
# $HeadURL: http://fgcz-svn.unizh.ch/repos/fgcz/testing/proteomics/R/protViz/R/PTM_MarkerFinder.R $
# $Id: PTM_MarkerFinder.R 6318 2014-04-02 13:02:51Z cpanse $
# $Date: 2014-04-02 15:02:51 +0200 (Wed, 02 Apr 2014) $
.PTM_MarkerFinder_writeMGF_Header <- function(mgfFilename){
FILE <- file(mgfFilename, "a")
writeLines(as.character(paste("# protViz packageVersion ", packageVersion('protViz'))), FILE)
writeLines(paste("# ",date(),sep=''), FILE)
close(FILE)
}
.PTM_MarkerFinder_writeMGF <- function(spec, mgfFilename, pattern=-1){
FILE <- file(mgfFilename, "a")
if (length(pattern)>0){
writeLines(as.character(paste("# PATTERN ", pattern, sep=' ')), FILE)
}
writeLines("BEGIN IONS", FILE)
writeLines(paste("TITLE=",as.character(spec$title),sep=''), FILE)
writeLines(paste("PEPMASS=",as.character(spec$pepmass),sep=''), FILE)
writeLines(paste("CHARGE=",as.character(spec$charge),'+',sep=''), FILE)
writeLines(paste("SCANS=",as.character(spec$scans),sep=''), FILE)
writeLines(paste("RTINSECONDS=",as.character(spec$rtinseconds),sep=''), FILE)
writeLines(as.character(paste(spec$mZ, spec$intensity, sep=' ')), FILE)
writeLines("END IONS", FILE)
writeLines("", FILE)
close(FILE)
}
.PTM_MarkerFinder_ <- function(peptideSequence, modification, modificationName){
result <- as.character()
n <- nchar(peptideSequence)
seq <- strsplit(peptideSequence,'')[[1]]
mod <- as.numeric((strsplit(modification,'')[[1]]))
for (i in 1:n){
result<-paste(result, seq[i], sep='')
if (mod[i+1] > 0){
result<-paste(result,'[', modificationName[mod[i+1]+1], ']', sep='')
}
}
return(result)
}
.PTM_MarkerFinder_no_peakplot<-function(data,
modification,
modificationName,
mZmarkerIons,
minNumberIons=2,
itol_ppm=10,
minMarkerIntensityRatio=5){
lapply(data, function(x){
idx <- findNN_(mZmarkerIons, x$mZ)
ppm.error <- 1e-06 * itol_ppm * x$mZ[idx]
b<-(abs(mZmarkerIons-x$mZ[idx]) < ppm.error)
sum.mZmarkerIons.intensity <- sum(x$intensity[idx[b]])
sum.intensity <- sum(x$intensity) - sum.mZmarkerIons.intensity
percent.mZmarkerIons <- round(100 * (sum.mZmarkerIons.intensity / (sum.mZmarkerIons.intensity + sum.intensity)),1)
if ((length(x$mZ[idx[b]]) >= minNumberIons) & percent.mZmarkerIons > minMarkerIntensityRatio ){
split.screen(c(2,1))
split.screen(c(1,4), 2)
screen(1)
plot(x$mZ, x$intensity, type='h', xlab='m/z',ylab='Intensity', xlim=c(0, max(x$mZ)), sub=paste(x$title))
points(x$mZ[idx[b]], x$intensity[idx[b]], pch=22, col='#6CC417AA', bg='#6CC417AA', cex=1)
legend("topright", paste(c('query', 'pepmass', 'charge'), c(x$id, x$pepmass, x$charge)), cex=0.75)
################################################################################
screen(6)
pie(x$intensity[idx[b]],
paste(round(100*x$intensity[idx[b]]/sum(x$intensity[idx[b]]),0),"% ",
round(x$mZ[idx[b]], 2), sep=""),
main=list("relative marker ions %", cex=0.8),
cex=0.75)
################################################################################
screen(5)
pie(c(sum.intensity, sum.mZmarkerIons.intensity),
c('Intensity', paste(percent.mZmarkerIons, '%',sep='')),
col=c('white', '#6CC417AA'),
main=list("% marker ion intensities", cex=0.8),
cex=0.75)
################################################################################
screen(3)
#plot(mZmarkerIons, (mZmarkerIons-data[[i]]$mZ[idx]),
plot(mZmarkerIons[b], 1e+06 * (mZmarkerIons[b] - x$mZ[idx[b]]) / x$mZ[idx[b]],
axes=TRUE, type='p', xlab='m/z', ylab='ppm error',log='')# ylim=c(-max(ppm.error),max(ppm.error)))
axis(3, mZmarkerIons[b], round(mZmarkerIons[b],1),las=2)
abline(h=0.0, col='grey')
################################################################################
####### P A T T E R N P L O T ##################################################
screen(4)
plot(x$mZ[idx[b]],
x$intensity[idx[b]],
col='#6CC417AA',
lwd=2,
type='h',
xlab='m/z',
ylab='Intensity',
axes=TRUE)
axis(3, x$mZ[idx[b]],round(x$mZ[idx[b]],1), las=2)
# readline("Press <Enter> to continue")
close.screen(all.screens = TRUE)
}
close.screen(all.screens = TRUE)
}
)
}
# mZmarker which ion you want to see
PTM_MarkerFinder <- function(data,
modification,
modificationName,
mZmarkerIons,
minNumberIons=2,
itol_ppm=10,
minMarkerIntensityRatio=5,
mgfFilename=-1,
PEAKPLOT=TRUE){
if(!PEAKPLOT){
.PTM_MarkerFinder_no_peakplot(data,
modification,
modificationName,
mZmarkerIons,
minNumberIons,
itol_ppm,
minMarkerIntensityRatio)
return (TRUE)
}
query.idx<-1:length(data)
query.to.scan<-as.integer(as.character(lapply(data, function(x){
if (length(x$scans) == 1){
return(x$scans)
}else{return (x$scans[1])}
})))
scan.to.query<-rep(-1, max(query.to.scan, na.rm=TRUE))
scan.to.query[query.to.scan[query.idx]]<-query.idx
if (mgfFilename != -1){
unlink(mgfFilename)
.PTM_MarkerFinder_writeMGF_Header(mgfFilename)
}
rr<-numeric()
for (i in 2:length(data)){
idx <- findNN_(mZmarkerIons, data[[i]]$mZ)
# ppm.itol.cutoff
ppm.error <- 1e-06 * itol_ppm * data[[i]]$mZ[idx]
b<-(abs(mZmarkerIons-data[[i]]$mZ[idx]) < ppm.error)
sum.mZmarkerIons.intensity <- sum(data[[i]]$intensity[idx[b]])
sum.intensity <- sum(data[[i]]$intensity) - sum.mZmarkerIons.intensity
percent.mZmarkerIons <- round(100 * (sum.mZmarkerIons.intensity / (sum.mZmarkerIons.intensity + sum.intensity)),1)
if ((length(data[[i]]$mZ[idx[b]]) >= minNumberIons)
& percent.mZmarkerIons > minMarkerIntensityRatio ){
if (mgfFilename != -1){
.PTM_MarkerFinder_writeMGF(data[[i]], mgfFilename,
pattern=paste(data[[i]]$mZ[idx[b]], data[[i]]$intensity[idx[b]]))
}
r <- cbind(scans=data[[i]]$scans,
mZ=data[[i]]$mZ[idx[b]],
markerIonMZ=mZmarkerIons[b],
markerIonIntensity=data[[i]]$intensity[idx[b]],
markerIonMzError=mZmarkerIons[b]-data[[i]]$mZ[idx[b]],
markerIonPpmError=1e+06 * (mZmarkerIons[b]-data[[i]]$mZ[idx[b]])/data[[i]]$mZ[idx[b]],
query=i,
pepmass=data[[i]]$pepmass,
peptideSequence=data[[i]]$peptideSequence,
modification=data[[i]]$modification
)
rr<-rbind(rr,r)
split.screen(c(3,1))
split.screen(c(1,5), 2)
screen(1)
####### P E A K P L O T ########################################################
if (!is.na(data[[i]]$mascotScore) ){
fi<-fragmentIon(sequence=data[[i]]$peptideSequence,
FUN=defaultIon,
modified=substr(data[[i]]$modification, 2, nchar(data[[i]]$modification)-1),
modification=modification)
fi.by<-as.data.frame(cbind(b=fi[[1]]$b, y=fi[[1]]$y), stringsAsFactors = TRUE)
peakplot(data[[i]]$peptideSequence,
spec=data[[i]], fi=fi.by, ion.axes=FALSE,
main=paste("scantype: HCD / peptide sequence: ",
.PTM_MarkerFinder_(data[[i]]$peptideSequence, data[[i]]$modification, modificationName), sep=''),
xlim=c(0,max(data[[i]]$mZ)))
}else{
par(cex=0.75)
plot(data[[i]]$mZ, data[[i]]$intensity, type='h',
xlab='m/z',ylab='Intensity',
main=paste("scantype: HCD"),
xlim=c(0, max(data[[i]]$mZ)),
sub=paste(i, data[[i]]$title))
}
points(data[[i]]$mZ[idx[b]],
data[[i]]$intensity[idx[b]],
pch=22,
col='#6CC417AA',
bg='#6CC417AA',
cex=1)
legend("topright", paste(c('protein', 'm/z', 'charge', 'ionScore', 'query#'),
c(data[[i]]$proteinInformation,
data[[i]]$pepmass, data[[i]]$charge, data[[i]]$mascotScore,i)), cex=0.75)
#c(data[[i]]$proteinInformation,
################################################################################
####### P E A K P L O T ########################################################
screen(3)
j<-1;
PLOTFLAG<-FALSE
while (j<5){
#k<-i+j
k <- scan.to.query[query.to.scan[i] + j]
if (k < 0){break;}
if (!is.na(data[[k]]$mascotScore) & abs(data[[k]]$pepmass - data[[i]]$pepmass) < 1){
if (mgfFilename != -1){
.PTM_MarkerFinder_writeMGF(data[[k]], mgfFilename)
}
fi <- fragmentIon(sequence=data[[k]]$peptideSequence,
FUN=defaultIon,
modified=substr(data[[k]]$modification, 2, nchar(data[[k]]$modification)-1),
modification=modification)
fi.cyz <- as.data.frame(cbind(y=fi[[1]]$y, c=fi[[1]]$c, z=fi[[1]]$z), stringsAsFactors = TRUE)
p <- peakplot(data[[k]]$peptideSequence, spec=data[[k]],
main=paste("scantype: ETD / peptide sequence: ",
.PTM_MarkerFinder_(data[[k]]$peptideSequence, data[[k]]$modification, modificationName), sep=''),
fi=fi.cyz,
itol=0.6,
xlim=c(0, max(data[[i]]$mZ)),
ion.axes=FALSE)
PLOTFLAG<-TRUE
legend("topleft", paste(c('protein', 'm/z', 'charge', 'ionScore', 'query#'),
c(data[[k]]$proteinInformation,
data[[k]]$pepmass, data[[k]]$charge, data[[k]]$mascotScore, k)), cex=0.75)
break;
}
j <- j+1
}
if (PLOTFLAG == FALSE){
plot(0,0,type='n', xlab='',ylab='', axes=FALSE)
text(0,0,"no peptide assignment", cex=3, col="lightgrey")
}
################################################################################
screen(8)
plot(0,0,axes=FALSE,type='n', xlab='',ylab='')
################################################################################
screen(7)
pie(data[[i]]$intensity[idx[b]],
paste(round(100*data[[i]]$intensity[idx[b]]/sum(data[[i]]$intensity[idx[b]]),0),"% ",
round(data[[i]]$mZ[idx[b]], 2), sep=""),
main=list("relative marker ions %", cex=0.8),
cex=0.75)
################################################################################
screen(6)
pie(c(sum.intensity, sum.mZmarkerIons.intensity),
c('Intensity', paste(percent.mZmarkerIons, '%',sep='')),
col=c('white', '#6CC417AA'),
main=list("% marker ion intensities", cex=0.8),
cex=0.75)
################################################################################
screen(4)
#plot(mZmarkerIons, (mZmarkerIons-data[[i]]$mZ[idx]),
plot(mZmarkerIons[b], 1e+06 * (mZmarkerIons[b]-data[[i]]$mZ[idx[b]])/data[[i]]$mZ[idx[b]],
axes=TRUE,type='p', xlab='m/z',ylab='ppm error',log='')# ylim=c(-max(ppm.error),max(ppm.error)))
axis(3, mZmarkerIons[b], round(mZmarkerIons[b],1),las=2)
abline(h=0.0, col='grey')
################################################################################
####### P A T T E R N P L O T ##################################################
screen(5)
plot(data[[i]]$mZ[idx[b]],
data[[i]]$intensity[idx[b]],
col='#6CC417AA',
lwd=2,
type='h',
xlab='m/z',
ylab='Intensity',
axes=TRUE)
axis(3,data[[i]]$mZ[idx[b]],round(data[[i]]$mZ[idx[b]],1), las=2)
# readline("Press <Enter> to continue")
close.screen(all.screens = TRUE)
}
}
close.screen(all.screens = TRUE)
# TODO(cp): make a S3 class
rr <- as.data.frame(rr, stringsAsFactors = TRUE)
rr$markerIonIntensity <- as.numeric(levels(rr$markerIonIntensity))[rr$markerIonIntensity]
rr$mZ <- as.numeric(levels(rr$mZ))[rr$mZ]
rr$pepmass <- as.numeric(levels(rr$pepmass))[rr$pepmass]
rr$markerIonMZ <- as.numeric(levels(rr$markerIonMZ))[rr$markerIonMZ]
rr$markerIonMzError <- as.numeric(levels(rr$markerIonMzError))[rr$markerIonMzError]
return(rr)
}
ptmmf <- setClass("ptmmf",
slots=c(
scans="character",
mZ="numeric",
markerIonMZ="numeric",
markerIonIntensity="numeric",
markerIonMzError="numeric",
markerIonPpmError="numeric",
relativeFragmentIntensity="numeric",
query ="character",
pepmass="numeric",
peptideSequence="character",
modification="character")
)
# code for screening multiple modifications
# by Paolo <paolo.nanni@fgcz.uzh.ch> and Christian <cp@fgcz.ethz.ch>, March 2013
#
# processed using R version 2.15.3 (2013-03-01)
# protViz packageVersion 0.1.41
# Linux fgcz-c-071 2.6.32-5-amd64 #1 SMP Mon Jan 16 16:22:28 UTC 2012 x86_64 GNU/Linux
# Intel(R) Xeon(R) CPU X5650 @ 2.67GHz
.PTM_MarkerFinder_summaryPlot<-function(s, config){
plot(0,0,xlab='', ylab='', axes=FALSE, type='n')
legend(-0.95,+0.5, paste(c(names(config), 'R', 'protViz.version', 'create date'),
c(config, R.version.string, as.character(packageVersion('protViz')), date()),
sep=": "), cex=2.0)
if (nrow(s) > 1){
op<-par(mfrow=c(4,5), mar=c(4,4,3,1));
dump <- lapply(split(s,s$query),
function(x){
plot(x$mZ, x$markerIonIntensity,
type='h',
col='lightblue',
cex=2,
ylab='intensity', xlab='m/z',
ylim=range(s$markerIonIntensity),
log='y',
main=paste("scan=",unique(x$scans),"/query=", unique(x$query), sep=''));
text(x$mZ, x$markerIonIntensity,round(x$mZ,1),col='red',cex=0.7)
})
par(op)
op<-par(mfrow=c(4,5), mar=c(1,1,3,1));
dump <- lapply(split(s,s$query),
function(x){ pie(x$markerIonIntensity,
cex=0.7,
round(x$markerIonMZ,2),
col='white',
main=paste("scan=",unique(x$scans),"/query=", unique(x$query), sep=''))})
par(op)
op<-par(mfrow=c(1,1), mar=c(5,5,5,5));
boxplot(s$markerIonIntensity ~ as.factor(s$markerIonMZ),
log='y',
main='summary marker ion intensity distribution',
xlab='markerIon m/z', ylab='log10 based marker ion intensity')
box()
#boxplot(markerIonPpmError~(scans),
# data=s,
# main='marker ions PPM error summary',
# xlab='scan number', ylab='ppm error')
#abline(h=0.0,col='grey')
par(mfrow=c(1,1));
barplot(tapply(s$markerIonIntensity, s$scans, FUN=sum),log='y',
main='sum marker ion intensities summary',
xlab='scan numbers',
ylab='sum marker ion intensities [log10]'
)
}
}
PTM_MarkerFinder_util<-function(dataFileName,
mZmarkerIons=c(126.0550, 138.0550, 144.0655, 168.0655, 186.0761, 204.0866),
minMarkerIntensityRatio=2,
minNumberIons=2,
itol_ppm=20,
write_csv=FALSE){
config<-list(
dataFileName=dataFileName,
MarkerIon=mZmarkerIons,
minMarkerIntensityRatio=minMarkerIntensityRatio,
minNumberIons=minNumberIons,
itol_ppm=itol_ppm)
datFilename <- paste(config$dataFileName, "RData", sep='.')
message(paste("load file", datFilename, "..."))
load(datFilename)
pdf(paste(config$dataFileName, "pdf", sep='.'), 19, 12)
s<-PTM_MarkerFinder(data=get(config$dataFileName),
modification=get(paste(config$dataFileName, "modification.mass", sep=".")),
modificationName=get(paste(config$dataFileName, "modification.description", sep = '.')),
minMarkerIntensityRatio=config$minMarkerIntensityRatio,
minNumberIons=config$minNumberIons,
itol_ppm=config$itol_ppm,
mgfFilename=paste(config$dataFileName,".mgf", sep=''),
mZmarkerIons=config$MarkerIon)
dev.off()
# overwrite the it again.
pdf(paste(config$dataFileName, "pdf", sep='.'), 19, 12)
.PTM_MarkerFinder_summaryPlot(s, config)
s<-PTM_MarkerFinder(data=get(config$dataFileName),
modification=get(paste(config$dataFileName, "modification.mass", sep=".")),
modificationName=get(paste(config$dataFileName, "modification.description", sep = '.')),
minMarkerIntensityRatio=config$minMarkerIntensityRatio,
minNumberIons=config$minNumberIons,
itol_ppm=config$itol_ppm,
mZmarkerIons=config$MarkerIon)
dev.off()
if ( nrow(s) > 1 & write_csv==TRUE ){
w<-reshape(s[, c(1,7,8,3,4,9,10)],
direction='wide',
timevar="markerIonMZ",
idvar=c('scans', 'query', 'pepmass', 'peptideSequence', 'modification'))
write.table(w, file=paste(config$dataFileName, "csv", sep='.'),
sep=',', row.names=FALSE, col.names=TRUE, quote=FALSE)
}
# return(s)
}
# ToDo adapt the parameter mZmarkerIon in all function
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.