This document will first explore differentially expressed genes in humans 4 hours after infection followed by the same question in mice.

1 Which genes are DE in human macrophages at 4 hours upon infection with L. major?

2 Gather annotation data

I want to perform a series of comparisons among the host cells: human and mouse. Thus I need to collect annotation data for both species and get the set of orthologs between them.

2.1 Start with the human annotation data

In the following block, I download the human annotations from biomart. In addition, I take a moment to recreate the transcript IDs as observed in the salmon count tables (yes, I know they are not actually count tables). Finally, I create a table which maps transcripts to genes, this will be used when we generate the expressionset so that we get gene expression levels from transcripts via the R package ‘tximport’.

hs_annot <- load_biomart_annotations()$annotation

## The biomart annotations file already exists, loading from it.

rownames(hs_annot) <- make.names(
  paste0(hs_annot[["ensembl_transcript_id"]], ".",
         hs_annot[["transcript_version"]]),
  unique=TRUE)
hs_tx_gene <- hs_annot[, c("ensembl_gene_id", "ensembl_transcript_id")]
hs_tx_gene[["id"]] <- rownames(hs_tx_gene)
hs_tx_gene <- hs_tx_gene[, c("id", "ensembl_gene_id")]
new_hs_annot <- hs_annot
rownames(new_hs_annot) <- make.names(hs_annot[["ensembl_gene_id"]], unique=TRUE)

2.2 Generate expressionsets

The question is reasonably self-contained. I want to compare the uninfected human samples against any samples which were infected for 4 hours. So let us first pull those samples and then poke at them a bit.

The following block creates an expressionset using all human-quantified samples. As mentioned previously, it uses the table of transcript<->gene mappings, and the biomart annotations.

Given this set of ~440 samples, it then drops the following:

All samples marked ‘skipped’.
All samples which are not from time ‘t4h’.

and resets the condition and batch factors to the ‘infection state’ metadatum and ‘study’, respectively.

sample_sheet <- "sample_sheets/leishmania_host_metasheet_20190401.xlsx"
hs_expt <- create_expt(sample_sheet,
                       file_column="hsapiensfile",
                       gene_info=new_hs_annot,
                       tx_gene_map=hs_tx_gene)

## Reading the sample metadata.

## The sample definitions comprises: 437 rows(samples) and 55 columns(metadata fields).

## Reading count tables.

## Using the transcript to gene mapping.

## Reading salmon data with tximport.

## Finished reading count tables.

## Matched 19629 annotations and counts.

## Bringing together the count matrix and gene information.

## The mapped IDs are not the rownames of your gene information, changing them now.

## Some annotations were lost in merging, setting them to 'undefined'.

hs_expt_noskipped <- subset_expt(hs_expt, subset="skipped!='yes'")

## There were 267, now there are 247 samples.

hs_t4h_expt <- subset_expt(hs_expt_noskipped, subset="expttime=='t4h'")

## There were 247, now there are 64 samples.

hs_t4h_expt <- set_expt_conditions(hs_t4h_expt, fact="infectstate")
hs_t4h_expt <- set_expt_batches(hs_t4h_expt, fact="study")
table(hs_t4h_expt$conditions)

## 
## bead   no stim  yes 
##    3   18   35    8

table(hs_t4h_expt$batches)

## 
## lps-timecourse       m-gm-csf           mbio 
##              8             39             17

hs_written <- write_expt(hs_t4h_expt, excel="excel/HsM0Lm4h_expt.xlsx")

## Writing the first sheet, containing a legend and some summary data.

## Writing the raw reads.

## Graphing the raw reads.

## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

## Writing the normalized reads.

## Graphing the normalized reads.

## Writing the median reads by factor.

## The factor bead has 3 rows.

## The factor no has 18 rows.

## The factor stim has 35 rows.

## The factor yes has 8 rows.

## Note: zip::zip() is deprecated, please use zip::zipr() instead

2.3 Examine t4h vs uninfected

Let us perform some generic metrics of the t4h human expressionset. As per usual, I plot the metrics first of the raw data; followed by the same metrics of log2(quantile(cpm(sva(filtered(data))))).

hs_t4h_plots <- sm(graph_metrics(hs_t4h_expt))

hs_t4h_norm <- normalize_expt(hs_t4h_expt, norm="quant", convert="cpm",
                              transform="log2", filter=TRUE, batch="svaseq")

## This function will replace the expt$expressionset slot with:

## log2(svaseq(cpm(quant(cbcb(data)))))

## It backs up the current data into a slot named:
##  expt$backup_expressionset. It will also save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
##  when invoking limma.  The appropriate libsize is the non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Warning in normalize_expt(hs_t4h_expt, norm = "quant", convert = "cpm", :
## Quantile normalization and sva do not always play well together.

## Step 1: performing count filter with option: cbcb

## Removing 7360 low-count genes (12269 remaining).

## Step 2: normalizing the data with quant.

## Using normalize.quantiles.robust due to a thread error in preprocessCore.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## transform_counts: Found 3822 values equal to 0, adding 1 to the matrix.

## Step 5: doing batch correction with svaseq.

## Note to self:  If you get an error like 'x contains missing values' The data has too many 0's and needs a stronger low-count filter applied.

## Passing off to all_adjusters.

## batch_counts: Before batch/surrogate estimation, 711279 entries are x>1: 90.6%.

## batch_counts: Before batch/surrogate estimation, 3822 entries are x==0: 0.487%.

## batch_counts: Before batch/surrogate estimation, 70115 entries are 0<x<1: 8.93%.

## The be method chose 12 surrogate variable(s).

## Attempting svaseq estimation with 12 surrogates.

## There are 1760 (0.224%) elements which are < 0 after batch correction.

hs_t4h_norm_plots <- sm(graph_metrics(hs_t4h_norm))

2.3.1 Print some of the plots

hs_t4h_plots$legend

hs_t4h_plots$libsize

hs_t4h_plots$boxplot

hs_t4h_norm_plots$pc_plot

2.4 Remove stimulated samples

I perhaps should have removed the stimulated samples sooner, but I was curious to see their effect on the distribution first.

hs_t4h_inf <- subset_expt(hs_t4h_expt, subset="condition!='stim'")

## There were 64, now there are 29 samples.

hs_t4h_inf_norm <- normalize_expt(hs_t4h_inf, transform="log2", convert="cpm",
                                  filter=TRUE, batch="svaseq")

## This function will replace the expt$expressionset slot with:

## log2(svaseq(cpm(cbcb(data))))

## It backs up the current data into a slot named:
##  expt$backup_expressionset. It will also save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
##  when invoking limma.  The appropriate libsize is the non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Leaving the data unnormalized.  This is necessary for DESeq, but
##  EdgeR/limma might benefit from normalization.  Good choices include quantile,
##  size-factor, tmm, etc.

## Step 1: performing count filter with option: cbcb

## Removing 7759 low-count genes (11870 remaining).

## Step 2: not normalizing the data.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## transform_counts: Found 3276 values equal to 0, adding 1 to the matrix.

## Step 5: doing batch correction with svaseq.

## Note to self:  If you get an error like 'x contains missing values' The data has too many 0's and needs a stronger low-count filter applied.

## Passing off to all_adjusters.

## batch_counts: Before batch/surrogate estimation, 319013 entries are x>1: 92.7%.

## batch_counts: Before batch/surrogate estimation, 3276 entries are x==0: 0.952%.

## batch_counts: Before batch/surrogate estimation, 21941 entries are 0<x<1: 6.37%.

## The be method chose 6 surrogate variable(s).

## Attempting svaseq estimation with 6 surrogates.

## There are 557 (0.162%) elements which are < 0 after batch correction.

hs_t4h_pca <- plot_pca(hs_t4h_inf_norm, plot_title="H. sapiens, L. major, t4h")
hs_t4h_pca$plot

keepers <- list("infection" = c("yes", "no"))
hs_t4h_de <- all_pairwise(hs_t4h_inf, model_batch="svaseq", filter=TRUE, force=TRUE)

## batch_counts: Before batch/surrogate estimation, 339179 entries are x>1: 98.5%.

## batch_counts: Before batch/surrogate estimation, 3276 entries are x==0: 0.952%.

## batch_counts: Before batch/surrogate estimation, 216 entries are 0<x<1: 0.0627%.

## The be method chose 5 surrogate variable(s).

## Attempting svaseq estimation with 5 surrogates.

## This function will replace the expt$expressionset slot with:

## log2(cpm(quant(cbcb(data))))

## It backs up the current data into a slot named:
##  expt$backup_expressionset. It will also save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
##  when invoking limma.  The appropriate libsize is the non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Not correcting the count-data for batch effects.  If batch is
##  included in EdgerR/limma's model, then this is probably wise; but in extreme
##  batch effects this is a good parameter to play with.

## Step 1: performing count filter with option: cbcb

## Removing 0 low-count genes (11870 remaining).

## Step 2: normalizing the data with quant.

## Using normalize.quantiles.robust due to a thread error in preprocessCore.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## transform_counts: Found 1964 values equal to 0, adding 1 to the matrix.

## Step 5: not doing batch correction.

## Plotting a PCA before surrogates/batch inclusion.

## Using svaseq to visualize before/after batch inclusion.

## Performing a test normalization with: raw

## This function will replace the expt$expressionset slot with:

## log2(svaseq(cpm(cbcb(data))))

## It backs up the current data into a slot named:
##  expt$backup_expressionset. It will also save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
##  when invoking limma.  The appropriate libsize is the non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Leaving the data unnormalized.  This is necessary for DESeq, but
##  EdgeR/limma might benefit from normalization.  Good choices include quantile,
##  size-factor, tmm, etc.

## Step 1: performing count filter with option: cbcb

## Removing 0 low-count genes (11870 remaining).

## Step 2: not normalizing the data.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## transform_counts: Found 3276 values equal to 0, adding 1 to the matrix.

## Step 5: doing batch correction with svaseq.

## Note to self:  If you get an error like 'x contains missing values' The data has too many 0's and needs a stronger low-count filter applied.

## Passing off to all_adjusters.

## batch_counts: Before batch/surrogate estimation, 319013 entries are x>1: 92.7%.

## batch_counts: Before batch/surrogate estimation, 3276 entries are x==0: 0.952%.

## batch_counts: Before batch/surrogate estimation, 21941 entries are 0<x<1: 6.37%.

## The be method chose 6 surrogate variable(s).

## Attempting svaseq estimation with 6 surrogates.

## There are 557 (0.162%) elements which are < 0 after batch correction.

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

hs_t4h_table <- combine_de_tables(hs_t4h_de, keepers=keepers,
                                  excel="excel/HsM0Lm4h_de_tables.xlsx")

## Deleting the file excel/HsM0Lm4h_de_tables.xlsx before writing the tables.

## Writing a legend of columns.

## Printing a pca plot before/after surrogates/batch estimation.

## Working on 1/1: infection which is: yes/no.

## Found table with yes_vs_no

## 20181210 a pthread error in normalize.quantiles leads me to robust.

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

## Adding venn plots for infection.

## Limma expression coefficients for infection; R^2: 0.92; equation: y = 0.996x - 0.0675

## Edger expression coefficients for infection; R^2: 0.92; equation: y = 1.08x - 0.728

## DESeq2 expression coefficients for infection; R^2: 0.92; equation: y = 1.08x - 0.701

## Writing summary information.

## Attempting to add the comparison plot to pairwise_summary at row: 23 and column: 1

## Performing save of the workbook.

hs_t4h_sig <- extract_significant_genes(hs_t4h_table, excel="excel/HsM0Lm4h_sig_tables.xlsx")

## Writing a legend of columns.

## The count is: 1 and the test is: limma.

## The count is: 2 and the test is: edger.

## The count is: 3 and the test is: deseq.

## The count is: 4 and the test is: ebseq.

## The count is: 5 and the test is: basic.

## Writing excel data according to limma for infection: 1/5.

## After (adj)p filter, the up genes table has 2976 genes.

## After (adj)p filter, the down genes table has 4109 genes.

## After fold change filter, the up genes table has 1081 genes.

## After fold change filter, the down genes table has 871 genes.

## Printing significant genes to the file: excel/HsM0Lm4h_sig_tables.xlsx

## 1/1: Creating significant table up_1limma_infection

## Writing excel data according to edger for infection: 1/5.

## After (adj)p filter, the up genes table has 3428 genes.

## After (adj)p filter, the down genes table has 3346 genes.

## After fold change filter, the up genes table has 1243 genes.

## After fold change filter, the down genes table has 1029 genes.

## Printing significant genes to the file: excel/HsM0Lm4h_sig_tables.xlsx

## 1/1: Creating significant table up_1edger_infection

## Writing excel data according to deseq for infection: 1/5.

## After (adj)p filter, the up genes table has 3600 genes.

## After (adj)p filter, the down genes table has 3712 genes.

## After fold change filter, the up genes table has 1226 genes.

## After fold change filter, the down genes table has 1056 genes.

## Printing significant genes to the file: excel/HsM0Lm4h_sig_tables.xlsx

## 1/1: Creating significant table up_1deseq_infection

## Writing excel data according to ebseq for infection: 1/5.

## After (adj)p filter, the up genes table has 1893 genes.

## After (adj)p filter, the down genes table has 1940 genes.

## After fold change filter, the up genes table has 895 genes.

## After fold change filter, the down genes table has 1031 genes.

## Printing significant genes to the file: excel/HsM0Lm4h_sig_tables.xlsx

## 1/1: Creating significant table up_1ebseq_infection

## Writing excel data according to basic for infection: 1/5.

## After (adj)p filter, the up genes table has 1335 genes.

## After (adj)p filter, the down genes table has 2590 genes.

## After fold change filter, the up genes table has 965 genes.

## After fold change filter, the down genes table has 1441 genes.

## Printing significant genes to the file: excel/HsM0Lm4h_sig_tables.xlsx

## 1/1: Creating significant table up_1basic_infection

## Adding significance bar plots.

3 Which genes are DE in mouse macrophages at 4 hours upon infection with L. major?

Most of this should be the same in process as what was performed for the human.

3.1 Gather annotation data

I want to perform a series of comparisons among the host cells: human and mouse. Thus I need to collect annotation data for both species and get the set of orthologs between them.

3.1.1 Start with the human annotation data

mm_annot <- load_biomart_annotations(species="mmusculus")$annotation

## The biomart annotations file already exists, loading from it.

rownames(mm_annot) <- make.names(
  paste0(mm_annot[["ensembl_transcript_id"]], ".",
         mm_annot[["transcript_version"]]),
  unique=TRUE)
mm_tx_gene <- mm_annot[, c("ensembl_gene_id", "ensembl_transcript_id")]
mm_tx_gene[["id"]] <- rownames(mm_tx_gene)
mm_tx_gene <- mm_tx_gene[, c("id", "ensembl_gene_id")]
new_mm_annot <- mm_annot
rownames(new_mm_annot) <- make.names(mm_annot[["ensembl_gene_id"]], unique=TRUE)

3.2 Generate expressionsets

mm_expt <- create_expt(sample_sheet,
                       file_column="mmusculusfile",
                       gene_info=new_mm_annot,
                       tx_gene_map=mm_tx_gene)

## Reading the sample metadata.

## The sample definitions comprises: 437 rows(samples) and 55 columns(metadata fields).

## Reading count tables.

## Using the transcript to gene mapping.

## Reading salmon data with tximport.

## Finished reading count tables.

## Matched 19660 annotations and counts.

## Bringing together the count matrix and gene information.

## The mapped IDs are not the rownames of your gene information, changing them now.

## Some annotations were lost in merging, setting them to 'undefined'.

mm_t4h_expt <- subset_expt(mm_expt, subset="expttime=='t4h'")

## There were 105, now there are 41 samples.

mm_t4h_expt <- set_expt_conditions(mm_t4h_expt, fact="infectstate")
table(mm_t4h_expt$conditions)

## 
##   no stim  yes 
##   11   24    6

table(mm_t4h_expt$batches)

## 
## undefined 
##        41

mm_written <- write_expt(mm_t4h_expt, excel="excel/MmM0Lm4h_expt.xlsx")

## Writing the first sheet, containing a legend and some summary data.

## Writing the raw reads.

## Graphing the raw reads.

## Writing the normalized reads.

## Graphing the normalized reads.

## Writing the median reads by factor.

## The factor no has 11 rows.

## The factor stim has 24 rows.

## The factor yes has 6 rows.

3.3 Examine t4h vs uninfected

