### R code from vignette source 'BRAIN-vignette.Rnw'

###################################################
### code chunk number 1: BRAIN-vignette.Rnw:14-15
###################################################
library(BRAIN)


###################################################
### code chunk number 2: BRAIN-vignette.Rnw:65-74
###################################################
#angiotensineII
angiotensineII <- list(C=50,H=71,N=13,O=12) 
monoisotopicMassAngiotensineII <- calculateMonoisotopicMass(aC = angiotensineII)
monoisotopicMassAngiotensineII

#human dynein heavy chain
humandynein <- list(C=23832,H=37816,N=6528,O=7031,S=170) 
monoisotopicMassHumandynein <- calculateMonoisotopicMass(aC = humandynein)
monoisotopicMassHumandynein


###################################################
### code chunk number 3: BRAIN-vignette.Rnw:89-96
###################################################
#angiotensineII
averageMassAngiotensineII <- calculateAverageMass(aC = angiotensineII)
averageMassAngiotensineII

#humandynein
averageMassHumandynein <- calculateAverageMass(aC = humandynein)
averageMassHumandynein


###################################################
### code chunk number 4: BRAIN-vignette.Rnw:102-110
###################################################
#angiotensineII
nrPeaksAngiotensineII <- calculateNrPeaks(aC = angiotensineII)
nrPeaksAngiotensineII


#human dynein heavy chain
nrPeaksHumandynein <- calculateNrPeaks(aC = humandynein)
nrPeaksHumandynein


###################################################
### code chunk number 5: BRAIN-vignette.Rnw:117-152
###################################################
#angiotensineII

#with default options
prob1 <- calculateIsotopicProbabilities(aC = angiotensineII)
print(length(prob1))

#with user defined number of requested aggregated isotopic variants
prob2 <- calculateIsotopicProbabilities(aC = angiotensineII, nrPeaks=20)
print(length(prob2))

#with user-defined coverage as stopping criterium
prob3 <- calculateIsotopicProbabilities(aC = angiotensineII, 
stopOption = "coverage", coverage = 0.99)
print(length(prob3))  

#with user-defined abundantEstim as stopping criterium
prob4 <- calculateIsotopicProbabilities(aC = angiotensineII, 
stopOption = "abundantEstim", abundantEstim = 10)
print(length(prob4))


#human dynein heavy chain
prob1 <- calculateIsotopicProbabilities(aC = humandynein)
print(length(prob1))

prob2 <- calculateIsotopicProbabilities(aC = humandynein, nrPeaks=300)
print(length(prob2))

prob3 <- calculateIsotopicProbabilities(aC = humandynein, 
stopOption = "coverage", coverage = 0.99)
print(length(prob3))

prob4 <- calculateIsotopicProbabilities(aC = humandynein, 
stopOption = "abundantEstim", abundantEstim = 150)
print(length(prob4))


###################################################
### code chunk number 6: BRAIN-vignette.Rnw:158-171
###################################################
#angiotensineII
headr <- expression(paste(C[50], H[71], N[13],O[12]))
#with default options
res <- useBRAIN(aC= angiotensineII) 
plot(res$masses,res$isoDistr,xlab="mass",ylab="abundances", type="h", 
xlim=c(min(res$masses)-1,max(res$masses)+1.5))
title(headr)
labelMono <- paste("mono-isotopic mass:", res$monoisotopicMass, "Da", sep=" ")
labelAvg <- paste("average mass:", res$avgMass, "Da", sep=" ")
text(x=1048.7, y=0.5, labelMono, col="purple")
text(x=1048.7, y=0.45, labelAvg, col="blue")
lines(x=rep(res$monoisotopicMass[1],2),y=c(0,res$isoDistr[1]), col = "purple")
lines(x=rep(res$avgMass,2),y=c(0,max(res$isoDistr)), col = "blue", lty=2)


###################################################
### code chunk number 7: BRAIN-vignette.Rnw:175-190
###################################################
#human dynein heavy chain
headr <- expression(paste(C[23832], H[37816], N[6528],O[7031], S[170]))
res <- useBRAIN(aC=humandynein, stopOption="coverage", coverage=0.99) 
plot(res$masses,res$isoDistr,xlab="mass",ylab="abundances", type="h", 
xlim=c(min(res$masses)-1,max(res$masses)+1))
title(headr)
labelMono <- paste("mono-isotopic mass: ", res$monoisotopicMass, "Da", sep="")
labelAvg <- paste("average mass: ", res$avgMass, "Da", sep="")
mostAbundant <- res$masses[which.max(res$isoDistr)]
labelAbundant <- paste("most abundant mass: ", mostAbundant, "Da", sep="")
text(x=533550, y=0.02, labelMono, col="purple")
text(x=533550, y=0.0175, labelAvg, col="blue")
text(x=533550, y=0.015, labelAbundant, col="red")
lines(x=rep(res$avgMass,2),y=c(0,max(res$isoDistr)), col = "blue", lty=2)
lines(x=rep(mostAbundant,2),y=c(0,max(res$isoDistr)), col = "red")


