M.tuberculosis 2017: Evaluating some MS/MS peptide search metrics.

1 Evaluate Previous Result: version: 20171205

1.1 TODO

1. Replace rnaseq data used with a more appropriate set. (Sauton medium)
2. Do a linear relationship for only the P/PE proteins vs. rpkm.
3. Learn how to quantify peptide abundance properly.

1.2 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.

## I pulled this function into hpgltools and cleaned it up I think.
## I also gave it a less stupid name: read_thermo_xlsx()
parse_difficult_excel <- function(xlsx_file, test_row=92150) {
  old_options <- options(java.parameters="-Xmx20G")
  message(paste0("Reading ", xlsx_file))
  result <- readxl::read_xlsx(path=xlsx_file, sheet=1, col_names=FALSE)
  message("Starting parse.")
  group_data <- list()
  bar <- utils::txtProgressBar(style=3)
  for (r in 1:nrow(result)) {
    row <- as.data.frame(result[r, ])
    row[, is.na(row)] <- ""
    pct_done <- r / nrow(result)
    setTxtProgressBar(bar, pct_done)
    ## The following 3 stanzas handle the creation of the levels of our data structure
    ## The first defines the protein group
    if (row[, 1] == "Checked") {
      group_colnames <- as.character(row)
      group_keepers <- !grepl(pattern="^$", x=group_colnames)
      group_keepers[1] <- FALSE
      group_colnames <- group_colnames[group_keepers]
      next
    }
    ## When the 2nd column is 'Checked', then this defines a new protein in the group.
    if (row[, 2] == "Checked") {
      protein_colnames <- as.character(row)
      protein_keepers <- !grepl(pattern="^$", x=protein_colnames)
      protein_keepers[2] <- FALSE
      protein_colnames <- protein_colnames[protein_keepers]
      next
    }
    ## When the 3rd column is 'Checked', then this starts a peptide definition
    if (row[, 3] == "Checked") {
      peptide_colnames <- as.character(row)
      peptide_keepers <- !grepl(pattern="^$", x=peptide_colnames)
      peptide_keepers[3] <- FALSE
      peptide_colnames <- peptide_colnames[peptide_keepers]
      next
    }
    ## Once the column names for the data are defined, we consider how to
    ## Fill in the actual data, the protein group is probably the least interesting.
    if (row[, 1] == FALSE | row[, 1] == TRUE) {
      group_information <- row[group_keepers]
      colnames(group_information) <- group_colnames
      group_information[["ID"]] <- sub(pattern="^.* GN=(\\w+) .*$",
                                       replacement="\\1",
                                       x=group_information[["Group Description"]])
      group_accession <- group_information[["Protein Group ID"]]
      group_list <- list(
        "summary" = group_information,
        "data" = list())
      group_data[[group_accession]] <- group_list
      next
    }
    ## When the 2nd column is FALSE, then this defined a protein in the group.
    ## The protein data structure is likely the most interesting.
    if (row[, 2] == FALSE | row[, 2] == TRUE) {
      protein_information <- row[protein_keepers]
      colnames(protein_information) <- protein_colnames
      protein_information[["ID"]] <- sub(pattern="^.* GN=(\\w+) .*$",
                                         replacement="\\1",
                                         x=protein_information[["Description"]])
      protein_accession <- protein_information[["Accession"]]
      protein_list <- list(
        "summary" = protein_information,
        "data" = data.frame())
      group_data[[group_accession]][["data"]][[protein_accession]] <- protein_list
      next
    }
    ## When the 3rd group is FALSE, then this adds a peptide.
    ## The peptide data structure is the most detailed, but probably not the most interesting.
    if (row[, 3] == FALSE | row[, 3] == TRUE) {
      peptide_information <- row[peptide_keepers]
      colnames(peptide_information) <- peptide_colnames
      current <- group_data[[group_accession]][["data"]][[protein_accession]][["data"]]
      new <- rbind(current, peptide_information)
      group_data[[group_accession]][["data"]][[protein_accession]][["data"]] <- new
      next
    }
  } ## End iterating over ever row of this unholy stupid data structure.
  close(bar)
  message("Finished parsing, reorganizing the protein data.")
  protein_df <- data.frame()
  peptide_df <- data.frame()
  protein_names <- c()
  message(paste0("Starting to iterate over ", length(group_data),  " groups."))
  bar <- utils::txtProgressBar(style=3)
  for (g in 1:length(group_data)) {
    pct_done <- g / length(group_data)
    setTxtProgressBar(bar, pct_done)
    group <- as.character(names(group_data)[g])
    protein_group <- group_data[[group]][["data"]]
    protein_accessions <- names(protein_group)
    for (p in 1:length(protein_accessions)) {
      protein <- protein_accessions[p]
      protein_names <- c(protein_names, protein)
      protein_summary <- group_data[[group]][["data"]][[protein]][["summary"]]
      protein_df <- rbind(protein_df, protein_summary)
      peptide_data <- group_data[[group]][["data"]][[protein]][["data"]]
      ## peptide_df <- rbind(peptide_df, peptide_data)
    }
  } ## End of the for loop
  close(bar)
  retlist <- list(
    "names" = protein_names,
    "group_data" = group_data,
    "protein_data" = protein_df,
    "peptide_data" = peptide_df)
  return(retlist)
}

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

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.

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.

