% \VignetteIndexEntry{Transformation and other Factors of the pairwise Mass Spectrometry peak-list Comparison Process}
% \VignetteDepends{}
% \VignetteKeywords{Mass Spectrometry}
% \VignettePackage{msbase} 
\documentclass[a4paper]{article}

\title{Transformation and other Factors of the pairwise Mass Spectrometry peak-list Comparison Process.(PRELIMINARY!!!)}
\author{Peter Martus and Witold Eryk Wolski}

\begin{document}

\maketitle

\setkeys{Gin}{width=0.6\textwidth}

\section*{Anova analysis of the PMF dataset}

We analyse here the 4 datasets (pmf.binary ,pmf.intensity, msms.binary, msms.intensity).
The datasets provide results of evaluating the sensitivity and specificity of the pairwise peak-list comparison
performed on an dataset of identified Tandem MS data (msms) and on an dataset of identified Peptide Mass Fingerprint spectra.
(publication submitted).

Load the results for the binary measures.

<<Readingandpreparing>>=
#binary based measures#
library(msbase)
data(pmf.binary)
pmf.binary <- pmf.binary[,c(2:6,9:10)]
pmf.binary$FPPAUC <- pmf.binary$FPPAUC*1000
pmf.binary$TPPAUC <- pmf.binary$TPPAUC*1000
pmf.binary$Theta <- as.factor(pmf.binary$Theta)
@
 
The minimal linear model (containing as few factors as possible) which sufficiently describes the outcome the specificity-PAUC (given small FP rates)
 of the experiment is given by.

<<anovamodel1PMFbinary>>=
tplm <- lm(TPPAUC ~ Measure + Theta + Length + Measure*Theta + Measure*Length + Theta*Length + Measure*Theta*Length,data=pmf.binary)
tplm <- lm(TPPAUC ~ Measure + Theta + Length + Measure*Theta , data=pmf.binary)
summary(tplm)$adj.r
anova(tplm)
@

To identify the best measure we compute the average S-PAUC for each CP having the same length , thata and measure factor.

<<tapplyTP>>=
with(pmf.binary,tapply(TPPAUC,list(Length=Length,Theta=Theta,Measure=Measure),mean))
@

The same model describe the specificity(Sp)-PAUC.

<<ANOVAFPPAUC>>= 
tplm <- lm(FPPAUC ~ Measure+Theta+Length+Measure:Theta+Measure:Length+Theta:Length+Measure:Theta:Length,data=pmf.binary)
anova(tplm)
tplm <- lm(FPPAUC ~ Measure*Theta+Measure:Length+Theta:Length+Measure:Theta:Length,data=pmf.binary)
anova(tplm)
summary(tplm)$adj.r
@
 
And again the average given this three factors is computed.

<<tapplyFPPAUCbinaryPMF>>=
with(pmf.binary,tapply(FPPAUC,list(Length=Length,Theta=Theta,Measure=Measure),mean))
@

<<>>=
boxplot(FPPAUC~Length*Measure,data=pmf.binary)
par(mar=c(6,3,3,3))
boxplot(TPPAUC~Length*Theta*Measure,data=pmf.binary,las=2,ylim=c(95,100))
abline(v=1:50,col="gray")
@ 
Looking at the output of the tapply function we identified the Fowlkes Mallows statistics as the best measure.

\subsubsection*{Fowlkes mallows}

<<>>=
bingow <- pmf.binary[pmf.binary$Measure=="gower",]
tplm <- lm(TPPAUC~Theta,data=bingow)
summary(tplm)$adj.r
anova(tplm)
@

Check if any sensible happens. 

<<fig=T>>=
par(mfrow=c(2,2))
boxplot(TPPAUC~Noncross,data=bingow,main="S-PAUC noncross")
boxplot(TPPAUC~Weight,data=bingow,main="S-PAUC weight")
boxplot(FPPAUC~Noncross,data=bingow,main="Sp-PAUC noncross")
boxplot(FPPAUC~Weight,data=bingow,main="Sp-PAUC weight")

@ 


Make some more boxplots on what the outcome of the comparison depends.

<<MeasureTheta,fig=TRUE>>=
boxplot(TPPAUC~Measure*Theta,data=pmf.binary,main="S-PAUC meausure",las=2)
@ 