###################################################
### code chunk number 8: BRAIN-vignette.Rnw:199-214
###################################################
#with user defined number of requested aggregated isotopic variants
res <- useBRAIN(aC = angiotensineII, nrPeaks = 20) 
plot(res$masses,res$isoDistr,xlab="mass",ylab="abundances", type="h", 
xlim=c(min(res$masses)-1,max(res$masses)+1))
title(headr)
#with user defined coverage as stopping criterium
res <- useBRAIN(aC = angiotensineII, stopOption = "coverage", coverage = 0.99) 
plot(res$masses,res$isoDistr,xlab="mass",ylab="abundances", type="h", 
xlim=c(min(res$masses)-1,max(res$masses)+1))
title(headr)
#with user defined abundantEstim as stopping criterium
res <- useBRAIN(aC = angiotensineII, stopOption = "abundantEstim", abundantEstim = 10) 
plot(res$masses,res$isoDistr,xlab="mass",ylab="abundances", type="h", 
xlim=c(min(res$masses)-1,max(res$masses)+1))
title(headr)


###################################################
### code chunk number 9: BRAIN-vignette.Rnw:232-235 (eval = FALSE)
###################################################
## tabUniprot <- read.table('data/uniprot.tab', sep="\t", header=TRUE) 
## tabUniprot$Sequence <- gsub(" ", "", tabUniprot$Sequence)
## howMany <- nrow(tabUniprot)


###################################################
### code chunk number 10: BRAIN-vignette.Rnw:245-248 (eval = FALSE)
###################################################
## AMINOS <- c("A", "R", "N", "D", "C", "E", "Q", "G", "H", "I", "L", "K", 
## "M", "F", "P", "S", "T", "W", "Y", "V")
## !is.na(sum(match(seqVector, AMINOS)))


###################################################
### code chunk number 11: BRAIN-vignette.Rnw:256-282 (eval = FALSE)
###################################################
## nrPeaks = 1000
## howMany <- nrow(tabUniprot)
## dfUniprot <- data.frame()
## for (idx in 1:howMany){
##   seq <- as.character(tabUniprot[idx,"Sequence"])
##   seqVector <- strsplit(seq, split="")[[1]]
##   if (!is.na(sum(match(seqVector, AMINOS)))){  
##   aC <- getAtomsFromSeq(seq)
##   res <- useBRAIN(aC = aC, nrPeaks = nrPeaks, stopOption = "abundantEstim", 
## abundantEstim = 10)
##   isoDistr <- res$isoDistr
##   masses <- res$masses
##   maxIdx <- which.max(isoDistr)
##   mostAbundantPeakMass <- masses[maxIdx]
##   monoMass <- res$monoisotopicMass  
##   dfAtomicComp <- data.frame(C=aC[1], H=aC[2], N=aC[3], O=aC[4], S=aC[5])
##   singleDfUniprot <- data.frame(dfAtomicComp, monoMass=monoMass, maxIdx=maxIdx, 
## mostAbundantPeakMass=mostAbundantPeakMass  
##   )
##   if (!is.na(monoMass)){ #for huge atomic configuration numerical problems may occur
##     dfUniprot <- rbind(dfUniprot, singleDfUniprot)
##   }
##   }
## }
## 
## write.table(unique(dfUniprot), "data/uniprotBRAIN.txt")


###################################################
### code chunk number 12: BRAIN-vignette.Rnw:293-295
###################################################
uniprotBRAIN <- read.table(system.file("extdata", "uniprotBRAIN.txt", package="BRAIN"))
nrow(uniprotBRAIN)


###################################################
### code chunk number 13: uniprot_10_5
###################################################
uniprot10to5 <- uniprotBRAIN[uniprotBRAIN$monoMass < 10^5,]
nrow(uniprot10to5)


###################################################
### code chunk number 14: BRAIN-vignette.Rnw:308-315
###################################################
library(lattice)
mm <- uniprot10to5$monoMass
mmMin <- floor(min(mm)) - 1
mmMax <- ceiling(max(mm)) + 1
sq <- seq(from=0, to=mmMax, by=5000)
bw <- bwplot(cut(monoMass, sq) ~ mostAbundantPeakMass, data=uniprot10to5, ylab="monoMass (Da)", xlab="mostAbundantPeakMass (Da)")
plot(bw)


###################################################
### code chunk number 15: BRAIN-vignette.Rnw:318-321
###################################################
bw2 <- bwplot(cut(monoMass, sq) ~ (mostAbundantPeakMass - monoMass), 
data=uniprot10to5,  ylab="monoMass (Da)", xlab="mostAbundantPeakMass - monoMass (Da)")
plot(bw2)


###################################################
### code chunk number 16: BRAIN-vignette.Rnw:331-333
###################################################
lmod <- lm(monoMass ~ mostAbundantPeakMass, data=uniprot10to5)
summary(lmod)


###################################################
### code chunk number 17: BRAIN-vignette.Rnw:352-359
###################################################
icptLm <- lmod$coefficients[1]
cffLm <- lmod$coefficients[2]
expected <- icptLm + cffLm * uniprot10to5$mostAbundantPeakMass
residuals <- uniprot10to5$monoMass - expected
hist <- hist(residuals, seq(floor(min(residuals)),ceiling(max(residuals)), by=0.1), 
main="", xlab="residuals of the linear model (Da)")
plot(hist)