I am also reasonably certain that I codified this logic in a separate function in hpgltools since writing this.

## Set some reasonable row names
protein_df <- less_difficult[["protein_data"]]
protein_names <- less_difficult[["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", "marked_as",
                  "num_peptides", "num_psms", "num_unique_pep", "num_prot_groups", "num_aas",
                  "mw", "pi", "mascot_score", "sequest_score", "num_peptides_mascot",
                  "num_peptides_sequest", "bp", "cc", "mf", "pfam_ids",
                  "entrez_gene_id", "ensembl_gene_id", "gene_symbol", "chromosome", "kegg",
                  "wikipathways", "found_trypsin_lysc_cf", "found_trypsin_lysc_wcl", "found_chymotrypsin_cf", "found_chymotrypsin_wcl",
                  "modifications", "id")
colnames(protein_df) <- new_colnames
## 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]])
}
factor_cols <- c("found_trypsin_lysc_cf", "found_trypsin_lysc_wcl",
                 "found_chymotrypsin_cf", "found_chymotrypsin_wcl")
for (col in factor_cols) {
  protein_df[[col]] <- as.factor(protein_df[[col]])
}
## 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[["found_chymotrypsin_cf"]] == "High" | protein_df[["found_chymotrypsin_wcl"]] == "High") &
  (protein_df[["found_trypsin_lysc_cf"]] != "High" & protein_df[["found_trypsin_lysc_wcl"]] != "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[["found_trypsin_lysc_cf"]] == "High" | protein_df[["found_trypsin_lysc_wcl"]] == "High") &
  (protein_df[["found_chymotrypsin_cf"]] != "High" & protein_df[["found_chymotrypsin_wcl"]] != "High")
protein_df[tryp_hits, "state"] <- "tryp_only"
no_hits <- protein_df[["found_trypsin_lysc_cf"]] == "Not Found" &
  protein_df[["found_trypsin_lysc_wcl"]] == "Not Found" &
  protein_df[["found_chymotrypsin_cf"]] == "Not Found" &
  protein_df[["found_chymotrypsin_wcl"]] == "Not Found"
protein_df[no_hits, "state"] <- "None"
protein_df[["state"]] <- as.factor(protein_df[["state"]])

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)
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.

coverage_wrt_pi

## Warning: Computation failed in `stat_density2d()`:
## bandwidths must be strictly positive

## 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.

coverage_wrt_mw

## Warning: Computation failed in `stat_density2d()`:
## bandwidths must be strictly positive

## 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

## Warning: Removed 96 rows containing non-finite values (stat_boxplot).

2 Bring together previous rnaseq and current peptides

This turned out to be a surprisingly difficult problem because the geneID<->annotation mappings were inconsistent. Therefore I downloaded the uniprot annotations for Mtb at:

http://www.uniprot.org/proteomes/UP000001584

The annotation data used by Keith appears to me to have come from either genbank, microbesonline, or ensembl. All of these match quite easily, but none of them match easily to uniprot.

2.1 Read in the limma data and merge the annotations.

Throughout the following steps I will print the dimensions of the resulting data so that I can see where I am losing information.

all_de <- read.table("limma_result.csv", header=TRUE, sep="\t")
colnames(all_de) <- c("ID", "logFC", "AveExpr", "adj.P.Val", "abbrev_id", "description", "gene_length", "type")
all_de <- merge(all_de, mtb_annotations, by.x="ID", by.y="ID")
dim(all_de)

## [1] 3995   21

2.2 Get separate PE proteins.

We will have occasion to look only at the P/PE proteins later. So lets separately grab that information now.

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

## Warning in merge.data.frame(pe_ids, mtb_annotations, by.x = "locus_tag", : column name
## 'locus_tag' is duplicated in the result

rownames(pe_annotations) <- pe_annotations[["ID"]]
colnames(pe_annotations) <- c("locus_tag", "seqnames", "start", "end", "width",
                              "strand", "source", "type", "score", "phase",
                              "locus_tag2", "gene", "description", "function")
dim(pe_annotations)

## [1] 161  14

2.3 Load uniprot data

There is in fact an existing interface for downloading annotation data from uniprot. It is terrible, because it seems to be missing important species and does not provide the mappings for important data. So, I downloaded the reference myself and wrote a parser for it.

if (file.exists("uniprot_annot.csv")) {
  message("Reusing a previous annotation database.")
  uniprot_annot <- read.csv(file="uniprot_annot.csv")
} else {
  uniprot_annot <- load_uniprot_annotations(file="reference/uniprot_3AUP000001584.txt.gz")
  write.csv(file="uniprot_annot.csv", uniprot_annot)
}

## Reusing a previous annotation database.

dim(uniprot_annot)

## [1] 3993   42

2.4 Perform merges

Now we have the relevant data sets, lets bring them together into what is hopefully a coherent whole.

## Start by merging Keith's DE and uniprot annotations.
uniprot_rna <- merge(all_de, uniprot_annot, by.x="ID", by.y="loci")
dim(uniprot_rna)

## [1] 3891   62

## Repeat for the set of 161 PE proteins.
uniprot_pe <- merge(pe_annotations, uniprot_annot, by.x="locus_tag", by.y="loci")
dim(uniprot_pe)

## [1] 157  55

## We appear to have lost 4 PE proteins in this process.
uniprot_pe_peptides <- merge(uniprot_annot, protein_df, by.x="primary_accessions", by.y="accession")
uniprot_pe_peptides <- merge(uniprot_pe_peptides, pe_ids, by.x="loci", by.y="locus_tag")
dim(uniprot_pe_peptides)

## [1] 40 80

uniprot_rna_protein <- merge(uniprot_rna, protein_df, by.x="primary_accessions", by.y="accession")
dim(protein_df)

## [1] 2671   39

dim(uniprot_rna_protein)

## [1] 2511  100

uniprot_rna_pe <- merge(uniprot_pe, uniprot_rna_protein, by="primary_id")

## Warning in merge.data.frame(uniprot_pe, uniprot_rna_protein, by = "primary_id"): column
## names 'description.x', 'description.y' are duplicated in the result

dim(uniprot_rna_pe)

## [1]  40 154

2.5 Use these combined data

Now that I finally got the various data sets to merge, I should be able to visualize how well they correspond to one another (not very).

But first, lets make some shorthand names for them so that this is a bit less obnoxious.

## This is exhausting, I will call these big and little from now on
big <- uniprot_rna_protein
rownames(big) <- make.names(big$ID, unique=TRUE)
ordered <- order(rownames(big))
big <- big[ordered, ]
little <- uniprot_rna_pe
rownames(little) <- make.names(little$ID, unique=TRUE)
ordered <- order(rownames(little))
little <- little[ordered, ]
shared_idx <- rownames(big) %in% rownames(little)
shared <- big[shared_idx, ]