mm_t4h_plots <- sm(graph_metrics(mm_t4h_expt))

mm_t4h_norm <- normalize_expt(mm_t4h_expt, norm="quant", convert="cpm",
                              transform="log2", filter=TRUE, batch="svaseq")

## This function will replace the expt$expressionset slot with:

## log2(svaseq(cpm(quant(cbcb(data)))))

## It backs up the current data into a slot named:
##  expt$backup_expressionset. It will also save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
##  when invoking limma.  The appropriate libsize is the non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Warning in normalize_expt(mm_t4h_expt, norm = "quant", convert = "cpm", :
## Quantile normalization and sva do not always play well together.

## Step 1: performing count filter with option: cbcb

## Removing 9350 low-count genes (10310 remaining).

## Step 2: normalizing the data with quant.

## Using normalize.quantiles.robust due to a thread error in preprocessCore.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## transform_counts: Found 38 values equal to 0, adding 1 to the matrix.

## Step 5: doing batch correction with svaseq.

## Note to self:  If you get an error like 'x contains missing values' The data has too many 0's and needs a stronger low-count filter applied.

## Passing off to all_adjusters.

## batch_counts: Before batch/surrogate estimation, 390888 entries are x>1: 92.5%.

## batch_counts: Before batch/surrogate estimation, 38 entries are x==0: 0.00899%.

## batch_counts: Before batch/surrogate estimation, 31784 entries are 0<x<1: 7.52%.

## The be method chose 8 surrogate variable(s).

## Attempting svaseq estimation with 8 surrogates.

## There are 807 (0.191%) elements which are < 0 after batch correction.

mm_t4h_norm_plots <- sm(graph_metrics(mm_t4h_norm))

3.3.1 Print some of the plots

mm_t4h_plots$legend

mm_t4h_plots$libsize

mm_t4h_plots$boxplot

mm_t4h_norm_plots$pc_plot

3.4 Perform de analyses

mm_t4h_inf_norm <- normalize_expt(mm_t4h_expt, transform="log2", convert="cpm",
                                  filter=TRUE, batch="svaseq")

## This function will replace the expt$expressionset slot with:

## log2(svaseq(cpm(cbcb(data))))

## It backs up the current data into a slot named:
##  expt$backup_expressionset. It will also save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
##  when invoking limma.  The appropriate libsize is the non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Leaving the data unnormalized.  This is necessary for DESeq, but
##  EdgeR/limma might benefit from normalization.  Good choices include quantile,
##  size-factor, tmm, etc.

## Step 1: performing count filter with option: cbcb

## Removing 9350 low-count genes (10310 remaining).

## Step 2: not normalizing the data.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## transform_counts: Found 3253 values equal to 0, adding 1 to the matrix.

## Step 5: doing batch correction with svaseq.

## Note to self:  If you get an error like 'x contains missing values' The data has too many 0's and needs a stronger low-count filter applied.

## Passing off to all_adjusters.

## batch_counts: Before batch/surrogate estimation, 385615 entries are x>1: 91.2%.

## batch_counts: Before batch/surrogate estimation, 3253 entries are x==0: 0.770%.

## batch_counts: Before batch/surrogate estimation, 33842 entries are 0<x<1: 8.01%.

## The be method chose 7 surrogate variable(s).

## Attempting svaseq estimation with 7 surrogates.

## There are 1186 (0.281%) elements which are < 0 after batch correction.

mm_t4h_pca <- plot_pca(mm_t4h_inf_norm, plot_title="M. musculus, L. major, t4h")
mm_t4h_pca$plot

mm_t4h_de <- all_pairwise(mm_t4h_expt, model_batch="svaseq", filter=TRUE, force=TRUE)

## batch_counts: Before batch/surrogate estimation, 417208 entries are x>1: 98.7%.

## batch_counts: Before batch/surrogate estimation, 3253 entries are x==0: 0.770%.

## batch_counts: Before batch/surrogate estimation, 392 entries are 0<x<1: 0.0927%.

## The be method chose 5 surrogate variable(s).

## Attempting svaseq estimation with 5 surrogates.

## This function will replace the expt$expressionset slot with:

## log2(cpm(quant(cbcb(data))))

## It backs up the current data into a slot named:
##  expt$backup_expressionset. It will also save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
##  when invoking limma.  The appropriate libsize is the non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Not correcting the count-data for batch effects.  If batch is
##  included in EdgerR/limma's model, then this is probably wise; but in extreme
##  batch effects this is a good parameter to play with.

## Step 1: performing count filter with option: cbcb

## Removing 0 low-count genes (10310 remaining).

## Step 2: normalizing the data with quant.

## Using normalize.quantiles.robust due to a thread error in preprocessCore.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## transform_counts: Found 38 values equal to 0, adding 1 to the matrix.

## Step 5: not doing batch correction.

## Plotting a PCA before surrogates/batch inclusion.

## Using svaseq to visualize before/after batch inclusion.

## Performing a test normalization with: raw

## This function will replace the expt$expressionset slot with:

## log2(svaseq(cpm(cbcb(data))))

## It backs up the current data into a slot named:
##  expt$backup_expressionset. It will also save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
##  when invoking limma.  The appropriate libsize is the non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Leaving the data unnormalized.  This is necessary for DESeq, but
##  EdgeR/limma might benefit from normalization.  Good choices include quantile,
##  size-factor, tmm, etc.

## Step 1: performing count filter with option: cbcb

## Removing 0 low-count genes (10310 remaining).

## Step 2: not normalizing the data.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## transform_counts: Found 3253 values equal to 0, adding 1 to the matrix.

## Step 5: doing batch correction with svaseq.

## Note to self:  If you get an error like 'x contains missing values' The data has too many 0's and needs a stronger low-count filter applied.

## Passing off to all_adjusters.

## batch_counts: Before batch/surrogate estimation, 385615 entries are x>1: 91.2%.

## batch_counts: Before batch/surrogate estimation, 3253 entries are x==0: 0.770%.

## batch_counts: Before batch/surrogate estimation, 33842 entries are 0<x<1: 8.01%.

## The be method chose 7 surrogate variable(s).

## Attempting svaseq estimation with 7 surrogates.

## There are 1186 (0.281%) elements which are < 0 after batch correction.

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

mm_t4h_table <- combine_de_tables(mm_t4h_de, keepers=keepers,
                                  excel="excel/MmM0Lm4h_de_tables.xlsx")

## Deleting the file excel/MmM0Lm4h_de_tables.xlsx before writing the tables.

## Writing a legend of columns.

## Printing a pca plot before/after surrogates/batch estimation.

## Working on 1/1: infection which is: yes/no.

## Found table with yes_vs_no

## 20181210 a pthread error in normalize.quantiles leads me to robust.

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

## Adding venn plots for infection.

## Limma expression coefficients for infection; R^2: 0.92; equation: y = 0.95x + 0.187

## Edger expression coefficients for infection; R^2: 0.91; equation: y = 1.01x - 0.0739

## DESeq2 expression coefficients for infection; R^2: 0.909; equation: y = 1.01x - 0.0692

## Writing summary information.

## Attempting to add the comparison plot to pairwise_summary at row: 23 and column: 1

## Performing save of the workbook.

mm_t4h_sig <- extract_significant_genes(mm_t4h_table,
                                        excel="excel/MmM0Lm4h_sig_tables.xlsx")

## Writing a legend of columns.

## The count is: 1 and the test is: limma.

## The count is: 2 and the test is: edger.

## The count is: 3 and the test is: deseq.

## The count is: 4 and the test is: ebseq.

## The count is: 5 and the test is: basic.

## Writing excel data according to limma for infection: 1/5.

## After (adj)p filter, the up genes table has 2234 genes.

## After (adj)p filter, the down genes table has 2524 genes.

## After fold change filter, the up genes table has 825 genes.

## After fold change filter, the down genes table has 895 genes.

## Printing significant genes to the file: excel/MmM0Lm4h_sig_tables.xlsx

## 1/1: Creating significant table up_1limma_infection

## Writing excel data according to edger for infection: 1/5.

## After (adj)p filter, the up genes table has 2153 genes.

## After (adj)p filter, the down genes table has 1967 genes.

## After fold change filter, the up genes table has 924 genes.

## After fold change filter, the down genes table has 911 genes.

## Printing significant genes to the file: excel/MmM0Lm4h_sig_tables.xlsx

## 1/1: Creating significant table up_1edger_infection

## Writing excel data according to deseq for infection: 1/5.

## After (adj)p filter, the up genes table has 2474 genes.

## After (adj)p filter, the down genes table has 2404 genes.

## After fold change filter, the up genes table has 929 genes.

## After fold change filter, the down genes table has 932 genes.

## Printing significant genes to the file: excel/MmM0Lm4h_sig_tables.xlsx

## 1/1: Creating significant table up_1deseq_infection

## Writing excel data according to ebseq for infection: 1/5.

## After (adj)p filter, the up genes table has 1149 genes.

## After (adj)p filter, the down genes table has 1310 genes.

## After fold change filter, the up genes table has 644 genes.

## After fold change filter, the down genes table has 710 genes.

## Printing significant genes to the file: excel/MmM0Lm4h_sig_tables.xlsx

## 1/1: Creating significant table up_1ebseq_infection

## Writing excel data according to basic for infection: 1/5.

## After (adj)p filter, the up genes table has 2 genes.

## After (adj)p filter, the down genes table has 3 genes.

## After fold change filter, the up genes table has 2 genes.

## After fold change filter, the down genes table has 3 genes.

## Printing significant genes to the file: excel/MmM0Lm4h_sig_tables.xlsx

## 1/1: Creating significant table up_1basic_infection

## Adding significance bar plots.

3.5 Compare this to the previous result.

Let us see if our human differential expression result is similar to that obtained in Table S2.

previous_hs <- readxl::read_excel("excel/inline-supplementary-material-5.xls", sheet=2)
previous_hs_lfc <- previous_hs[, c("ID", "Fold change")]
neg_idx <- previous_hs_lfc[[2]] < 0
previous_hs_lfc[neg_idx, 2] <- -1 * (1 / previous_hs_lfc[neg_idx, 2])
previous_hs_lfc[[2]] <- log2(previous_hs_lfc[[2]])

merged <- merge(previous_hs_lfc, hs_t4h_table$data[[1]], by.x="ID", by.y="row.names")
cor.test(merged[["limma_logfc"]], merged[["Fold change"]])

## 
##  Pearson's product-moment correlation
## 
## data:  merged[["limma_logfc"]] and merged[["Fold change"]]
## t = 127, df = 5082, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8645 0.8778
## sample estimates:
##    cor 
## 0.8713

3.6 Compare the previous mouse and these mouse results.

previous_mm <- readxl::read_excel("excel/12864_2015_2237_MOESM3_ESM.xls", sheet=2, skip=1)
previous_mm_lfc <- previous_mm[, c("ID", "Fold change")]
neg_idx <- previous_mm_lfc[[2]] < 0
previous_mm_lfc[neg_idx, 2] <- -1 * (1 / previous_mm_lfc[neg_idx, 2])
previous_mm_lfc[[2]] <- log2(previous_mm_lfc[[2]])

merged <- merge(previous_mm_lfc, mm_t4h_table$data[[1]], by.x="ID", by.y="row.names")
cor.test(merged[["limma_logfc"]], merged[["Fold change"]])

## 
##  Pearson's product-moment correlation
## 
## data:  merged[["limma_logfc"]] and merged[["Fold change"]]
## t = 223, df = 5718, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9445 0.9498
## sample estimates:
##    cor 
## 0.9472

4 What genes are shared in the mouse and human data?

This is one method of addressing Najibs big question #1c. In this method, I am taking the tables for the human analysis and mouse analysis separately, then merging them using a table of orthologs.

Side note: a different way of addressing this question resides in 20190220_host_comparisons.Rmd. In this competing method, the table of orthologs is used in the beginning to make a single set of IDs for the human and mouse genes, then perform the differential expression analysis.

4.1 Extract human mouse orthologs

My load_biomart_orthologs() function should provide this mapping gene ID table.

## The defaults of this function are suitable for mouse/human queries.
mm_hs_ortho <- load_biomart_orthologs()$all_linked_genes

## Unable to perform useMart, perhaps the host/mart is incorrect: dec2016.archive.ensembl.org ENSEMBL_MART_ENSEMBL.

## The available first_marts are:

## ENSEMBL_MART_ENSEMBLENSEMBL_MART_MOUSEENSEMBL_MART_SNPENSEMBL_MART_FUNCGENENSEMBL_MART_VEGA

## Trying the first one.

mm_table <- mm_t4h_table$data[[1]]
hs_table <- hs_t4h_table$data[[1]]

mm_table <- merge(mm_hs_ortho, mm_table, by.x="mmusculus", by.y="row.names", all.y=TRUE)
hs_table <- merge(mm_hs_ortho, hs_table, by.x="hsapiens", by.y="row.names", all.y=TRUE)
both_table_hs <- merge(hs_table, mm_table, by.x="hsapiens", by.y="hsapiens")
both_table_mm <- merge(hs_table, mm_table, by.x="mmusculus", by.y="mmusculus")

cor.test(both_table_hs[["limma_logfc.x"]], both_table_hs[["limma_logfc.y"]])

## 
##  Pearson's product-moment correlation
## 
## data:  both_table_hs[["limma_logfc.x"]] and both_table_hs[["limma_logfc.y"]]
## t = 34, df = 13699, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2656 0.2964
## sample estimates:
##    cor 
## 0.2811

cor.test(both_table_mm[["limma_logfc.x"]], both_table_mm[["limma_logfc.y"]])

## 
##  Pearson's product-moment correlation
## 
## data:  both_table_mm[["limma_logfc.x"]] and both_table_mm[["limma_logfc.y"]]
## t = 34, df = 21859, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2135 0.2387
## sample estimates:
##    cor 
## 0.2261

tt <- plot_scatter(both_table_hs[, c("limma_logfc.x", "limma_logfc.y")])
tt

I believe these both_table_hs and both_table_mm tables are good candidates for the set of genes which are shared across the human and mouse samples, from the perspective of the human.

Separately

Now lets write some of these out.

shared_hsmm_up_idx <- both_table_hs[["deseq_logfc.x"]] > 1 &
  both_table_hs[["deseq_logfc.y"]] > 1 &
  both_table_hs[["deseq_adjp.x"]] <= 0.05 &
  both_table_hs[["deseq_adjp.y"]] <= 0.05
shared_hsmm_up <- both_table_hs[shared_hsmm_up_idx, ]
written <- write_xls(data=shared_hsmm_up,
                     excel="excel/HsM0Lm4h_vs_MmM0Lm4h_shared_up_hs.xlsx")

## Saving to: excel/HsM0Lm4h_vs_MmM0Lm4h_shared_up_hs.xlsx

shared_mmhs_up_idx <- both_table_mm[["deseq_logfc.x"]] > 1 &
  both_table_mm[["deseq_logfc.y"]] > 1 &
  both_table_mm[["deseq_adjp.x"]] <= 0.05 &
  both_table_mm[["deseq_adjp.y"]] <= 0.05
shared_mmhs_up <- both_table_mm[shared_mmhs_up_idx, ]
written <- write_xls(data=shared_mmhs_up,
                     excel="excel/HsM0Lm4h_vs_MmM0Lm4h_shared_up_mm.xlsx")

## Saving to: excel/HsM0Lm4h_vs_MmM0Lm4h_shared_up_mm.xlsx

shared_hsmm_down_idx <- both_table_hs[["deseq_logfc.x"]] < -1 &
  both_table_hs[["deseq_logfc.y"]] < -1 &
  both_table_hs[["deseq_adjp.x"]] <= 0.05 &
  both_table_hs[["deseq_adjp.y"]] <= 0.05
shared_hsmm_down <- both_table_hs[shared_hsmm_down_idx, ]
written <- write_xls(data=shared_hsmm_down,
                     excel="excel/HsM0Lm4h_vs_MmM0Lm4h_shared_down_hs.xlsx")

## Saving to: excel/HsM0Lm4h_vs_MmM0Lm4h_shared_down_hs.xlsx

shared_mmhs_down_idx <- both_table_mm[["deseq_logfc.x"]] < -1 &
  both_table_mm[["deseq_logfc.y"]] < -1 &
  both_table_mm[["deseq_adjp.x"]] <= 0.05 &
  both_table_mm[["deseq_adjp.y"]] <= 0.05