<<fig=TRUE>>=
par(mfcol=c(2,3))
boxplot(TPPAUC~Measure*Theta,data=pmf.binary,main="S-PAUC meausure",las=2)
boxplot(FPPAUC~Measure,data=pmf.binary,main="Sp-PAUC meausure")
boxplot(TPPAUC~Theta,data=pmf.binary,main="S-PAUC theta")
boxplot(FPPAUC~Theta,data=pmf.binary,main="Sp-PAUC theta")
boxplot(TPPAUC~Length,data=pmf.binary,main="S-PAUC length")
boxplot(FPPAUC~Length,data=pmf.binary,main="Sp-PAUC length")
@ 

\section*{PMF data - intensity based maesures}

Load the evaluation results.

<<LOADANDPREPARE>>=
data(pmf.intensity)
pmf.intensity$Measure<-factor(pmf.intensity$Measure,levels=c("canberra","simindex","manhattan","euclidean","dotprod","cov","soai"))
pmf.intensity$Scale<-factor(pmf.intensity$Scale,levels=c("T","N","S","Z"))
pmf.intensity$Trans<-factor(pmf.intensity$Trans,levels=c("N","S","L","R"))
pmf.intensity<-pmf.intensity[,c(2:8,11:12)]
pmf.intensity$FPPAUC <- pmf.intensity$FPPAUC*1000
pmf.intensity$TPPAUC <- pmf.intensity$TPPAUC*1000
pmf.intensity$Theta <- as.factor(pmf.intensity$Theta)
@ 

<<>>=
par(mar=c(8,2,2,2))
boxplot(FPPAUC~Length+Measure+Scale+Theta,data=pmf.intensity,las=2,cex.axis=0.5)
abline(v=1:250,col="gray")

tmpb<- pmf.intensity[pmf.intensity$Measure=="manhattan" & pmf.intensity$Scale=="T",]
boxplot(FPPAUC~Length+Theta,data=tmpb,las=2,cex.axis=0.5)

@ 

The minimal model explaining as much as possible variance is:
<<>>=
intlm<-lm(FPPAUC ~ Measure+Scale+Theta+Length+Measure:Scale+Measure:Theta+Measure:Length,data=pmf.intensity)
anova(intlm)
summary(intlm)$adj.r

@ 

<<IntensityPMFBoxplot,fig=TRUE>>=

#pdf(file="IntensityPMFBoxplot.pdf",width=7,height=6)
nf <- layout(matrix(c(1,2,3,4),2,2,byrow=TRUE), c(1.5,2.4),c(1,3), TRUE)
par(mar=c(1,3.5,3,1))
boxplot(TPPAUC~Theta+Measure,data=pmf.binary,las=2,axes=F,ylim=c(97.5,98.5),cex.axis=0.9)
axis(2,cex.axis=0.9,col=3)
box(,col=3)
abline(v=c(3.5,6.5,9.5),col="gray70")
abline(h=max(pmf.intensity$TPPAUC),lty=2,col="gray50")
abline(h=max(pmf.binary$TPPAUC),lty=4,col="gray50")
mtext("S-PAUC",side=2,line=2,cex=1)
mtext("Measure + Theta (Bin.)",side=3,line=1)
mtext("A",side=3,line=1,at=-3,cex=1.5)

par(mar=c(1,3,3,1))
boxplot(TPPAUC ~ Scale+Measure,data=pmf.intensity,las=2,axes=F,ylim=c(97.5,98.5),cex.axis=0.9)
abline(v=c(4.5,8.5,12.5,16.5,20.5,24.5),col="gray70")
abline(v=max(pmf.binary$TPPAUC))
axis(2,cex.axis=0.9,col=3)
box(,col=3)
abline(h=max(pmf.intensity$TPPAUC),lty=2,col="gray50")
abline(h=max(pmf.binary$TPPAUC),lty=4,col="gray50")
mtext("S-PAUC",side=2,line=2,cex=1)
mtext("Measure+Scaling (Int.)",side=3,line=1)
mtext("B",side=3,line=1,at=-1,cex=1.5)

par(mar=c(6,3.5,0,1))
plot(1,1,xlim=c(0,13),ylim=c(0,100),type="n",axes=F,ylab="",xlab="",main="")
#abline(h=c(97.5,98.5),col=3)
rect(-1,97.5,14,98.5,col="lightgreen",border="lightgreen")
par(new=TRUE)
boxplot(TPPAUC~Theta+Measure,data=pmf.binary,las=2,ylim=c(0,100),ylab="",cex.axis=0.9)
mtext("S-PAUC",side=2,line=2)
abline(v=c(3.5,6.5,9.5),col="gray70")

