Evaluating some MS/MS peptide search metrics.

1 Evaluate Previous Result: version: 20170912

1.1 Reading in the existing result

As per Yan’s suggestion, I am using the existing results, found in the file: “2017_0811BrikenMultiConsensus.xlsx”

Unfortunately, this table is in a particularly strange format: it is a triply nested table of tables of tables in one worksheet in Excel. At one level, this is a really novel and interesting use of Excel, but on every other level it is sub-optimal. The primary reason it is sub-optimal is that the column headings end up being repeated once for every new sub-table (874 of the inner tables plus however many middle tables), which makes it nearly unparseable.

I am displeased I need to do this level of data wrangling.

The following function is my first attempt at handling this parsing puzzle.

For the interested, the structure is:

• all_protein_groups is a list of two elements: group_id => (“summary”, “data”)
• for every protein group discovered, one of these two element lists is created.
  • The summary slot gets the one-row summary for the protein group from the excel.
• The data slot is a list of groups: protein => (“summary”, “data”)
  • Once again, the summary is the single row from excel.
• Again, the data slot is a list of peptides: peptide => data.frame(peptides)
  • Each peptide data frame contains all of the inner-most rows describing the observed peptides.

parse_difficult_excel <- function(xlsx_file) {
  old_options <- options(java.parameters="-Xmx20G")
  message(paste0("Reading ", xlsx_file))
  result <- readxl::read_xlsx(path=xlsx_file, sheet=1)
  message("Starting parse.")
  group_data <- list()

  ## This is the toplevel heading for the protein groups.
  ## Dropping the first column because it is just used by excel to define the top-level table.
  outer_colnames <- colnames(result)
  outer_keepers <- !grepl(pattern="^X__", x=outer_colnames)
  outer_keepers[1] <- FALSE ## Get rid of the one that says 'Checked'
  outer_colnames <- outer_colnames[outer_keepers]

  for (r in 1:nrow(result)) {
    ## if ((1000 %% r) == 0) {
    ##   message(paste0("Completed ", r, " rows."))
    ## }
    row <- as.data.frame(result[r, ])
    row[, is.na(row)] <- ""

    ## When the first column is false, the next row will start a new protein group.
    if (row[, 1] == FALSE) {
      group_information <- row[outer_keepers]
      group_id <- group_information[, 1]
      group_data[[group_id]] <- list(
        "summary" = group_information,
        "data" = list())
      next
    }

    ## When the 2nd column is 'Checked', then this defines a new protein in the group.
    ## dropping the first 2 columns as they are used to define the middle-table.
    if (row[, 2] == "Checked") {
      group_keepers <- !grepl(pattern="X__", x=names(row))
      group_keepers <- names(row)[group_keepers]
      group_colnames <- as.character(row[, group_keepers])
      next
    }

    ## When the 2nd column is FALSE, then this defined a protein in the group.
    if (row[, 2] == FALSE) {
      group_row <- row[, group_keepers]
      colnames(group_row) <- group_colnames
      group_accession <- group_row[["Accession"]]
      protein_list <- list(
        "summary" = group_row,
        "data" = data.frame())
      group_data[[group_id]][["data"]][[group_accession]] <- protein_list
      next
    }

    ## When the 3rd column is 'Checked', then this starts a peptide definition
    if (row[, 3] == "Checked") {
      ## dropping the first 3 columns, they define the inner table.
      peptide_colnames <- row
      group_data[[group_id]][["data"]][[group_accession]][["data"]] <- data.frame()
      next
    }

    ## When the 3rd group is FALSE, then this adds a peptide.
    if (row[, 3] == FALSE) {
      peptide_data <- row
      colnames(peptide_data) <- peptide_colnames
      current <- group_data[[group_id]][["data"]][[group_accession]][["data"]]
      new <- rbind(current, peptide_data)
      group_data[[group_id]][["data"]][[group_accession]][["data"]] <- new
    }
  } ## End iterating over ever row of this unholy stupid data structure.

  return(group_data)
}

##xlsx_file <- "previous_results/2017_0811BrikenMultiConsensus.xlsx"
xlsx_file <- "previous_results/2017_1101BrikenMulti_PE_PPE_Complete.xlsx"
less_difficult <- parse_difficult_excel(xlsx_file)

## Reading previous_results/2017_1101BrikenMulti_PE_PPE_Complete.xlsx

## Starting parse.

The function above is not immediately helpful, as it only wrangles the crazy excel into a less-crazy data structure, but if we wish to graph the data or play with it, we need to extract the elements of interest. For the moment those are primarily the peptide fragments and whether they were found in the trypsin digestions, chymotrypsin digestions, or both.