shared_mmhs_down <- both_table_mm[shared_mmhs_down_idx, ]
written <- write_xls(data=shared_mmhs_down,
                     excel="excel/HsM0Lm4h_vs_MmM0Lm4h_shared_down_mm.xlsx")

## Saving to: excel/HsM0Lm4h_vs_MmM0Lm4h_shared_down_mm.xlsx

4.2 Compare the above with Table S7 from the Laura and Cecilia paper.

Table S7 has a set of genes from human which are also up-regulated in mouse upon infection.

fig_s7_up <- readxl::read_excel("excel/inline-supplementary-material-7.xls", sheet=3)
fig_s7_hs_up <- unique(fig_s7_up[[5]])

fig_s7_down <- readxl::read_excel("excel/inline-supplementary-material-7.xls", sheet=4)
fig_s7_hs_down <- unique(fig_s7_down[[5]])

both_up_idx <- both_table_hs[["limma_logfc.x"]] >= 0.8 &
  both_table_hs[["limma_logfc.y"]] >= 0.8 &
  both_table_hs[["limma_adjp.x"]] <= 0.1 &
  both_table_hs[["limma_adjp.y"]] <= 0.1
both_up_ids <- both_table_hs[both_up_idx, "hsapiens"]
up_venn <- Vennerable::Venn(Sets=list("figs7" = fig_s7_hs_up, "tables" = both_up_ids))
Vennerable::plot(up_venn)

both_down_idx <- both_table_hs[["limma_logfc.x"]] <= -0.8 &
  both_table_hs[["limma_logfc.y"]] <= -0.8 &
  both_table_hs[["limma_adjp.x"]] <= 0.1 &
  both_table_hs[["limma_adjp.y"]] <= 0.1
both_down_ids <- both_table_hs[both_down_idx, "hsapiens"]
down_venn <- Vennerable::Venn(Sets=list("figs7" = fig_s7_hs_down, "tables" = both_down_ids))
Vennerable::plot(down_venn)

5 Next question: Which 4 hour genes are shared with CIDEIM?

“Which 4 hour DE genes are shared with the CIDEIM panamensis infected human macrophages?”

So, once again there are two ways of approaching this question:

Examine the question of infected vs. uninfected for the two data sets separately. Then query the results and see what comes out.
Examine this as a single question using the union of both data sets.

Since I already have a putative answer for the human 4 hour macrophage data, The simplest method is to just look at the cideim data and do #1 above. So let us do that first.

cideim_macr <- subset_expt(hs_expt, subset="lab=='tmrc'&hostcelltype=='macrophage'")

## There were 267, now there are 23 samples.

cideim_macr <- set_expt_conditions(cideim_macr, fact="infectstate")
cideim_norm <- normalize_expt(cideim_macr, transform="log2", convert="cpm",
                              norm="quant", filter="simple", batch="svaseq")

## This function will replace the expt$expressionset slot with:

## log2(svaseq(cpm(quant(simple(data)))))

## It backs up the current data into a slot named:
##  expt$backup_expressionset. It will also save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
##  when invoking limma.  The appropriate libsize is the non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Warning in normalize_expt(cideim_macr, transform = "log2", convert =
## "cpm", : Quantile normalization and sva do not always play well together.

## Step 1: performing count filter with option: simple

## Removing 2127 low-count genes (17502 remaining).

## Step 2: normalizing the data with quant.

## Using normalize.quantiles.robust due to a thread error in preprocessCore.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## transform_counts: Found 33659 values equal to 0, adding 1 to the matrix.

## Step 5: doing batch correction with svaseq.

## Note to self:  If you get an error like 'x contains missing values' The data has too many 0's and needs a stronger low-count filter applied.

## Passing off to all_adjusters.

## batch_counts: Before batch/surrogate estimation, 278638 entries are x>1: 69.2%.

## batch_counts: Before batch/surrogate estimation, 33659 entries are x==0: 8.36%.

## batch_counts: Before batch/surrogate estimation, 90249 entries are 0<x<1: 22.4%.

## The be method chose 5 surrogate variable(s).

## Attempting svaseq estimation with 5 surrogates.

## There are 13850 (3.44%) elements which are < 0 after batch correction.

plot_pca(cideim_norm)$plot

cideim_de <- all_pairwise(cideim_macr, filter=TRUE, parallel=FALSE, force=TRUE)

## This function will replace the expt$expressionset slot with:

## log2(cpm(quant(cbcb(data))))

## It backs up the current data into a slot named:
##  expt$backup_expressionset. It will also save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
##  when invoking limma.  The appropriate libsize is the non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Not correcting the count-data for batch effects.  If batch is
##  included in EdgerR/limma's model, then this is probably wise; but in extreme
##  batch effects this is a good parameter to play with.

## Step 1: performing count filter with option: cbcb

## Removing 0 low-count genes (12499 remaining).

## Step 2: normalizing the data with quant.

## Using normalize.quantiles.robust due to a thread error in preprocessCore.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## transform_counts: Found 547 values equal to 0, adding 1 to the matrix.

## Step 5: not doing batch correction.

## Plotting a PCA before surrogates/batch inclusion.

## Using limma's removeBatchEffect to visualize with(out) batch inclusion.

## Starting basic_pairwise().

## Starting basic pairwise comparison.

## Leaving the data alone, regardless of normalization state.

## Basic step 0/3: Transforming data.

## Basic step 1/3: Creating median and variance tables.

## Basic step 2/3: Performing 3 comparisons.

## Basic step 3/3: Creating faux DE Tables.

## Basic: Returning tables.

## Starting deseq_pairwise().

## Starting DESeq2 pairwise comparisons.

## About to round the data, this is a pretty terrible thing to do. But if you, like me, want to see what happens when you put non-standard data into deseq, then here you go.

## Warning in choose_binom_dataset(input, force = force): This data was
## inappropriately forced into integers.

## The condition+batch model failed. Does your experimental design support both condition and batch? Using only a conditional model.

## Choosing the non-intercept containing model.

## DESeq2 step 1/5: Including batch and condition in the deseq model.

## converting counts to integer mode

## DESeq2 step 2/5: Estimate size factors.

## DESeq2 step 3/5: Estimate dispersions.

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## Using a parametric fitting seems to have worked.

## DESeq2 step 4/5: nbinomWaldTest.

## Starting ebseq_pairwise().

## The data should be suitable for EdgeR/DESeq/EBSeq. If they freak out, check the state of the count table and ensure that it is in integer counts.

## Starting EBSeq pairwise subset.

## Choosing the non-intercept containing model.

## Starting EBTest of no vs. yes.

## Copying ppee values as ajusted p-values until I figure out how to deal with them.

## Starting edger_pairwise().

## Starting edgeR pairwise comparisons.

## About to round the data, this is a pretty terrible thing to do. But if you, like me, want to see what happens when you put non-standard data into deseq, then here you go.

## Warning in choose_binom_dataset(input, force = force): This data was
## inappropriately forced into integers.

## The condition+batch model failed. Does your experimental design support both condition and batch? Using only a conditional model.

## Choosing the non-intercept containing model.

## EdgeR step 1/9: Importing and normalizing data.

## EdgeR step 2/9: Estimating the common dispersion.

## EdgeR step 3/9: Estimating dispersion across genes.

## EdgeR step 4/9: Estimating GLM Common dispersion.

## EdgeR step 5/9: Estimating GLM Trended dispersion.

## EdgeR step 6/9: Estimating GLM Tagged dispersion.

## EdgeR step 7/9: Running glmFit, switch to glmQLFit by changing the argument 'edger_test'.

## EdgeR step 8/9: Making pairwise contrasts.

## Starting limma_pairwise().

## Starting limma pairwise comparison.

## Leaving the data alone, regardless of normalization state.

## libsize was not specified, this parameter has profound effects on limma's result.

## Using the libsize from expt$best_libsize.

## Limma step 1/6: choosing model.

## The condition+batch model failed. Does your experimental design support both condition and batch? Using only a conditional model.

## Choosing the non-intercept containing model.

## Limma step 2/6: running limma::voom(), switch with the argument 'which_voom'.

## Using normalize.method=quantile for voom.

## Limma step 3/6: running lmFit with method: ls.

## Limma step 4/6: making and fitting contrasts with no intercept. (~ 0 + factors)

## Limma step 5/6: Running eBayes with robust=FALSE and trend=FALSE.

## Limma step 6/6: Writing limma outputs.

## Limma step 6/6: 1/1: Creating table: yes_vs_no.  Adjust=BH

## Limma step 6/6: 1/2: Creating table: no.  Adjust=BH

## Limma step 6/6: 2/2: Creating table: yes.  Adjust=BH

## Comparing analyses.

cideim_table <- combine_de_tables(cideim_de, keepers=keepers)

## Writing a legend of columns.

## Working on 1/1: infection which is: yes/no.

## Found table with yes_vs_no

## 20181210 a pthread error in normalize.quantiles leads me to robust.

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

cideim_merged <- merge(hs_t4h_table$data[[1]], cideim_table$data[[1]], by="row.names")
cor.test(cideim_merged[["deseq_logfc.x"]], cideim_merged[["deseq_logfc.y"]])

## 
##  Pearson's product-moment correlation
## 
## data:  cideim_merged[["deseq_logfc.x"]] and cideim_merged[["deseq_logfc.y"]]
## t = 31, df = 11570, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2574 0.2911
## sample estimates:
##    cor 
## 0.2743

plot_scatter(cideim_merged[, c("limma_logfc.x", "limma_logfc.y")])

5.1 Try again as one big expressionset

hs_expt_noskipped <- subset_expt(hs_expt, subset="skipped!='yes'")

## There were 267, now there are 247 samples.

hs_t4h_expt <- subset_expt(hs_expt_noskipped, subset="expttime=='t4h'")

## There were 247, now there are 64 samples.

hs_t4h_expt <- set_expt_conditions(hs_t4h_expt, fact="infectstate")
cideim_t4h <- subset_expt(
  hs_expt,
  subset="(skipped!='yes'&expttime=='t4h')|(lab=='tmrc'&hostcelltype=='macrophage')")

## There were 267, now there are 87 samples.

cideim_t4h <- set_expt_conditions(cideim_t4h, fact="infectstate")

cideim_t4h_norm <- normalize_expt(cideim_t4h, filter=TRUE, norm="quant",
                                  convert="cpm", transform="log2")

## This function will replace the expt$expressionset slot with:

## log2(cpm(quant(cbcb(data))))

## It backs up the current data into a slot named:
##  expt$backup_expressionset. It will also save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
##  when invoking limma.  The appropriate libsize is the non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Not correcting the count-data for batch effects.  If batch is
##  included in EdgerR/limma's model, then this is probably wise; but in extreme
##  batch effects this is a good parameter to play with.

## Step 1: performing count filter with option: cbcb

## Removing 6648 low-count genes (12981 remaining).

## Step 2: normalizing the data with quant.

## Using normalize.quantiles.robust due to a thread error in preprocessCore.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## transform_counts: Found 1937 values equal to 0, adding 1 to the matrix.

## Step 5: not doing batch correction.

plot_pca(cideim_t4h_norm)$plot

cideim_t4h_normbatch <- normalize_expt(cideim_t4h, filter=TRUE, norm="quant",
                                       convert="cpm", batch="svaseq", transform="log2")

## This function will replace the expt$expressionset slot with:

## log2(svaseq(cpm(quant(cbcb(data)))))

## It backs up the current data into a slot named:
##  expt$backup_expressionset. It will also save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
##  when invoking limma.  The appropriate libsize is the non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Warning in normalize_expt(cideim_t4h, filter = TRUE, norm = "quant",
## convert = "cpm", : Quantile normalization and sva do not always play well
## together.

## Step 1: performing count filter with option: cbcb

## Removing 6648 low-count genes (12981 remaining).

## Step 2: normalizing the data with quant.

## Using normalize.quantiles.robust due to a thread error in preprocessCore.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## transform_counts: Found 1937 values equal to 0, adding 1 to the matrix.

## Step 5: doing batch correction with svaseq.

## Note to self:  If you get an error like 'x contains missing values' The data has too many 0's and needs a stronger low-count filter applied.

## Passing off to all_adjusters.

## batch_counts: Before batch/surrogate estimation, 997436 entries are x>1: 88.3%.

## batch_counts: Before batch/surrogate estimation, 1937 entries are x==0: 0.172%.

## batch_counts: Before batch/surrogate estimation, 129974 entries are 0<x<1: 11.5%.

## The be method chose 14 surrogate variable(s).

## Attempting svaseq estimation with 14 surrogates.

## There are 3170 (0.281%) elements which are < 0 after batch correction.

plot_pca(cideim_t4h_normbatch)$plot

cideim_t4h_de <- all_pairwise(cideim_t4h, model_batch="svaseq",
                              filter=TRUE, force=TRUE, parallel=FALSE)

## batch_counts: Before batch/surrogate estimation, 1089986 entries are x>1: 96.5%.

## batch_counts: Before batch/surrogate estimation, 25706 entries are x==0: 2.28%.

## batch_counts: Before batch/surrogate estimation, 1303 entries are 0<x<1: 0.115%.

## The be method chose 10 surrogate variable(s).

## Attempting svaseq estimation with 10 surrogates.

## This function will replace the expt$expressionset slot with:

## log2(cpm(quant(cbcb(data))))

## It backs up the current data into a slot named:
##  expt$backup_expressionset. It will also save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
##  when invoking limma.  The appropriate libsize is the non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Not correcting the count-data for batch effects.  If batch is
##  included in EdgerR/limma's model, then this is probably wise; but in extreme
##  batch effects this is a good parameter to play with.

## Step 1: performing count filter with option: cbcb

## Removing 0 low-count genes (12981 remaining).

## Step 2: normalizing the data with quant.

## Using normalize.quantiles.robust due to a thread error in preprocessCore.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## transform_counts: Found 1937 values equal to 0, adding 1 to the matrix.

## Step 5: not doing batch correction.

## Plotting a PCA before surrogates/batch inclusion.

## Using svaseq to visualize before/after batch inclusion.

## Performing a test normalization with: raw

## This function will replace the expt$expressionset slot with:

## log2(svaseq(cpm(cbcb(data))))

## It backs up the current data into a slot named:
##  expt$backup_expressionset. It will also save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
##  when invoking limma.  The appropriate libsize is the non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Leaving the data unnormalized.  This is necessary for DESeq, but
##  EdgeR/limma might benefit from normalization.  Good choices include quantile,
##  size-factor, tmm, etc.

## Step 1: performing count filter with option: cbcb

## Removing 0 low-count genes (12981 remaining).

## Step 2: not normalizing the data.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## transform_counts: Found 25706 values equal to 0, adding 1 to the matrix.

## Step 5: doing batch correction with svaseq.

## Note to self:  If you get an error like 'x contains missing values' The data has too many 0's and needs a stronger low-count filter applied.

## Passing off to all_adjusters.

## batch_counts: Before batch/surrogate estimation, 988589 entries are x>1: 87.5%.

## batch_counts: Before batch/surrogate estimation, 25706 entries are x==0: 2.28%.

## batch_counts: Before batch/surrogate estimation, 115052 entries are 0<x<1: 10.2%.

## The be method chose 11 surrogate variable(s).

## Attempting svaseq estimation with 11 surrogates.

## There are 5446 (0.482%) elements which are < 0 after batch correction.

## Starting basic_pairwise().

## Starting basic pairwise comparison.

## Leaving the data alone, regardless of normalization state.

## Basic step 0/3: Transforming data.

## Basic step 1/3: Creating median and variance tables.

## Basic step 2/3: Performing 10 comparisons.

## Basic step 3/3: Creating faux DE Tables.

## Basic: Returning tables.

## Starting deseq_pairwise().

## Starting DESeq2 pairwise comparisons.

## About to round the data, this is a pretty terrible thing to do. But if you, like me, want to see what happens when you put non-standard data into deseq, then here you go.

## Warning in choose_binom_dataset(input, force = force): This data was
## inappropriately forced into integers.

## Including batch estimates from sva/ruv/pca in the model.

## Choosing the non-intercept containing model.

## DESeq2 step 1/5: Including a matrix of batch estimates in the deseq model.

## converting counts to integer mode

## DESeq2 step 2/5: Estimate size factors.

## DESeq2 step 3/5: Estimate dispersions.

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## Using a parametric fitting seems to have worked.

## DESeq2 step 4/5: nbinomWaldTest.

## Starting ebseq_pairwise().

## The data should be suitable for EdgeR/DESeq/EBSeq. If they freak out, check the state of the count table and ensure that it is in integer counts.

