M.tuberculosis 2017: Testing out Proteomics analyses.

1 Simulated digestion of PE proteins: 20171205

I did a little reading of the paper: “Structural basis of the PE-PPE protein interaction in Mycobacterium tuberculosis” and quickly realized that when Dr. Briken was asking for the set of PE containing genes, he was not in fact asking for the set of all genes with ‘PE’ in them, but in fact a set of annotated genes.

I therefore extracted those gene IDs from the gene annotations and left them in the file ‘annotated_pe_genes.txt’. In my previous test, I did a string search of all the peptides for ones containing ‘PE’ or ‘PPE’. This time I will just use the annotated peptides.

Let us therefore read peptide sequences into memory and extract only these 161 peptides. Then perform some digestions and see what happens.

1.1 Reading the peptides

This is performed via the R Biostrings package.

## Start out reading in all peptides.
mtb_peptides <- Biostrings::readAAStringSet(filepath="reference/mtb_cds.fasta")
summary(mtb_peptides)

##      Length       Class        Mode 
##        4008 AAStringSet          S4

## Get the IDs of the PE proteins.
pe_ids <- read.table("reference/annotated_pe_genes.txt", header=FALSE)[["V1"]]
## Subset all peptides for just these.
mtb_pe <- mtb_peptides[pe_ids, ]
summary(mtb_pe)

##      Length       Class        Mode 
##         161 AAStringSet          S4

pe_fasta <- Biostrings::writeXStringSet(x=mtb_pe, filepath="reference/mtb_pe.fasta")

1.2 Look for idealized digestion products

The cleaver package provides easy sample digestions of arbitrary peptide sequences. In addition, the BRAIN package provides center-mass calculations for the resulting peptide fragments. These may therefore be plotted as example spectra. I will copy the plotting implementation from the cleaver documentation in the following code block.

At this point, it is useful to note some information provide by Yan: “I’d like to be able to look at the sizes of expected trypsin and chymotrypsin digests and compare the unidentified spectra in the raw data to see if they match up. That will give us a rough idea on whether our problem is more in digestion or MSMS data itself. … Would you be able to do in silico digestion of all PE proteins by trypsin, cymotrypsin, Glu-C and Asp-N? Assuming we use the same LCMS method, our best bet is to identify peptides in the range of 700-3000 Da…”

1.3 Perform test digestions

The above function should probably be changed a little to provide some better visualization, but for the moment let us just see what happens.

1.3.1 Trypsin

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

cleavage_histogram <- function(pep_sequences, enzyme="trypsin",
                               start=600, end=1500, color="black") {
  products <- cleaver::cleave(pep_sequences, enzym=enzyme)
  prod_df <- as.data.frame(products)
  prod_df <- as.tbl(prod_df[, c("group_name", "value")])
  colnames(prod_df) <- c("group_name", "sequence")

  gather_masses <- function(sequence) {
    atoms <- try(BRAIN::getAtomsFromSeq(sequence), silent=TRUE)
    if (class(atoms) != "try-error") {
      d <- BRAIN::useBRAIN(atoms)
      ret <- round(d[["avgMass"]])
    } else {
      ret <- 0
    }
    return(ret)
  }

  new_df <- prod_df %>% rowwise() %>% mutate(mass=gather_masses(sequence))

  plot <- ggplot2::ggplot(data=new_df, ggplot2::aes_string(x="mass")) +
    ggplot2::geom_histogram(binwidth=1, colour=color) +
    ggplot2::scale_x_continuous(limits=c(start, end))

  retlist <- list(
    "plot" = plot,
    "masses" = new_df)
  return(retlist)
}

tt <- pp(file="images/tryp_intensity_600_1500.png")

## Going to write the image to: images/tryp_intensity_600_1500.png when dev.off() is called.

tryp_low <- sm(plot_cleaved(pep_sequences=mtb_pe, enzyme="trypsin",
                            start=600, end=1500))
dev.off()

## X11cairo 
##        2

tt <- pp(file="images/tryp_intensity_1500_3000.png")

## Going to write the image to: images/tryp_intensity_1500_3000.png when dev.off() is called.

tryp_mid <- sm(plot_cleaved(pep_sequences=mtb_pe, enzyme="trypsin",
                            start=1500, end=3000))
dev.off()

## X11cairo 
##        2

tt <- pp(file="images/tryp_intensity_1500_3000.png")

## Going to write the image to: images/tryp_intensity_1500_3000.png when dev.off() is called.

tryp_high <- sm(plot_cleaved(pep_sequences=mtb_pe, enzyme="trypsin",
                             start=3000, end=5000))