1.2 Extract peptides

The following block of code handles that task. It just iterates over the new data structure and pulls out the peptides.

protein_df <- data.frame()
protein_names <- c()
for (group in names(less_difficult)) {
  group <- as.character(group)
  protein_group <- less_difficult[[group]][["data"]]
  protein_accessions <- names(protein_group)
  for (protein in protein_accessions) {
    protein_names <- c(protein_names, protein)
    protein_summary <- less_difficult[[group]][["data"]][[protein]][["summary"]]
    protein_df <- rbind(protein_df, protein_summary)
    peptide_df <- less_difficult[[group]][["data"]][[protein]][["data"]]
  }
}
head(protein_df)

##    Checked Protein FDR Confidence: Mascot Protein FDR Confidence: Sequest HT
## 1    FALSE                           High                               High
## 2    FALSE                           High                               High
## 3    FALSE                           High                               High
## 4    FALSE                           High                               High
## 5    FALSE                           High                               High
## 6    FALSE                           High                               High
##           Master Accession
## 1 Master Protein    P95029
## 2 Master Protein    I6X8D2
## 3 Master Protein    P9WGY7
## 4 Master Protein    P9WPE7
## 5 Master Protein    P9WGY9
## 6 Master Protein    P9WK07
##                                                                                                                                         Description
## 1                                                 3-oxoacyl-ACP synthase OS=Mycobacterium tuberculosis (strain ATCC 25618 / H37Rv) GN=fas PE=1 SV=3
## 2                                                  Polyketide synthase OS=Mycobacterium tuberculosis (strain ATCC 25618 / H37Rv) GN=pks13 PE=1 SV=1
## 3                             DNA-directed RNA polymerase subunit beta' OS=Mycobacterium tuberculosis (strain ATCC 25618 / H37Rv) GN=rpoC PE=1 SV=1
## 4                                                 60 kDa chaperonin 2 OS=Mycobacterium tuberculosis (strain ATCC 25618 / H37Rv) GN=groEL2 PE=1 SV=1
## 5                              DNA-directed RNA polymerase subunit beta OS=Mycobacterium tuberculosis (strain ATCC 25618 / H37Rv) GN=rpoB PE=1 SV=1
## 6 5-methyltetrahydropteroyltriglutamate--homocysteine methyltransferase OS=Mycobacterium tuberculosis (strain ATCC 25618 / H37Rv) GN=metE PE=1 SV=1
##   Exp. q-value: Mascot Exp. q-value: Sequest HT Coverage [%] Contaminant
## 1                0.000                    0.000           57       FALSE
## 2                0.000                    0.000           74       FALSE
## 3                0.000                    0.000           76       FALSE
## 4                0.000                    0.000           93       FALSE
## 5                0.000                    0.000           75       FALSE
## 6                0.000                    0.000           81       FALSE

This data frame now contains everything we are likely to want to visualize.

1.3 Cleaning the data frame

Unfortunately, this data structure also suffers from a couple weaknesses:

1. The protein IDs are not unique (because of the groups)
2. The column names from Excel are insane (punctuation, special characters, and really long)
3. Excel left the numeric columns of data as characters, not numbers!

Therefore, I will take a moment to fix these problems before trying any plots.

## Set some reasonable row names
rownames(protein_df) <- make.names(protein_names, unique=TRUE)
## Set the column names to something legible.
## Each line denotes my mental organizational unit
new_colnames <- c("mascot_fdr", "sequest_fdr", "master", "accession", "description",
                  "mascot_q", "sequest_q", "coverage", "contaminant",
                  "num_peptides", "num_psms", "num_unique_pep", "num_prot_groups", "num_aas",
                  "mw", "pi", "mascot_score", "sequest_score",
                  "num_peptides_mascot", "num_peptides_sequest",
                  "chymo_found_1", "chymo_found_2", "tryp_found_1", "tryp_found_2",
                  "modifications")
colnames(protein_df) <- new_colnames

## Error in names(x) <- value: 'names' attribute [25] must be the same length as the vector [11]

## Now make sure the columns which should be numeric, are numeric.
numeric_cols <- c("mascot_q", "sequest_q", "coverage", "num_peptides", "num_psms",
                  "num_unique_pep", "num_prot_groups", "num_aas", "mw", "pi",
                  "mascot_score", "sequest_score",
                  "num_peptides_mascot", "num_peptides_sequest")
for (col in numeric_cols) {
  protein_df[[col]] <- as.numeric(protein_df[[col]])
}

## Error in `[[<-.data.frame`(`*tmp*`, col, value = numeric(0)): replacement has 0 rows, data has 2633

