M.tuberculosis 2017: Relative PE protein expression.

index.html

1 A query from Volker

The following is some text from an email from Volker:

“The last proteomics analysis found a total of 39 PE/PPE proteins in Mtb (Trypsin and Chemotrypsin digestions). Could you please find out using the RNAseq data of the wild-type Mtb transciptome which of the 161 PE/PPE proteins are actually expressed and at what relative level to the highest expressed PE/PPE protein? I thought that this would be very helpful in guiding our proteomics optimization.”

If I interpret this correctly, I should cross reference some of Keith’s results against the set of annotated PPE proteins and make a quick table providing a percentage of each against the largest value.

For the differential expression data, I used the file: mtb-invitro-manuscript/output/2016-11-23/invitro_3167c/invitro_3167c_Mtb-3167c_diff_expr_genes_all.csv

which I think contains the limma result for all samples/genes.

all_de <- read.table("limma_result.csv", header=TRUE, sep="\t")

pe_ids <- read.table("reference/annotated_pe_genes.txt")
colnames(pe_ids) <- "gene"
pe_annotations <- merge(pe_ids, mtb_annotations, by.x="gene", by.y="locus_tag")

## Error in fix.by(by.y, y): 'by' must specify a uniquely valid column

chosen_colnames <- c("rv_id", "chr", "start", "end", "width", "strand",
                     "source", "type", "score", "phase", "id", "gene_id",
                     "description", "function")
colnames(pe_annotations) <- chosen_colnames
rownames(pe_annotations) <- pe_annotations[["id"]]

tmp_de <- all_de
colnames(tmp_de) <- c("id", "logfc", "ave_expr", "adjp", "abbrev_id", "description", "gene_length", "type")
colnames(mtb_annotations) <- chosen_colnames
all_de <- merge(tmp_de, mtb_annotations, by="id")
colnames(all_de) <- c("id","logfc", "ave_expr", "adjp", "abbrev_id", "description", "gene_length", "type",
                      "rv_id", "start", "end", "width", "strand", "source", "type", "na", "anotherna",
                      "mtb_id", "another_abbrev_id", "another_description", "html_string")
all_de <- all_de[, c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15)]

pe_de <- merge(pe_annotations, all_de, by="id")
colnames(pe_de) <- c("id", "rv_id", "chromosome", "start", "end", "width", "strand",
                     "source", "type", "score", "phase", "gene_id", "description", "function",
                     "logfc", "ave_expr", "adjp", "another_gene", "another_description", "gene_length",
                     "another_type", "another_id", "another_start", "another_end", "another_width",
                     "another_strand", "another_source", "another_type")
pe_de <- pe_de[, c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17)]
## pe_de_annot <- merge(pe_de, mtb_annotations, by="locus_tag")
## The contrast is written as wt - delta, so add back the logFC I think
## Hopefully I did not reverse this in my head.
pe_de$ave_minus_log <- pe_de$ave_expr - pe_de$logfc

## Error in pe_de$ave_expr - pe_de$logfc: non-numeric argument to binary operator

all_de$ave_minus_log <- all_de$ave_expr - all_de$logfc
pe_de$pseudo_rpkm <- (pe_de$ave_minus_log / pe_de$width) * 1000
all_de$pseudo_rpkm <- (all_de$ave_minus_log / all_de$width) * 1000

maximum_expression <- max(pe_de[["pseudo_rpkm"]])

## Warning in max(pe_de[["pseudo_rpkm"]]): no non-missing arguments to max; returning -Inf

maximum_all <- max(all_de$pseudo_rpkm)

## Warning in max(all_de$pseudo_rpkm): no non-missing arguments to max; returning -Inf

pe_de[["exprs_vs_max"]] <- pe_de[["pseudo_rpkm"]] / maximum_expression
all_de$exprs_vs_max <- all_de$pseudo_rpkm / maximum_all

table_order <- order(pe_de[["exprs_vs_max"]], decreasing=TRUE)
pe_de <- pe_de[table_order, ]

knitr::kable(head(pe_de, n=50))

id rv_id chromosome start end width strand source type score phase gene_id description function logfc ave_expr adjp pseudo_rpkm exprs_vs_max — —— ———– —— —- —— ——- ——- —– —— —— ——– ———— ——— —— ——— —– ———— ————-

write.csv(x=pe_de, file="pe_relative_expression.csv")

library(hpgltools)
plot_histogram(pe_de$exprs_vs_max, bins=20)

## Error in FUN(X[[i]], ...): only defined on a data frame with all numeric variables