par(mar=c(6,3,0,1))
plot(1,1,xlim=c(0,25),ylim=c(0,100),type="n",axes=F,ylab="",xlab="",main="")
rect(-4,97.5,30,98.5,col="lightgreen",border="lightgreen")
par(new=TRUE)
boxplot(TPPAUC ~ Scale+Measure,data=pmf.intensity,las=2,type="n",ylim=c(0,100),cex.axis=0.9)
mtext("S-PAUC",side=2,line=2,cex=1)
abline(v=c(4.5,8.5,12.5,16.5,20.5,24.5),col="gray70")
#dev.off()

@ 

Now we tabulate the scores according to the identified factors.

<<>>=
with(pmf.intensity,tapply(FPPAUC,list(Length=Length,Theta=Theta,Scale=Scale,Measure=Measure),mean,data=pmf.intensity))
@ 

The same procedure with the second score and again the same model can be used.

<<>>=
intlm <-lm(TPPAUC ~ Measure+Scale+Theta+Length+Measure*Scale+Measure*Theta+Measure*Length,data=pmf.intensity)
anova(intlm)
summary(intlm)$adj.r
@ 

We tabulated the scores accroding to the identified factors.

<<>>=
with(pmf.intensity,tapply(TPPAUC,list(Length=Length,Theta=Theta,Scale=Scale,Measure=Measure),mean))
@
 
By analysing the table we identify the dotproduct measure computed on vector norm scaled data as performing best.
Also the euclidean and manhatten distances perfrom well but only with a sensible parameter choice.

<<>>=

intsoai <- pmf.intensity[(pmf.intensity$Measure == "soai" | pmf.intensity$Measure == "dotprod") & pmf.intensity$Scale == "N",]
intdp <- pmf.intensity[pmf.intensity$Measure == "soai" & pmf.intensity$Scale == "S",]

lmdp <-lm(FPPAUC~Length+Theta+Trans+Noncross+Weight,data=intdp)
anova(lmdp)
summary(lmdp)$adj.r
@ 

By the anova anlysis above we idnetify the Transformation as the only parameter influencing the DP measure

<<TransformIntensity,fig=TRUE>>=

#pdf(file="TransformationIntenisty.pdf",width=5,height=5)
par(mfrow=c(1,2))
par(mar=c(6,4,1,1))
boxplot(FPPAUC~Trans+Measure,data=intsoai,main="",las=2)
mtext("Specificity-PAUC",side=2,line=3)
mtext("A",side=3,line=-1,at=-3,cex=2)
abline(v=4.5,col="gray")
boxplot(TPPAUC~Trans+Measure,data=intsoai,main="",las=2)
mtext("Sensitivity-PAUC",side=2,line=3)
mtext("B",side=3,line=-1,at=-3,cex=2)
abline(v=4.5,col="gray")
#dev.off()

@ 

The boxplot shows that the log transformation gives highest sensitivities given small FP-rates as well highest specificities given high sensitivities.

\section{Analysing the MS/MS dataset}


The MSMS dataset we use to examine if the conclusion drawn anylysing the PMF dataset can be generalized.



\subsection{Binary measure}

<<ReadTHERESULTS>>=
data(msms.binary)
msms.binary <- msms.binary[,c(2:6,9:10)]
msms.binary$FPPAUC <- msms.binary$FPPAUC*1000
msms.binary$TPPAUC <- msms.binary$TPPAUC*1000
msms.binary$Theta <- as.factor(msms.binary$Theta)
tmpA<-msms.binary$FPPAUC[msms.binary$Noncross=="A"]
tmpO<-msms.binary$FPPAUC[msms.binary$Noncross=="O"]
sum(tmpA>tmpO) #35
sum(tmpA<tmpO) #13
tmpA<-msms.binary$TPPAUC[msms.binary$Noncross=="A"]
tmpO<-msms.binary$TPPAUC[msms.binary$Noncross=="O"]
sum(tmpA>tmpO) # 43
sum(tmpA<tmpO) #3

@ 


<<ReadINTENMSMS>>=
data(msms.intensity)
msms.intensity <- msms.intensity[,c(2:8,11:12)]
pmf.intensity$Measure <- factor(pmf.intensity$Measure,levels=c("canberra","simindex","manhattan","euclidean","dotprod","cov","soai"))
pmf.intensity$Scale <- factor(pmf.intensity$Scale,levels=c("T","N","S","Z"))
pmf.intensity$Trans <- factor(pmf.intensity$Trans,levels=c("N","S","L","R"))
msms.intensity$FPPAUC <- msms.intensity$FPPAUC*1000
msms.intensity$TPPAUC <- msms.intensity$TPPAUC*1000
msms.intensity$Theta <- as.factor(pmf.intensity$Theta)