## Finally, set a 'state' column describing the 4 columns of chymotrypsin/trypsin digestions.
protein_df[["state"]] <- "ambiguous"

## If something is found in one chymo trypsin, but in neither trypsin, set it to chymo_only.
chymo_hits <- (protein_df[["chymo_found_1"]] == "High" | protein_df[["chymo_found_2"]] == "High") &
  (protein_df[["tryp_found_1"]] != "High" & protein_df[["tryp_found_2"]] != "High")
protein_df[chymo_hits, "state"] <- "chymo_only"
## Conversely, if found in trypsin but not chymo, set it to tryp_only.
tryp_hits <- (protein_df[["tryp_found_1"]] == "High" | protein_df[["tryp_found_2"]] == "High") &
  (protein_df[["chymo_found_1"]] != "High" & protein_df[["chymo_found_2"]] != "High")
protein_df[tryp_hits, "state"] <- "tryp_only"

chymo_df <- protein_df[chymo_hits, ]
tryp_df <- protein_df[tryp_hits, ]

1.4 Make a couple plots

Now that the data has been wrangled, we can use the various R plotters to look at the data relatively quickly and easily.

library(ggplot2)
coverage_wrt_pi <- ggplot(protein_df, aes_string(x="pi", y="coverage", colour="state")) +
  geom_point(alpha=0.6, size=2) +
  geom_density_2d()

coverage_wrt_mw <- ggplot(protein_df, aes_string(x="mw", y="coverage", colour="state")) +
  geom_point(alpha=0.6, size=2) +
  geom_density_2d()

protein_df$logmascot <- log2(protein_df$mascot_score)

## Error in log2(protein_df$mascot_score): non-numeric argument to mathematical function

scores_by_state <- ggplot(protein_df, aes_string(x="state", y="logmascot")) +
  geom_boxplot(aes_string(fill="state"))

1.4.1 Coverage vs. pI

Yan suggested that peptides with low-pIs would be less covered, especially in the chymotrypsin samples.

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

## NULL

coverage_wrt_pi

## Error in FUN(X[[i]], ...): object 'coverage' not found

dev.off()

## png
##   2

## To me, it looks like there might be something to that idea.

1.4.2 Coverage vs. Molecular Weight

Maybe there is a trend of coverage and molecular weight by sample-type.

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

## NULL

coverage_wrt_mw

## Error in FUN(X[[i]], ...): object 'mw' not found

dev.off()

## png
##   2

## It looks to me that the trypsin is slightly better than chymotrypsin.
## It is not clear in the excel, but I presume this is in kDa.

1.4.3 Mascot scores

A quick boxplot showing mascot scores by what was observed.

scores_by_state

## Error in FUN(X[[i]], ...): object 'logmascot' not found

1.5 Evaluate what proteins were found by mascot for the chymo1 digestion

keepers <- protein_df[["chymo_found_1"]] == "High"
chymo1_df <- protein_df[keepers, ]
dim(chymo1_df)

## [1]  0 12

## unfortunately, this does not tell me exactly which spectra were observed to
## support these 274 protein.  However, I should be able to cross reference the
## proteins I find against these.
chymo_accessions <- chymo1_df[["accession"]]

get_xref <- function(accession) {
  request_url <- paste0("http://www.uniprot.org/uniprot/", accession, ".xml")
  data <- xml2::read_xml(request_url)
  data_list <- xml2::as_list(data)
  gene_names <- data_list[["entry"]][["gene"]]
  names <- reshape2::melt(gene_names)
  names <- names[["value"]]
  return(names)
}

acc_df <- data.frame()
for (acc in chymo_accessions) {
  hits <- try(get_xref(acc), silent=TRUE)
  if (class(hits) != "try-error") {
    for (h in hits) {
      row <- c(acc, h)
      acc_df <- rbind(acc_df, row)
    }
  }
}
colnames(acc_df) <- c("uniprot", "alias")

## Error in names(x) <- value: 'names' attribute [2] must be the same length as the vector [0]

head(acc_df)

## data frame with 0 columns and 0 rows

## This now contains only those IDs which were in the chymo digestion.

dim(chymo1_df)

## [1]  0 12

## So now lets get the set of all RvIDs
mtb_ids <- read.table("reference/mtb_ids.txt")
## And extract from the accessions list those which match
mtb_acc_df <- merge(mtb_ids, acc_df, by.x="V1", by.y="alias")

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

dim(mtb_acc_df)

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

## So we lost 7 (274 to 267), I can accept that.
colnames(mtb_acc_df) <- c("rv_id", "uniprot")