## Starting EBSeq pairwise subset.

## Choosing the non-intercept containing model.

## Starting EBTest of bead vs. no.

## Copying ppee values as ajusted p-values until I figure out how to deal with them.

## Starting EBTest of bead vs. stim.

## Copying ppee values as ajusted p-values until I figure out how to deal with them.

## Starting EBTest of bead vs. yes.

## Copying ppee values as ajusted p-values until I figure out how to deal with them.

## Starting EBTest of no vs. stim.

## Copying ppee values as ajusted p-values until I figure out how to deal with them.

## Starting EBTest of no vs. yes.

## Copying ppee values as ajusted p-values until I figure out how to deal with them.

## Starting EBTest of stim vs. yes.

## Copying ppee values as ajusted p-values until I figure out how to deal with them.

## Starting edger_pairwise().

## Starting edgeR pairwise comparisons.

## About to round the data, this is a pretty terrible thing to do. But if you, like me, want to see what happens when you put non-standard data into deseq, then here you go.

## Warning in choose_binom_dataset(input, force = force): This data was
## inappropriately forced into integers.

## Including batch estimates from sva/ruv/pca in the model.

## Choosing the non-intercept containing model.

## EdgeR step 1/9: Importing and normalizing data.

## EdgeR step 2/9: Estimating the common dispersion.

## EdgeR step 3/9: Estimating dispersion across genes.

## EdgeR step 4/9: Estimating GLM Common dispersion.

## EdgeR step 5/9: Estimating GLM Trended dispersion.

## EdgeR step 6/9: Estimating GLM Tagged dispersion.

## EdgeR step 7/9: Running glmFit, switch to glmQLFit by changing the argument 'edger_test'.

## EdgeR step 8/9: Making pairwise contrasts.

## Starting limma_pairwise().

## Starting limma pairwise comparison.

## Leaving the data alone, regardless of normalization state.

## libsize was not specified, this parameter has profound effects on limma's result.

## Using the libsize from expt$best_libsize.

## Limma step 1/6: choosing model.

## Including batch estimates from sva/ruv/pca in the model.

## Choosing the non-intercept containing model.

## Limma step 2/6: running limma::voom(), switch with the argument 'which_voom'.

## Using normalize.method=quantile for voom.

## Limma step 3/6: running lmFit with method: ls.

## Limma step 4/6: making and fitting contrasts with no intercept. (~ 0 + factors)

## Limma step 5/6: Running eBayes with robust=FALSE and trend=FALSE.

## Limma step 6/6: Writing limma outputs.

## Limma step 6/6: 1/6: Creating table: no_vs_bead.  Adjust=BH

## Limma step 6/6: 2/6: Creating table: stim_vs_bead.  Adjust=BH

## Limma step 6/6: 3/6: Creating table: yes_vs_bead.  Adjust=BH

## Limma step 6/6: 4/6: Creating table: stim_vs_no.  Adjust=BH

## Limma step 6/6: 5/6: Creating table: yes_vs_no.  Adjust=BH

## Limma step 6/6: 6/6: Creating table: yes_vs_stim.  Adjust=BH

## Limma step 6/6: 1/4: Creating table: bead.  Adjust=BH

## Limma step 6/6: 2/4: Creating table: no.  Adjust=BH

## Limma step 6/6: 3/4: Creating table: stim.  Adjust=BH

## Limma step 6/6: 4/4: Creating table: yes.  Adjust=BH

## Comparing analyses.

cideim_t4h_tables <- combine_de_tables(cideim_t4h_de, keepers=keepers,
                                       excel="excel/HsM0Lm4h_vs_HsM0Lp_de_tables.xlsx")

## Writing a legend of columns.

## Printing a pca plot before/after surrogates/batch estimation.

## Working on 1/1: infection which is: yes/no.

## Found table with yes_vs_no

## 20181210 a pthread error in normalize.quantiles leads me to robust.

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

## Adding venn plots for infection.

## Limma expression coefficients for infection; R^2: 0.977; equation: y = 0.983x + 0.0397

## Edger expression coefficients for infection; R^2: 0.975; equation: y = 0.981x + 0.154

## DESeq2 expression coefficients for infection; R^2: 0.975; equation: y = 0.981x + 0.173

## Writing summary information.

## Attempting to add the comparison plot to pairwise_summary at row: 23 and column: 1

## Performing save of the workbook.

cideim_t4h_sig <- combine_de_tables(cideim_t4h_tables,
                                    excel="excel/HsM0Lm4h_vs_HsM0Lp_sig_tables.xlsx")

## Not adding plots, limma had an error.

## Not adding plots, deseq had an error.

## Not adding plots, edger had an error.

## Not adding plots, basic had an error.

## Writing a legend of columns.

## Error in combine_de_tables(cideim_t4h_tables, excel = "excel/HsM0Lm4h_vs_HsM0Lp_sig_tables.xlsx"): None of the DE tools appear to have worked.

6 Next, overlap between macrophages and neutrophils

All of the neutrophil data is in mouse, apparently. This will make it more difficult, perhaps impossible to get an accurate answer.

So instead, look at infection vs. uninfected in mouse and then compare to the earliest Sacks’ timepoints in neutrophils.

subset <- "(hostcelltype=='PMN'&host=='mus_musculus'&expttime=='t12h') |
           (hostcelltype=='macrophage'&host=='mus_musculus')"
neut_macr_mus <- subset_expt(mm_expt, subset=subset)

## There were 105, now there are 80 samples.

neut_macr_mus <- subset_expt(neut_macr_mus, subset="infectstate!='stim'")

## There were 80, now there are 56 samples.

neut_macr_mus <- set_expt_conditions(neut_macr_mus, fact="infectstate")
neut_macr_mus <- set_expt_batches(neut_macr_mus, fact="hostcelltype")
neut_macr_mus_norm <- normalize_expt(neut_macr_mus, convert="cpm",
                                     norm="quant", filter=TRUE)

## This function will replace the expt$expressionset slot with:

## cpm(quant(cbcb(data)))

## It backs up the current data into a slot named:
##  expt$backup_expressionset. It will also save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
##  when invoking limma.  The appropriate libsize is the non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Leaving the data in its current base format, keep in mind that
##  some metrics are easier to see when the data is log2 transformed, but
##  EdgeR/DESeq do not accept transformed data.

## Not correcting the count-data for batch effects.  If batch is
##  included in EdgerR/limma's model, then this is probably wise; but in extreme
##  batch effects this is a good parameter to play with.

## Step 1: performing count filter with option: cbcb

## Removing 0 low-count genes (19660 remaining).

## Step 2: normalizing the data with quant.

## Using normalize.quantiles.robust due to a thread error in preprocessCore.

## Step 3: converting the data with cpm.

## Step 4: not transforming the data.

## Step 5: not doing batch correction.

plot_pca(neut_macr_mus_norm)$plot

neut_macr_mus_normbatch <- normalize_expt(neut_macr_mus, convert="cpm",
                                          norm="quant", filter=TRUE, batch="svaseq")

## This function will replace the expt$expressionset slot with:

## svaseq(cpm(quant(cbcb(data))))

## It backs up the current data into a slot named:
##  expt$backup_expressionset. It will also save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
##  when invoking limma.  The appropriate libsize is the non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Leaving the data in its current base format, keep in mind that
##  some metrics are easier to see when the data is log2 transformed, but
##  EdgeR/DESeq do not accept transformed data.

## Warning in normalize_expt(neut_macr_mus, convert = "cpm", norm = "quant", :
## Quantile normalization and sva do not always play well together.

## Step 1: performing count filter with option: cbcb

## Removing 0 low-count genes (19660 remaining).

## Step 2: normalizing the data with quant.

## Using normalize.quantiles.robust due to a thread error in preprocessCore.

## Step 3: converting the data with cpm.

## Step 4: not transforming the data.

## Step 5: doing batch correction with svaseq.

## Note to self:  If you get an error like 'x contains missing values' The data has too many 0's and needs a stronger low-count filter applied.

## Passing off to all_adjusters.

## batch_counts: Before batch/surrogate estimation, 515332 entries are x>1: 46.8%.

## batch_counts: Before batch/surrogate estimation, 323375 entries are x==0: 29.4%.

## batch_counts: Before batch/surrogate estimation, 262253 entries are 0<x<1: 23.8%.

## The be method chose 10 surrogate variable(s).

## Attempting svaseq estimation with 10 surrogates.

## There are 30709 (2.79%) elements which are < 0 after batch correction.

plot_pca(neut_macr_mus_normbatch)$plot

neut_macr_mus_filt <- normalize_expt(neut_macr_mus, filter="simple")

## This function will replace the expt$expressionset slot with:

## simple(data)

## It backs up the current data into a slot named:
##  expt$backup_expressionset. It will also save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
##  when invoking limma.  The appropriate libsize is the non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Leaving the data in its current base format, keep in mind that
##  some metrics are easier to see when the data is log2 transformed, but
##  EdgeR/DESeq do not accept transformed data.

## Leaving the data unconverted.  It is often advisable to cpm/rpkm
##  the data to normalize for sampling differences, keep in mind though that rpkm
##  has some annoying biases, and voom() by default does a cpm (though hpgl_voom()
##  will try to detect this).

## Leaving the data unnormalized.  This is necessary for DESeq, but
##  EdgeR/limma might benefit from normalization.  Good choices include quantile,
##  size-factor, tmm, etc.

## Not correcting the count-data for batch effects.  If batch is
##  included in EdgerR/limma's model, then this is probably wise; but in extreme
##  batch effects this is a good parameter to play with.

## Step 1: performing count filter with option: simple

## Removing 3175 low-count genes (16485 remaining).

## Step 2: not normalizing the data.

## Step 3: not converting the data.

## Step 4: not transforming the data.

## Step 5: not doing batch correction.

neut_macr_mus_de <- all_pairwise(neut_macr_mus_filt, parallel=FALSE,
                                 force=TRUE, model_batch="svaseq")

## batch_counts: Before batch/surrogate estimation, 616059 entries are x>1: 66.7%.

## batch_counts: Before batch/surrogate estimation, 269846 entries are x==0: 29.2%.

## batch_counts: Before batch/surrogate estimation, 2347 entries are 0<x<1: 0.254%.

## The be method chose 3 surrogate variable(s).

## Attempting svaseq estimation with 3 surrogates.

## This function will replace the expt$expressionset slot with:

## log2(cpm(quant(cbcb(data))))

## It backs up the current data into a slot named:
##  expt$backup_expressionset. It will also save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
##  when invoking limma.  The appropriate libsize is the non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Not correcting the count-data for batch effects.  If batch is
##  included in EdgerR/limma's model, then this is probably wise; but in extreme
##  batch effects this is a good parameter to play with.

## Step 1: performing count filter with option: cbcb

## Removing 0 low-count genes (16485 remaining).

## Step 2: normalizing the data with quant.

## Using normalize.quantiles.robust due to a thread error in preprocessCore.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## transform_counts: Found 117451 values equal to 0, adding 1 to the matrix.

## Step 5: not doing batch correction.

## Plotting a PCA before surrogates/batch inclusion.

## Using svaseq to visualize before/after batch inclusion.

## Performing a test normalization with: raw

## This function will replace the expt$expressionset slot with:

## log2(svaseq(cpm(cbcb(data))))

## It backs up the current data into a slot named:
##  expt$backup_expressionset. It will also save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
##  when invoking limma.  The appropriate libsize is the non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Leaving the data unnormalized.  This is necessary for DESeq, but
##  EdgeR/limma might benefit from normalization.  Good choices include quantile,
##  size-factor, tmm, etc.

## Step 1: performing count filter with option: cbcb

## Removing 0 low-count genes (16485 remaining).

## Step 2: not normalizing the data.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## transform_counts: Found 269846 values equal to 0, adding 1 to the matrix.

## Step 5: doing batch correction with svaseq.

## Note to self:  If you get an error like 'x contains missing values' The data has too many 0's and needs a stronger low-count filter applied.

## Passing off to all_adjusters.

## batch_counts: Before batch/surrogate estimation, 511083 entries are x>1: 55.4%.

## batch_counts: Before batch/surrogate estimation, 269846 entries are x==0: 29.2%.

## batch_counts: Before batch/surrogate estimation, 142231 entries are 0<x<1: 15.4%.

## The be method chose 7 surrogate variable(s).

## Attempting svaseq estimation with 7 surrogates.

## There are 61691 (6.68%) elements which are < 0 after batch correction.

## Starting basic_pairwise().

## Starting basic pairwise comparison.

## Leaving the data alone, regardless of normalization state.

## Basic step 0/3: Transforming data.

## Basic step 1/3: Creating median and variance tables.

## Basic step 2/3: Performing 3 comparisons.

## Basic step 3/3: Creating faux DE Tables.

## Basic: Returning tables.

## Starting deseq_pairwise().

## Starting DESeq2 pairwise comparisons.

## About to round the data, this is a pretty terrible thing to do. But if you, like me, want to see what happens when you put non-standard data into deseq, then here you go.

## Warning in choose_binom_dataset(input, force = force): This data was
## inappropriately forced into integers.

## Including batch estimates from sva/ruv/pca in the model.

## Choosing the non-intercept containing model.

## DESeq2 step 1/5: Including a matrix of batch estimates in the deseq model.

## converting counts to integer mode

## DESeq2 step 2/5: Estimate size factors.

## DESeq2 step 3/5: Estimate dispersions.

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## Using a parametric fitting seems to have worked.

## DESeq2 step 4/5: nbinomWaldTest.

## Starting ebseq_pairwise().

## The data should be suitable for EdgeR/DESeq/EBSeq. If they freak out, check the state of the count table and ensure that it is in integer counts.

## Starting EBSeq pairwise subset.

## Choosing the non-intercept containing model.

## Starting EBTest of no vs. yes.

## Copying ppee values as ajusted p-values until I figure out how to deal with them.

## Starting edger_pairwise().

## Starting edgeR pairwise comparisons.

## About to round the data, this is a pretty terrible thing to do. But if you, like me, want to see what happens when you put non-standard data into deseq, then here you go.

## Warning in choose_binom_dataset(input, force = force): This data was
## inappropriately forced into integers.

## Including batch estimates from sva/ruv/pca in the model.

## Choosing the non-intercept containing model.

## EdgeR step 1/9: Importing and normalizing data.

## EdgeR step 2/9: Estimating the common dispersion.

## EdgeR step 3/9: Estimating dispersion across genes.

## EdgeR step 4/9: Estimating GLM Common dispersion.

## EdgeR step 5/9: Estimating GLM Trended dispersion.

## EdgeR step 6/9: Estimating GLM Tagged dispersion.

## EdgeR step 7/9: Running glmFit, switch to glmQLFit by changing the argument 'edger_test'.

## EdgeR step 8/9: Making pairwise contrasts.

## Starting limma_pairwise().

## Starting limma pairwise comparison.

## Leaving the data alone, regardless of normalization state.

## libsize was not specified, this parameter has profound effects on limma's result.

## Using the libsize from expt$best_libsize.

## Limma step 1/6: choosing model.

## Including batch estimates from sva/ruv/pca in the model.

## Choosing the non-intercept containing model.

## Limma step 2/6: running limma::voom(), switch with the argument 'which_voom'.

## Using normalize.method=quantile for voom.

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Limma step 3/6: running lmFit with method: ls.

## Limma step 4/6: making and fitting contrasts with no intercept. (~ 0 + factors)

## Limma step 5/6: Running eBayes with robust=FALSE and trend=FALSE.

## Limma step 6/6: Writing limma outputs.

## Limma step 6/6: 1/1: Creating table: yes_vs_no.  Adjust=BH

## Limma step 6/6: 1/2: Creating table: no.  Adjust=BH

## Limma step 6/6: 2/2: Creating table: yes.  Adjust=BH

## Comparing analyses.

neut_macr_mus_table <- combine_de_tables(neut_macr_mus_de, keepers=keepers,
                                         excel="excel/MmM0Lm4h_vs_MmPMNLm12h_de_tables.xlsx")

## Writing a legend of columns.

## Printing a pca plot before/after surrogates/batch estimation.

## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

## Working on 1/1: infection which is: yes/no.

## Found table with yes_vs_no

## 20181210 a pthread error in normalize.quantiles leads me to robust.

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

## Adding venn plots for infection.

## Limma expression coefficients for infection; R^2: 0.992; equation: y = 1.01x - 0.066

## Edger expression coefficients for infection; R^2: 0.99; equation: y = 1.01x - 0.12