tmpA<-msms.intensity$FPPAUC[msms.intensity$Weight=="W" & msms.intensity$Measure!="dotprod"]
tmpO<-msms.intensity$FPPAUC[msms.intensity$Weight=="M" & msms.intensity$Measure!="dotprod"]
sum(tmpA>tmpO) #35
sum(tmpA<tmpO) #13
tmpA<-msms.intensity$TPPAUC[msms.intensity$Weight=="W" & msms.intensity$Measure!="dotprod"]
tmpO<-msms.intensity$TPPAUC[msms.intensity$Weight=="M" & msms.intensity$Measure!="dotprod"]
sum(tmpA>tmpO) # 43
sum(tmpA<tmpO) #3


tmpA<-msms.intensity$FPPAUC[msms.intensity$Noncross=="A"]
tmpO<-msms.intensity$FPPAUC[msms.intensity$Noncross=="O"]
sum(tmpA>tmpO) #35
sum(tmpA<tmpO) #13
tmpA<-msms.intensity$TPPAUC[msms.intensity$Noncross=="A"]
tmpO<-msms.intensity$TPPAUC[msms.intensity$Noncross=="O"]
sum(tmpA>tmpO) # 43
sum(tmpA<tmpO) #3

@ 

The descriptive plot.

<<fig=TRUE>>=

#pdf(file="MSMSres.pdf",width=5,height=5)
idfm <- rep("other CP",length(msms.binary$Measure))
idfm[which(msms.binary$Measure=="fm")] <- "fm"
par(mfrow=c(1,2))
par(mar=c(7,4,3,1))
#boxplot(msms.binary$TPPAUC~as.factor(idfm),main="Sp-PAUC")
boxplot(TPPAUC~as.factor(idfm),data=msms.binary,las=2,main="")
mtext("binary",side=3,line=1)
mtext("sensitivity-PAUC",side=2,line=2.5,cex=1.2)
mtext("A",side=3,line=1,at=0,cex=2)
par(mar=c(7,4,3,1))
ind <- rep("other CP",length(msms.intensity$Measure))
ind[msms.intensity$Measure == "dotprod" & msms.intensity$Scale == "N" & msms.intensity$Trans=="L"]<-"dot-product"
#par(mfrow=c(1,2))
#boxplot(TPPAUC~as.factor(ind),data=msms.intensity,main="Sp-PAUC - MS/MS")
boxplot(TPPAUC~as.factor(ind),data=msms.intensity,las=2,main="")
mtext("intensity",side=3,line=1)
mtext("B",side=3,line=1,at=0,cex=2)
mtext("sensitivity-PAUC",side=2,line=2.5,cex=1.2)
#dev.off()

@ 

In case of the MS/MS data there is not such a unambiguous superiority of the fowlkes mallows statistics compared to other measures.

\subsubsection*{Some further analysis}

<<OLDMODEL>>=
tplm <- lm(TPPAUC ~ Measure+Theta+Length+Measure*Theta+Measure*Length+Theta*Length+Measure*Theta*Length,data=msms.binary)
anova(tplm)
summary(tplm)$adj.rq
@ 
We see that the old model is not performing as well as in case of the PMF data.
A model which includes computing the noncrossing matching performs much better.

<<SUMMARYX>>=
tplm <- lm(TPPAUC~Measure*Theta*Length*Noncross,data=msms.binary)
summary(tplm)$adj
anova(tplm)
@ 

Tabulate the data according to identified factors.

<<>>=
with(msms.binary,tapply(TPPAUC,list(Length=Length,Theta=Theta,Noncross=Noncross,Measure=Measure),mean))
@ 
We see that the binary measures computed with noncrossing matching perform better than they associates without. Furthermore Huberts Gamma computed with $M_{00}=0$ can be recognized as the best performing measure.

The same analysis is repeated for the Sp-PAUC.

<<SUMMARYX>>=
boxplot(FPPAUC~Measure,data=msms.binary)
tplm <- lm(FPPAUC~Measure*Theta*Length*Noncross,data=msms.binary)
summary(tplm)$adj
anova(tplm)
@ 

Tabulate the data according to identified factors.

<<>>=
with(msms.binary,tapply(TPPAUC,list(Length=Length,Theta=Theta,Noncross=Noncross,Measure=Measure),mean))
@ 

Again hg performs best.

Further we analysed how the HG measure depends on factors like e.g. Weighting of mass measurement error, noncrossing matching and theta.