## Error in colnames(mtb_acc_df) <- c("rv_id", "uniprot"): object 'mtb_acc_df' not found

chymo_translated <- merge(chymo1_df, mtb_acc_df, by.x="accession", by.y="uniprot", all.x=TRUE)

## Error in as.data.frame(y): object 'mtb_acc_df' not found

## Ok, so now I can translate between the various IDs.

head(chymo1_df)

##  [1]                                    Checked
##  [3] Protein FDR Confidence: Mascot     Protein FDR Confidence: Sequest HT
##  [5] Master                             Accession
##  [7] Description                        Exp. q-value: Mascot
##  [9] Exp. q-value: Sequest HT           Coverage [%]
## [11] Contaminant                        state
## <0 rows> (or 0-length row.names)

head(mtb_annotations)

##    seqnames start  end width strand                 source type score phase       db_xref
## 1 NC_000962     1 1524  1524      + EMBL/GenBank/SwissProt gene    NA     1 GeneID:885041
## 2 NC_000962  2052 3260  1209      + EMBL/GenBank/SwissProt gene    NA     1 GeneID:887092
## 3 NC_000962  3280 4437  1158      + EMBL/GenBank/SwissProt gene    NA     1 GeneID:887089
## 4 NC_000962  4434 4997   564      + EMBL/GenBank/SwissProt gene    NA     1 GeneID:887088
## 5 NC_000962  5240 7267  2028      + EMBL/GenBank/SwissProt gene    NA     1 GeneID:887081
## 6 NC_000962  7302 9818  2517      + EMBL/GenBank/SwissProt gene    NA     1 GeneID:887105
##                                                       experiment gene locus_tag
## 1 DESCRIPTION:Mutation analysis, gene expression[PMID: 10375628] dnaA    Rv0001
## 2                                                           <NA> dnaN    Rv0002
## 3                                                           <NA> recF    Rv0003
## 4                                                           <NA> <NA>    Rv0004
## 5                    DESCRIPTION:Mutation analysis[PMID:8031045] gyrB    Rv0005
## 6                    DESCRIPTION:Mutation analysis[PMID:8031045] gyrA    Rv0006
##   gene_synonym pseudogene
## 1         <NA>       <NA>
## 2         <NA>       <NA>
## 3         <NA>       <NA>
## 4         <NA>       <NA>
## 5         <NA>       <NA>
## 6         <NA>       <NA>

chymo_t2 <- merge(chymo1_df, mtb_annotations, by.x=)


chymo2_df <- chymo1_df
chymo2_df[["gn"]] <- gsub(x=chymo2_df[["description"]], pattern="^.* GN=(.+) PE=.*$", replacement="\\1")
chymo2_translated <- merge(chymo2_df, mtb_annotations, by.x="gn", by.y="gene")

pander::pander(sessionInfo())

R version 3.4.3 (2017-11-30)

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

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

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

other attached packages: ggplot2(v.2.2.1) and hpgltools(v.2017.10)

loaded via a namespace (and not attached): Rcpp(v.0.12.15), cellranger(v.1.1.0), pillar(v.1.2.1), compiler(v.3.4.3), plyr(v.1.8.4), base64enc(v.0.1-3), iterators(v.1.0.9), tools(v.3.4.3), digest(v.0.6.15), memoise(v.1.1.0), evaluate(v.0.10.1), tibble(v.1.4.2), gtable(v.0.2.0), rlang(v.0.2.0), foreach(v.1.4.4), commonmark(v.1.4), yaml(v.2.1.16), parallel(v.3.4.3), withr(v.2.1.1), stringr(v.1.3.0), knitr(v.1.20), roxygen2(v.6.0.1), xml2(v.1.2.0), devtools(v.1.13.5), rprojroot(v.1.3-2), grid(v.3.4.3), data.table(v.1.10.4-3), Biobase(v.2.38.0), R6(v.2.2.2), readxl(v.1.0.0), rmarkdown(v.1.8), pander(v.0.6.1), magrittr(v.1.5), scales(v.0.5.0), backports(v.1.1.2), codetools(v.0.2-15), htmltools(v.0.3.6), BiocGenerics(v.0.24.0), colorspace(v.1.3-2), stringi(v.1.1.6), lazyeval(v.0.2.1) and munsell(v.0.4.3)

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 04051c4dfae6fbfe706f11b1ea8cb454e1cafa2c

## R> packrat::restore()

## This is hpgltools commit: Sun Feb 25 23:26:32 2018 -0500: 04051c4dfae6fbfe706f11b1ea8cb454e1cafa2c

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

## Saving to 02_eval_result-v20170912.rda.xz

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