dev.off()

## X11cairo 
##        2

1.3.2 Again as histograms rather than approximate intensity

tt <- pp(file="images/pe_tryp_600_1500.png")

## Going to write the image to: images/pe_tryp_600_1500.png when dev.off() is called.

cleavage_histogram(pep_sequences=mtb_pe, color="red")$plot

## Warning: Removed 1492 rows containing non-finite values (stat_bin).

dev.off()

## X11cairo 
##        2

tt <- pp(file="images/pe_tryp_1500_3000.png")

## Going to write the image to: images/pe_tryp_1500_3000.png when dev.off() is called.

cleavage_histogram(pep_sequences=mtb_pe, start=1500, end=3000, color="red")$plot

## Warning: Removed 1461 rows containing non-finite values (stat_bin).

dev.off()

## X11cairo 
##        2

tt <- pp(file="images/pe_tryp_3000_5000.png")

## Going to write the image to: images/pe_tryp_3000_5000.png when dev.off() is called.

cleavage_histogram(pep_sequences=mtb_pe, start=3000, end=5000, color="red")$plot

## Warning: Removed 1618 rows containing non-finite values (stat_bin).

dev.off()

## X11cairo 
##        2

1.3.3 ChymoTrypsin high specificity

tt <- pp(file="images/chymo_intensity_600_1500.png")

## Going to write the image to: images/chymo_intensity_600_1500.png when dev.off() is called.

chy_low <- sm(plot_cleaved(pep_sequences=mtb_pe, enzyme="chymotrypsin-high",))
dev.off()

## X11cairo 
##        2

tt <- pp(file="images/chymo_intensity_1500_3000.png")

## Going to write the image to: images/chymo_intensity_1500_3000.png when dev.off() is called.

chy_mid <-  sm(plot_cleaved(pep_sequences=mtb_pe, enzyme="chymotrypsin-high",
                            start=2000, end=3500))
dev.off()

## X11cairo 
##        2

tt <- pp(file="images/chymo_intensity_3000_5000.png")

## Going to write the image to: images/chymo_intensity_3000_5000.png when dev.off() is called.

chy_high <- sm(plot_cleaved(pep_sequences=mtb_pe, enzyme="chymotrypsin-high",
                            start=3500, end=5000))
dev.off()

## X11cairo 
##        2

tt <- pp(file="images/pe_chymo_600_1500.png")

## Going to write the image to: images/pe_chymo_600_1500.png when dev.off() is called.

cleavage_histogram(pep_sequences=mtb_pe, enzyme="chymotrypsin-high", color="blue")$plot

## Warning: Removed 2945 rows containing non-finite values (stat_bin).

dev.off()

## X11cairo 
##        2

tt <- pp(file="images/pe_chymo_1500_3000.png")

## Going to write the image to: images/pe_chymo_1500_3000.png when dev.off() is called.

cleavage_histogram(pep_sequences=mtb_pe, enzyme="chymotrypsin-high", start=1500, end=3000, color="blue")$plot

## Warning: Removed 3785 rows containing non-finite values (stat_bin).

dev.off()

## X11cairo 
##        2

tt <- pp(file="images/pe_chymo_3000_5000.png")

## Going to write the image to: images/pe_chymo_3000_5000.png when dev.off() is called.

cleavage_histogram(pep_sequences=mtb_pe, enzyme="chymotrypsin-high", start=3000, end=5000, color="blue")$plot

## Warning: Removed 4267 rows containing non-finite values (stat_bin).

dev.off()

## X11cairo 
##        2

1.3.4 asp-n endopeptidase

aspn_low <- sm(plot_cleaved(pep_sequences=mtb_pe, enzyme="asp-n endopeptidase"))

aspn_mid <- sm(plot_cleaved(pep_sequences=mtb_pe, enzyme="asp-n endopeptidase",
                             start=2000, end=3500))

aspn_high <- sm(plot_cleaved(pep_sequences=mtb_pe, enzyme="asp-n endopeptidase",
                             start=3500, end=5000))

cleavage_histogram(pep_sequences=mtb_pe, enzyme="asp-n endopeptidase", color="darkgreen")$plot

## Warning: Removed 2105 rows containing non-finite values (stat_bin).

dev.off()

## X11cairo 
##        2

tt <- pp(file="images/pe_aspn_1500_3000.png")

## Going to write the image to: images/pe_aspn_1500_3000.png when dev.off() is called.

cleavage_histogram(pep_sequences=mtb_pe, enzyme="asp-n endopeptidase", start=1500, end=3000, color="darkgreen")$plot