<<MORBOX,fig=TRUE>>=
msms.binaryhg<-msms.binary[msms.binary$Measure=="hg",]
par(mfcol=c(2,3))
boxplot(TPPAUC ~ Weight+Noncross,data=msms.binaryhg,main="S-PAUC weight")
boxplot(FPPAUC ~ Weight,data=msms.binaryhg,main="Sp-PAUC weight")
boxplot(TPPAUC ~ Noncross,data=msms.binaryhg,main="S-PAUC noncross")
boxplot(FPPAUC ~ Noncross,data=msms.binaryhg,main="Sp-PAUC noncross")
boxplot(TPPAUC ~ Theta,data=msms.binaryhg,main="S-PAUC theta")
boxplot(FPPAUC ~ Theta,data=msms.binaryhg,main="Sp-PAUC theta")
@ 

The boxplots reveal that weighting of match accuracy and resolving unambigous matches by computing the noncrossing matching increases the performance of the HG measure. Furthermore the optimal choice of the theta is 2.


\subsection{MS/MS-Intensity based measures}

Again we first prove if the result obtained for the PMF data can be generalized to the MS/MS data.

<<EXPLOT2,fig=TRUE>>=
ind <- rep("other measure",length(msms.intensity$Measure))
ind[msms.intensity$Measure == "dotprod" & msms.intensity$Scale == "N" & msms.intensity$Trans=="L"]<-"dotprod ; N(log(I))"
par(mfrow=c(1,2))
boxplot(TPPAUC~as.factor(ind),data=msms.intensity,main="S-PAUC - MS/MS")
boxplot(FPPAUC~as.factor(ind),data=msms.intensity,main="Sp-PAUC - MS/MS")
@ 

This time the observation done using the PMF data can be generalized to the MS/MS data.
Because we were interested to identify the other measures which can be used to classify the data we tabulated the scores
according to theta, scaling, length and the measures.

<<MODELFPPAUCMSMS>>=
intlm <- lm(FPPAUC ~ Measure+Scale+Theta+Length+Measure:Scale+Measure:Theta+Measure:Length,data=msms.intensity)
anova(intlm)
summary(intlm)$adj.r
@

<<TABULATEINTMSMS>>=
with(msms.intensity,tapply(FPPAUC,list(Length=Length,Theta=Theta,Scale=Scale,Measure=Measure),mean))
@ 

<<MODELTPPAUCMSMS>>=
intlm <-lm(TPPAUC ~ Measure+Scale+Theta+Length+Measure*Scale+Measure*Theta+Measure*Length,data=msms.intensity)
anova(intlm)
summary(intlm)$adj.r
@

<<TabulateTheMEASURES>>=
with(msms.intensity,tapply(TPPAUC,list(Length=Length,Theta=Theta,Scale=Scale,Measure=Measure),mean))
@ 

The spectral angle measure is the highest scoring one. None of the other measures is able to obtain similar scores. This may be due to the fact that the database search of the MS/MS data is performed using the normalized crosscorelation which has a very similar mathematical property than the spectra angle.

Finally we analyse how factors like intensity transformation, weighting of mass measurement accuracy and computing the noncrossing matching influences the performance of the spectral angle.

<<EXTRACTDPMSMS>>=
intdp <- msms.intensity[(msms.intensity$Measure == "euclidean") & msms.intensity$Scale == "S" & msms.intensity$Trans=="L",]
boxplot(TPPAUC~Weight*Noncross*Measure*Trans,data=intdp,main="S-PAUC weight",las=2)
par(mar=c(8,5,2,2))
boxplot(FPPAUC~Noncross,data=intdp,main="S-PAUC weight",las=2)

boxplot(TPPAUC~Weight,data=intdp,main="S-PAUC weight",las=2)

@ 
<<>>=

boxplot(TPPAUC~Weight*Noncross*Theta*Measure,data=intdp,main="S-PAUC weight",las=2)
@ 
<<ANALYSEDPMSMS>>=
lmdp<-lm(FPPAUC~Weight+Noncross+Weight:Noncross,data=intdp)
anova(lmdp)
summary(lmdp)$adj.r
@ 

<<BOXPLOTDP,fig=TRUE>>=
par(mfcol=c(2,3))
boxplot(TPPAUC~Trans,data=intdp,main="S-PAUC trans")
boxplot(FPPAUC~Trans,data=intdp,main="Sp-PAUC trans")

boxplot(FPPAUC~Weight,data=intdp,main="Sp-PAUC weight")
boxplot(FPPAUC~Noncross,data=intdp,main="S-PAUC noncross")
boxplot(TPPAUC~Noncross,data=intdp,main="Sp-PAUC noncross")
@ 

\end{document}
