Testing out Proteomics analyses.

1 Examples version: 20170905

I am trying out the examples from ‘Using R and Bioconductor for Proteomics Data Analysis.’ Most of the following is taken verbatim from that.

2 Setting up

Make sure I have the various required packages.

library("hpgltools")
tt <- sm(library("RforProteomics"))
tt <- sm(library("mzR"))
tt <- sm(library("msdata"))

3 Query 1 from Volker

What is the size distribution of all Mtb proteins containing PP and PPE? More specifically, how many of them are smaller than 30 kDa?

Of these proteins, how many have trypsin digestion sites?

3.1 Method

I have a database of the translated CDS sequences from MTb. Lets query it.

pattern <- "PE"
cds <- "mtb_cds.fasta"
mtb_cds <- Biostrings::readAAStringSet(filepath=cds)

pattern <- BiocGenerics::sapply(pattern, toString)
mtb_string <- BiocGenerics::sapply(mtb_cds, toString)
## Perform the search, this returns a list of hits by position and name of pattern.
mtb_pp <- AhoCorasickTrie::AhoCorasickSearch(keywords=pattern,
                                             text=mtb_string,
                                             alphabet="ascii")
## Melt this into a somewhat usable data frame.
pp_df <- reshape2::melt(data=mtb_pp, level=1)
## Drop the columns I don't want
pp_df <- pp_df[, c("L1", "value")]
## Drop the rows which have useless information
drop <- pp_df[["value"]] == "PP"
pp_df <- pp_df[!drop, ]
## Rename the columns to something useful
colnames(pp_df) <- c("ID", "position")
## Look at the first 20.
head(pp_df, n=20)

##        ID position
## 1  Rv0001       PE
## 2  Rv0001      333
## 3  Rv0002       PE
## 4  Rv0002      133
## 5  Rv0003       PE
## 6  Rv0003      119
## 7  Rv0003       PE
## 8  Rv0003      210
## 9  Rv0003       PE
## 10 Rv0003      221
## 11 Rv0005       PE
## 12 Rv0005      578
## 13 Rv0006       PE
## 14 Rv0006       42
## 15 Rv0006       PE
## 16 Rv0006      472
## 17 Rv0006       PE
## 18 Rv0006      585
## 19 Rv0006       PE
## 20 Rv0006      604

## Lets check that this is correct, Rv0007 has 43, 44, 61, 62, so that should be easy to check.
## Here are the first 80 AA of Rv0007        v  <- pos. 43     v <- pos. 61 v <- 74
## VTAPNEPGALSKGDGPNADGLVDRGGAHRAATGPGRIPDAGDPPPWQRAATRQSQAGHRQPPPVSHPEGRPTNPPAAADA
## Yay!

## Ok, so the actual question was, how many PP genes are there.
length(unique(pp_df[["ID"]]))

## [1] 2292

pp_hits <- unique(pp_df$ID)

3.2 Repeat for PPE

## What about PPE?
pattern <- "PPE"
pattern <- BiocGenerics::sapply(pattern, toString)
mtb_ppe <- AhoCorasickTrie::AhoCorasickSearch(keywords=pattern,
                                             text=mtb_string,
                                             alphabet="ascii")
ppe_df <- reshape2::melt(data=mtb_ppe, level=1)
ppe_df <- ppe_df[, c("L1", "value")]
drop <- ppe_df[["value"]] == "PPE"
ppe_df <- ppe_df[!drop, ]
colnames(ppe_df) <- c("ID", "position")
head(ppe_df, n=20)

##         ID position
## 2   Rv0001      332
## 4  Rv0014c      281
## 6  Rv0018c      390
## 8  Rv0020c      203
## 10 Rv0020c      341
## 12 Rv0043c      217
## 14 Rv0046c      348
## 16  Rv0050       46
## 18  Rv0050      630
## 20  Rv0073      225
## 22  Rv0096        4
## 24  Rv0111      452
## 26  Rv0148      206
## 28  Rv0185       47
## 30  Rv0186      531
## 32 Rv0198c      624
## 34 Rv0227c      417
## 36 Rv0256c       10
## 38 Rv0275c      224
## 40  Rv0280        8

length(unique(ppe_df[["ID"]]))

## [1] 280

3.3 Find overlap of trypsin products with these

tryp_hits <- cleaver::cleavageRanges(mtb_cds, enzym="trypsin")
tryp_df <- as.data.frame(tryp_hits)
tryp_names <- unique(tryp_df$group_name)

3.4 And the likely molecular weights

mtb_ppe <- mtb_cds[ppe_df$ID, ]
tryp_products <- cleaver::cleave(mtb_ppe, enzym="trypsin")
tryp_products_df <- as.data.frame(tryp_products)
tryp_products_df$mw <- Peptides::mw(tryp_products_df$value)
dim(tryp_products_df)

## [1] 11631     4

## What products are there from 10,000 Daltons to 30,000?
hist(tryp_products_df$mw, breaks=1000, xlim=c(10000, 30000), ylim=c(0, 15))

