### 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)