2.6 Make pseudo_rpkm

Our gene-size-normalized expression values have not yet been generated. Do that here.

## 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.
big$ave_minus_log <- big$AveExpr - big$logFC
big$pseudo_rpkm <- (big$ave_minus_log / big$width) * 1000
maximum_rpkm <- max(big$pseudo_rpkm)
big$express_vs_max <- big$pseudo_rpkm / maximum_rpkm
## Since we constant ordered the dataframes above, this is valid.
little$pseudo_rpkm <- big$pseudo_rpkm[shared_idx]
little$express_vs_max <- big$express_vs_max[shared_idx]

write.csv(x=little, file="pe_relative_expression_annotation-20171205.csv")

2.7 Finally plot

Now all the pieces are finally in place, lets look at relationships between our available metrics of protein abundance and rna abundance.

2.8 Histograms of all expression vs. P/PE

library(hpgltools)
library(ggplot2)
test <- list(
  "all_peptides" = big$express_vs_max,
  "pe_peptides" = little$express_vs_max)
multi <- plot_multihistogram(test)

## Used Bon Ferroni corrected t test(s) between columns.

multi$plot

2.9 RPKM vs PSMS

psms_coverage_peptides <- big[, c("num_psms", "coverage", "num_peptides")]

rpkm_vs_psms <- big[, c("pseudo_rpkm", "num_psms")]
little_rpkm_vs_psms <- little[, c("pseudo_rpkm", "num_psms")]
rpkm_vs_psms$num_psms <- log2(rpkm_vs_psms$num_psms + 1)
little_rpkm_vs_psms$num_psms <- log2(little_rpkm_vs_psms$num_psms + 1)
rpkm_psms_linear <- plot_linear_scatter(rpkm_vs_psms)

## Used Bon Ferroni corrected t test(s) between columns.

rpkm_psms_linear$scatter

rpkm_psms_linear$cor

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 5.8, df = 2500, p-value = 9e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.07575 0.15296
## sample estimates:
##    cor 
## 0.1145

rpkm_psms_linear$lm_rsq

## [1] 0.007207

big_little <- ggplot(data=rpkm_vs_psms, aes_string(x="pseudo_rpkm", y="num_psms")) +
  geom_point()
big_little

big_little <- big_little +
  geom_point(data=little_rpkm_vs_psms, aes_string(x="pseudo_rpkm", y="num_psms"), colour="red")
favorites <- little_rpkm_vs_psms[["num_psms"]] > 3
labeled <- little_rpkm_vs_psms[favorites, ]
big_little +
  ggrepel::geom_text_repel(data=labeled,
                           arrow=arrow(length=unit(0.01, "npc")),
                           nudge_x=-4.0, nudge_y=-0.2,
                           ggplot2::aes(x=pseudo_rpkm, y=num_psms, label=rownames(labeled)))

labeled

##         pseudo_rpkm num_psms
## Rv0285       28.621    5.426
## Rv0442c       3.750    3.585
## Rv1195       34.594    4.807
## Rv1196        9.911    5.931
## Rv1386       20.836    6.895
## Rv1387        6.264    3.907
## Rv1768        3.798    3.700
## Rv2162c       5.764    4.170
## Rv2352c       5.946    3.807
## Rv2430c      15.085    3.907
## Rv2431c      27.463    4.954
## Rv3020c      14.002    4.907
## Rv3477       34.108    7.011
## Rv3478        9.491    4.807

2.10 rpkm vs peptides

rpkm_vs_peptides <- big[, c("pseudo_rpkm", "num_peptides")]
rpkm_peptides_linear <- plot_linear_scatter(rpkm_vs_peptides)

## Used Bon Ferroni corrected t test(s) between columns.

rpkm_peptides_linear$scatter

rpkm_peptides_linear$cor

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = -4.8, df = 2500, p-value = 2e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.13318 -0.05564
## sample estimates:
##      cor 
## -0.09455

rpkm_psms_linear$lm_rsq