## DESeq2 expression coefficients for infection; R^2: 0.991; equation: y = 0.995x - 0.0119

## Writing summary information.

## Attempting to add the comparison plot to pairwise_summary at row: 23 and column: 1

## Performing save of the workbook.

neut_macr_mus_sig <- extract_significant_genes(neut_macr_mus_table,
                                         excel="excel/MmM0Lm4h_vs_MmPMNLm12h_sig_tables.xlsx")

## Writing a legend of columns.

## The count is: 1 and the test is: limma.

## The count is: 2 and the test is: edger.

## The count is: 3 and the test is: deseq.

## The count is: 4 and the test is: ebseq.

## The count is: 5 and the test is: basic.

## Writing excel data according to limma for infection: 1/5.

## After (adj)p filter, the up genes table has 68 genes.

## After (adj)p filter, the down genes table has 37 genes.

## After fold change filter, the up genes table has 50 genes.

## After fold change filter, the down genes table has 31 genes.

## Printing significant genes to the file: excel/MmM0Lm4h_vs_MmPMNLm12h_sig_tables.xlsx

## 1/1: Creating significant table up_1limma_infection

## Writing excel data according to edger for infection: 1/5.

## After (adj)p filter, the up genes table has 350 genes.

## After (adj)p filter, the down genes table has 155 genes.

## After fold change filter, the up genes table has 193 genes.

## After fold change filter, the down genes table has 102 genes.

## Printing significant genes to the file: excel/MmM0Lm4h_vs_MmPMNLm12h_sig_tables.xlsx

## 1/1: Creating significant table up_1edger_infection

## Writing excel data according to deseq for infection: 1/5.

## After (adj)p filter, the up genes table has 287 genes.

## After (adj)p filter, the down genes table has 299 genes.

## After fold change filter, the up genes table has 163 genes.

## After fold change filter, the down genes table has 144 genes.

## Printing significant genes to the file: excel/MmM0Lm4h_vs_MmPMNLm12h_sig_tables.xlsx

## 1/1: Creating significant table up_1deseq_infection

## Writing excel data according to ebseq for infection: 1/5.

## After (adj)p filter, the up genes table has 288 genes.

## After (adj)p filter, the down genes table has 684 genes.

## After fold change filter, the up genes table has 173 genes.

## After fold change filter, the down genes table has 564 genes.

## Printing significant genes to the file: excel/MmM0Lm4h_vs_MmPMNLm12h_sig_tables.xlsx

## 1/1: Creating significant table up_1ebseq_infection

## Writing excel data according to basic for infection: 1/5.

## After (adj)p filter, the up genes table has 6 genes.

## After (adj)p filter, the down genes table has 0 genes.

## After fold change filter, the up genes table has 6 genes.

## After fold change filter, the down genes table has 0 genes.

## Printing significant genes to the file: excel/MmM0Lm4h_vs_MmPMNLm12h_sig_tables.xlsx

## 1/1: Creating significant table up_1basic_infection

## Adding significance bar plots.

I think this handles questions a through e?

pander::pander(sessionInfo())

R version 3.6.0 beta (2019-04-11 r76379)

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: parallel, stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: edgeR(v.3.24.3), foreach(v.1.4.4), ruv(v.0.9.7), hpgltools(v.1.0), Biobase(v.2.42.0) and BiocGenerics(v.0.28.0)

loaded via a namespace (and not attached): tidyselect(v.0.2.5), lme4(v.1.1-21), htmlwidgets(v.1.3), RSQLite(v.2.1.1), AnnotationDbi(v.1.44.0), grid(v.3.6.0), BiocParallel(v.1.16.6), Rtsne(v.0.15), devtools(v.2.0.1), munsell(v.0.5.0), codetools(v.0.2-16), preprocessCore(v.1.45.0), withr(v.2.1.2), colorspace(v.1.4-1), GOSemSim(v.2.8.0), knitr(v.1.22), rstudioapi(v.0.10), stats4(v.3.6.0), Vennerable(v.3.1.0.9000), robustbase(v.0.93-4), DOSE(v.3.8.2), labeling(v.0.3), urltools(v.1.7.2), tximport(v.1.10.1), GenomeInfoDbData(v.1.2.0), polyclip(v.1.10-0), bit64(v.0.9-7), farver(v.1.1.0), rprojroot(v.1.3-2), xfun(v.0.6), R6(v.2.4.0), doParallel(v.1.0.14), GenomeInfoDb(v.1.18.2), locfit(v.1.5-9.1), bitops(v.1.0-6), fgsea(v.1.8.0), gridGraphics(v.0.3-0), DelayedArray(v.0.8.0), assertthat(v.0.2.1), scales(v.1.0.0), ggraph(v.1.0.2), nnet(v.7.3-12), enrichplot(v.1.2.0), gtable(v.0.3.0), sva(v.3.30.1), processx(v.3.3.0), rlang(v.0.3.3), genefilter(v.1.64.0), splines(v.3.6.0), rtracklayer(v.1.42.2), lazyeval(v.0.2.2), acepack(v.1.4.1), checkmate(v.1.9.1), europepmc(v.0.3), yaml(v.2.2.0), reshape2(v.1.4.3), GenomicFeatures(v.1.34.7), backports(v.1.1.3), qvalue(v.2.14.1), Hmisc(v.4.2-0), RBGL(v.1.58.2), clusterProfiler(v.3.10.1), tools(v.3.6.0), usethis(v.1.4.0), ggplotify(v.0.0.3), ggplot2(v.3.1.0), gplots(v.3.0.1.1), RColorBrewer(v.1.1-2), blockmodeling(v.0.3.4), sessioninfo(v.1.1.1), ggridges(v.0.5.1), Rcpp(v.1.0.1), plyr(v.1.8.4), base64enc(v.0.1-3), progress(v.1.2.0), zlibbioc(v.1.28.0), purrr(v.0.3.2), RCurl(v.1.95-4.12), ps(v.1.3.0), prettyunits(v.1.0.2), rpart(v.4.1-13), viridis(v.0.5.1), cowplot(v.0.9.4), S4Vectors(v.0.20.1), SummarizedExperiment(v.1.12.0), ggrepel(v.0.8.0), cluster(v.2.0.8), colorRamps(v.2.3), fs(v.1.2.7), variancePartition(v.1.12.3), magrittr(v.1.5), data.table(v.1.12.0), DO.db(v.2.9), openxlsx(v.4.1.0), triebeard(v.0.3.0), packrat(v.0.5.0), matrixStats(v.0.54.0), pkgload(v.1.0.2), hms(v.0.4.2), evaluate(v.0.13), xtable(v.1.8-3), pbkrtest(v.0.4-7), XML(v.3.98-1.19), readxl(v.1.3.1), IRanges(v.2.16.0), gridExtra(v.2.3), testthat(v.2.0.1), compiler(v.3.6.0), biomaRt(v.2.38.0), tibble(v.2.1.1), KernSmooth(v.2.23-15), crayon(v.1.3.4), minqa(v.1.2.4), htmltools(v.0.3.6), mgcv(v.1.8-28), corpcor(v.1.6.9), snow(v.0.4-3), Formula(v.1.2-3), geneplotter(v.1.60.0), tidyr(v.0.8.3), DBI(v.1.0.0), tweenr(v.1.0.1), MASS(v.7.3-51.3), boot(v.1.3-20), Matrix(v.1.2-17), readr(v.1.3.1), cli(v.1.1.0), quadprog(v.1.5-5), gdata(v.2.18.0), igraph(v.1.2.4), GenomicRanges(v.1.34.0), pkgconfig(v.2.0.2), registry(v.0.5-1), rvcheck(v.0.1.3), GenomicAlignments(v.1.18.1), foreign(v.0.8-71), xml2(v.1.2.0), annotate(v.1.60.1), rngtools(v.1.3.1), pkgmaker(v.0.27), XVector(v.0.22.0), bibtex(v.0.4.2), doRNG(v.1.7.1), EBSeq(v.1.22.1), stringr(v.1.4.0), callr(v.3.2.0), digest(v.0.6.18), graph(v.1.60.0), Biostrings(v.2.50.2), cellranger(v.1.1.0), rmarkdown(v.1.12), fastmatch(v.1.1-0), htmlTable(v.1.13.1), directlabels(v.2018.05.22), curl(v.3.3), Rsamtools(v.1.34.1), gtools(v.3.8.1), nloptr(v.1.2.1), nlme(v.3.1-137), jsonlite(v.1.6), desc(v.1.2.0), viridisLite(v.0.3.0), limma(v.3.38.3), pillar(v.1.3.1), lattice(v.0.20-38), DEoptimR(v.1.0-8), httr(v.1.4.0), pkgbuild(v.1.0.3), survival(v.2.44-1.1), GO.db(v.3.7.0), glue(v.1.3.1), remotes(v.2.0.2), zip(v.2.0.1), UpSetR(v.1.3.3), iterators(v.1.0.10), pander(v.0.6.3), bit(v.1.1-14), ggforce(v.0.2.1), stringi(v.1.4.3), blob(v.1.1.1), DESeq2(v.1.22.2), doSNOW(v.1.0.16), latticeExtra(v.0.6-28), caTools(v.1.17.1.2), memoise(v.1.1.0) and dplyr(v.0.8.0.1)