2 TODO from 20171129

We had a recent meeting with Volker and Dallas in which a few fun requests were made:

1. Take the existing rnaseq based metric of abundance PPE and cross reference against the potential metrics of the ms/ms data. (#PSM suggested by Yan, peptides from Nathan), perhaps rename last to something like pseudorpkm.
2. Provide a plot showing PSMs/peptides vs abundance along with rnaseq.
3. A flow diagram describing:
   • Starting with raw dia data through identifying PPE, acquire putative relative peptide abundance
   • Methodologies used in openSWATH analysis.

peptide_file <- "02_eval_result.Rmd"
tmp <- try(sm(loadme(filename=paste0(gsub(pattern="\\.Rmd", replace="", x=peptide_file), "-v", ver, ".rda.xz"))))

## The data structure which combines the chymotrypsin digested, found peptides, and the mtb annotation
## data is: chymo_translated.
head(chymo2_translated)
colnames(chymo2_translated)
head(all_de)
## The column I can use to look for correlations is 'rv_id'.

peptides_and_rna <- merge(x=chymo2_translated, y=all_de, by.x="id", by.y="locus_tag")
dim(peptides_and_rna)
## As per our conversation and emails with Yan and Nathan, the likely peptide columns are:
## coverage and num_psms, but I will try num_peptides too just to see.
##rownames(peptides_and_rna) <- peptides_and_rna[["rv_id"]]

rna_vs_peptides <- peptides_and_rna[, c("pseudo_rpkm", "num_peptides")]
rvp <- plot_linear_scatter(df=rna_vs_peptides, cormethod="pearson")
rvp$scatter
rvp$correlation
rvp$lm_rsq

## It looks like if you ignore psms with < 10, then the correlation is quite nice.
rna_vs_psms <- peptides_and_rna[, c("pseudo_rpkm", "num_psms")]
rvps <- plot_linear_scatter(df=rna_vs_psms, cormethod="pearson")
rvps$scatter
rvps$correlation
rvps$lm_rsq

rna_vs_coverage <- peptides_and_rna[, c("pseudo_rpkm", "coverage")]
rvc <- plot_linear_scatter(df=rna_vs_coverage, cormethod="spearman")
rvc$scatter
rvc$correlation
rvc$lm_rsq

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: hpgltools(v.2018.03)

loaded via a namespace (and not attached): SummarizedExperiment(v.1.8.1), pander(v.0.6.1), lattice(v.0.20-35), colorspace(v.1.3-2), htmltools(v.0.3.6), stats4(v.3.4.4), base64enc(v.0.1-3), rtracklayer(v.1.38.3), yaml(v.2.1.18), XML(v.3.98-1.10), rlang(v.0.2.0), pillar(v.1.2.1), withr(v.2.1.2), DBI(v.0.8), BiocParallel(v.1.12.0), BiocGenerics(v.0.24.0), matrixStats(v.0.53.1), GenomeInfoDbData(v.1.0.0), foreach(v.1.4.4), plyr(v.1.8.4), stringr(v.1.3.0), zlibbioc(v.1.24.0), Biostrings(v.2.46.0), munsell(v.0.4.3), commonmark(v.1.4), gtable(v.0.2.0), devtools(v.1.13.5), codetools(v.0.2-15), memoise(v.1.1.0), evaluate(v.0.10.1), labeling(v.0.3), Biobase(v.2.38.0), knitr(v.1.20), IRanges(v.2.12.0), GenomeInfoDb(v.1.14.0), parallel(v.3.4.4), highr(v.0.6), Rcpp(v.0.12.16), scales(v.0.5.0), backports(v.1.1.2), DelayedArray(v.0.4.1), S4Vectors(v.0.16.0), XVector(v.0.18.0), Rsamtools(v.1.30.0), ggplot2(v.2.2.1), RMySQL(v.0.10.14), digest(v.0.6.15), stringi(v.1.1.7), GenomicRanges(v.1.30.3), grid(v.3.4.4), rprojroot(v.1.3-2), tools(v.3.4.4), bitops(v.1.0-6), magrittr(v.1.5), lazyeval(v.0.2.1), RCurl(v.1.95-4.10), tibble(v.1.4.2), Matrix(v.1.2-12), data.table(v.1.10.4-3), xml2(v.1.2.0), rmarkdown(v.1.9), roxygen2(v.6.0.1), iterators(v.1.0.9), R6(v.2.2.2), GenomicAlignments(v.1.14.1) and compiler(v.3.4.4)

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 02_relative_pe_expression-v20170912.rda.xz

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