## [1] 0.007207

little_rpkm_vs_peptides <- little[, c("pseudo_rpkm", "num_peptides")]
big_little <- ggplot(data=rpkm_vs_peptides, aes_string(x="pseudo_rpkm", y="num_peptides")) +
  geom_point()
big_little

big_little <- big_little +
  geom_point(data=little_rpkm_vs_peptides, aes_string(x="pseudo_rpkm", y="num_peptides"), colour="red")
big_little

rpkm_vs_coverage <- big[, c("pseudo_rpkm", "coverage")]
rpkm_coverage_linear <- plot_linear_scatter(rpkm_vs_coverage)

## Used Bon Ferroni corrected t test(s) between columns.

rpkm_coverage_linear$scatter

rpkm_coverage_linear$cor

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 21, df = 2500, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3532 0.4198
## sample estimates:
##   cor 
## 0.387

rpkm_coverage_linear$lm_rsq

## [1] 0.1618

little_rpkm_vs_coverage <- little[, c("pseudo_rpkm", "coverage")]
big_little <- ggplot(data=rpkm_vs_coverage, aes_string(x="pseudo_rpkm", y="coverage")) +
  geom_point()
big_little

big_little <- big_little +
  geom_point(data=little_rpkm_vs_coverage, aes_string(x="pseudo_rpkm", y="coverage"), colour="red")
big_little

3 PView vs rpkm

I recently spent some time playing with pview and believe it might provide a reasonable proxy for peptide abundance for our data. Lets see if that is true.

I first will need to copy the pview results from scratch/proteomics to here. I copied a test run of pview to preprocessing/pview_20171207.csv

pview_result <- read.csv("preprocessing/pview_20171207.csv", header=TRUE, sep="\t")
small_df <- pview_result[, c("cond.01.median", "cond.02.median")]
rownames(small_df) <- gsub(pattern="\\[mtb_cds\\] ", replacement="", x=pview_result[["protein.group"]])
na_positions <- is.na(small_df)
small_df[na_positions] <- 0
small_log_df <- log2(small_df + 1)

test_df <- small_df
nonzero_rows <- test_df[["cond.02.median"]] != 0
test_df <- test_df[nonzero_rows, ]
test_min <- min(test_df[["cond.02.median"]])
test_df[["cond.02.median"]] <- test_df[["cond.02.median"]] - (test_min - 1)
test_df[["l2"]] <- log2(test_df[["cond.02.median"]])

big_pview <- merge(big, test_df, by="row.names")
rownames(big_pview) <- big_pview[["Row.names"]]
rpkm_vs_pview <- big_pview[, c("AveExpr", "l2")]

rpkm_pview_linear <- plot_linear_scatter(rpkm_vs_pview, cormethod="spearman")
rpkm_pview_linear$scatter
rpkm_pview_linear$cor
rpkm_pview_linear$lm_rsq

little_rpkm_vs_coverage <- little[, c("pseudo_rpkm", "coverage")]
big_little <- ggplot(data=rpkm_vs_coverage, aes_string(x="pseudo_rpkm", y="coverage")) +
  geom_point()
big_little
big_little <- big_little +
  geom_point(data=little_rpkm_vs_coverage, aes_string(x="pseudo_rpkm", y="coverage"), colour="red")
big_little

if (!isTRUE(get0("skip_load"))) {
  pander::pander(sessionInfo())
  message(paste0("This is hpgltools commit: ", get_git_commit()))
  this_save <- paste0(gsub(pattern="\\.Rmd", replace="", x=rmd_file), "-v", ver, ".rda.xz")
  message(paste0("Saving to ", this_save))
  tmp <- sm(saveme(filename=this_save))
}

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

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

## > git reset 7de4503f6bb5724c28cce24af5dbee22bb1c0cae

## R> packrat::restore()

## This is hpgltools commit: Thu Apr 12 22:08:53 2018 -0400: 7de4503f6bb5724c28cce24af5dbee22bb1c0cae

## Saving to 02_eval_result_201711-v20171205.rda.xz