## Warning: Removed 2003 rows containing non-finite values (stat_bin).

dev.off()

## X11cairo 
##        2

tt <- pp(file="images/pe_aspn_3000_5000.png")

## Going to write the image to: images/pe_aspn_3000_5000.png when dev.off() is called.

cleavage_histogram(pep_sequences=mtb_pe, enzyme="asp-n endopeptidase", start=3000, end=5000, color="darkgreen")$plot

## Warning: Removed 2259 rows containing non-finite values (stat_bin).

dev.off()

## X11cairo 
##        2

tt <- pp(file="images/pe_glu_1500_3000.png")

## Going to write the image to: images/pe_glu_1500_3000.png when dev.off() is called.

cleavage_histogram(pep_sequences=mtb_pe, enzyme="glutamyl endopeptidase", color="purple")$plot

## Warning: Removed 1333 rows containing non-finite values (stat_bin).

dev.off()

## X11cairo 
##        2

tt <- pp(file="images/pe_glu_1500_3000.png")

## Going to write the image to: images/pe_glu_1500_3000.png when dev.off() is called.

cleavage_histogram(pep_sequences=mtb_pe, enzyme="glutamyl endopeptidase", start=1500, end=3000, color="darkgreen")$plot

## Warning: Removed 1388 rows containing non-finite values (stat_bin).

dev.off()

## X11cairo 
##        2

tt <- pp(file="images/pe_glu_3000_5000.png")

## Going to write the image to: images/pe_glu_3000_5000.png when dev.off() is called.

cleavage_histogram(pep_sequences=mtb_pe, enzyme="glutamyl endopeptidase", start=3000, end=5000, color="darkgreen")$plot

## Warning: Removed 1479 rows containing non-finite values (stat_bin).

dev.off()

## X11cairo 
##        2

1.3.5 glu-c endopeptidase

glu_low <- sm(plot_cleaved(pep_sequences=mtb_pe, enzyme="glutamyl endopeptidase",
                           start=700, end=2000))

glu_high <- sm(plot_cleaved(pep_sequences=mtb_pe, enzyme="glutamyl endopeptidase",
                            start=2000, end=3000))

1.4 Compare simulated PE vs all

all_tryp <- cleavage_histogram(pep_sequences=mtb_pe)$masses
all_tryp$enzyme <- "tryp"
all_chy <- cleavage_histogram(pep_sequences=mtb_pe, enzyme="chymotrypsin-high")$masses
all_chy$enzyme <- "chymo"
all_aspn <- cleavage_histogram(pep_sequences=mtb_pe, enzyme="asp-n endopeptidase")$masses
all_aspn$enzyme <- "aspn"
all_glu <- cleavage_histogram(pep_sequences=mtb_pe, enzyme="glutamyl endopeptidase")$masses
all_glu$enzyme <- "glu"

all_df <- all_tryp[, c("enzyme", "mass")]
all_df <- rbind(all_df, all_chy[, c("enzyme","mass")])
all_df <- rbind(all_df, all_aspn[, c("enzyme", "mass")])
all_df <- rbind(all_df, all_glu[, c("enzyme", "mass")])
all_df <- as.data.frame(all_df)
all_df$enzyme <- as.factor(all_df$enzyme)
all_df$l2mass <- log2(all_df$mass + 1)

comparison <- ggplot(data=all_df, aes_string(x="mass", group="enzyme", colour="enzyme")) +
  geom_density() +
  scale_x_continuous(limits=c(0,10000))

## Error in ggplot(data = all_df, aes_string(x = "mass", group = "enzyme", : could not find function "ggplot"

pp(file="images/density_enzymes.png")

## Going to write the image to: images/density_enzymes.png when dev.off() is called.

## NULL

comparison

## Error in eval(expr, envir, enclos): object 'comparison' not found

dev.off()

## X11cairo 
##        2

pe_distribution <- tryp_low[["sizes"]]
all_trypsin_distribution <- plot_cleaved(pep_sequences=mtb_peptides, enzyme="trypsin", start=2000, end=2010)

all_chymo_distribution <- plot_cleaved(pep_sequences=mtb_peptides, enzyme="chymotrypsin-high", start=2000, end=2010)

pe_names <- unique(pe_distribution[["cds"]])