message("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 0ebb3165f07d676a83da460824f337251efbcf69

## This is hpgltools commit: Fri Apr 12 15:03:48 2019 -0400: 0ebb3165f07d676a83da460824f337251efbcf69

## message(paste0("Saving to ", savefile))
## tmp <- sm(saveme(filename=savefile))

LS0tCnRpdGxlOiAiMjAxOTA0MTcgQW4gYXR0ZW1wdCB0byBhbnN3ZXIgb25lIG9mIHRoZSBiaWcgcXVlc3Rpb25zIGZyb20gTmFqaWIuIgphdXRob3I6ICJhdGIgYWJlbGV3QGdtYWlsLmNvbSIKZGF0ZTogImByIFN5cy5EYXRlKClgIgpvdXRwdXQ6CiAgaHRtbF9kb2N1bWVudDoKICAgIGNvZGVfZG93bmxvYWQ6IHRydWUKICAgIGNvZGVfZm9sZGluZzogc2hvdwogICAgZmlnX2NhcHRpb246IHRydWUKICAgIGZpZ19oZWlnaHQ6IDcKICAgIGZpZ193aWR0aDogNwogICAgaGlnaGxpZ2h0OiB0YW5nbwogICAga2VlcF9tZDogZmFsc2UKICAgIG1vZGU6IHNlbGZjb250YWluZWQKICAgIG51bWJlcl9zZWN0aW9uczogdHJ1ZQogICAgc2VsZl9jb250YWluZWQ6IHRydWUKICAgIHRoZW1lOiByZWFkYWJsZQogICAgdG9jOiB0cnVlCiAgICB0b2NfZmxvYXQ6CiAgICAgIGNvbGxhcHNlZDogZmFsc2UKICAgICAgc21vb3RoX3Njcm9sbDogZmFsc2UKICBybWRmb3JtYXRzOjpyZWFkdGhlZG93bjoKICAgIGNvZGVfZG93bmxvYWQ6IHRydWUKICAgIGNvZGVfZm9sZGluZzogc2hvdwogICAgZGZfcHJpbnQ6IHBhZ2VkCiAgICBmaWdfY2FwdGlvbjogdHJ1ZQogICAgZmlnX2hlaWdodDogNwogICAgZmlnX3dpZHRoOiA3CiAgICBoaWdobGlnaHQ6IHRhbmdvCiAgICB3aWR0aDogMzAwCiAgICBrZWVwX21kOiBmYWxzZQogICAgbW9kZTogc2VsZmNvbnRhaW5lZAogICAgdG9jX2Zsb2F0OiB0cnVlCiAgQmlvY1N0eWxlOjpodG1sX2RvY3VtZW50OgogICAgY29kZV9kb3dubG9hZDogdHJ1ZQogICAgY29kZV9mb2xkaW5nOiBzaG93CiAgICBmaWdfY2FwdGlvbjogdHJ1ZQogICAgZmlnX2hlaWdodDogNwogICAgZmlnX3dpZHRoOiA3CiAgICBoaWdobGlnaHQ6IHRhbmdvCiAgICBrZWVwX21kOiBmYWxzZQogICAgbW9kZTogc2VsZmNvbnRhaW5lZAogICAgdG9jX2Zsb2F0OiB0cnVlCi0tLQoKPHN0eWxlIHR5cGU9InRleHQvY3NzIj4KYm9keSwgdGQgewogIGZvbnQtc2l6ZTogMTZweDsKfQpjb2RlLnJ7CiAgZm9udC1zaXplOiAxNnB4Owp9CnByZSB7CiBmb250LXNpemU6IDE2cHgKfQo8L3N0eWxlPgoKYGBge3Igb3B0aW9ucywgaW5jbHVkZT1GQUxTRX0KbGlicmFyeShocGdsdG9vbHMpCnR0IDwtIHNtKGRldnRvb2xzOjpsb2FkX2FsbCgifi9ocGdsdG9vbHMiKSkKa25pdHI6Om9wdHNfa25pdCRzZXQocHJvZ3Jlc3M9VFJVRSwKICAgICAgICAgICAgICAgICAgICAgdmVyYm9zZT1UUlVFLAogICAgICAgICAgICAgICAgICAgICB3aWR0aD0xMjAsCiAgICAgICAgICAgICAgICAgICAgIGVjaG89VFJVRSkKa25pdHI6Om9wdHNfY2h1bmskc2V0KGVycm9yPVRSVUUsCiAgICAgICAgICAgICAgICAgICAgICBmaWcud2lkdGg9OCwKICAgICAgICAgICAgICAgICAgICAgIGZpZy5oZWlnaHQ9OCwKICAgICAgICAgICAgICAgICAgICAgIGRwaT05NikKb2xkX29wdGlvbnMgPC0gb3B0aW9ucyhkaWdpdHM9NCwKICAgICAgICAgICAgICAgICAgICAgICBzdHJpbmdzQXNGYWN0b3JzPUZBTFNFLAogICAgICAgICAgICAgICAgICAgICAgIGtuaXRyLmR1cGxpY2F0ZS5sYWJlbD0iYWxsb3ciKQpnZ3Bsb3QyOjp0aGVtZV9zZXQoZ2dwbG90Mjo6dGhlbWVfYncoYmFzZV9zaXplPTEyKSkKdmVyIDwtICIyMDE5MDQxNyIKcnVuZGF0ZSA8LSBmb3JtYXQoU3lzLkRhdGUoKSwgZm9ybWF0PSIlWSVtJWQiKQpybWRfZmlsZSA8LSAiMjAxOTA0MTdfSHNNbV9NMExtNGguUm1kIgpgYGAKClRoaXMgZG9jdW1lbnQgd2lsbCBmaXJzdCBleHBsb3JlIGRpZmZlcmVudGlhbGx5IGV4cHJlc3NlZCBnZW5lcyBpbiBodW1hbnMgNApob3VycyBhZnRlciBpbmZlY3Rpb24gZm9sbG93ZWQgYnkgdGhlIHNhbWUgcXVlc3Rpb24gaW4gbWljZS4KCiMgV2hpY2ggZ2VuZXMgYXJlIERFIGluIGh1bWFuIG1hY3JvcGhhZ2VzIGF0IDQgaG91cnMgdXBvbiBpbmZlY3Rpb24gd2l0aCBMLiBtYWpvcj8KCiMgR2F0aGVyIGFubm90YXRpb24gZGF0YQoKSSB3YW50IHRvIHBlcmZvcm0gYSBzZXJpZXMgb2YgY29tcGFyaXNvbnMgYW1vbmcgdGhlIGhvc3QgY2VsbHM6IGh1bWFuIGFuZCBtb3VzZS4KVGh1cyBJIG5lZWQgdG8gY29sbGVjdCBhbm5vdGF0aW9uIGRhdGEgZm9yIGJvdGggc3BlY2llcyBhbmQgZ2V0IHRoZSBzZXQgb2YKb3J0aG9sb2dzIGJldHdlZW4gdGhlbS4KCiMjIFN0YXJ0IHdpdGggdGhlIGh1bWFuIGFubm90YXRpb24gZGF0YQoKSW4gdGhlIGZvbGxvd2luZyBibG9jaywgSSBkb3dubG9hZCB0aGUgaHVtYW4gYW5ub3RhdGlvbnMgZnJvbSBiaW9tYXJ0LiBJbgphZGRpdGlvbiwgSSB0YWtlIGEgbW9tZW50IHRvIHJlY3JlYXRlIHRoZSB0cmFuc2NyaXB0IElEcyBhcyBvYnNlcnZlZCBpbiB0aGUKc2FsbW9uIGNvdW50IHRhYmxlcyAoeWVzLCBJIGtub3cgdGhleSBhcmUgbm90IGFjdHVhbGx5IGNvdW50IHRhYmxlcykuICBGaW5hbGx5LApJIGNyZWF0ZSBhIHRhYmxlIHdoaWNoIG1hcHMgdHJhbnNjcmlwdHMgdG8gZ2VuZXMsIHRoaXMgd2lsbCBiZSB1c2VkIHdoZW4gd2UKZ2VuZXJhdGUgdGhlIGV4cHJlc3Npb25zZXQgc28gdGhhdCB3ZSBnZXQgZ2VuZSBleHByZXNzaW9uIGxldmVscyBmcm9tCnRyYW5zY3JpcHRzIHZpYSB0aGUgUiBwYWNrYWdlICd0eGltcG9ydCcuCgpgYGB7ciBodW1hbl9hbm5vdGF0aW9uc30KaHNfYW5ub3QgPC0gbG9hZF9iaW9tYXJ0X2Fubm90YXRpb25zKCkkYW5ub3RhdGlvbgpyb3duYW1lcyhoc19hbm5vdCkgPC0gbWFrZS5uYW1lcygKICBwYXN0ZTAoaHNfYW5ub3RbWyJlbnNlbWJsX3RyYW5zY3JpcHRfaWQiXV0sICIuIiwKICAgICAgICAgaHNfYW5ub3RbWyJ0cmFuc2NyaXB0X3ZlcnNpb24iXV0pLAogIHVuaXF1ZT1UUlVFKQpoc190eF9nZW5lIDwtIGhzX2Fubm90WywgYygiZW5zZW1ibF9nZW5lX2lkIiwgImVuc2VtYmxfdHJhbnNjcmlwdF9pZCIpXQpoc190eF9nZW5lW1siaWQiXV0gPC0gcm93bmFtZXMoaHNfdHhfZ2VuZSkKaHNfdHhfZ2VuZSA8LSBoc190eF9nZW5lWywgYygiaWQiLCAiZW5zZW1ibF9nZW5lX2lkIildCm5ld19oc19hbm5vdCA8LSBoc19hbm5vdApyb3duYW1lcyhuZXdfaHNfYW5ub3QpIDwtIG1ha2UubmFtZXMoaHNfYW5ub3RbWyJlbnNlbWJsX2dlbmVfaWQiXV0sIHVuaXF1ZT1UUlVFKQpgYGAKCiMjIEdlbmVyYXRlIGV4cHJlc3Npb25zZXRzCgpUaGUgcXVlc3Rpb24gaXMgcmVhc29uYWJseSBzZWxmLWNvbnRhaW5lZC4gIEkgd2FudCB0byBjb21wYXJlIHRoZSB1bmluZmVjdGVkCmh1bWFuIHNhbXBsZXMgYWdhaW5zdCBhbnkgc2FtcGxlcyB3aGljaCB3ZXJlIGluZmVjdGVkIGZvciA0IGhvdXJzLgpTbyBsZXQgdXMgZmlyc3QgcHVsbCB0aG9zZSBzYW1wbGVzIGFuZCB0aGVuIHBva2UgYXQgdGhlbSBhIGJpdC4KClRoZSBmb2xsb3dpbmcgYmxvY2sgY3JlYXRlcyBhbiBleHByZXNzaW9uc2V0IHVzaW5nIGFsbCBodW1hbi1xdWFudGlmaWVkCnNhbXBsZXMuIEFzIG1lbnRpb25lZCBwcmV2aW91c2x5LCBpdCB1c2VzIHRoZSB0YWJsZSBvZiB0cmFuc2NyaXB0PC0+Z2VuZQptYXBwaW5ncywgYW5kIHRoZSBiaW9tYXJ0IGFubm90YXRpb25zLgoKR2l2ZW4gdGhpcyBzZXQgb2YgfjQ0MCBzYW1wbGVzLCBpdCB0aGVuIGRyb3BzIHRoZSBmb2xsb3dpbmc6CgoxLiAgQWxsIHNhbXBsZXMgbWFya2VkICdza2lwcGVkJy4KMi4gIEFsbCBzYW1wbGVzIHdoaWNoIGFyZSBub3QgZnJvbSB0aW1lICd0NGgnLgoKYW5kIHJlc2V0cyB0aGUgY29uZGl0aW9uIGFuZCBiYXRjaCBmYWN0b3JzIHRvIHRoZSAnaW5mZWN0aW9uIHN0YXRlJyBtZXRhZGF0dW0KYW5kICdzdHVkeScsIHJlc3BlY3RpdmVseS4KCmBgYHtyIGV4cHRzfQpzYW1wbGVfc2hlZXQgPC0gInNhbXBsZV9zaGVldHMvbGVpc2htYW5pYV9ob3N0X21ldGFzaGVldF8yMDE5MDQwMS54bHN4Igpoc19leHB0IDwtIGNyZWF0ZV9leHB0KHNhbXBsZV9zaGVldCwKICAgICAgICAgICAgICAgICAgICAgICBmaWxlX2NvbHVtbj0iaHNhcGllbnNmaWxlIiwKICAgICAgICAgICAgICAgICAgICAgICBnZW5lX2luZm89bmV3X2hzX2Fubm90LAogICAgICAgICAgICAgICAgICAgICAgIHR4X2dlbmVfbWFwPWhzX3R4X2dlbmUpCgpoc19leHB0X25vc2tpcHBlZCA8LSBzdWJzZXRfZXhwdChoc19leHB0LCBzdWJzZXQ9InNraXBwZWQhPSd5ZXMnIikKaHNfdDRoX2V4cHQgPC0gc3Vic2V0X2V4cHQoaHNfZXhwdF9ub3NraXBwZWQsIHN1YnNldD0iZXhwdHRpbWU9PSd0NGgnIikKaHNfdDRoX2V4cHQgPC0gc2V0X2V4cHRfY29uZGl0aW9ucyhoc190NGhfZXhwdCwgZmFjdD0iaW5mZWN0c3RhdGUiKQpoc190NGhfZXhwdCA8LSBzZXRfZXhwdF9iYXRjaGVzKGhzX3Q0aF9leHB0LCBmYWN0PSJzdHVkeSIpCnRhYmxlKGhzX3Q0aF9leHB0JGNvbmRpdGlvbnMpCnRhYmxlKGhzX3Q0aF9leHB0JGJhdGNoZXMpCmhzX3dyaXR0ZW4gPC0gd3JpdGVfZXhwdChoc190NGhfZXhwdCwgZXhjZWw9ImV4Y2VsL0hzTTBMbTRoX2V4cHQueGxzeCIpCmBgYAoKIyMgRXhhbWluZSB0NGggdnMgdW5pbmZlY3RlZAoKTGV0IHVzIHBlcmZvcm0gc29tZSBnZW5lcmljIG1ldHJpY3Mgb2YgdGhlIHQ0aCBodW1hbiBleHByZXNzaW9uc2V0LiAgQXMgcGVyCnVzdWFsLCBJIHBsb3QgdGhlIG1ldHJpY3MgZmlyc3Qgb2YgdGhlIHJhdyBkYXRhOyBmb2xsb3dlZCBieSB0aGUgc2FtZSBtZXRyaWNzIG9mCmxvZzIocXVhbnRpbGUoY3BtKHN2YShmaWx0ZXJlZChkYXRhKSkpKSkuCgpgYGB7ciBoc19leGFtaW5lX3Q0aCwgZmlnLnNob3c9J2hpZGUnfQpoc190NGhfcGxvdHMgPC0gc20oZ3JhcGhfbWV0cmljcyhoc190NGhfZXhwdCkpCgpoc190NGhfbm9ybSA8LSBub3JtYWxpemVfZXhwdChoc190NGhfZXhwdCwgbm9ybT0icXVhbnQiLCBjb252ZXJ0PSJjcG0iLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICB0cmFuc2Zvcm09ImxvZzIiLCBmaWx0ZXI9VFJVRSwgYmF0Y2g9InN2YXNlcSIpCmhzX3Q0aF9ub3JtX3Bsb3RzIDwtIHNtKGdyYXBoX21ldHJpY3MoaHNfdDRoX25vcm0pKQpgYGAKCiMjIyBQcmludCBzb21lIG9mIHRoZSBwbG90cwoKYGBge3IgaHNfZXhhbWluZV9wbG90c30KaHNfdDRoX3Bsb3RzJGxlZ2VuZApoc190NGhfcGxvdHMkbGlic2l6ZQpoc190NGhfcGxvdHMkYm94cGxvdAoKaHNfdDRoX25vcm1fcGxvdHMkcGNfcGxvdApgYGAKCiMjIFJlbW92ZSBzdGltdWxhdGVkIHNhbXBsZXMKCkkgcGVyaGFwcyBzaG91bGQgaGF2ZSByZW1vdmVkIHRoZSBzdGltdWxhdGVkIHNhbXBsZXMgc29vbmVyLCBidXQgSSB3YXMgY3VyaW91cwp0byBzZWUgdGhlaXIgZWZmZWN0IG9uIHRoZSBkaXN0cmlidXRpb24gZmlyc3QuCgpgYGB7ciBoc190NGhfbm91bnN0aW19CmhzX3Q0aF9pbmYgPC0gc3Vic2V0X2V4cHQoaHNfdDRoX2V4cHQsIHN1YnNldD0iY29uZGl0aW9uIT0nc3RpbSciKQpoc190NGhfaW5mX25vcm0gPC0gbm9ybWFsaXplX2V4cHQoaHNfdDRoX2luZiwgdHJhbnNmb3JtPSJsb2cyIiwgY29udmVydD0iY3BtIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGZpbHRlcj1UUlVFLCBiYXRjaD0ic3Zhc2VxIikKCmhzX3Q0aF9wY2EgPC0gcGxvdF9wY2EoaHNfdDRoX2luZl9ub3JtLCBwbG90X3RpdGxlPSJILiBzYXBpZW5zLCBMLiBtYWpvciwgdDRoIikKaHNfdDRoX3BjYSRwbG90CgprZWVwZXJzIDwtIGxpc3QoImluZmVjdGlvbiIgPSBjKCJ5ZXMiLCAibm8iKSkKaHNfdDRoX2RlIDwtIGFsbF9wYWlyd2lzZShoc190NGhfaW5mLCBtb2RlbF9iYXRjaD0ic3Zhc2VxIiwgZmlsdGVyPVRSVUUsIGZvcmNlPVRSVUUpCmhzX3Q0aF90YWJsZSA8LSBjb21iaW5lX2RlX3RhYmxlcyhoc190NGhfZGUsIGtlZXBlcnM9a2VlcGVycywKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGV4Y2VsPSJleGNlbC9Ic00wTG00aF9kZV90YWJsZXMueGxzeCIpCmhzX3Q0aF9zaWcgPC0gZXh0cmFjdF9zaWduaWZpY2FudF9nZW5lcyhoc190NGhfdGFibGUsIGV4Y2VsPSJleGNlbC9Ic00wTG00aF9zaWdfdGFibGVzLnhsc3giKQpgYGAKCiMgV2hpY2ggZ2VuZXMgYXJlIERFIGluIG1vdXNlIG1hY3JvcGhhZ2VzIGF0IDQgaG91cnMgdXBvbiBpbmZlY3Rpb24gd2l0aCBMLiBtYWpvcj8KCk1vc3Qgb2YgdGhpcyBzaG91bGQgYmUgdGhlIHNhbWUgaW4gcHJvY2VzcyBhcyB3aGF0IHdhcyBwZXJmb3JtZWQgZm9yIHRoZSBodW1hbi4KCiMjIEdhdGhlciBhbm5vdGF0aW9uIGRhdGEKCkkgd2FudCB0byBwZXJmb3JtIGEgc2VyaWVzIG9mIGNvbXBhcmlzb25zIGFtb25nIHRoZSBob3N0IGNlbGxzOiBodW1hbiBhbmQgbW91c2UuClRodXMgSSBuZWVkIHRvIGNvbGxlY3QgYW5ub3RhdGlvbiBkYXRhIGZvciBib3RoIHNwZWNpZXMgYW5kIGdldCB0aGUgc2V0IG9mCm9ydGhvbG9ncyBiZXR3ZWVuIHRoZW0uCgojIyMgU3RhcnQgd2l0aCB0aGUgaHVtYW4gYW5ub3RhdGlvbiBkYXRhCgpgYGB7ciBtb3VzZV9hbm5vdGF0aW9uc30KbW1fYW5ub3QgPC0gbG9hZF9iaW9tYXJ0X2Fubm90YXRpb25zKHNwZWNpZXM9Im1tdXNjdWx1cyIpJGFubm90YXRpb24Kcm93bmFtZXMobW1fYW5ub3QpIDwtIG1ha2UubmFtZXMoCiAgcGFzdGUwKG1tX2Fubm90W1siZW5zZW1ibF90cmFuc2NyaXB0X2lkIl1dLCAiLiIsCiAgICAgICAgIG1tX2Fubm90W1sidHJhbnNjcmlwdF92ZXJzaW9uIl1dKSwKICB1bmlxdWU9VFJVRSkKbW1fdHhfZ2VuZSA8LSBtbV9hbm5vdFssIGMoImVuc2VtYmxfZ2VuZV9pZCIsICJlbnNlbWJsX3RyYW5zY3JpcHRfaWQiKV0KbW1fdHhfZ2VuZVtbImlkIl1dIDwtIHJvd25hbWVzKG1tX3R4X2dlbmUpCm1tX3R4X2dlbmUgPC0gbW1fdHhfZ2VuZVssIGMoImlkIiwgImVuc2VtYmxfZ2VuZV9pZCIpXQpuZXdfbW1fYW5ub3QgPC0gbW1fYW5ub3QKcm93bmFtZXMobmV3X21tX2Fubm90KSA8LSBtYWtlLm5hbWVzKG1tX2Fubm90W1siZW5zZW1ibF9nZW5lX2lkIl1dLCB1bmlxdWU9VFJVRSkKYGBgCgojIyBHZW5lcmF0ZSBleHByZXNzaW9uc2V0cwoKVGhlIHF1ZXN0aW9uIGlzIHJlYXNvbmFibHkgc2VsZi1jb250YWluZWQuICBJIHdhbnQgdG8gY29tcGFyZSB0aGUgdW5pbmZlY3RlZApodW1hbiBzYW1wbGVzIGFnYWluc3QgYW55IHNhbXBsZXMgd2hpY2ggd2VyZSBpbmZlY3RlZCBmb3IgNCBob3Vycy4KU28gbGV0IHVzIGZpcnN0IHB1bGwgdGhvc2Ugc2FtcGxlcyBhbmQgdGhlbiBwb2tlIGF0IHRoZW0gYSBiaXQuCgpgYGB7ciBtb3VzZV9leHB0c30KbW1fZXhwdCA8LSBjcmVhdGVfZXhwdChzYW1wbGVfc2hlZXQsCiAgICAgICAgICAgICAgICAgICAgICAgZmlsZV9jb2x1bW49Im1tdXNjdWx1c2ZpbGUiLAogICAgICAgICAgICAgICAgICAgICAgIGdlbmVfaW5mbz1uZXdfbW1fYW5ub3QsCiAgICAgICAgICAgICAgICAgICAgICAgdHhfZ2VuZV9tYXA9bW1fdHhfZ2VuZSkKbW1fdDRoX2V4cHQgPC0gc3Vic2V0X2V4cHQobW1fZXhwdCwgc3Vic2V0PSJleHB0dGltZT09J3Q0aCciKQptbV90NGhfZXhwdCA8LSBzZXRfZXhwdF9jb25kaXRpb25zKG1tX3Q0aF9leHB0LCBmYWN0PSJpbmZlY3RzdGF0ZSIpCnRhYmxlKG1tX3Q0aF9leHB0JGNvbmRpdGlvbnMpCnRhYmxlKG1tX3Q0aF9leHB0JGJhdGNoZXMpCm1tX3dyaXR0ZW4gPC0gd3JpdGVfZXhwdChtbV90NGhfZXhwdCwgZXhjZWw9ImV4Y2VsL01tTTBMbTRoX2V4cHQueGxzeCIpCmBgYAoKIyMgRXhhbWluZSB0NGggdnMgdW5pbmZlY3RlZAoKYGBge3IgZXhhbWluZV90NGgsIGZpZy5zaG93PSdoaWRlJ30KbW1fdDRoX3Bsb3RzIDwtIHNtKGdyYXBoX21ldHJpY3MobW1fdDRoX2V4cHQpKQptbV90NGhfbm9ybSA8LSBub3JtYWxpemVfZXhwdChtbV90NGhfZXhwdCwgbm9ybT0icXVhbnQiLCBjb252ZXJ0PSJjcG0iLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICB0cmFuc2Zvcm09ImxvZzIiLCBmaWx0ZXI9VFJVRSwgYmF0Y2g9InN2YXNlcSIpCm1tX3Q0aF9ub3JtX3Bsb3RzIDwtIHNtKGdyYXBoX21ldHJpY3MobW1fdDRoX25vcm0pKQpgYGAKCiMjIyBQcmludCBzb21lIG9mIHRoZSBwbG90cwoKYGBge3IgbW1fZXhhbWluZV9wbG90c30KbW1fdDRoX3Bsb3RzJGxlZ2VuZAptbV90NGhfcGxvdHMkbGlic2l6ZQptbV90NGhfcGxvdHMkYm94cGxvdAoKbW1fdDRoX25vcm1fcGxvdHMkcGNfcGxvdApgYGAKCiMjIFBlcmZvcm0gZGUgYW5hbHlzZXMKCmBgYHtyIG1tX3Q0aF9kZX0KbW1fdDRoX2luZl9ub3JtIDwtIG5vcm1hbGl6ZV9leHB0KG1tX3Q0aF9leHB0LCB0cmFuc2Zvcm09ImxvZzIiLCBjb252ZXJ0PSJjcG0iLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZmlsdGVyPVRSVUUsIGJhdGNoPSJzdmFzZXEiKQoKbW1fdDRoX3BjYSA8LSBwbG90X3BjYShtbV90NGhfaW5mX25vcm0sIHBsb3RfdGl0bGU9Ik0uIG11c2N1bHVzLCBMLiBtYWpvciwgdDRoIikKbW1fdDRoX3BjYSRwbG90CgptbV90NGhfZGUgPC0gYWxsX3BhaXJ3aXNlKG1tX3Q0aF9leHB0LCBtb2RlbF9iYXRjaD0ic3Zhc2VxIiwgZmlsdGVyPVRSVUUsIGZvcmNlPVRSVUUpCm1tX3Q0aF90YWJsZSA8LSBjb21iaW5lX2RlX3RhYmxlcyhtbV90NGhfZGUsIGtlZXBlcnM9a2VlcGVycywKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGV4Y2VsPSJleGNlbC9NbU0wTG00aF9kZV90YWJsZXMueGxzeCIpCm1tX3Q0aF9zaWcgPC0gZXh0cmFjdF9zaWduaWZpY2FudF9nZW5lcyhtbV90NGhfdGFibGUsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBleGNlbD0iZXhjZWwvTW1NMExtNGhfc2lnX3RhYmxlcy54bHN4IikKYGBgCgojIyBDb21wYXJlIHRoaXMgdG8gdGhlIHByZXZpb3VzIHJlc3VsdC4KCkxldCB1cyBzZWUgaWYgb3VyIGh1bWFuIGRpZmZlcmVudGlhbCBleHByZXNzaW9uIHJlc3VsdCBpcyBzaW1pbGFyIHRvIHRoYXQKb2J0YWluZWQgaW4gVGFibGUgUzIuCgpgYGB7ciBjb21wYXJlX3ByZXZpb3VzX2hzfQpwcmV2aW91c19ocyA8LSByZWFkeGw6OnJlYWRfZXhjZWwoImV4Y2VsL2lubGluZS1zdXBwbGVtZW50YXJ5LW1hdGVyaWFsLTUueGxzIiwgc2hlZXQ9MikKcHJldmlvdXNfaHNfbGZjIDwtIHByZXZpb3VzX2hzWywgYygiSUQiLCAiRm9sZCBjaGFuZ2UiKV0KbmVnX2lkeCA8LSBwcmV2aW91c19oc19sZmNbWzJdXSA8IDAKcHJldmlvdXNfaHNfbGZjW25lZ19pZHgsIDJdIDwtIC0xICogKDEgLyBwcmV2aW91c19oc19sZmNbbmVnX2lkeCwgMl0pCnByZXZpb3VzX2hzX2xmY1tbMl1dIDwtIGxvZzIocHJldmlvdXNfaHNfbGZjW1syXV0pCgptZXJnZWQgPC0gbWVyZ2UocHJldmlvdXNfaHNfbGZjLCBoc190NGhfdGFibGUkZGF0YVtbMV1dLCBieS54PSJJRCIsIGJ5Lnk9InJvdy5uYW1lcyIpCmNvci50ZXN0KG1lcmdlZFtbImxpbW1hX2xvZ2ZjIl1dLCBtZXJnZWRbWyJGb2xkIGNoYW5nZSJdXSkKYGBgCgojIyBDb21wYXJlIHRoZSBwcmV2aW91cyBtb3VzZSBhbmQgdGhlc2UgbW91c2UgcmVzdWx0cy4KCmBgYHtyIGNvbXBhcmVfcHJldmlvdXNfbW19CnByZXZpb3VzX21tIDwtIHJlYWR4bDo6cmVhZF9leGNlbCgiZXhjZWwvMTI4NjRfMjAxNV8yMjM3X01PRVNNM19FU00ueGxzIiwgc2hlZXQ9Miwgc2tpcD0xKQpwcmV2aW91c19tbV9sZmMgPC0gcHJldmlvdXNfbW1bLCBjKCJJRCIsICJGb2xkIGNoYW5nZSIpXQpuZWdfaWR4IDwtIHByZXZpb3VzX21tX2xmY1tbMl1dIDwgMApwcmV2aW91c19tbV9sZmNbbmVnX2lkeCwgMl0gPC0gLTEgKiAoMSAvIHByZXZpb3VzX21tX2xmY1tuZWdfaWR4LCAyXSkKcHJldmlvdXNfbW1fbGZjW1syXV0gPC0gbG9nMihwcmV2aW91c19tbV9sZmNbWzJdXSkKCm1lcmdlZCA8LSBtZXJnZShwcmV2aW91c19tbV9sZmMsIG1tX3Q0aF90YWJsZSRkYXRhW1sxXV0sIGJ5Lng9IklEIiwgYnkueT0icm93Lm5hbWVzIikKY29yLnRlc3QobWVyZ2VkW1sibGltbWFfbG9nZmMiXV0sIG1lcmdlZFtbIkZvbGQgY2hhbmdlIl1dKQpgYGAKCiMgV2hhdCBnZW5lcyBhcmUgc2hhcmVkIGluIHRoZSBtb3VzZSBhbmQgaHVtYW4gZGF0YT8KClRoaXMgaXMgb25lIG1ldGhvZCBvZiBhZGRyZXNzaW5nIE5hamlicyBiaWcgcXVlc3Rpb24gIzFjLiAgSW4gdGhpcyBtZXRob2QsIEkgYW0KdGFraW5nIHRoZSB0YWJsZXMgZm9yIHRoZSBodW1hbiBhbmFseXNpcyBhbmQgbW91c2UgYW5hbHlzaXMgc2VwYXJhdGVseSwgdGhlbgptZXJnaW5nIHRoZW0gdXNpbmcgYSB0YWJsZSBvZiBvcnRob2xvZ3MuCgpTaWRlIG5vdGU6IGEgZGlmZmVyZW50IHdheSBvZiBhZGRyZXNzaW5nIHRoaXMgcXVlc3Rpb24gcmVzaWRlcyBpbgoyMDE5MDIyMF9ob3N0X2NvbXBhcmlzb25zLlJtZC4gSW4gdGhpcyBjb21wZXRpbmcgbWV0aG9kLCB0aGUgdGFibGUgb2Ygb3J0aG9sb2dzCmlzIHVzZWQgaW4gdGhlIGJlZ2lubmluZyB0byBtYWtlIGEgc2luZ2xlIHNldCBvZiBJRHMgZm9yIHRoZSBodW1hbiBhbmQgbW91c2UKZ2VuZXMsIHRoZW4gcGVyZm9ybSB0aGUgZGlmZmVyZW50aWFsIGV4cHJlc3Npb24gYW5hbHlzaXMuCgojIyBFeHRyYWN0IGh1bWFuIG1vdXNlIG9ydGhvbG9ncwoKTXkgbG9hZF9iaW9tYXJ0X29ydGhvbG9ncygpIGZ1bmN0aW9uIHNob3VsZCBwcm92aWRlIHRoaXMgbWFwcGluZyBnZW5lIElEIHRhYmxlLgoKYGBge3IgbW1faHNfb3J0aG9sb2d9CiMjIFRoZSBkZWZhdWx0cyBvZiB0aGlzIGZ1bmN0aW9uIGFyZSBzdWl0YWJsZSBmb3IgbW91c2UvaHVtYW4gcXVlcmllcy4KbW1faHNfb3J0aG8gPC0gbG9hZF9iaW9tYXJ0X29ydGhvbG9ncygpJGFsbF9saW5rZWRfZ2VuZXMKCm1tX3RhYmxlIDwtIG1tX3Q0aF90YWJsZSRkYXRhW1sxXV0KaHNfdGFibGUgPC0gaHNfdDRoX3RhYmxlJGRhdGFbWzFdXQoKbW1fdGFibGUgPC0gbWVyZ2UobW1faHNfb3J0aG8sIG1tX3RhYmxlLCBieS54PSJtbXVzY3VsdXMiLCBieS55PSJyb3cubmFtZXMiLCBhbGwueT1UUlVFKQpoc190YWJsZSA8LSBtZXJnZShtbV9oc19vcnRobywgaHNfdGFibGUsIGJ5Lng9ImhzYXBpZW5zIiwgYnkueT0icm93Lm5hbWVzIiwgYWxsLnk9VFJVRSkKYm90aF90YWJsZV9ocyA8LSBtZXJnZShoc190YWJsZSwgbW1fdGFibGUsIGJ5Lng9ImhzYXBpZW5zIiwgYnkueT0iaHNhcGllbnMiKQpib3RoX3RhYmxlX21tIDwtIG1lcmdlKGhzX3RhYmxlLCBtbV90YWJsZSwgYnkueD0ibW11c2N1bHVzIiwgYnkueT0ibW11c2N1bHVzIikKCmNvci50ZXN0KGJvdGhfdGFibGVfaHNbWyJsaW1tYV9sb2dmYy54Il1dLCBib3RoX3RhYmxlX2hzW1sibGltbWFfbG9nZmMueSJdXSkKY29yLnRlc3QoYm90aF90YWJsZV9tbVtbImxpbW1hX2xvZ2ZjLngiXV0sIGJvdGhfdGFibGVfbW1bWyJsaW1tYV9sb2dmYy55Il1dKQp0dCA8LSBwbG90X3NjYXR0ZXIoYm90aF90YWJsZV9oc1ssIGMoImxpbW1hX2xvZ2ZjLngiLCAibGltbWFfbG9nZmMueSIpXSkKdHQKYGBgCgpJIGJlbGlldmUgdGhlc2UgYm90aF90YWJsZV9ocyBhbmQgYm90aF90YWJsZV9tbSB0YWJsZXMgYXJlIGdvb2QgY2FuZGlkYXRlcyBmb3IKdGhlIHNldCBvZiBnZW5lcyB3aGljaCBhcmUgc2hhcmVkIGFjcm9zcyB0aGUgaHVtYW4gYW5kIG1vdXNlIHNhbXBsZXMsIGZyb20gdGhlCnBlcnNwZWN0aXZlIG9mIHRoZSBodW1hbi4KClNlcGFyYXRlbHkKCk5vdyBsZXRzIHdyaXRlIHNvbWUgb2YgdGhlc2Ugb3V0LgoKYGBge3Igd3JpdGVfc2hhcmVkX3RhYmxlc30Kc2hhcmVkX2hzbW1fdXBfaWR4IDwtIGJvdGhfdGFibGVfaHNbWyJkZXNlcV9sb2dmYy54Il1dID4gMSAmCiAgYm90aF90YWJsZV9oc1tbImRlc2VxX2xvZ2ZjLnkiXV0gPiAxICYKICBib3RoX3RhYmxlX2hzW1siZGVzZXFfYWRqcC54Il1dIDw9IDAuMDUgJgogIGJvdGhfdGFibGVfaHNbWyJkZXNlcV9hZGpwLnkiXV0gPD0gMC4wNQpzaGFyZWRfaHNtbV91cCA8LSBib3RoX3RhYmxlX2hzW3NoYXJlZF9oc21tX3VwX2lkeCwgXQp3cml0dGVuIDwtIHdyaXRlX3hscyhkYXRhPXNoYXJlZF9oc21tX3VwLAogICAgICAgICAgICAgICAgICAgICBleGNlbD0iZXhjZWwvSHNNMExtNGhfdnNfTW1NMExtNGhfc2hhcmVkX3VwX2hzLnhsc3giKQpzaGFyZWRfbW1oc191cF9pZHggPC0gYm90aF90YWJsZV9tbVtbImRlc2VxX2xvZ2ZjLngiXV0gPiAxICYKICBib3RoX3RhYmxlX21tW1siZGVzZXFfbG9nZmMueSJdXSA+IDEgJgogIGJvdGhfdGFibGVfbW1bWyJkZXNlcV9hZGpwLngiXV0gPD0gMC4wNSAmCiAgYm90aF90YWJsZV9tbVtbImRlc2VxX2FkanAueSJdXSA8PSAwLjA1CnNoYXJlZF9tbWhzX3VwIDwtIGJvdGhfdGFibGVfbW1bc2hhcmVkX21taHNfdXBfaWR4LCBdCndyaXR0ZW4gPC0gd3JpdGVfeGxzKGRhdGE9c2hhcmVkX21taHNfdXAsCiAgICAgICAgICAgICAgICAgICAgIGV4Y2VsPSJleGNlbC9Ic00wTG00aF92c19NbU0wTG00aF9zaGFyZWRfdXBfbW0ueGxzeCIpCgpzaGFyZWRfaHNtbV9kb3duX2lkeCA8LSBib3RoX3RhYmxlX2hzW1siZGVzZXFfbG9nZmMueCJdXSA8IC0xICYKICBib3RoX3RhYmxlX2hzW1siZGVzZXFfbG9nZmMueSJdXSA8IC0xICYKICBib3RoX3RhYmxlX2hzW1siZGVzZXFfYWRqcC54Il1dIDw9IDAuMDUgJgogIGJvdGhfdGFibGVfaHNbWyJkZXNlcV9hZGpwLnkiXV0gPD0gMC4wNQpzaGFyZWRfaHNtbV9kb3duIDwtIGJvdGhfdGFibGVfaHNbc2hhcmVkX2hzbW1fZG93bl9pZHgsIF0Kd3JpdHRlbiA8LSB3cml0ZV94bHMoZGF0YT1zaGFyZWRfaHNtbV9kb3duLAogICAgICAgICAgICAgICAgICAgICBleGNlbD0iZXhjZWwvSHNNMExtNGhfdnNfTW1NMExtNGhfc2hhcmVkX2Rvd25faHMueGxzeCIpCnNoYXJlZF9tbWhzX2Rvd25faWR4IDwtIGJvdGhfdGFibGVfbW1bWyJkZXNlcV9sb2dmYy54Il1dIDwgLTEgJgogIGJvdGhfdGFibGVfbW1bWyJkZXNlcV9sb2dmYy55Il1dIDwgLTEgJgogIGJvdGhfdGFibGVfbW1bWyJkZXNlcV9hZGpwLngiXV0gPD0gMC4wNSAmCiAgYm90aF90YWJsZV9tbVtbImRlc2VxX2FkanAueSJdXSA8PSAwLjA1CnNoYXJlZF9tbWhzX2Rvd24gPC0gYm90aF90YWJsZV9tbVtzaGFyZWRfbW1oc19kb3duX2lkeCwgXQp3cml0dGVuIDwtIHdyaXRlX3hscyhkYXRhPXNoYXJlZF9tbWhzX2Rvd24sCiAgICAgICAgICAgICAgICAgICAgIGV4Y2VsPSJleGNlbC9Ic00wTG00aF92c19NbU0wTG00aF9zaGFyZWRfZG93bl9tbS54bHN4IikKYGBgCgojIyBDb21wYXJlIHRoZSBhYm92ZSB3aXRoIFRhYmxlIFM3IGZyb20gdGhlIExhdXJhIGFuZCBDZWNpbGlhIHBhcGVyLgoKVGFibGUgUzcgaGFzIGEgc2V0IG9mIGdlbmVzIGZyb20gaHVtYW4gd2hpY2ggYXJlIGFsc28gdXAtcmVndWxhdGVkIGluIG1vdXNlIHVwb24KaW5mZWN0aW9uLgoKYGBge3Igb3Blbl9maWdzN30KZmlnX3M3X3VwIDwtIHJlYWR4bDo6cmVhZF9leGNlbCgiZXhjZWwvaW5saW5lLXN1cHBsZW1lbnRhcnktbWF0ZXJpYWwtNy54bHMiLCBzaGVldD0zKQpmaWdfczdfaHNfdXAgPC0gdW5pcXVlKGZpZ19zN191cFtbNV1dKQoKZmlnX3M3X2Rvd24gPC0gcmVhZHhsOjpyZWFkX2V4Y2VsKCJleGNlbC9pbmxpbmUtc3VwcGxlbWVudGFyeS1tYXRlcmlhbC03LnhscyIsIHNoZWV0PTQpCmZpZ19zN19oc19kb3duIDwtIHVuaXF1ZShmaWdfczdfZG93bltbNV1dKQoKYm90aF91cF9pZHggPC0gYm90aF90YWJsZV9oc1tbImxpbW1hX2xvZ2ZjLngiXV0gPj0gMC44ICYKICBib3RoX3RhYmxlX2hzW1sibGltbWFfbG9nZmMueSJdXSA+PSAwLjggJgogIGJvdGhfdGFibGVfaHNbWyJsaW1tYV9hZGpwLngiXV0gPD0gMC4xICYKICBib3RoX3RhYmxlX2hzW1sibGltbWFfYWRqcC55Il1dIDw9IDAuMQpib3RoX3VwX2lkcyA8LSBib3RoX3RhYmxlX2hzW2JvdGhfdXBfaWR4LCAiaHNhcGllbnMiXQp1cF92ZW5uIDwtIFZlbm5lcmFibGU6OlZlbm4oU2V0cz1saXN0KCJmaWdzNyIgPSBmaWdfczdfaHNfdXAsICJ0YWJsZXMiID0gYm90aF91cF9pZHMpKQpWZW5uZXJhYmxlOjpwbG90KHVwX3Zlbm4pCgpib3RoX2Rvd25faWR4IDwtIGJvdGhfdGFibGVfaHNbWyJsaW1tYV9sb2dmYy54Il1dIDw9IC0wLjggJgogIGJvdGhfdGFibGVfaHNbWyJsaW1tYV9sb2dmYy55Il1dIDw9IC0wLjggJgogIGJvdGhfdGFibGVfaHNbWyJsaW1tYV9hZGpwLngiXV0gPD0gMC4xICYKICBib3RoX3RhYmxlX2hzW1sibGltbWFfYWRqcC55Il1dIDw9IDAuMQpib3RoX2Rvd25faWRzIDwtIGJvdGhfdGFibGVfaHNbYm90aF9kb3duX2lkeCwgImhzYXBpZW5zIl0KZG93bl92ZW5uIDwtIFZlbm5lcmFibGU6OlZlbm4oU2V0cz1saXN0KCJmaWdzNyIgPSBmaWdfczdfaHNfZG93biwgInRhYmxlcyIgPSBib3RoX2Rvd25faWRzKSkKVmVubmVyYWJsZTo6cGxvdChkb3duX3Zlbm4pCmBgYAoKIyBOZXh0IHF1ZXN0aW9uOiBXaGljaCA0IGhvdXIgZ2VuZXMgYXJlIHNoYXJlZCB3aXRoIENJREVJTT8KCiJXaGljaCA0IGhvdXIgREUgZ2VuZXMgYXJlIHNoYXJlZCB3aXRoIHRoZSBDSURFSU0gcGFuYW1lbnNpcyBpbmZlY3RlZCBodW1hbgptYWNyb3BoYWdlcz8iCgpTbywgb25jZSBhZ2FpbiB0aGVyZSBhcmUgdHdvIHdheXMgb2YgYXBwcm9hY2hpbmcgdGhpcyBxdWVzdGlvbjoKCjEuICBFeGFtaW5lIHRoZSBxdWVzdGlvbiBvZiBpbmZlY3RlZCB2cy4gdW5pbmZlY3RlZCBmb3IgdGhlIHR3byBkYXRhIHNldHMKICAgIHNlcGFyYXRlbHkuICBUaGVuIHF1ZXJ5IHRoZSByZXN1bHRzIGFuZCBzZWUgd2hhdCBjb21lcyBvdXQuCjIuICBFeGFtaW5lIHRoaXMgYXMgYSBzaW5nbGUgcXVlc3Rpb24gdXNpbmcgdGhlIHVuaW9uIG9mIGJvdGggZGF0YSBzZXRzLgoKU2luY2UgSSBhbHJlYWR5IGhhdmUgYSBwdXRhdGl2ZSBhbnN3ZXIgZm9yIHRoZSBodW1hbiA0IGhvdXIgbWFjcm9waGFnZSBkYXRhLCBUaGUKc2ltcGxlc3QgbWV0aG9kIGlzIHRvIGp1c3QgbG9vayBhdCB0aGUgY2lkZWltIGRhdGEgYW5kIGRvICMxIGFib3ZlLiAgU28gbGV0IHVzCmRvIHRoYXQgZmlyc3QuCgpgYGB7ciBjaWRlaW1fc2VwYXJhdGV9CmNpZGVpbV9tYWNyIDwtIHN1YnNldF9leHB0KGhzX2V4cHQsIHN1YnNldD0ibGFiPT0ndG1yYycmaG9zdGNlbGx0eXBlPT0nbWFjcm9waGFnZSciKQpjaWRlaW1fbWFjciA8LSBzZXRfZXhwdF9jb25kaXRpb25zKGNpZGVpbV9tYWNyLCBmYWN0PSJpbmZlY3RzdGF0ZSIpCmNpZGVpbV9ub3JtIDwtIG5vcm1hbGl6ZV9leHB0KGNpZGVpbV9tYWNyLCB0cmFuc2Zvcm09ImxvZzIiLCBjb252ZXJ0PSJjcG0iLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICBub3JtPSJxdWFudCIsIGZpbHRlcj0ic2ltcGxlIiwgYmF0Y2g9InN2YXNlcSIpCgpwbG90X3BjYShjaWRlaW1fbm9ybSkkcGxvdAoKY2lkZWltX2RlIDwtIGFsbF9wYWlyd2lzZShjaWRlaW1fbWFjciwgZmlsdGVyPVRSVUUsIHBhcmFsbGVsPUZBTFNFLCBmb3JjZT1UUlVFKQpjaWRlaW1fdGFibGUgPC0gY29tYmluZV9kZV90YWJsZXMoY2lkZWltX2RlLCBrZWVwZXJzPWtlZXBlcnMpCgpjaWRlaW1fbWVyZ2VkIDwtIG1lcmdlKGhzX3Q0aF90YWJsZSRkYXRhW1sxXV0sIGNpZGVpbV90YWJsZSRkYXRhW1sxXV0sIGJ5PSJyb3cubmFtZXMiKQpjb3IudGVzdChjaWRlaW1fbWVyZ2VkW1siZGVzZXFfbG9nZmMueCJdXSwgY2lkZWltX21lcmdlZFtbImRlc2VxX2xvZ2ZjLnkiXV0pCgpwbG90X3NjYXR0ZXIoY2lkZWltX21lcmdlZFssIGMoImxpbW1hX2xvZ2ZjLngiLCAibGltbWFfbG9nZmMueSIpXSkKYGBgCgojIyBUcnkgYWdhaW4gYXMgb25lIGJpZyBleHByZXNzaW9uc2V0CgpgYGB7ciBjaWRlaW1fYmlnX2V4cHJzfQpoc19leHB0X25vc2tpcHBlZCA8LSBzdWJzZXRfZXhwdChoc19leHB0LCBzdWJzZXQ9InNraXBwZWQhPSd5ZXMnIikKaHNfdDRoX2V4cHQgPC0gc3Vic2V0X2V4cHQoaHNfZXhwdF9ub3NraXBwZWQsIHN1YnNldD0iZXhwdHRpbWU9PSd0NGgnIikKaHNfdDRoX2V4cHQgPC0gc2V0X2V4cHRfY29uZGl0aW9ucyhoc190NGhfZXhwdCwgZmFjdD0iaW5mZWN0c3RhdGUiKQpjaWRlaW1fdDRoIDwtIHN1YnNldF9leHB0KAogIGhzX2V4cHQsCiAgc3Vic2V0PSIoc2tpcHBlZCE9J3llcycmZXhwdHRpbWU9PSd0NGgnKXwobGFiPT0ndG1yYycmaG9zdGNlbGx0eXBlPT0nbWFjcm9waGFnZScpIikKY2lkZWltX3Q0aCA8LSBzZXRfZXhwdF9jb25kaXRpb25zKGNpZGVpbV90NGgsIGZhY3Q9ImluZmVjdHN0YXRlIikKCmNpZGVpbV90NGhfbm9ybSA8LSBub3JtYWxpemVfZXhwdChjaWRlaW1fdDRoLCBmaWx0ZXI9VFJVRSwgbm9ybT0icXVhbnQiLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgY29udmVydD0iY3BtIiwgdHJhbnNmb3JtPSJsb2cyIikKcGxvdF9wY2EoY2lkZWltX3Q0aF9ub3JtKSRwbG90CmNpZGVpbV90NGhfbm9ybWJhdGNoIDwtIG5vcm1hbGl6ZV9leHB0KGNpZGVpbV90NGgsIGZpbHRlcj1UUlVFLCBub3JtPSJxdWFudCIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGNvbnZlcnQ9ImNwbSIsIGJhdGNoPSJzdmFzZXEiLCB0cmFuc2Zvcm09ImxvZzIiKQpwbG90X3BjYShjaWRlaW1fdDRoX25vcm1iYXRjaCkkcGxvdAoKY2lkZWltX3Q0aF9kZSA8LSBhbGxfcGFpcndpc2UoY2lkZWltX3Q0aCwgbW9kZWxfYmF0Y2g9InN2YXNlcSIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGZpbHRlcj1UUlVFLCBmb3JjZT1UUlVFLCBwYXJhbGxlbD1GQUxTRSkKY2lkZWltX3Q0aF90YWJsZXMgPC0gY29tYmluZV9kZV90YWJsZXMoY2lkZWltX3Q0aF9kZSwga2VlcGVycz1rZWVwZXJzLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBleGNlbD0iZXhjZWwvSHNNMExtNGhfdnNfSHNNMExwX2RlX3RhYmxlcy54bHN4IikKY2lkZWltX3Q0aF9zaWcgPC0gY29tYmluZV9kZV90YWJsZXMoY2lkZWltX3Q0aF90YWJsZXMsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGV4Y2VsPSJleGNlbC9Ic00wTG00aF92c19Ic00wTHBfc2lnX3RhYmxlcy54bHN4IikKYGBgCgojIE5leHQsIG92ZXJsYXAgYmV0d2VlbiBtYWNyb3BoYWdlcyBhbmQgbmV1dHJvcGhpbHMKCkFsbCBvZiB0aGUgbmV1dHJvcGhpbCBkYXRhIGlzIGluIG1vdXNlLCBhcHBhcmVudGx5LiAgVGhpcyB3aWxsIG1ha2UgaXQgbW9yZQpkaWZmaWN1bHQsIHBlcmhhcHMgaW1wb3NzaWJsZSB0byBnZXQgYW4gYWNjdXJhdGUgYW5zd2VyLgoKU28gaW5zdGVhZCwgbG9vayBhdCBpbmZlY3Rpb24gdnMuIHVuaW5mZWN0ZWQgaW4gbW91c2UgYW5kIHRoZW4gY29tcGFyZSB0byB0aGUKZWFybGllc3QgU2Fja3MnIHRpbWVwb2ludHMgaW4gbmV1dHJvcGhpbHMuCgpgYGB7ciBuZXV0cm9waGlsX21hY3JfbW91c2V9CnN1YnNldCA8LSAiKGhvc3RjZWxsdHlwZT09J1BNTicmaG9zdD09J211c19tdXNjdWx1cycmZXhwdHRpbWU9PSd0MTJoJykgfAogICAgICAgICAgIChob3N0Y2VsbHR5cGU9PSdtYWNyb3BoYWdlJyZob3N0PT0nbXVzX211c2N1bHVzJykiCm5ldXRfbWFjcl9tdXMgPC0gc3Vic2V0X2V4cHQobW1fZXhwdCwgc3Vic2V0PXN1YnNldCkKbmV1dF9tYWNyX211cyA8LSBzdWJzZXRfZXhwdChuZXV0X21hY3JfbXVzLCBzdWJzZXQ9ImluZmVjdHN0YXRlIT0nc3RpbSciKQpuZXV0X21hY3JfbXVzIDwtIHNldF9leHB0X2NvbmRpdGlvbnMobmV1dF9tYWNyX211cywgZmFjdD0iaW5mZWN0c3RhdGUiKQpuZXV0X21hY3JfbXVzIDwtIHNldF9leHB0X2JhdGNoZXMobmV1dF9tYWNyX211cywgZmFjdD0iaG9zdGNlbGx0eXBlIikKbmV1dF9tYWNyX211c19ub3JtIDwtIG5vcm1hbGl6ZV9leHB0KG5ldXRfbWFjcl9tdXMsIGNvbnZlcnQ9ImNwbSIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBub3JtPSJxdWFudCIsIGZpbHRlcj1UUlVFKQpwbG90X3BjYShuZXV0X21hY3JfbXVzX25vcm0pJHBsb3QKCm5ldXRfbWFjcl9tdXNfbm9ybWJhdGNoIDwtIG5vcm1hbGl6ZV9leHB0KG5ldXRfbWFjcl9tdXMsIGNvbnZlcnQ9ImNwbSIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIG5vcm09InF1YW50IiwgZmlsdGVyPVRSVUUsIGJhdGNoPSJzdmFzZXEiKQpwbG90X3BjYShuZXV0X21hY3JfbXVzX25vcm1iYXRjaCkkcGxvdAoKbmV1dF9tYWNyX211c19maWx0IDwtIG5vcm1hbGl6ZV9leHB0KG5ldXRfbWFjcl9tdXMsIGZpbHRlcj0ic2ltcGxlIikKbmV1dF9tYWNyX211c19kZSA8LSBhbGxfcGFpcndpc2UobmV1dF9tYWNyX211c19maWx0LCBwYXJhbGxlbD1GQUxTRSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZm9yY2U9VFJVRSwgbW9kZWxfYmF0Y2g9InN2YXNlcSIpCm5ldXRfbWFjcl9tdXNfdGFibGUgPC0gY29tYmluZV9kZV90YWJsZXMobmV1dF9tYWNyX211c19kZSwga2VlcGVycz1rZWVwZXJzLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGV4Y2VsPSJleGNlbC9NbU0wTG00aF92c19NbVBNTkxtMTJoX2RlX3RhYmxlcy54bHN4IikKbmV1dF9tYWNyX211c19zaWcgPC0gZXh0cmFjdF9zaWduaWZpY2FudF9nZW5lcyhuZXV0X21hY3JfbXVzX3RhYmxlLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGV4Y2VsPSJleGNlbC9NbU0wTG00aF92c19NbVBNTkxtMTJoX3NpZ190YWJsZXMueGxzeCIpCmBgYAoKSSB0aGluayB0aGlzIGhhbmRsZXMgcXVlc3Rpb25zIGEgdGhyb3VnaCBlPwoKYGBge3Igc2F2ZW1lfQpwYW5kZXI6OnBhbmRlcihzZXNzaW9uSW5mbygpKQptZXNzYWdlKCJUaGlzIGlzIGhwZ2x0b29scyBjb21taXQ6ICIsIGdldF9naXRfY29tbWl0KCkpCiMjIG1lc3NhZ2UocGFzdGUwKCJTYXZpbmcgdG8gIiwgc2F2ZWZpbGUpKQojIyB0bXAgPC0gc20oc2F2ZW1lKGZpbGVuYW1lPXNhdmVmaWxlKSkKYGBgCg==

20190417 An attempt to answer one of the big questions from Najib.

atb abelew@gmail.com

2019-04-19