FASTA p1875 db10 / WU160118 / R444589 - Postprocessing / Summary
2017-11-02, 2017-11-03, 2017-11-04, 2017-11-06, 2017-11-08, 2018-01-30, 2018-02-01, 2018-12-13
Abstract
This vignette contains code snippets to display results of the NestLink project (p1875). Briefly, it analyzes the amno acid frequencies of the in-silico composed and measured flycodes.
Package
NestLink 1.20.0
Contents
1 Load requiered R packages
library(NestLink)2 Display ESP
ESP_Prediction was generated using an application https://genepattern.broadinstitute.org (???).
library(ggplot2)
ESP <- rbind(getFC(), getNB())## see ?NestLink and browseVignettes('NestLink') for documentation## loading from cache## see ?NestLink and browseVignettes('NestLink') for documentation## loading from cachep <- ggplot(ESP, aes(x = ESP_Prediction, fill = cond)) +
geom_histogram(bins = 50, alpha = 0.4, position="identity") +
labs(x = "detectability in LC-MS (ESP prediction)") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))
print(p)
ESP <- rbind(getFC(), NB.unambiguous(getNB()))## see ?NestLink and browseVignettes('NestLink') for documentation
## loading from cache## see ?NestLink and browseVignettes('NestLink') for documentation## loading from cachep <- ggplot(ESP, aes(x = ESP_Prediction, fill = cond)) +
geom_histogram(bins = 50, alpha = 0.4, position="identity") +
labs(x = "detectability in LC-MS (ESP prediction)") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))
print(p)
ESP <- rbind(getFC(), NB.unique(getNB()))## see ?NestLink and browseVignettes('NestLink') for documentation
## loading from cache## see ?NestLink and browseVignettes('NestLink') for documentation## loading from cachep <- ggplot(ESP, aes(x = ESP_Prediction, fill = cond)) +
geom_histogram(bins = 50, alpha = 0.4, position="identity") +
labs(x = "detectability in LC-MS (ESP prediction)") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))
print(p)
ESP <- rbind(getNB(), NB.unambiguous(getNB()), NB.unique(getNB()))## see ?NestLink and browseVignettes('NestLink') for documentation
## loading from cache## see ?NestLink and browseVignettes('NestLink') for documentation## loading from cache## see ?NestLink and browseVignettes('NestLink') for documentation## loading from cachep <- ggplot(ESP, aes(x = ESP_Prediction, fill = cond)) +
geom_histogram(bins = 50, alpha = 0.4, position="identity") +
labs(x = "detectability in LC-MS (ESP prediction)") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))
print(p)
table(ESP$cond)##
## NB NB.unambiguous NB.unique
## 27681 5832 77053 Charge State Frequency
3.1 Load workunit 160118
# library(ExperimentHub)
# eh <- ExperimentHub();
# load(query(eh, c("NestLink", "WU160118.RData"))[[1]])
load(getExperimentHubFilename("WU160118.RData"))
WU <- WU160118filtering
PATTERN <- "^GS[ASTNQDEFVLYWGP]{7}(WR|WLTVR|WQEGGR|WLR|WQSR)$"
idx <- grepl(PATTERN, WU$pep_seq)
WU <- WU[idx & WU$pep_score > 25,]determine charge state frequency through counting
WU$chargeInt <- as.integer(substr(WU$query_charge, 0, 1))
table(WU$chargeInt)##
## 2 3 4
## 33392 157 1in percent
round(100 * (table(WU$chargeInt) / nrow(WU)), 1)##
## 2 3 4
## 99.5 0.5 0.0as histograms
library(scales)
ggplot(WU, aes(x = moverz * chargeInt, fill = (query_charge),
colour = (query_charge))) +
facet_wrap(~ datfilename, ncol = 2) +
geom_histogram(binwidth= 10, alpha=.3, position="identity") +
labs(title = "Precursor mass to charge frequency plot",
subtitle = "Plotting frequency of precursor masses for each charge state",
x = "Precursor mass [neutral mass]",
y = "Frequency [counts]",
fill = "Charge State",
colour = "Charge State") +
scale_x_continuous(breaks = pretty_breaks(8)) +
theme_light()
We confirmed this prediction by experimental data and found that 99.9 percent of flycode precursor ions correspond to doubly charge species (data not shown). The omission of positively charged residues is also critical in order to render trypsin a site-specific protease.
3.2 FlyCode Mass Distribution
WU$suffix <- substr(WU$pep_seq, 10, 100)
ggplot(WU, aes(x = moverz * chargeInt, fill = suffix, colour = suffix)) +
geom_histogram(binwidth= 10, alpha=.2, position="identity") +
#facet_wrap(~ substr(pep_seq, 10, 100)) +
theme_light()
ggplot(WU, aes(x = moverz * chargeInt, fill = suffix)) +
geom_histogram(binwidth= 10, alpha=.9, position="identity") +
facet_wrap(~ substr(pep_seq, 10, 100)) +
theme_light()
3.3 Mascot Ion Score Distribution
ggplot(WU, aes(x = pep_score, fill = query_charge, colour = query_charge)) +
geom_histogram(binwidth= 2, alpha=.5, position="identity") +
facet_wrap(~ substr(pep_seq, 10, 100)) +
theme_light()
4 Compare Real and in-silico FlyCodes
Reads the flycodes from p1875 db10 fasta file
# FC <- read.table(system.file("extdata/FC.tryptic", package = "NestLink"),
# header = TRUE)
FC <- getFC()## see ?NestLink and browseVignettes('NestLink') for documentation## loading from cacheFC$peptide <- (as.character(FC$peptide))
idx <- grep (PATTERN, FC$peptide)
FC$cond <- "FC"
FC$pim <- parentIonMass(FC$peptide)
FC <- FC[nchar(FC$peptide) >2, ]
FC$ssrc <- sapply(FC$peptide, ssrc)
FC$peptideLength <- nchar(as.character(FC$peptide))
FC <- FC[idx,]
head(FC) ## peptide ESP_Prediction cond pim ssrc peptideLength
## 120 GSAAAAADSWLTVR 0.75450 FC 1375.696 27.80445 14
## 121 GSAAAAATDWLTVR 0.76422 FC 1389.712 29.00445 14
## 122 GSAAAAATGWLTVR 0.65522 FC 1331.707 28.60445 14
## 123 GSAAAAATVWLR 0.65496 FC 1173.637 29.10445 12
## 124 GSAAAAAYEWLTVR 0.72754 FC 1465.743 33.10445 14
## 125 GSAAAADAAWQEGGR 0.53588 FC 1417.645 11.70445 154.1 AA design frequency
aa_pool_x7 <- c(rep('A', 18), rep('S', 6), rep('T', 12), rep('N', 1),
rep('Q', 1), rep('D', 11), rep('E', 11), rep('V', 12), rep('L', 2),
rep('F', 1), rep('Y', 4), rep('W', 1), rep('G', 8), rep('P', 12))
aa_pool_x7.post <- c(rep('WR', 20), rep('WLTVR', 20), rep('WQEGGR', 20),
rep('WLR', 20), rep('WQSR', 20))FC.pep_seq.freq <- table(unlist(strsplit(substr(FC$peptide, 3, 9), "")))
FC.freq <- round(FC.pep_seq.freq / sum(FC.pep_seq.freq) * 100, 3)
FC.pep_seq.freq.post <- table(unlist((substr(FC$peptide, 10, 100))))
FC.freq.post <- round(FC.pep_seq.freq.post / sum(FC.pep_seq.freq.post) * 100, 3)
WU.pep_seq.freq <- table(unlist(strsplit(substr(WU$pep_seq, 3, 9), "")))
WU.freq <- round(WU.pep_seq.freq / sum(WU.pep_seq.freq) * 100,3)
WU.pep_seq.freq.post <- table(unlist(substr(WU$pep_seq, 10, 100)))
WU.freq.post <- round(WU.pep_seq.freq.post / sum(WU.pep_seq.freq.post) *100,3)table(aa_pool_x7)## aa_pool_x7
## A D E F G L N P Q S T V W Y
## 18 11 11 1 8 2 1 12 1 6 12 12 1 4FC.freq##
## A D E F G L N P Q S T
## 16.027 13.452 10.097 0.872 6.918 1.658 0.618 10.662 1.142 5.637 10.663
## V W Y
## 16.533 1.245 4.478WU.freq##
## A D E F G L N P Q S T
## 15.434 16.983 11.939 0.532 7.077 1.165 0.784 11.813 1.484 5.748 10.532
## V W Y
## 12.645 0.773 3.090AAfreq1 <- data.frame(aa = table(aa_pool_x7))
names(AAfreq1) <- c('AA', 'freq')
AAfreq1 <-rbind(AAfreq1, df<-data.frame(AA=c('C','H','I','M','K','R'), freq=rep(0,6)))
AAfreq1$cond <- 'theoretical frequency'
AAfreq2 <- data.frame(aa=FC.freq)
names(AAfreq2) <- c('AA', 'freq')
AAfreq2$cond <- "db10.fasta"
AAfreq3 <- data.frame(aa=WU.freq)
names(AAfreq3) <- c('AA', 'freq')
AAfreq3$cond <- "WU"
#FC.all <- read.table(system.file("extdata/FC.tryptic", package = "NestLink"),
# header = TRUE)
FC.all <- getFC()## see ?NestLink and browseVignettes('NestLink') for documentation## loading from cachePATTERN.all <- "^GS[A-Z]{7}(WR|WLTVR|WQEGGR|WLR|WQSR)$"
FC.all <- as.character(FC.all$peptide)[grep(PATTERN.all, as.character(FC.all$peptide))]
FC.all.pep_seq.freq <- table(unlist(strsplit(substr(FC.all, 3, 9), "")))
FC.all.freq <- round(FC.all.pep_seq.freq / sum(FC.all.pep_seq.freq) * 100, 3)
FC.all.freq[20] <- 0
names(FC.all.freq) <- c(names(FC.all.freq)[1:19], "K")
AAfreq4 <- data.frame(AA=names(FC.all.freq), freq=as.numeric(FC.all.freq))
AAfreq4$cond <- "sequenced frequency"
AAfreq <- do.call('rbind', list(AAfreq1, AAfreq2, AAfreq3, AAfreq4))
AAfreq$freq <- as.numeric(AAfreq$freq)ggplot(data=AAfreq, aes(x=AA, y=freq, fill=cond)) +
geom_bar(stat="identity", position="dodge") +
labs(title = "amino acid frequency plot") +
theme_light()## Warning: Removed 5 rows containing missing values or values outside the scale range
## (`geom_bar()`).