distribution <- all_trypsin_distribution
plot_spectra <- function(distribution, subset=pe_names) {
  all_sizes <- distribution[["sizes"]]
  all_sizes[["subset"]] <- 0
  subset_found <- all_sizes[["cds"]] %in% subset
  all_sizes[subset_found, "subset"] <- 1
  all_sizes[["average_mass"]] <- as.numeric(all_sizes[["average_mass"]])
  all_sizes[["highest_likelihood"]] <- as.numeric(all_sizes[["highest_likelihood"]])

  a_plot <- ggplot2::ggplot(data=all_sizes,
                            ggplot2::aes_string(x="average_mass",
                                                y="highest_likelihood",
                                                colour="as.factor(subset)")) +
    ggplot2::geom_point(alpha=0.3, size=2, stat="identity") +
    ggplot2::scale_color_manual(name="as.factor(subset)",
                                values=c("0"="darkblue", "1"="darkred"))
  return(a_plot)
}

trypsin_dist_plot <- plot_spectra(all_trypsin_distribution)
pp(file="images/tryp_dist_plot.png")

## Going to write the image to: images/tryp_dist_plot.png when dev.off() is called.

## NULL

trypsin_dist_plot
dev.off()

## X11cairo 
##        2

chymotrypsin_dist_plot <- plot_spectra(all_chymo_distribution)

pp(file="images/chymo_dist_plot.png")

## Going to write the image to: images/chymo_dist_plot.png when dev.off() is called.

## NULL

chymotrypsin_dist_plot
dev.off()

## X11cairo 
##        2

## I think this demonstrates nicely how much more efficient chymotrypsin is for pe proteins.

pander::pander(sessionInfo())

R version 3.4.4 (2018-03-15)

**Platform:** x86_64-pc-linux-gnu (64-bit)

locale: LC_CTYPE=en_US.utf8, LC_NUMERIC=C, LC_TIME=en_US.utf8, LC_COLLATE=en_US.utf8, LC_MONETARY=en_US.utf8, LC_MESSAGES=en_US.utf8, LC_PAPER=en_US.utf8, LC_NAME=C, LC_ADDRESS=C, LC_TELEPHONE=C, LC_MEASUREMENT=en_US.utf8 and LC_IDENTIFICATION=C

attached base packages: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: bindrcpp(v.0.2), dplyr(v.0.7.4) and hpgltools(v.2018.03)

loaded via a namespace (and not attached): Rcpp(v.0.12.16), bindr(v.0.1.1), pillar(v.1.2.1), compiler(v.3.4.4), plyr(v.1.8.4), XVector(v.0.18.0), base64enc(v.0.1-3), iterators(v.1.0.9), tools(v.3.4.4), zlibbioc(v.1.24.0), digest(v.0.6.15), evaluate(v.0.10.1), lattice(v.0.20-35), memoise(v.1.1.0), tibble(v.1.4.2), gtable(v.0.2.0), pkgconfig(v.2.0.1), rlang(v.0.2.0), foreach(v.1.4.4), commonmark(v.1.4), yaml(v.2.1.18), parallel(v.3.4.4), withr(v.2.1.2), stringr(v.1.3.0), knitr(v.1.20), roxygen2(v.6.0.1), xml2(v.1.2.0), Biostrings(v.2.46.0), S4Vectors(v.0.16.0), devtools(v.1.13.5), IRanges(v.2.12.0), rprojroot(v.1.3-2), stats4(v.3.4.4), grid(v.3.4.4), glue(v.1.2.0), data.table(v.1.10.4-3), Biobase(v.2.38.0), PolynomF(v.1.0-1), R6(v.2.2.2), rmarkdown(v.1.9), pander(v.0.6.1), ggplot2(v.2.2.1), cleaver(v.1.16.0), magrittr(v.1.5), backports(v.1.1.2), htmltools(v.0.3.6), scales(v.0.5.0), codetools(v.0.2-15), BiocGenerics(v.0.24.0), assertthat(v.0.2.0), colorspace(v.1.3-2), labeling(v.0.3), stringi(v.1.1.7), lazyeval(v.0.2.1), munsell(v.0.4.3) and BRAIN(v.1.24.0)

message(paste0("This is hpgltools commit: ", get_git_commit()))

## If you wish to reproduce this exact build of hpgltools, invoke the following:

## > git clone http://github.com/abelew/hpgltools.git

## > git reset 2a0661d6e37f8a3d8831eb3bbd6347c0d9c4b3b7

## R> packrat::restore()

## This is hpgltools commit: Thu Mar 29 16:59:07 2018 -0400: 2a0661d6e37f8a3d8831eb3bbd6347c0d9c4b3b7

this_save <- paste0(gsub(pattern="\\.Rmd", replace="", x=rmd_file), "-v", ver, ".rda.xz")
message(paste0("Saving to ", this_save))

## Saving to 03_simulated_pe_digestions-v20171205.rda.xz

tmp <- sm(saveme(filename=this_save))