## Call these 'intermediate'
intermediate <- tryp_products_df
## Drop anything bigger than 30,000 Daltons
keepers <- intermediate$mw <= 30000
intermediate <- intermediate[keepers, ]
## Drop anything smaller than 10,000 Daltons
keepers <- intermediate$mw > 10000
intermediate <- intermediate[keepers, ]
unique(intermediate$group_name)

##  [1] "Rv0096"  "Rv0256c" "Rv0280"  "Rv0286"  "Rv0355c" "Rv0442c" "Rv0453"  "Rv0878c"
##  [9] "Rv0915c" "Rv1039c" "Rv1040c" "Rv1135c" "Rv1196"  "Rv1233c" "Rv1361c" "Rv1387" 
## [17] "Rv1548c" "Rv1753c" "Rv1789"  "Rv1808"  "Rv1917c" "Rv1918c" "Rv2264c" "Rv2352c"
## [25] "Rv2444c" "Rv2768c" "Rv2892c" "Rv3018c" "Rv3136"  "Rv3159c" "Rv3343c" "Rv3347c"
## [33] "Rv3350c" "Rv3478"  "Rv3494c" "Rv3533c" "Rv3558"  "Rv3873"

3.5 Look for overlap of high-specificity chymotrypsin cleavages

chtryp_hits <- cleaver::cleavageRanges(mtb_cds, enzym="chymotrypsin-high")
chtryp_df <- as.data.frame(chtryp_hits)
chtryp_names <- unique(chtryp_df$group_name)

3.6 The likely molecular weights for high-specificity chymotrypsin

In this case we look for C->[FYW] not before P.

mtb_ppe <- mtb_cds[ppe_df$ID, ]
tryp_products <- cleaver::cleave(mtb_ppe, enzym="chymotrypsin-high")
tryp_products_df <- as.data.frame(tryp_products)
tryp_products_df$mw <- Peptides::mw(tryp_products_df$value)
dim(tryp_products_df)

## [1] 10629     4

## What products are there from 10,000 Daltons to 30,000?
hist(tryp_products_df$mw, breaks=1000, xlim=c(10000, 30000), ylim=c(0, 15))

## Call these 'intermediate'
intermediate <- tryp_products_df
## Drop anything bigger than 30,000 Daltons
keepers <- intermediate$mw <= 30000
intermediate <- intermediate[keepers, ]
## Drop anything smaller than 10,000 Daltons
keepers <- intermediate$mw > 10000
intermediate <- intermediate[keepers, ]
unique(intermediate$group_name)

##  [1] "Rv0018c" "Rv0227c" "Rv0336"  "Rv0339c" "Rv0425c" "Rv0494"  "Rv0515"  "Rv0668" 
##  [9] "Rv0904c" "Rv1020"  "Rv1186c" "Rv1409"  "Rv1796"  "Rv1887"  "Rv1917c" "Rv2139" 
## [17] "Rv2207"  "Rv2264c" "Rv2352c" "Rv2444c" "Rv2578c" "Rv2847c" "Rv2921c" "Rv2931" 
## [25] "Rv2934"  "Rv2941"  "Rv3144c" "Rv3209"  "Rv3245c" "Rv3296"  "Rv3333c" "Rv3663c"
## [33] "Rv3721c" "Rv3723"  "Rv3763"  "Rv3823c" "Rv3835"  "Rv3873"  "Rv3886c" "Rv3903c"
## [41] "Rv3910"

3.7 Finally, the low-specificity chymotrypsin digestion products

cltryp_hits <- cleaver::cleavageRanges(mtb_cds, enzym="chymotrypsin-low")
cltryp_df <- as.data.frame(cltryp_hits)
cltryp_names <- unique(cltryp_df$group_name)

## All of the pp containing peptides have trypsin digestion products.
sum(pp_hits %in% tryp_names)

## [1] 2292

3.8 And Chymotrypsin-low size distribution

And here we look for C->[FYWML] not before P.

mtb_ppe <- mtb_cds[ppe_df$ID, ]
tryp_products <- cleaver::cleave(mtb_ppe, enzym="chymotrypsin-low")
tryp_products_df <- as.data.frame(tryp_products)
tryp_products_df$mw <- Peptides::mw(tryp_products_df$value)
dim(tryp_products_df)

## [1] 25847     4

## What products are there from 10,000 Daltons to 30,000?
hist(tryp_products_df$mw, breaks=1000, xlim=c(1000, 7000), ylim=c(0, 15))

## Call these 'intermediate'
intermediate <- tryp_products_df
## Drop anything bigger than 30,000 Daltons
keepers <- intermediate$mw <= 30000
intermediate <- intermediate[keepers, ]
## Drop anything smaller than 10,000 Daltons
keepers <- intermediate$mw > 10000
intermediate <- intermediate[keepers, ]
unique(intermediate$group_name)

## character(0)