TMRC3 202411: Differential Expression analyses, Tumaco only.

1 Changelog

• 202406: Added an explicit comparison of different model constructions using our most variable cell type, the neutrophils.
• 202406: Working entirely out of the container now, separated GSE/GSEA analyses, added a full treatment with clusterProfiler; I am not currently writing the cp results out as xlsx files until/unless someone expresses interest in them.
• 202309: Disabled GSVA analyses until/unless we get permission to include the mSigDB 7.5.1 release. I will simplify the filenames so that one may easily drop in a downloaded copy of the data and run those blocks. Until then, I guess you (fictitious reader) will have to trust me when I say those blocks all work? (Also, GSVA was moved to a separate document)
• 202309: Moved all gene set enrichment analyses to 04lrt_gsea_gsva.Rmd
• 202309 next day: Moving gene set enrichment back because it adds too much complexity to save/reload the DE results for gProfiler and friends.
• Still hunting for messed up colors, changed input data to match new version.

2 Introduction

The various differential expression analyses of the data generated in tmrc3_datasets will occur in this document. Most of the actual work is via the function ‘all_pairwise()’; the word ‘all’ in the name does a lot of work; it is responsible for performing all possible pairwise contrasts using all possible methods for which I have sufficient understanding to be able to write a reasonably robust pairwise function. Currently this is limited to:

• DESeq2 (Love, Huber, and Anders (2014)): Our ‘default’
• edgeR (McCarthy, Chen, and Smyth (2012)): shares a close conceptual lineage with DESeq2 I think.
• limma (Ritchie et al. (2015)): along with voom this provides a nicely robust set of tools.
• EBseq (Leng et al. (2013)): I think it is not as robust as the previous entries, but I like using it because it is an almost purely bayesian method and as such provides a different perspective on any dataset.
• Noiseq (Tarazona et al. (2015)): I noticed this method relatively recently and was sufficiently intrigued that I threw a method together using it. The authors appear to me to be looking to understand a lot of the questions on which I spend a lot of time.
• Dream (Hoffman and Roussos (2020)): I mostly like this because it uses variancePartition, which I think is a really nice toy when trying to understand what is going on in a dataset.
• basic is my own, explicitly uninformed analysis. It is my ‘negative control’ method because, if something agrees entirely with it, then I know that all the fancy math and statistics performed by that method worked out just the same as some doofus (me) just log2 subtracting the expression values. It is not quite that basic, but pretty close.

The first 3 methods allow one to add surrogate variable estimates to the model when performing the differential expression analyses. Noiseq and Dream handle surrogates using their own heuristics, EBSeq is inimicable to that kind of model, and I explicitly chose to not make that possible for basic. With that in mind, in most instances I usually deal with surrogates/batches in one of a few ways:

1. If the data is absurdly pretty, do nothing (pretty much only for well-controlled bacterial data).
2. Add a known batch factor to the model (the default).
3. Try to ensure the data is suitable and invoke sva
   1. to acquire estimates and add them to the model.
4. If the data has a known batch factor and it is particularly pathological, use the combat implementation in sva. As a general rule I do not like this option because it is data destructive.

The last two options are handled via a function named ‘all_adjusters’ in hpgltools which is responsible for ensuring that the data is sane for the assumptions made by each method and invokes each method (hopefully) properly. It returns both modified counts and model estimates when possible and has implementations for a fair number of methods in this realm. sva is my favorite by a pretty big margin, though I do sometimes use RUV (Risso et al. (2014)) and of course, in writing this document I stumbled into another interesting contender: (Molania et al. (2023)) all_adjusters() also has implementations of every example/method I got out of the papers for sva (e.g. ssva/fsva), isva, smartsva, and some others.

2.1 Define contrasts for DE analyses

Each of the following lists describes the set of contrasts that I think are interesting for the various ways one might consider the TMRC3 dataset. The variables are named according to the assumed data with which they will be used, thus tc_cf_contrasts is expected to be used for the Tumaco+Cali data and provide a series of cure/fail comparisons which (to the extent possible) across both locations. In every case, the name of the list element will be used as the contrast name, and will thus be seen as the sheet name in the output xlsx file(s); the two pieces of the character vector value are the numerator and denominator of the associated contrast.

• Our primary question: fail/cure: Any excel file written using this contrast will get a single worksheet comparing fail/cure.
• Compare fail/cure for each visit: This takes a more granular view of the previous contrast. If one is so-inclined, one could compare results from the following contrast against the previous and following contrast to learn about the dynamics of the healing (or not) process.
• All samples by visit: This is effectively the opposite of the previous and compares all samples of visit x against visit y.
• Visit 1 vs everything else: When I first did the previous set of contrasts I quickly realized that visits 2 and 3 are relatively similar and that it may be possible to gain a little power and learn a little more by combining them.
• Directly compare celltypes: We have three clinical cell types in the data and the differences among them are quite interesting.
• Ethnicities: We also have three ethnic groups in the data, though there are some wacky confounded variables when considering them through the lense of cure/fail; so any results comparing them should be treated with caution.
• Powerless visits+celltype+cf: This is a last-minute addition requested by Maria Adelaida. I assume it was suggested by a reviewer, though I do not recall seeing anything in the reviews which made this request. The number of samples we have in the data just barely supports these contrasts, and given the strength of all the various surrogates, I would be somewhat reluctant to trust any genes deemed DE in them without some other evidence. It should be noted that this is the intellectual counterpoint to the critique from a different reviewer, that artifically merging factors like this is problematic (I personally tend to agree with the later argument more than the former with the caveat that the added complexity (with respect to what is actually typed by the person (me)) can be a problem. Thus I tend to do the thing which is explicitly less statistically correct (but I can also show pretty definitively that the results are very nearly identical) in order to make it easier to show that no mistakes were made. E.g. tension between ‘correctness’ and ‘robustness’.

t_cf_contrast <- list(
  "outcome" = c("tumaco_failure", "tumaco_cure"))
cf_contrast <- list(
  "outcome" = c("failure", "cure"))
visitcf_contrasts <- list(
  "v1cf" = c("v1_failure", "v1_cure"),
  "v2cf" = c("v2_failure", "v2_cure"),
  "v3cf" = c("v3_failure", "v3_cure"))
visit_contrasts <- list(
  "v2v1" = c("c2", "c1"),
  "v3v1" = c("c3", "c1"),
  "v3v2" = c("c3", "c2"))
visit_v1later <- list(
  "later_vs_first" = c("later", "first"))
celltypes <- list(
  "eo_mono" = c("eosinophils", "monocytes"),
  "ne_mono" = c("neutrophils", "monocytes"),
  "eo_ne" = c("eosinophils", "neutrophils"))
ethnicity_contrasts <- list(
  "mestizo_indigenous" = c("mestiza", "indigena"),
  "mestizo_afrocol" = c("mestiza", "afrocol"),
  "indigenous_afrocol" = c("indigena", "afrocol"))
visittype_contrasts_mono <- list(
  "v2v1_mono_cure" = c("monocytes_2_cure", "monocytes_1_cure"),
  "v2v1_mono_failure" = c("monocytes_2_failure", "monocytes_1_failure"),
  "v3v1_mono_cure" = c("monocytes_3_cure", "monocytes_1_cure"),
  "v3v1_mono_failure" = c("monocytes_3_failure", "monocytes_1_failure"))
visittype_contrasts_eo <- list(
  "v2v1_eo_cure" = c("eosinophils_2_cure", "eosinophils_1_cure"),
  "v2v1_eo_failure" = c("eosinophils_2_failure", "eosinophils_1_failure"),
  "v3v1_eo_cure" = c("eosinophils_3_cure", "eosinophils_1_cure"),
  "v3v1_eo_failure" = c("eosinophils_3_failure", "eosinophils_1_failure"))
visittype_contrasts_ne <- list(
  "v2v1_ne_cure" = c("neutrophils_2_cure", "neutrophils_1_cure"),
  "v2v1_ne_failure" = c("neutrophils_2_failure", "neutrophils_1_failure"),
  "v3v1_ne_cure" = c("neutrophils_3_cure", "neutrophils_1_cure"),
  "v3v1_ne_failure" = c("neutrophils_3_failure", "neutrophils_1_failure"))
visittype_contrasts <- c(visittype_contrasts_mono,
                         visittype_contrasts_eo,
                         visittype_contrasts_ne)

2.2 Gene Set Enrichment

Previously, the gene set enrichment followed each DE analysis during this document. I am adding a series of explicitly GSEA analyses in this most recent iteration, in these I will pass the full DE table and check the distribution of logFC values against the genes in each category as opposed to the simpler over-enrichment of the high/low DE genes.

However, I moved those analyses to a separate document in the hopes of improving their organization.

3 Only Tumaco samples

Start over, this time with only the samples from Tumaco. We currently are assuming these will prove to be the only analyses used for final interpretation. This is primarily because we have insufficient failed treatment samples from Cali. There is one disadvantage when using these samples: they had to travel further than the samples taken in Cali and there is significant variance observed between the two locations and we cannot discern its source. In the worst case scenario (one which I think unlikely), the variance is caused by degraded RNA during transit. We do know that the samples were well-stored in RNALater and frozen/etc, so I am inclined to discount that possibility. I think a more compelling difference lies in the different population demographics observed in the two locations. Actually, now that I have typed these sentences out, I think I can semi-test this hypothesis by looking at the set of DE genes between the two locations and compare that result to the Tumaco (and/or Cali) ethnicity comparison which is most representative of the ethnicity differences between them. If I get it into my head to try this, I will need to load the DE tables from the 03differential_expression_both.Rmd document; so I am most likely to try it out in the 07var_coef document, which was mostly written by Theresa and is already examining some similar questions.

3.1 All samples

Start by considering all Tumaco cell types. Note that in this case we only use SVA, primarily because I am not certain what would be an appropriate batch factor, perhaps visit?

t_cf_clinical_de_sva <- all_pairwise(t_clinical, model_batch = "svaseq",
                                     parallel = parallel, filter = TRUE,
                                     methods = methods)

## 
##    cure failure 
##      67      56

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

## Setting 17331 low elements to zero.

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

t_cf_clinical_de_sva

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

## The logFC agreement among the methods follows:

##                 falr_vs_cr
## limma_vs_deseq      0.8063
## limma_vs_edger      0.8177
## limma_vs_basic      0.8704
## limma_vs_noiseq    -0.7287
## deseq_vs_edger      0.9845
## deseq_vs_basic      0.8242
## deseq_vs_noiseq    -0.7318
## edger_vs_basic      0.8285
## edger_vs_noiseq    -0.7370
## basic_vs_noiseq    -0.7978

t_cf_clinical_table_sva <- combine_de_tables(
  t_cf_clinical_de_sva, keepers = cf_contrast,
  excel = glue("{cf_prefix}/All_Samples/t_clinical_cf_table_sva-v{ver}.xlsx"))
t_cf_clinical_table_sva

## A set of combined differential expression results.

##             table deseq_sigup deseq_sigdown edger_sigup edger_sigdown
## 1 failure_vs_cure          94           183         103           159
##   limma_sigup limma_sigdown
## 1          50            38

## `geom_line()`: Each group consists of only one observation.
## i Do you need to adjust the group aesthetic?

## Plot describing unique/shared genes in a differential expression table.

t_cf_clinical_table_sva[["plots"]][["outcome"]][["deseq_ma_plots"]]

t_cf_clinical_sig_sva <- extract_significant_genes(
  t_cf_clinical_table_sva,
  excel = glue("{cf_prefix}/All_Samples/t_clinical_cf_sig_sva-v{ver}.xlsx"))
t_cf_clinical_sig_sva

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##         limma_up limma_down edger_up edger_down deseq_up deseq_down basic_up
## outcome       50         38      103        159       94        183       16
##         basic_down
## outcome          4

dim(t_cf_clinical_sig_sva$deseq$ups[[1]])

## [1] 94 59

dim(t_cf_clinical_sig_sva$deseq$downs[[1]])

## [1] 183  59

Repeat without the biopsies.

t_cf_clinicalnb_de_sva <- all_pairwise(t_clinical_nobiop, model_batch = "svaseq",
                                       parallel = parallel, filter = TRUE,
                                       methods = methods)

## 
##    cure failure 
##      58      51

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

## Setting 9605 low elements to zero.

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

t_cf_clinicalnb_de_sva

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

## The logFC agreement among the methods follows:

##                 falr_vs_cr
## limma_vs_deseq      0.8463
## limma_vs_edger      0.8506
## limma_vs_basic      0.8571
## limma_vs_noiseq    -0.7532
## deseq_vs_edger      0.9964
## deseq_vs_basic      0.8266
## deseq_vs_noiseq    -0.7746
## edger_vs_basic      0.8367
## edger_vs_noiseq    -0.7818
## basic_vs_noiseq    -0.8524

t_cf_clinicalnb_table_sva <- combine_de_tables(
  t_cf_clinicalnb_de_sva, keepers = cf_contrast,
  excel = glue("{cf_prefix}/All_Samples/t_clinical_nobiop_cf_table_sva-v{ver}.xlsx"))
t_cf_clinicalnb_table_sva

## A set of combined differential expression results.

##             table deseq_sigup deseq_sigdown edger_sigup edger_sigdown
## 1 failure_vs_cure         140            75         142            67
##   limma_sigup limma_sigdown
## 1          54            46

## `geom_line()`: Each group consists of only one observation.
## i Do you need to adjust the group aesthetic?

## Plot describing unique/shared genes in a differential expression table.

t_cf_clinicalnb_table_sva[["plots"]][["outcome"]][["deseq_ma_plots"]]

t_cf_clinicalnb_sig_sva <- extract_significant_genes(
  t_cf_clinicalnb_table_sva,
  excel = glue("{cf_prefix}/All_Samples/t_clinical_nobiop_cf_sig_sva-v{ver}.xlsx"))
t_cf_clinicalnb_sig_sva

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##         limma_up limma_down edger_up edger_down deseq_up deseq_down basic_up
## outcome       54         46      142         67      140         75       29
##         basic_down
## outcome          7

dim(t_cf_clinicalnb_sig_sva$deseq$ups[[1]])

## [1] 140  59

dim(t_cf_clinicalnb_sig_sva$deseq$downs[[1]])

## [1] 75 59

As the data structure’s name suggests, the above comparison seeks to learn if there are fail/cure differences discernable across all clinical celltypes in samples taken in Tumaco.

The set of steps taken in this previous block will be essentially repeated for every set of contrasts and way of mixing/matching the data and follows the path:

1. Run all_pairwise to run deseq and friends using surogate estimates provided by sva when appropriate/possible. This creates an unwieldy datastructure containing the results from all methods and all contrasts as a series of nested lists.
2. Mash them together with combine_de_tables, use the ‘keepers’ argument to define the desired numerators/denominators, and write the tables to the file provided in the ‘excel’ argument.
3. Yank out the ‘significant’ genes and send them to a separate excel document. In all cases, ‘significant’ is the set with a |log2FC| >= 1.0 and adjusted p-value <= 0.05. This reminds me, one of the reviewers mentioned a set of international guidelines for significant genes, I thought I basically know what I am doing, but this caught me completely unaware. If anyone ever reads this (no one will, let us be honest) I would love to know. The closest thing I found is: (Chung et al. (2021)), but I do not think it really addresses this idea (I have not yet read it carefully).

These datastructures are all exposed to various functions in hpgltools which allow one to poke/compare them; I am not a fan of Excel, but I think the xlsx documents creates are pretty decent, too.

4 Visit comparisons

Later in this document I do a bunch of visit/cf comparisons. In this block I want to explicitly only compare v1 to other visits. This is something I did quite a lot in the 2019 datasets, but never actually moved to this document.

tv1_vs_later <- all_pairwise(t_v1vs, model_batch = "svaseq",
                             parallel = parallel, filter = TRUE,
                             methods = methods)

## 
## first later 
##    40    69

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

## Setting 9640 low elements to zero.

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

tv1_vs_later

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

## The logFC agreement among the methods follows:

##                 ltr_vs_frs
## limma_vs_deseq      0.8394
## limma_vs_edger      0.8457
## limma_vs_basic      0.8133
## limma_vs_noiseq    -0.7094
## deseq_vs_edger      0.9983
## deseq_vs_basic      0.7946
## deseq_vs_noiseq    -0.7474
## edger_vs_basic      0.7983
## edger_vs_noiseq    -0.7524
## basic_vs_noiseq    -0.8068

tv1_vs_later_table <- combine_de_tables(
  tv1_vs_later, keepers = visit_v1later,
  excel = glue("{xlsx_prefix}/DE_Visits/tv1_vs_later_tables-v{ver}.xlsx"))
tv1_vs_later_table

## A set of combined differential expression results.

##            table deseq_sigup deseq_sigdown edger_sigup edger_sigdown
## 1 later_vs_first          24             7          22             7
##   limma_sigup limma_sigdown
## 1          23             7

## `geom_line()`: Each group consists of only one observation.
## i Do you need to adjust the group aesthetic?

## Plot describing unique/shared genes in a differential expression table.

tv1_vs_later_sig <- extract_significant_genes(
  tv1_vs_later_table,
  excel = glue("{xlsx_prefix}/DE_Visits/tv1_vs_later_sig-v{ver}.xlsx"))
tv1_vs_later_sig

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##                limma_up limma_down edger_up edger_down deseq_up deseq_down
## later_vs_first       23          7       22          7       24          7
##                basic_up basic_down
## later_vs_first        0          0

5 Sex comparison

t_sex <- subset_expt(tc_sex, subset = "clinic == 'tumaco'")

## subset_expt(): There were 184, now there are 123 samples.

t_sex

## A modified expressionSet containing 19952  and 123 sample. There are 164 metadata columns and 15 annotation columns.
## The primary condition is comprised of:
## female, male.
## Its current state is: raw(data).

t_sex_de <- all_pairwise(t_sex, model_batch = "svaseq", methods = methods,
                         parallel = parallel, filter = TRUE)

## 
## female   male 
##     22    101

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

## Setting 17311 low elements to zero.

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

t_sex_de

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

## The logFC agreement among the methods follows:

##                 mal_vs_fml
## limma_vs_deseq      0.8596
## limma_vs_edger      0.8663
## limma_vs_basic      0.9481
## limma_vs_noiseq    -0.7351
## deseq_vs_edger      0.9909
## deseq_vs_basic      0.8703
## deseq_vs_noiseq    -0.7263
## edger_vs_basic      0.8748
## edger_vs_noiseq    -0.7283
## basic_vs_noiseq    -0.7498

t_sex_table <- combine_de_tables(
  t_sex_de, excel = glue("{xlsx_prefix}/Gene_Set_Enrichment/t_sex_table-v{ver}.xlsx"))
t_sex_table

## A set of combined differential expression results.

##            table deseq_sigup deseq_sigdown edger_sigup edger_sigdown
## 1 male_vs_female         129            96         116            95
##   limma_sigup limma_sigdown
## 1          54            74

## `geom_line()`: Each group consists of only one observation.
## i Do you need to adjust the group aesthetic?

## Plot describing unique/shared genes in a differential expression table.

t_sex_sig <- extract_significant_genes(
  t_sex_table, excel = glue("{xlsx_prefix}/Gene_Set_Enrichment/t_sex_sig-v{ver}.xlsx"))
t_sex_sig

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##                limma_up limma_down edger_up edger_down deseq_up deseq_down
## male_vs_female       54         74      116         95      129         96
##                basic_up basic_down
## male_vs_female       15         10

tc_sex_cure <- subset_expt(tc_sex, subset = "finaloutcome=='cure'")

## subset_expt(): There were 184, now there are 122 samples.

t_sex_cure <- subset_expt(tc_sex_cure, subset = "clinic == 'tumaco'")

## subset_expt(): There were 122, now there are 67 samples.

t_sex_cure_de <- all_pairwise(t_sex_cure, model_batch = "svaseq",
                              parallel = parallel, filter = TRUE,
                              methods = methods)

## 
## female   male 
##     13     54

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

## Setting 8998 low elements to zero.

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

t_sex_cure_de

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

## The logFC agreement among the methods follows:

##                 mal_vs_fml
## limma_vs_deseq      0.7804
## limma_vs_edger      0.8380
## limma_vs_basic      0.9284
## limma_vs_noiseq    -0.7515
## deseq_vs_edger      0.9294
## deseq_vs_basic      0.7995
## deseq_vs_noiseq    -0.6524
## edger_vs_basic      0.8474
## edger_vs_noiseq    -0.6952
## basic_vs_noiseq    -0.7882

t_sex_cure_table <- combine_de_tables(
  t_sex_cure_de, excel = glue("{xlsx_prefix}/DE_Sex/t_sex_cure_table-v{ver}.xlsx"))
t_sex_cure_table

## A set of combined differential expression results.

##            table deseq_sigup deseq_sigdown edger_sigup edger_sigdown
## 1 male_vs_female         176           134         162           143
##   limma_sigup limma_sigdown
## 1          64           108

## `geom_line()`: Each group consists of only one observation.
## i Do you need to adjust the group aesthetic?

## Plot describing unique/shared genes in a differential expression table.

t_sex_cure_sig <- extract_significant_genes(
  t_sex_cure_table, excel = glue("{xlsx_prefix}/DE_Sex/t_sex_cure_sig-v{ver}.xlsx"))
t_sex_cure_sig

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##                limma_up limma_down edger_up edger_down deseq_up deseq_down
## male_vs_female       64        108      162        143      176        134
##                basic_up basic_down
## male_vs_female       13          5

6 Ethnicity comparisons

t_ethnicity_de <- all_pairwise(t_etnia_expt, model_batch = "svaseq",
                               parallel = parallel, filter = TRUE,
                               methods = methods)

## 
##  afrocol indigena  mestiza 
##       76       19       28

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

## Setting 15870 low elements to zero.

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

t_ethnicity_de

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

t_ethnicity_table <- combine_de_tables(
  t_ethnicity_de, keepers = ethnicity_contrasts,
  excel = glue("{xlsx_prefix}/DE_Ethnicity/t_ethnicity_table-v{ver}.xlsx"))
t_ethnicity_table

## A set of combined differential expression results.

##                 table deseq_sigup deseq_sigdown edger_sigup edger_sigdown
## 1 mestiza_vs_indigena          83            97          67           108
## 2  mestiza_vs_afrocol          57            92          52            96
## 3 indigena_vs_afrocol         165           236         187           216
##   limma_sigup limma_sigdown
## 1          58            56
## 2          42            53
## 3         165           147

## Plot describing unique/shared genes in a differential expression table.

t_ethnicity_sig <- extract_significant_genes(
  t_ethnicity_table, according_to = "deseq",
  excel = glue("{xlsx_prefix}/DE_Ethnicity/t_ethnicity_sig-v{ver}.xlsx"))
t_ethnicity_sig

## A set of genes deemed significant according to deseq.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##                    deseq_up deseq_down
## mestizo_indigenous       83         97
## mestizo_afrocol          57         92
## indigenous_afrocol      165        236

6.0.1 Separate the Tumaco data by visit

One of the most compelling ideas in the data is the opportunity to find genes in the first visit which may help predict the likelihood that a person will respond well to treatment. The following block will therefore look at cure/fail from Tumaco at visit 1.

6.0.1.1 Cure/Fail, Tumaco Visit 1

t_cf_clinical_v1_de_sva <- all_pairwise(tv1_samples, model_batch = "svaseq",
                                        parallel = parallel, filter = TRUE,
                                        methods = methods)

## 
##    cure failure 
##      30      24

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

## Setting 7655 low elements to zero.

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

t_cf_clinical_v1_de_sva

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

## The logFC agreement among the methods follows:

##                 falr_vs_cr
## limma_vs_deseq      0.7398
## limma_vs_edger      0.7829
## limma_vs_basic      0.6886
## limma_vs_noiseq    -0.5605
## deseq_vs_edger      0.9537
## deseq_vs_basic      0.6955
## deseq_vs_noiseq    -0.6068
## edger_vs_basic      0.7228
## edger_vs_noiseq    -0.6303
## basic_vs_noiseq    -0.7772

t_cf_clinical_v1_table_sva <- combine_de_tables(
  t_cf_clinical_v1_de_sva, keepers = cf_contrast,
  excel = glue("{cf_prefix}/Visits/t_clinical_v1_cf_table_sva-v{ver}.xlsx"))
t_cf_clinical_v1_table_sva

## A set of combined differential expression results.

##             table deseq_sigup deseq_sigdown edger_sigup edger_sigdown
## 1 failure_vs_cure          27            75          28            55
##   limma_sigup limma_sigdown
## 1           3             3

## `geom_line()`: Each group consists of only one observation.
## i Do you need to adjust the group aesthetic?

## Plot describing unique/shared genes in a differential expression table.

t_cf_clinical_v1_sig_sva <- extract_significant_genes(
  t_cf_clinical_v1_table_sva,
  excel = glue("{cf_prefix}/Visits/t_clinical_v1_cf_sig_sva-v{ver}.xlsx"))
t_cf_clinical_v1_sig_sva

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##         limma_up limma_down edger_up edger_down deseq_up deseq_down basic_up
## outcome        3          3       28         55       27         75        0
##         basic_down
## outcome          0

dim(t_cf_clinical_v1_sig_sva$deseq$ups[[1]])

## [1] 27 59

dim(t_cf_clinical_v1_sig_sva$deseq$downs[[1]])

## [1] 75 59

6.0.1.2 Cure/Fail, Tumaco Visit 2

The visit 2 and visit 3 samples are interesting because they provide an opportunity to see if we can observe changes in response in the middle and end of treatment…

t_cf_clinical_v2_de_sva <- all_pairwise(tv2_samples, model_batch = "svaseq",
                                        parallel = parallel, filter = TRUE,
                                        methods = methods)

## 
##    cure failure 
##      20      15

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

## Setting 2857 low elements to zero.

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

t_cf_clinical_v2_de_sva

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

## The logFC agreement among the methods follows:

##                 falr_vs_cr
## limma_vs_deseq      0.8138
## limma_vs_edger      0.8162
## limma_vs_basic      0.7404
## limma_vs_noiseq    -0.6360
## deseq_vs_edger      0.9986
## deseq_vs_basic      0.7689
## deseq_vs_noiseq    -0.7586
## edger_vs_basic      0.7701
## edger_vs_noiseq    -0.7602
## basic_vs_noiseq    -0.8078

t_cf_clinical_v2_table_sva <- combine_de_tables(
  t_cf_clinical_v2_de_sva, keepers = cf_contrast,
  excel = glue("{cf_prefix}/Visits/t_clinical_v2_cf_table_sva-v{ver}.xlsx"))
t_cf_clinical_v2_table_sva

## A set of combined differential expression results.

##             table deseq_sigup deseq_sigdown edger_sigup edger_sigdown
## 1 failure_vs_cure          51            15          50            11
##   limma_sigup limma_sigdown
## 1           0             0

## `geom_line()`: Each group consists of only one observation.
## i Do you need to adjust the group aesthetic?

## Plot describing unique/shared genes in a differential expression table.

t_cf_clinical_v2_sig_sva <- extract_significant_genes(
  t_cf_clinical_v2_table_sva,
  excel = glue("{cf_prefix}/Visits/t_clinical_v2_cf_sig_sva-v{ver}.xlsx"))
t_cf_clinical_v2_sig_sva

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##         limma_up limma_down edger_up edger_down deseq_up deseq_down basic_up
## outcome        0          0       50         11       51         15        0
##         basic_down
## outcome          0

dim(t_cf_clinical_v2_sig_sva$deseq$ups[[1]])

## [1] 51 59

dim(t_cf_clinical_v2_sig_sva$deseq$downs[[1]])

## [1] 15 59

6.0.2 By cell type

Now let us switch our view to each individual cell type collected. The hope here is that we will be able to learn some cell-specific differences in the response for people who did(not) respond well.

6.0.2.1 Cure/Fail, Biopsies

t_cf_biopsy_de_sva <- all_pairwise(t_biopsies, model_batch = "svaseq",
                                   parallel = parallel, filter = TRUE,
                                   methods = methods)

## 
##    tumaco_cure tumaco_failure 
##              9              5

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

## Setting 146 low elements to zero.

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

t_cf_biopsy_de_sva

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

## The logFC agreement among the methods follows:

##                 tmc_flr___
## limma_vs_deseq      0.7927
## limma_vs_edger      0.8628
## limma_vs_basic      0.8497
## limma_vs_noiseq    -0.7628
## deseq_vs_edger      0.9516
## deseq_vs_basic      0.8164
## deseq_vs_noiseq    -0.7927
## edger_vs_basic      0.8809
## edger_vs_noiseq    -0.8535
## basic_vs_noiseq    -0.8819

t_cf_biopsy_table_sva <- combine_de_tables(
  t_cf_biopsy_de_sva, keepers = t_cf_contrast,
  excel = glue("{cf_prefix}/Biopsies/t_biopsy_cf_table_sva-v{ver}.xlsx"))
t_cf_biopsy_table_sva

## A set of combined differential expression results.

##                           table deseq_sigup deseq_sigdown edger_sigup
## 1 tumaco_failure_vs_tumaco_cure          17            11          19
##   edger_sigdown limma_sigup limma_sigdown
## 1            15           0             0

## `geom_line()`: Each group consists of only one observation.
## i Do you need to adjust the group aesthetic?

## Plot describing unique/shared genes in a differential expression table.

t_cf_biopsy_sig_sva <- extract_significant_genes(
  t_cf_biopsy_table_sva,
  excel = glue("{cf_prefix}/Biopsies/t_cf_biopsy_sig_sva-v{ver}.xlsx"))
t_cf_biopsy_sig_sva

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##         limma_up limma_down edger_up edger_down deseq_up deseq_down basic_up
## outcome        0          0       19         15       17         11        0
##         basic_down
## outcome          0

dim(t_cf_biopsy_sig_sva$deseq$ups[[1]])

## [1] 17 59

dim(t_cf_biopsy_sig_sva$deseq$downs[[1]])

## [1] 11 59

6.0.2.2 Cure/Fail, Monocytes

Same question, but this time looking at monocytes. In addition, this comparison was done twice, once using SVA and once using visit as a batch factor.

t_cf_monocyte_de_sva <- all_pairwise(t_monocytes, model_batch = "svaseq",
                                     parallel = parallel, filter = TRUE,
                                     methods = methods)

## 
##    tumaco_cure tumaco_failure 
##             21             21

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

## Setting 736 low elements to zero.

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

t_cf_monocyte_de_sva

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

## The logFC agreement among the methods follows:

##                 tmc_flr___
## limma_vs_deseq      0.8614
## limma_vs_edger      0.8663
## limma_vs_basic      0.9210
## limma_vs_noiseq    -0.8064
## deseq_vs_edger      0.9989
## deseq_vs_basic      0.8506
## deseq_vs_noiseq    -0.7947
## edger_vs_basic      0.8560
## edger_vs_noiseq    -0.8011
## basic_vs_noiseq    -0.8949

t_cf_monocyte_table_sva <- combine_de_tables(
  t_cf_monocyte_de_sva, keepers = t_cf_contrast,
  excel = glue("{cf_prefix}/Monocytes/t_monocyte_cf_table_sva-v{ver}.xlsx"))
t_cf_monocyte_table_sva

## A set of combined differential expression results.

##                           table deseq_sigup deseq_sigdown edger_sigup
## 1 tumaco_failure_vs_tumaco_cure          60            52          56
##   edger_sigdown limma_sigup limma_sigdown
## 1            51          11            34

## `geom_line()`: Each group consists of only one observation.
## i Do you need to adjust the group aesthetic?

## Plot describing unique/shared genes in a differential expression table.

head(t_cf_monocyte_table_sva[["data"]][["outcome"]][["deseq_logfc"]])

## [1]  0.33760 -0.07193  0.09665 -0.09082 -0.13500  0.23270

t_cf_monocyte_sig_sva <- extract_significant_genes(
  t_cf_monocyte_table_sva,
  excel = glue("{cf_prefix}/Monocytes/t_monocyte_cf_sig_sva-v{ver}.xlsx"))
t_cf_monocyte_sig_sva

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##         limma_up limma_down edger_up edger_down deseq_up deseq_down basic_up
## outcome       11         34       56         51       60         52        8
##         basic_down
## outcome         21

dim(t_cf_monocyte_sig_sva$deseq$ups[[1]])

## [1] 60 59

dim(t_cf_monocyte_sig_sva$deseq$downs[[1]])

## [1] 52 59

t_cf_monocyte_de_batchvisit <- all_pairwise(t_monocytes, model_batch = TRUE,
                                            parallel = parallel, filter = TRUE,
                                            methods = methods)

## 
##    tumaco_cure tumaco_failure 
##             21             21 
## 
##  3  2  1 
## 13 13 16

t_cf_monocyte_de_batchvisit

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: batch in model/limma.

## The primary analysis performed 10 comparisons.

## The logFC agreement among the methods follows:

##                 tmc_flr___
## limma_vs_deseq      0.8120
## limma_vs_edger      0.8150
## limma_vs_basic      0.9509
## limma_vs_noiseq    -0.8326
## deseq_vs_edger      0.9998
## deseq_vs_basic      0.8505
## deseq_vs_noiseq    -0.7807
## edger_vs_basic      0.8540
## edger_vs_noiseq    -0.7852
## basic_vs_noiseq    -0.8949

t_cf_monocyte_table_batchvisit <- combine_de_tables(
  t_cf_monocyte_de_batchvisit, keepers = t_cf_contrast,
  excel = glue("{cf_prefix}/Monocytes/t_monocyte_cf_table_batchvisit-v{ver}.xlsx"))
t_cf_monocyte_table_batchvisit

## A set of combined differential expression results.

##                           table deseq_sigup deseq_sigdown edger_sigup
## 1 tumaco_failure_vs_tumaco_cure          43            93          47
##   edger_sigdown limma_sigup limma_sigdown
## 1           105           6            13

## `geom_line()`: Each group consists of only one observation.
## i Do you need to adjust the group aesthetic?

## Plot describing unique/shared genes in a differential expression table.

t_cf_monocyte_sig_batchvisit <- extract_significant_genes(
  t_cf_monocyte_table_batchvisit,
  excel = glue("{cf_prefix}/Monocytes/t_monocyte_cf_sig_batchvisit-v{ver}.xlsx"))
t_cf_monocyte_sig_batchvisit

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##         limma_up limma_down edger_up edger_down deseq_up deseq_down basic_up
## outcome        6         13       47        105       43         93        8
##         basic_down
## outcome         21

dim(t_cf_monocyte_sig_batchvisit$deseq$ups[[1]])

## [1] 43 59

dim(t_cf_monocyte_sig_batchvisit$deseq$downs[[1]])

## [1] 93 59

6.0.2.2.1 Repeat test of different models: Monocytes

Either above or below this section I have a nearly identical block which seeks to demonstrate the similarities/difference observed between my preferred/simplified model vs. a more explicitly correct and complex model. If the trend holds from what we observed with the eosinophils and neutrophils, I would expect to see that the results are marginally ‘better’ (as defined by the strength of the perceived interleukin response and raw number of ‘significant’ genes); but I remain worried that this will prove a more brittle and error-prone analysis.

## The original pairwise invocation with sva:
##t_cf_monocyte_de_sva <- all_pairwise(t_monocyte, model_batch = "svaseq",
##                                     filter = TRUE, parallel = FALSE,
##                                     methods = methods)
test_monocytes <- normalize_expt(t_monocytes, filter = "simple")

## Removing 2619 low-count genes (17333 remaining).

test_mono_design <- pData(test_monocytes)
test_formula <- as.formula("~ finaloutcome + visitnumber")
test_model <- model.matrix(test_formula, data = test_mono_design)
null_formula <- as.formula("~ visitnumber")
null_model <- model.matrix(null_formula, data = test_mono_design)

linear_mtrx <- exprs(test_monocytes)
l2_mtrx <- log2(linear_mtrx + 1)
chosen_surrogates <- sva::num.sv(dat = l2_mtrx, mod = test_model)
chosen_surrogates

## [1] 3

surrogate_result <- sva::svaseq(
  dat = linear_mtrx, n.sv = chosen_surrogates, mod = test_model, mod0 = null_model)

## Number of significant surrogate variables is:  3 
## Iteration (out of 5 ):1  2  3  4  5

model_adjust <- as.matrix(surrogate_result[["sv"]])

colnames(model_adjust) <- c("SV1", "SV2", "SV3")
rownames(model_adjust) <- rownames(pData(test_monocytes))
longer_model <- as.formula("~ finaloutcome + visitnumber + SV1 + SV2 + SV3")
mono_design_svs <- cbind(test_mono_design, model_adjust)
summarized <- DESeq2::DESeqDataSetFromMatrix(countData = linear_mtrx,
                                             colData = mono_design_svs,
                                             design = longer_model)

## converting counts to integer mode

deseq_run <- DESeq2::DESeq(summarized)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

deseq_table <- as.data.frame(DESeq2::results(object = deseq_run,
                                             contrast = c("finaloutcome", "failure", "cure"),
                                             format = "DataFrame"))

big_table <- t_cf_monocyte_table_sva[["data"]][["outcome"]]
only_deseq <- big_table[, c("deseq_logfc", "deseq_adjp")]
merged <- merge(deseq_table, only_deseq, by = "row.names")
rownames(merged) <- merged[["Row.names"]]
merged[["Row.names"]] <- NULL

cor_value <- cor.test(merged[["log2FoldChange"]], merged[["deseq_logfc"]])
cor_value

## 
##  Pearson's product-moment correlation
## 
## data:  merged[["log2FoldChange"]] and merged[["deseq_logfc"]]
## t = 300, df = 10860, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9426 0.9466
## sample estimates:
##    cor 
## 0.9446

logfc_plotter <- plot_linear_scatter(merged[, c("log2FoldChange", "deseq_logfc")])
logfc_plot <- logfc_plotter[["scatter"]] +
  xlab("DESeq2 log2FC: Visit explicitly in model") +
  ylab("DESeq2 log2FC: Default pairwise comparison") +
  ggtitle(glue("Comparing results from models: {prettyNum(cor_value[['estimate']])} (pearson)"))
pp(file = "figures/compare_cf_and_visit_in_model_monocyte_logfc.svg")
logfc_plot
dev.off()

## png 
##   2

logfc_plot

cor_value <- cor.test(merged[["padj"]], merged[["deseq_adjp"]], method = "spearman")

## Warning in cor.test.default(merged[["padj"]], merged[["deseq_adjp"]], method =
## "spearman"): Cannot compute exact p-value with ties

cor_value

## 
##  Spearman's rank correlation rho
## 
## data:  merged[["padj"]] and merged[["deseq_adjp"]]
## S = 2.9e+10, p-value <2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##    rho 
## 0.8642

adjp_plotter <- plot_linear_scatter(merged[, c("padj", "deseq_adjp")])
adjp_plot <- adjp_plotter[["scatter"]] +
  xlab("DESeq2 adjp: Visit explicitly in model") +
  ylab("DESeq2 adjp: Default pairwise comparison") +
  ggtitle(glue("Comparing results from models: {prettyNum(cor_value[['estimate']])} (spearman)"))
pp(file = "images/compare_cf_and_visit_in_model_monocyte_adjp.svg")
adjp_plot
dev.off()

## png 
##   2

adjp_plot

previous_sig_idx <- big_table[["deseq_adjp"]] <= 0.05 & abs(big_table[["deseq_logfc"]] >= 1.0)
summary(previous_sig_idx)

##    Mode   FALSE    TRUE 
## logical   10802      60

previous_genes <- rownames(big_table)[previous_sig_idx]

new_sig_idx <- abs(deseq_table[["log2FoldChange"]]) >= 1.0 & deseq_table[["padj"]] < 0.05
new_genes <- rownames(deseq_table)[new_sig_idx]
na_idx <- is.na(new_genes)
new_genes <- new_genes[!na_idx]

Vennerable::Venn(list("previous" = previous_genes, "new" = new_genes))

## A Venn object on 2 sets named
## previous,new 
##  00  10  01  11 
##   0  28 156  32

test_new <- simple_gprofiler(new_genes)
test_new

## A set of ontologies produced by gprofiler using 188
## genes against the hsapiens annotations and significance cutoff 0.05.
## There are: 
## 13 MF
## 36 BP
## 3 KEGG
## 4 REAC
## 2 WP
## 3 TF
## 0 MIRNA
## 0 HPA
## 0 CORUM
## 0 HP hits.

test_old <- simple_gprofiler(previous_genes)
test_old

## A set of ontologies produced by gprofiler using 60
## genes against the hsapiens annotations and significance cutoff 0.05.
## There are: 
## 0 MF
## 44 BP
## 3 KEGG
## 0 REAC
## 2 WP
## 0 TF
## 0 MIRNA
## 0 HPA
## 0 CORUM
## 0 HP hits.

new_annotated <- merge(fData(t_monocytes), deseq_table, by = "row.names")
rownames(new_annotated) <- new_annotated[["Row.names"]]
new_annotated[["Row.names"]] <- NULL
write_xlsx(data = new_annotated, excel = "excel/monocyte_visit_in_model_sva_cf_new.xlsx")

## write_xlsx() wrote excel/monocyte_visit_in_model_sva_cf_new.xlsx.
## The cursor is on sheet first, row: 17336 column: 23.

old_annotated <- merge(fData(t_eosinophils), big_table, by = "row.names")
rownames(old_annotated) <- old_annotated[["Row.names"]]
old_annotated[["Row.names"]] <- NULL
write_xlsx(data = old_annotated, excel = "excel/monocyte_visit_in_model_sva_cf_old.xlsx")

## write_xlsx() wrote excel/monocyte_visit_in_model_sva_cf_old.xlsx.
## The cursor is on sheet first, row: 10865 column: 76.

The previous couple of blocks I think are not approaching this problem correctly. In our discussions with Neal, he suggested that a mixed linear model is the most appropriate method for this type of query. I think I can perform that with limma, via voom. Let us try and see what happens. After doing some reading, I think the most appropriate way to perform this is to use dream() from varianceParition, which is cool because I really like it.

I have already written a skeleton function ‘dream_pairwise()’ as a sibling to my other *_pairwise() functions. I think that with some minor modifications (or maybe none at all, when I wrote it I was thinking about fun models that variancePartition supports) it can accept the mixed linear model of interest.

mixed_model <- "~ 0 + finaloutcome + visitnumber + (1|donor)"
mixed_fstring <- as.formula(mixed_model)
get_formula_factors(mixed_model)

## Getting factors from: ~ 0 + finaloutcome + visitnumber + (1|donor).

## $type
## [1] "cellmeans"
## 
## $interaction
## [1] FALSE
## 
## $mixed
## [1] TRUE
## 
## $mixers
## [1] "1"     "donor"
## 
## $cellmeans_intercept
## [1] "1"
## 
## $factors
## [1] "finaloutcome" "visitnumber"  "donor"

mixed_eosinophil_de <- dream_pairwise(t_eosinophils, alt_model = mixed_fstring)

## Starting limma/varpart pairwise comparison.

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

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Error in 1 | donor: operations are possible only for numeric, logical or complex types

mixed_monocyte_de <- dream_pairwise(t_monocytes, alt_model = mixed_fstring)

## Starting limma/varpart pairwise comparison.

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

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Error in 1 | donor: operations are possible only for numeric, logical or complex types

mixed_neutrophil_de <- dream_pairwise(t_neutrophils, alt_model = mixed_fstring)

## Starting limma/varpart pairwise comparison.

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

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Error in 1 | donor: operations are possible only for numeric, logical or complex types

mixed_model_fv <- "~ 0 + finaloutcome + visitnumber"
mixed_fstring_fv <- as.formula(mixed_model_fv)
get_formula_factors(mixed_model_fv)

## Getting factors from: ~ 0 + finaloutcome + visitnumber.

## $type
## [1] "cellmeans"
## 
## $interaction
## [1] FALSE
## 
## $mixed
## [1] FALSE
## 
## $cellmeans_intercept
## [1] "0"
## 
## $factors
## [1] "finaloutcome" "visitnumber"

mixed_eosinophil_fv_de <- dream_pairwise(t_eosinophils, alt_model = mixed_fstring_fv)

## Starting limma/varpart pairwise comparison.

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

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Attempting voomWithDreamWeights.

## Error in exprObj$E: $ operator is invalid for atomic vectors

mixed_monocyte_fv_de <- dream_pairwise(t_monocytes, alt_model = mixed_fstring_fv)

## Starting limma/varpart pairwise comparison.

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

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Attempting voomWithDreamWeights.

## Error in exprObj$E: $ operator is invalid for atomic vectors

mixed_neutrophil_fv_de <- dream_pairwise(t_neutrophils, alt_model = mixed_fstring_fv)

## Starting limma/varpart pairwise comparison.

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

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Attempting voomWithDreamWeights.

## Error in exprObj$E: $ operator is invalid for atomic vectors

Now compare the results of pre/post donor

monocyte_visit_with_donor <- mixed_monocyte_de$all_tables$contrasts[[1]]

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

monocyte_visit_without_donor <- mixed_monocyte_fv_de$all_tables$contrasts[[1]]

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

donor_aucc <- calculate_aucc(monocyte_visit_with_donor, monocyte_visit_without_donor,
                             px = "adj.P.Val", py = "adj.P.Val",
                             lx = "logFC", ly = "logFC")

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

with_donor_genes <- abs(monocyte_visit_with_donor[["logFC"]]) >= 1.0 &
  monocyte_visit_with_donor[["P.Value"]] <= 0.05

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

without_donor_genes <- abs(monocyte_visit_without_donor[["logFC"]]) >= 1.0 &
  monocyte_visit_with_donor[["P.Value"]] <= 0.05

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

donor_genes <- rownames(monocyte_visit_with_donor)[with_donor_genes]

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'rownames': object 'monocyte_visit_with_donor' not found

visit_genes <- rownames(monocyte_visit_with_donor)[without_donor_genes]

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'rownames': object 'monocyte_visit_with_donor' not found

venn_lst <- list(
  "with_donor" = donor_genes,
  "with_visit" = visit_genes)

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

Vennerable::Venn(venn_lst)

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

Compare back to deseq with SVA and with SVA+visit

deseq_aucc <- calculate_aucc(merged, monocyte_visit_without_donor,
                             px = "deseq_adjp", py = "P.Value",
                             lx = "deseq_logfc", ly = "logFC")

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

deseq_aucc

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

deseq_genes_idx <- abs(merged[["deseq_logfc"]]) >= 1.0 &
  merged[["deseq_adjp"]] <= 0.05
without_donor_genes_idx <- abs(monocyte_visit_without_donor[["logFC"]]) >= 1.0 &
  monocyte_visit_with_donor[["P.Value"]] <= 0.05

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

deseq_genes <- rownames(merged)[deseq_genes_idx]
visit_genes <- rownames(monocyte_visit_with_donor)[without_donor_genes_idx]

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'rownames': object 'monocyte_visit_with_donor' not found

venn_lst <- list(
  "with_donor" = deseq_genes,
  "with_visit" = visit_genes)

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

Vennerable::Venn(venn_lst)

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

And with visit

deseq_aucc <- calculate_aucc(merged, monocyte_visit_without_donor,
                             px = "log2FoldChange", py = "padj",
                             lx = "adj.P.Val", ly = "logFC")

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

deseq_aucc

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

deseq_genes_idx <- abs(merged[["log2FoldChange"]]) >= 1.0 &
  merged[["padj"]] <= 0.05
without_donor_genes_idx <- abs(monocyte_visit_without_donor[["logFC"]]) >= 1.0 &
  monocyte_visit_with_donor[["P.Value"]] <= 0.05

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

deseq_genes <- rownames(merged)[deseq_genes_idx]
visit_genes <- rownames(monocyte_visit_with_donor)[without_donor_genes_idx]

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'rownames': object 'monocyte_visit_with_donor' not found

venn_lst <- list(
  "with_donor" = deseq_genes,
  "with_visit" = visit_genes)

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

Vennerable::Venn(venn_lst)

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

mixed_model_fd <- "~ 0 + finaloutcome + (1|donor)"
mixed_fstring_fd <- as.formula(mixed_model_fd)
get_formula_factors(mixed_model_fd)

## Getting factors from: ~ 0 + finaloutcome + (1|donor).

## $type
## [1] "cellmeans"
## 
## $interaction
## [1] FALSE
## 
## $mixed
## [1] TRUE
## 
## $mixers
## [1] "1"     "donor"
## 
## $cellmeans_intercept
## [1] "1"
## 
## $factors
## [1] "finaloutcome" "donor"

mixed_eosinophil_fd_de <- dream_pairwise(t_eosinophils, alt_model = mixed_fstring_fd)

## Starting limma/varpart pairwise comparison.

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

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Error in 1 | donor: operations are possible only for numeric, logical or complex types

mixed_monocyte_fd_de <- dream_pairwise(t_monocytes, alt_model = mixed_fstring_fd)

## Starting limma/varpart pairwise comparison.

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

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Error in 1 | donor: operations are possible only for numeric, logical or complex types

mixed_neutrophil_fd_de <- dream_pairwise(t_neutrophils, alt_model = mixed_fstring_fd)

## Starting limma/varpart pairwise comparison.

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

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Error in 1 | donor: operations are possible only for numeric, logical or complex types

6.1 Compare monocytes

dream_result <- mixed_monocyte_de[["all_tables"]][["contrasts"]][[1]]

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

big_table <- t_cf_monocyte_table_sva[["data"]][["outcome"]]
merged <- merge(big_table, dream_result, by = "row.names")

## Error in h(simpleError(msg, call)): error in evaluating the argument 'y' in selecting a method for function 'merge': object 'dream_result' not found

rownames(merged) <- merged[["Row.names"]]
merged[["Row.names"]] <- NULL
cor_value <- cor.test(merged[["logFC"]], merged[["deseq_logfc"]])

## Error in cor.test.default(merged[["logFC"]], merged[["deseq_logfc"]]): 'x' must be a numeric vector

cor_value

## 
##  Spearman's rank correlation rho
## 
## data:  merged[["padj"]] and merged[["deseq_adjp"]]
## S = 2.9e+10, p-value <2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##    rho 
## 0.8642

t_cf_monocyte_de_sva[["dream"]] <- mixed_monocyte_de

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

test <- combine_de_tables(t_cf_monocyte_de_sva, excel = "/tmp/test_combined.xlsx")

## Deleting the file /tmp/test_combined.xlsx before writing the tables.

test_aucc <- calculate_aucc(big_table, tbl2 = dream_result, px = "deseq_adjp", py = "adj.P.Val", lx = "deseq_logfc", ly = "logFC")

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

logfc_plotter <- plot_linear_scatter(merged[, c("logFC", "deseq_logfc")])

## Error in `[.data.frame`(merged, , c("logFC", "deseq_logfc")): undefined columns selected

logfc_plot <- logfc_plotter[["scatter"]] +
  xlab("Dream log2FC with (1|donor) and visit in model") +
  ylab("DESeq2 log2FC: Default pairwise comparison") +
  ggtitle(glue("Comparing results from models: {prettyNum(cor_value[['estimate']])} (pearson)"))
pp(file = "figures/compare_cf_and_visit_in_model_monocyte_logfc.svg")
logfc_plot
dev.off()

## png 
##   2

logfc_plot

previous_sig_idx <- merged[["deseq_adjp"]] <= 0.05 & abs(merged[["deseq_logfc"]] >= 1.0)
summary(previous_sig_idx)

##    Mode   FALSE    TRUE 
## logical   10802      60

previous_genes <- rownames(merged)[previous_sig_idx]

new_sig_idx <- abs(merged[["logFC"]]) >= 1.0 & merged[["P.Value"]] < 0.05

## Error in abs(merged[["logFC"]]): non-numeric argument to mathematical function

summary(new_sig_idx)

##    Mode   FALSE    TRUE    NA's 
## logical   16603     188     542

new_genes <- rownames(merged)[new_sig_idx]
na_idx <- is.na(new_genes)
new_genes <- new_genes[!na_idx]

annot <- fData(t_monocytes)
compare <- Vennerable::Venn(list("previous" = previous_genes, "new" = new_genes))
shared_genes <- compare@IntersectionSets[["11"]]
name_idx <- rownames(annot) %in% shared_genes
annot[name_idx, ]

##  [1] ensembl_gene_id       ensembl_transcript_id version              
##  [4] transcript_version    description           gene_biotype         
##  [7] cds_length            chromosome_name       strand               
## [10] start_position        end_position          hgnc_symbol          
## [13] uniprot_gn_symbol     transcript            mean_cds_len         
## <0 rows> (or 0-length row.names)

7 TODO

1. Repeat this process, clean it up for: monocytes/neutrophils
2. Repeat this process and compare the results when using models which are:
   1. ~ finaloutcome + visitnumber + (1|donor)
   2. ~ finaloutcome + visitnumber
3. Compare the results from limma for a,b
4. Extract the set of ‘significant’ genes via logFC/pvalue for all of the above and see the shared/unique genes.

8 Test some code to do this automagically

I think I implemented some changes to my toys which make them able to set up appropriate models and contrasts for arbitrarily chosen models instead of only my simplified version. Let us test that here!

tt = normalize_expt(t_monocytes, transform = "log2", convert = "cpm",
                    filter = TRUE, batch = "svaseq")

t_monocyte_cf_visit_sva_de <- all_pairwise(t_monocytes, model_batch = "svaseq",
                                           alt_model = "~ 0 + finaloutcome + visitnumber",
                                           null_model = "~ 0 + visitnumber", filter = TRUE,
                                           parallel = FALSE, methods = methods)

8.0.1 Individual visits, Monocytes

Now focus in on the monocyte samples on a per-visit basis.

8.0.1.1 Visit 1

t_cf_monocyte_v1_de_sva <- all_pairwise(tv1_monocytes, model_batch = "svaseq",
                                        parallel = parallel, filter = TRUE,
                                        methods = methods)

## 
##    tumaco_cure tumaco_failure 
##              8              8

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

## Setting 190 low elements to zero.

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

t_cf_monocyte_v1_de_sva

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

## The logFC agreement among the methods follows:

##                 tmc_flr___
## limma_vs_deseq      0.8902
## limma_vs_edger      0.8905
## limma_vs_basic      0.9440
## limma_vs_noiseq    -0.8547
## deseq_vs_edger      0.9999
## deseq_vs_basic      0.8866
## deseq_vs_noiseq    -0.8760
## edger_vs_basic      0.8870
## edger_vs_noiseq    -0.8766
## basic_vs_noiseq    -0.9062

t_cf_monocyte_v1_table_sva <- combine_de_tables(
  t_cf_monocyte_v1_de_sva, keepers = t_cf_contrast,
  excel = glue("{cf_prefix}/Monocytes/t_monocyte_v1_cf_table_sva-v{ver}.xlsx"))
t_cf_monocyte_v1_table_sva

## A set of combined differential expression results.

##                           table deseq_sigup deseq_sigdown edger_sigup
## 1 tumaco_failure_vs_tumaco_cure          14            52          15
##   edger_sigdown limma_sigup limma_sigdown
## 1            57           0             0

## `geom_line()`: Each group consists of only one observation.
## i Do you need to adjust the group aesthetic?

## Plot describing unique/shared genes in a differential expression table.

t_cf_monocyte_v1_sig_sva <- extract_significant_genes(
  t_cf_monocyte_v1_table_sva,
  excel = glue("{cf_prefix}/Monocytes/t_monocyte_v1_cf_sig_sva-v{ver}.xlsx"))
t_cf_monocyte_v1_sig_sva

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##         limma_up limma_down edger_up edger_down deseq_up deseq_down basic_up
## outcome        0          0       15         57       14         52        0
##         basic_down
## outcome          0

dim(t_cf_monocyte_v1_sig_sva$deseq$ups[[1]])

## [1] 14 59

dim(t_cf_monocyte_v1_sig_sva$deseq$downs[[1]])

## [1] 52 59

8.0.1.2 Monocytes: Compare sva to batch-in-model

sva_aucc <- calculate_aucc(t_cf_monocyte_table_sva[["data"]][[1]],
                           tbl2 = t_cf_monocyte_table_batchvisit[["data"]][[1]],
                           py = "deseq_adjp", ly = "deseq_logfc")
sva_aucc

## These two tables have an aucc value of: 0.694200173169544 and correlation:

## 
##  Pearson's product-moment correlation
## 
## data:  tbl[[lx]] and tbl[[ly]]
## t = 182, df = 10860, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8633 0.8726
## sample estimates:
##    cor 
## 0.8681

shared_ids <- rownames(t_cf_monocyte_table_sva[["data"]][[1]]) %in%
  rownames(t_cf_monocyte_table_batchvisit[["data"]][[1]])
first <- t_cf_monocyte_table_sva[["data"]][[1]][shared_ids, ]
second <- t_cf_monocyte_table_batchvisit[["data"]][[1]][rownames(first), ]
cor.test(first[["deseq_logfc"]], second[["deseq_logfc"]])

## 
##  Pearson's product-moment correlation
## 
## data:  first[["deseq_logfc"]] and second[["deseq_logfc"]]
## t = 182, df = 10860, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8633 0.8726
## sample estimates:
##    cor 
## 0.8681

8.0.2 Neutrophil samples

Switch context to the Neutrophils, once again repeat the analysis using SVA and visit as a batch factor.

t_cf_neutrophil_de_sva <- all_pairwise(t_neutrophils, model_batch = "svaseq",
                                       parallel = parallel, filter = TRUE,
                                       methods = methods)

## 
##    tumaco_cure tumaco_failure 
##             20             21

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

## Setting 754 low elements to zero.

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

t_cf_neutrophil_de_sva

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

## The logFC agreement among the methods follows:

##                 tmc_flr___
## limma_vs_deseq      0.8742
## limma_vs_edger      0.8784
## limma_vs_basic      0.9321
## limma_vs_noiseq    -0.8237
## deseq_vs_edger      0.9994
## deseq_vs_basic      0.8754
## deseq_vs_noiseq    -0.8196
## edger_vs_basic      0.8812
## edger_vs_noiseq    -0.8268
## basic_vs_noiseq    -0.8903

t_cf_neutrophil_table_sva <- combine_de_tables(
  t_cf_neutrophil_de_sva, keepers = t_cf_contrast,
  excel = glue("{cf_prefix}/Neutrophils/t_neutrophil_cf_table_sva-v{ver}.xlsx"))
t_cf_neutrophil_table_sva

## A set of combined differential expression results.

##                           table deseq_sigup deseq_sigdown edger_sigup
## 1 tumaco_failure_vs_tumaco_cure         130            30         120
##   edger_sigdown limma_sigup limma_sigdown
## 1            27          12            12

## `geom_line()`: Each group consists of only one observation.
## i Do you need to adjust the group aesthetic?

## Plot describing unique/shared genes in a differential expression table.

t_cf_neutrophil_sig_sva <- extract_significant_genes(
  t_cf_neutrophil_table_sva,
  excel = glue("{cf_prefix}/Neutrophils/t_neutrophil_cf_sig_sva-v{ver}.xlsx"))
t_cf_neutrophil_sig_sva

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##         limma_up limma_down edger_up edger_down deseq_up deseq_down basic_up
## outcome       12         12      120         27      130         30        4
##         basic_down
## outcome          2

dim(t_cf_neutrophil_sig_sva$deseq$ups[[1]])

## [1] 130  59

dim(t_cf_neutrophil_sig_sva$deseq$downs[[1]])

## [1] 30 59

t_cf_neutrophil_de_batchvisit <- all_pairwise(t_neutrophils, model_batch = TRUE,
                                              parallel = parallel, filter = TRUE,
                                              methods = methods)

## 
##    tumaco_cure tumaco_failure 
##             20             21 
## 
##  3  2  1 
## 12 13 16

t_cf_neutrophil_de_batchvisit

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: batch in model/limma.

## The primary analysis performed 10 comparisons.

## The logFC agreement among the methods follows:

##                 tmc_flr___
## limma_vs_deseq      0.8380
## limma_vs_edger      0.8401
## limma_vs_basic      0.9658
## limma_vs_noiseq    -0.8544
## deseq_vs_edger      0.9999
## deseq_vs_basic      0.8644
## deseq_vs_noiseq    -0.8164
## edger_vs_basic      0.8671
## edger_vs_noiseq    -0.8202
## basic_vs_noiseq    -0.8903

t_cf_neutrophil_table_batchvisit <- combine_de_tables(
  t_cf_neutrophil_de_batchvisit, keepers = t_cf_contrast,
  excel = glue("{cf_prefix}/Neutrophils/t_neutrophil_cf_table_batchvisit-v{ver}.xlsx"))
t_cf_neutrophil_table_batchvisit

## A set of combined differential expression results.

##                           table deseq_sigup deseq_sigdown edger_sigup
## 1 tumaco_failure_vs_tumaco_cure          92            47         101
##   edger_sigdown limma_sigup limma_sigdown
## 1            44           3             1

## `geom_line()`: Each group consists of only one observation.
## i Do you need to adjust the group aesthetic?

## Plot describing unique/shared genes in a differential expression table.

t_cf_neutrophil_sig_batchvisit <- extract_significant_genes(
  t_cf_neutrophil_table_batchvisit,
  excel = glue("{cf_prefix}/Neutrophils/t_neutrophil_cf_sig_batchvisit-v{ver}.xlsx"))
t_cf_neutrophil_sig_batchvisit

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##         limma_up limma_down edger_up edger_down deseq_up deseq_down basic_up
## outcome        3          1      101         44       92         47        4
##         basic_down
## outcome          2

dim(t_cf_neutrophil_sig_batchvisit$deseq$ups[[1]])

## [1] 92 59

dim(t_cf_neutrophil_sig_batchvisit$deseq$downs[[1]])

## [1] 47 59

Let us add a new block in which we test a concern: if we explicitly add visit to the model (with sva, potentially without too), will that change the results we observe? My assumption is that it should change the results very minimally; but we should make absolutely certain that this is true. The neutrophils are the place to test this first because they have some of the most variance observed in the data.

Therefore I want to have an instance of the pairwise contrast that has a model of ~ finaloutcome + visitnumber + SVs where the SVs come from an invocation of sva which also has finaloutcome + visitnumber before the null model.

In theory, all_pairwise() is able to do this via the argument alt_model, but it may be safer to do it manually in order to absolutely ensure that nothing unintended happens.

## The original pairwise invocation with sva:
## t_cf_neutrophil_de_sva <- all_pairwise(t_neutrophils, model_batch = "svaseq",
##                                        parallel = parallel, filter = TRUE,
##                                        methods = methods)
test_neutrophils <- normalize_expt(t_neutrophils, filter = "simple")

## Removing 2652 low-count genes (17300 remaining).

test_neut_design <- pData(test_neutrophils)
test_formula <- as.formula("~ 0 + finaloutcome + visitnumber")
test_model <- model.matrix(test_formula, data = test_neut_design)
## Note to self: double-check that the following line is correct.
null_formula <- as.formula("~ 0 + visitnumber")
## null_model <- test_model[, c(1, 2)]
null_model <- model.matrix(null_formula, data = test_neut_design)

linear_mtrx <- exprs(test_neutrophils)
l2_mtrx <- log2(linear_mtrx + 1)
chosen_surrogates <- sva::num.sv(dat = l2_mtrx, mod = test_model)

surrogate_result <- sva::svaseq(
  dat = linear_mtrx, n.sv = chosen_surrogates, mod = test_model, mod0 = null_model)

## Number of significant surrogate variables is:  4 
## Iteration (out of 5 ):1  2  3  4  5

model_adjust <- as.matrix(surrogate_result[["sv"]])

## I don't think the following is actually required, but it is weird to just have this
## unnamed matrix hangingout.
## Set the columns to the SV#s
colnames(model_adjust) <- c("SV1", "SV2", "SV3", "SV4")
## Set the rows the sample IDs
rownames(model_adjust) <- rownames(pData(test_neutrophils))

longer_model <- as.formula("~ finaloutcome + visitnumber + SV1 + SV2 + SV3 + SV4")
neut_design_svs <- cbind(test_neut_design, model_adjust)
summarized <- DESeq2::DESeqDataSetFromMatrix(countData = linear_mtrx,
                                             colData = neut_design_svs,
                                             design = longer_model)

## converting counts to integer mode

deseq_run <- DESeq2::DESeq(summarized)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

deseq_table <- as.data.frame(DESeq2::results(object = deseq_run,
                                             contrast = c("finaloutcome", "failure", "cure"),
                                             format = "DataFrame"))

## We should be able to directly compare this to the the deseq columns from the above
## data structure named: t_cf_neutrophil_table_sva

big_table <- t_cf_neutrophil_table_sva[["data"]][["outcome"]]
only_deseq <- big_table[, c("deseq_logfc", "deseq_adjp")]
merged <- merge(deseq_table, only_deseq, by = "row.names")
rownames(merged) <- merged[["Row.names"]]
merged[["Row.names"]] <- NULL

cor_value <- cor.test(merged[["log2FoldChange"]], merged[["deseq_logfc"]])
cor_value

## 
##  Pearson's product-moment correlation
## 
## data:  merged[["log2FoldChange"]] and merged[["deseq_logfc"]]
## t = 393, df = 9099, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9706 0.9729
## sample estimates:
##    cor 
## 0.9718

logfc_plotter <- plot_linear_scatter(merged[, c("log2FoldChange", "deseq_logfc")])
logfc_plot <- logfc_plotter[["scatter"]] +
  xlab("DESeq2 log2FC: Visit explicitly in model") +
  ylab("DESeq2 log2FC: Default pairwise comparison") +
  ggtitle(glue("Comparing results from models: {prettyNum(cor_value[['estimate']])} (pearson)"))
pp(file = "figures/compare_cf_and_visit_in_model_neutrophil_logfc.svg")
logfc_plot
dev.off()

## png 
##   2

logfc_plot

cor_value <- cor.test(merged[["padj"]], merged[["deseq_adjp"]], method = "spearman")

## Warning in cor.test.default(merged[["padj"]], merged[["deseq_adjp"]], method =
## "spearman"): Cannot compute exact p-value with ties

cor_value

## 
##  Spearman's rank correlation rho
## 
## data:  merged[["padj"]] and merged[["deseq_adjp"]]
## S = 1e+10, p-value <2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##    rho 
## 0.9202

adjp_plotter <- plot_linear_scatter(merged[, c("padj", "deseq_adjp")])
adjp_plot <- adjp_plotter[["scatter"]] +
  xlab("DESeq2 adjp: Visit explicitly in model") +
  ylab("DESeq2 adjp: Default pairwise comparison") +
  ggtitle(glue("Comparing results from models: {prettyNum(cor_value[['estimate']])} (spearman)"))
pp(file = "images/compare_cf_and_visit_in_model_neutrophil_adjp.svg")
adjp_plot
dev.off()

## png 
##   2

adjp_plot

previous_sig_idx <- big_table[["deseq_adjp"]] <= 0.05 & abs(big_table[["deseq_logfc"]] >= 1.0)
summary(previous_sig_idx)

##    Mode   FALSE    TRUE 
## logical    8971     130

previous_genes <- rownames(big_table)[previous_sig_idx]

new_sig_idx <- abs(deseq_table[["log2FoldChange"]]) >= 1.0 & deseq_table[["padj"]] < 0.05
new_genes <- rownames(deseq_table)[new_sig_idx]
na_idx <- is.na(new_genes)
new_genes <- new_genes[!na_idx]

Vennerable::Venn(list("previous" = previous_genes, "new" = new_genes))

## A Venn object on 2 sets named
## previous,new 
## 00 10 01 11 
##  0 51 92 79

test_new <- simple_gprofiler(new_genes)
test_new

## A set of ontologies produced by gprofiler using 171
## genes against the hsapiens annotations and significance cutoff 0.05.
## There are: 
## 1 MF
## 12 BP
## 0 KEGG
## 2 REAC
## 0 WP
## 3 TF
## 2 MIRNA
## 0 HPA
## 0 CORUM
## 0 HP hits.

test_old <- simple_gprofiler(previous_genes)
test_old

## A set of ontologies produced by gprofiler using 130
## genes against the hsapiens annotations and significance cutoff 0.05.
## There are: 
## 4 MF
## 67 BP
## 0 KEGG
## 5 REAC
## 2 WP
## 57 TF
## 0 MIRNA
## 0 HPA
## 0 CORUM
## 0 HP hits.

new_annotated <- merge(fData(t_neutrophils), deseq_table, by = "row.names")
rownames(new_annotated) <- new_annotated[["Row.names"]]
new_annotated[["Row.names"]] <- NULL
write_xlsx(data = new_annotated, excel = "excel/neutrophil_visit_in_model_sva_cf_new.xlsx")

## write_xlsx() wrote excel/neutrophil_visit_in_model_sva_cf_new.xlsx.
## The cursor is on sheet first, row: 17303 column: 23.

old_annotated <- merge(fData(t_neutrophils), big_table, by = "row.names")
rownames(old_annotated) <- old_annotated[["Row.names"]]
old_annotated[["Row.names"]] <- NULL
write_xlsx(data = old_annotated, excel = "excel/neutrophil_visit_in_model_sva_cf_old.xlsx")

## write_xlsx() wrote excel/neutrophil_visit_in_model_sva_cf_old.xlsx.
## The cursor is on sheet first, row: 9104 column: 76.

8.0.2.1 Neutrophils by visit

When I did this with the monocytes, I split it up into multiple blocks for each visit. This time I am just going to run them all together.

visitcf_factor <- paste0("v", pData(t_neutrophils)[["visitnumber"]],
                         pData(t_neutrophils)[["finaloutcome"]])
t_neutrophil_visitcf <- set_expt_conditions(t_neutrophils, fact=visitcf_factor)

## The numbers of samples by condition are:

## 
##    v1cure v1failure    v2cure v2failure    v3cure v3failure 
##         8         8         7         6         5         7

t_cf_neutrophil_visits_de_sva <- all_pairwise(t_neutrophil_visitcf, model_batch = "svaseq",
                                              parallel = parallel, filter = TRUE,
                                              methods = methods)

## 
##    v1cure v1failure    v2cure v2failure    v3cure v3failure 
##         8         8         7         6         5         7

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

## Setting 685 low elements to zero.

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

t_cf_neutrophil_visits_de_sva

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

t_cf_neutrophil_visits_table_sva <- combine_de_tables(
  t_cf_neutrophil_visits_de_sva, keepers = visitcf_contrasts,
  excel = glue("{cf_prefix}/Neutrophils/t_neutrophil_visitcf_table_sva-v{ver}.xlsx"))

## The keepers has no elements in the coefficients.

## Here are the keepers: v1_failure, v1_cure, v2_failure, v2_cure, v3_failure, v3_cure

## Here are the coefficients: v1failure, v1cure, v2cure, v1cure, v2failure, v1cure, v3cure, v1cure, v3failure, v1cure, v2cure, v1failure, v2failure, v1failure, v3cure, v1failure, v3failure, v1failure, v2failure, v2cure, v3cure, v2cure, v3failure, v2cure, v3cure, v2failure, v3failure, v2failure, v3failure, v3cure

## Error in extract_keepers(extracted, keepers, table_names, all_coefficients, : Unable to find the set of contrasts to keep, fix this and try again.

t_cf_neutrophil_visits_table_sva

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

t_cf_neutrophil_visits_sig_sva <- extract_significant_genes(
  t_cf_neutrophil_visits_table_sva,
  excel = glue("{cf_prefix}/Neutrophils/t_neutrophil_visitcf_sig_sva-v{ver}.xlsx"))

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

t_cf_neutrophil_visits_sig_sva

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

dim(t_cf_neutrophil_visits_sig_sva$deseq$ups[[1]])

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

dim(t_cf_neutrophil_visits_sig_sva$deseq$downs[[1]])

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

Now V1

t_cf_neutrophil_v1_de_sva <- all_pairwise(tv1_neutrophils, model_batch = "svaseq",
                                          parallel = parallel, filter = TRUE,
                                          methods = methods)

## 
##    tumaco_cure tumaco_failure 
##              8              8

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

## Setting 145 low elements to zero.

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

t_cf_neutrophil_v1_de_sva

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

## The logFC agreement among the methods follows:

##                 tmc_flr___
## limma_vs_deseq      0.8180
## limma_vs_edger      0.8319
## limma_vs_basic      0.8909
## limma_vs_noiseq    -0.8084
## deseq_vs_edger      0.9946
## deseq_vs_basic      0.8627
## deseq_vs_noiseq    -0.8288
## edger_vs_basic      0.8792
## edger_vs_noiseq    -0.8474
## basic_vs_noiseq    -0.9014

t_cf_neutrophil_v1_table_sva <- combine_de_tables(
  t_cf_neutrophil_v1_de_sva, keepers = t_cf_contrast,
  excel = glue("{cf_prefix}/Neutrophils/t_neutrophil_v1_cf_table_sva-v{ver}.xlsx"))
t_cf_neutrophil_v1_table_sva

## A set of combined differential expression results.

##                           table deseq_sigup deseq_sigdown edger_sigup
## 1 tumaco_failure_vs_tumaco_cure           5             8           5
##   edger_sigdown limma_sigup limma_sigdown
## 1            11           0             0

## `geom_line()`: Each group consists of only one observation.
## i Do you need to adjust the group aesthetic?

## Plot describing unique/shared genes in a differential expression table.

t_cf_neutrophil_v1_sig_sva <- extract_significant_genes(
  t_cf_neutrophil_v1_table_sva,
  excel = glue("{cf_prefix}/Neutrophils/t_neutrophil_v1_cf_sig_sva-v{ver}.xlsx"))
t_cf_neutrophil_v1_sig_sva

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##         limma_up limma_down edger_up edger_down deseq_up deseq_down basic_up
## outcome        0          0        5         11        5          8        0
##         basic_down
## outcome          0

dim(t_cf_neutrophil_v1_sig_sva$deseq$ups[[1]])

## [1]  5 59

dim(t_cf_neutrophil_v1_sig_sva$deseq$downs[[1]])

## [1]  8 59

t_cf_neutrophil_v2_de_sva <- all_pairwise(tv2_neutrophils, model_batch = "svaseq",
                                          parallel = parallel, filter = TRUE,
                                          methods = methods)

## 
##    tumaco_cure tumaco_failure 
##              7              6

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

## Setting 78 low elements to zero.

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

t_cf_neutrophil_v2_de_sva

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

## The logFC agreement among the methods follows:

##                 tmc_flr___
## limma_vs_deseq      0.8893
## limma_vs_edger      0.8880
## limma_vs_basic      0.9485
## limma_vs_noiseq    -0.8394
## deseq_vs_edger      0.9986
## deseq_vs_basic      0.9021
## deseq_vs_noiseq    -0.9024
## edger_vs_basic      0.9026
## edger_vs_noiseq    -0.9009
## basic_vs_noiseq    -0.8742

t_cf_neutrophil_v2_table_sva <- combine_de_tables(
  t_cf_neutrophil_v2_de_sva,
  keepers = t_cf_contrast,
  excel = glue("{cf_prefix}/Neutrophils/t_neutrophil_v2_cf_table_sva-v{ver}.xlsx"))
t_cf_neutrophil_v2_table_sva

## A set of combined differential expression results.

##                           table deseq_sigup deseq_sigdown edger_sigup
## 1 tumaco_failure_vs_tumaco_cure           9             3          20
##   edger_sigdown limma_sigup limma_sigdown
## 1             6           0             0

## `geom_line()`: Each group consists of only one observation.
## i Do you need to adjust the group aesthetic?

## Plot describing unique/shared genes in a differential expression table.

t_cf_neutrophil_v2_sig_sva <- extract_significant_genes(
  t_cf_neutrophil_v2_table_sva,
  excel = glue("{cf_prefix}/Neutrophils/t_neutrophil_v2_cf_sig_sva-v{ver}.xlsx"))
t_cf_neutrophil_v2_sig_sva

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##         limma_up limma_down edger_up edger_down deseq_up deseq_down basic_up
## outcome        0          0       20          6        9          3        0
##         basic_down
## outcome          0

dim(t_cf_neutrophil_v2_sig_sva$deseq$ups[[1]])

## [1]  9 59

dim(t_cf_neutrophil_v2_sig_sva$deseq$downs[[1]])

## [1]  3 59

t_cf_neutrophil_v3_de_sva <- all_pairwise(tv3_neutrophils, model_batch = "svaseq",
                                          parallel = parallel, filter = TRUE,
                                          methods = methods)

## 
##    tumaco_cure tumaco_failure 
##              5              7

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

## Setting 83 low elements to zero.

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

t_cf_neutrophil_v3_de_sva

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

## The logFC agreement among the methods follows:

##                 tmc_flr___
## limma_vs_deseq      0.8849
## limma_vs_edger      0.8859
## limma_vs_basic      0.7952
## limma_vs_noiseq    -0.7095
## deseq_vs_edger      0.9993
## deseq_vs_basic      0.7528
## deseq_vs_noiseq    -0.7574
## edger_vs_basic      0.7514
## edger_vs_noiseq    -0.7563
## basic_vs_noiseq    -0.8892

t_cf_neutrophil_v3_table_sva <- combine_de_tables(
  t_cf_neutrophil_v3_de_sva, keepers = t_cf_contrast,
  excel = glue("{cf_prefix}/Neutrophils/t_neutrophil_v3_cf_table_sva-v{ver}.xlsx"))
t_cf_neutrophil_v3_table_sva

## A set of combined differential expression results.

##                           table deseq_sigup deseq_sigdown edger_sigup
## 1 tumaco_failure_vs_tumaco_cure           5             1           5
##   edger_sigdown limma_sigup limma_sigdown
## 1             1           0             0

## `geom_line()`: Each group consists of only one observation.
## i Do you need to adjust the group aesthetic?

## Plot describing unique/shared genes in a differential expression table.

t_cf_neutrophil_v3_sig_sva <- extract_significant_genes(
  t_cf_neutrophil_v3_table_sva,
  excel = glue("{cf_prefix}/Neutrophils/t_neutrophil_v3_cf_sig_sva-v{ver}.xlsx"))
t_cf_neutrophil_v3_sig_sva

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##         limma_up limma_down edger_up edger_down deseq_up deseq_down basic_up
## outcome        0          0        5          1        5          1        0
##         basic_down
## outcome          0

dim(t_cf_neutrophil_v3_sig_sva$deseq$ups[[1]])

## [1]  5 59

dim(t_cf_neutrophil_v3_sig_sva$deseq$downs[[1]])

## [1]  1 59

8.0.2.2 Neutrophils: Compare sva to batch-in-model

sva_aucc <- calculate_aucc(t_cf_neutrophil_table_sva[["data"]][[1]],
                           tbl2 = t_cf_neutrophil_table_batchvisit[["data"]][[1]],
                           py = "deseq_adjp", ly = "deseq_logfc")
sva_aucc

## These two tables have an aucc value of: 0.673209505652166 and correlation:

## 
##  Pearson's product-moment correlation
## 
## data:  tbl[[lx]] and tbl[[ly]]
## t = 209, df = 9099, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9060 0.9131
## sample estimates:
##    cor 
## 0.9096

shared_ids <- rownames(t_cf_neutrophil_table_sva[["data"]][[1]]) %in%
  rownames(t_cf_neutrophil_table_batchvisit[["data"]][[1]])
first <- t_cf_neutrophil_table_sva[["data"]][[1]][shared_ids, ]
second <- t_cf_neutrophil_table_batchvisit[["data"]][[1]][rownames(first), ]
cor.test(first[["deseq_logfc"]], second[["deseq_logfc"]])

## 
##  Pearson's product-moment correlation
## 
## data:  first[["deseq_logfc"]] and second[["deseq_logfc"]]
## t = 209, df = 9099, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9060 0.9131
## sample estimates:
##    cor 
## 0.9096

8.0.3 Eosinophils

This time, with feeling! Repeating the same set of tasks with the eosinophil samples.

t_cf_eosinophil_de_sva <- all_pairwise(t_eosinophils, model_batch = "svaseq",
                                       parallel = parallel, filter = TRUE,
                                       methods = methods)

## 
##    tumaco_cure tumaco_failure 
##             17              9

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

## Setting 327 low elements to zero.

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

t_cf_eosinophil_de_sva

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

## The logFC agreement among the methods follows:

##                 tmc_flr___
## limma_vs_deseq      0.9099
## limma_vs_edger      0.9174
## limma_vs_basic      0.8756
## limma_vs_noiseq    -0.7909
## deseq_vs_edger      0.9973
## deseq_vs_basic      0.8488
## deseq_vs_noiseq    -0.7864
## edger_vs_basic      0.8546
## edger_vs_noiseq    -0.7964
## basic_vs_noiseq    -0.8713

t_cf_eosinophil_table_sva <- combine_de_tables(
  t_cf_eosinophil_de_sva, keepers = t_cf_contrast,
  excel = glue("{cf_prefix}/Eosinophils/t_eosinophil_cf_table_sva-v{ver}.xlsx"))
t_cf_eosinophil_table_sva

## A set of combined differential expression results.

##                           table deseq_sigup deseq_sigdown edger_sigup
## 1 tumaco_failure_vs_tumaco_cure         116            75         112
##   edger_sigdown limma_sigup limma_sigdown
## 1            63          57            34

## `geom_line()`: Each group consists of only one observation.
## i Do you need to adjust the group aesthetic?

## Plot describing unique/shared genes in a differential expression table.

t_cf_eosinophil_sig_sva <- extract_significant_genes(
  t_cf_eosinophil_table_sva,
  excel = glue("{cf_prefix}/Eosinophils/t_eosinophil_cf_sig_sva-v{ver}.xlsx"))
t_cf_eosinophil_sig_sva

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##         limma_up limma_down edger_up edger_down deseq_up deseq_down basic_up
## outcome       57         34      112         63      116         75        0
##         basic_down
## outcome          0

dim(t_cf_eosinophil_sig_sva$deseq$ups[[1]])

## [1] 116  59

dim(t_cf_eosinophil_sig_sva$deseq$downs[[1]])

## [1] 75 59

knitr::kable(head(t_cf_eosinophil_sig_sva$deseq$ups[[1]]))

	ensembl_gene_id	ensembl_transcript_id	version	transcript_version	description	gene_biotype	cds_length	chromosome_name	strand	start_position	end_position	hgnc_symbol	uniprot_gn_symbol	transcript	mean_cds_len	deseq_logfc	deseq_adjp	edger_logfc	edger_adjp	limma_logfc	limma_adjp	noiseq_logfc	noiseq_adjp	basic_num	basic_den	basic_numvar	basic_denvar	basic_logfc	basic_t	basic_p	basic_adjp	deseq_basemean	deseq_lfcse	deseq_stat	deseq_p	deseq_num	deseq_den	edger_logcpm	edger_lr	edger_p	limma_ave	limma_t	limma_b	limma_p	noiseq_num	noiseq_den	noiseq_theta	noiseq_prob	noiseq_p	limma_adjp_fdr	deseq_adjp_fdr	edger_adjp_fdr	basic_adjp_fdr	noiseq_adjp_fdr	lfc_meta	lfc_var	lfc_varbymed	p_meta	p_var
ENSG00000198178	ENSG00000198178	ENST00000537530	10	1	C-type lectin domain family 4 member C [Source:HGNC Symbol;Acc:HGNC:13258]	protein_coding	267	12	-	7729415	7751605	CLEC4C	CLEC4C	ENSG00000198178.1	510.5	5.537	2e-03	5.152	0.1934	4.225	0.0191	-1.4675	1	2.0400	-2.764	8.784	8.946	4.804	3.920	0.0012	0.1398	193.80	1.3030	4.249	0e+00	9.617	4.0806	2.2170	5.953	0.0147	-1.4190	4.4420	0.1256	0.0002	3.580	9.899	-1.445	0.9473	0.0527	0.0191	0.0020	0.1934	0.1398	0.8357	4.982	4.880e-03	9.796e-04	4.956e-03	7.107e-05
ENSG00000187569	ENSG00000187569	ENST00000345088	3	3	developmental pluripotency associated 3 [Source:HGNC Symbol;Acc:HGNC:19199]	protein_coding	480	12	+	7711433	7717559	DPPA3	DPPA3	ENSG00000187569.3	480	5.448	7e-03	4.721	0.0547	3.504	0.0470	-0.6819	1	-0.6872	-4.301	7.263	2.839	3.614	3.662	0.0035	0.1989	23.10	1.4300	3.810	1e-04	5.866	0.4183	-0.5913	9.779	0.0018	-3.3720	3.6990	-1.6550	0.0011	2.088	3.349	-1.177	0.9103	0.0897	0.0470	0.0069	0.0547	0.1989	0.8357	4.443	4.098e-01	9.225e-02	9.866e-04	6.647e-07
ENSG00000136235	ENSG00000136235	ENST00000479625	16	1	glycoprotein nmb [Source:HGNC Symbol;Acc:HGNC:4462]	protein_coding	undefined	7	+	23235967	23275108	GPNMB	GPNMB	ENSG00000136235.1	1447.5	5.410	4e-04	5.374	0.0001	3.867	0.1754	-1.2238	1	-1.1190	-3.695	12.101	2.665	2.576	2.102	0.0621	0.4987	53.03	1.1380	4.752	0e+00	6.906	1.4965	0.4580	25.540	0.0000	-3.2370	2.5210	-3.2570	0.0184	2.100	4.905	-1.087	0.9259	0.0741	0.1754	0.0004	0.0001	0.4987	0.8357	4.798	5.245e-02	1.093e-02	6.124e-03	1.125e-04
ENSG00000089012	ENSG00000089012	ENST00000497407	14	2	signal regulatory protein gamma [Source:HGNC Symbol;Acc:HGNC:15757]	protein_coding	undefined	20	-	1629152	1657779	SIRPG	SIRPG	ENSG00000089012.2	880.8	4.040	0e+00	4.018	0.0000	1.598	0.6625	-2.7349	1	0.8317	-1.681	13.876	1.355	2.513	1.974	0.0805	0.5427	272.50	0.7310	5.526	0e+00	7.574	3.5336	2.7060	32.750	0.0000	-1.1480	0.9807	-5.0020	0.3360	2.385	15.874	-2.097	0.9883	0.0117	0.6624	0.0000	0.0000	0.5427	0.8357	3.148	1.579e+00	5.016e-01	1.120e-01	3.763e-02
ENSG00000089127	ENSG00000089127	ENST00000540589	13	2	2’-5’-oligoadenylate synthetase 1 [Source:HGNC Symbol;Acc:HGNC:8086]	protein_coding	68	12	+	112906783	112933222	OAS1	OAS1	ENSG00000089127.2	682.8	3.933	0e+00	3.943	0.0000	3.339	0.0596	-2.1220	1	1.9510	-1.301	8.237	1.116	3.252	3.284	0.0092	0.2669	184.60	0.5478	7.180	0e+00	7.841	3.9081	2.1560	44.580	0.0000	-0.4596	3.4950	-1.3000	0.0018	2.522	10.978	-2.205	0.9849	0.0151	0.0596	0.0000	0.0000	0.2669	0.8357	3.722	1.794e-02	4.820e-03	5.900e-04	1.044e-06
ENSG00000137959	ENSG00000137959	ENST00000450498	16	1	interferon induced protein 44 like [Source:HGNC Symbol;Acc:HGNC:17817]	protein_coding	699	1	+	78619922	78646145	IFI44L	IFI44L	ENSG00000137959.1	783.333333333333	3.828	0e+00	3.831	0.0000	3.443	0.0123	-3.6168	0	5.5560	1.793	6.584	3.318	3.763	3.909	0.0020	0.1645	1932.00	0.5401	7.087	0e+00	11.304	7.4755	5.4900	57.400	0.0000	3.0090	4.7380	1.6850	0.0001	7.837	96.140	-3.316	1.0000	0.0000	0.0123	0.0000	0.0000	0.1645	0.0000	3.691	2.660e-03	7.207e-04	2.392e-05	1.716e-09

knitr::kable(head(t_cf_eosinophil_sig_sva$deseq$downs[[1]]))

	ensembl_gene_id	ensembl_transcript_id	version	transcript_version	description	gene_biotype	cds_length	chromosome_name	strand	start_position	end_position	hgnc_symbol	uniprot_gn_symbol	transcript	mean_cds_len	deseq_logfc	deseq_adjp	edger_logfc	edger_adjp	limma_logfc	limma_adjp	noiseq_logfc	noiseq_adjp	basic_num	basic_den	basic_numvar	basic_denvar	basic_logfc	basic_t	basic_p	basic_adjp	deseq_basemean	deseq_lfcse	deseq_stat	deseq_p	deseq_num	deseq_den	edger_logcpm	edger_lr	edger_p	limma_ave	limma_t	limma_b	limma_p	noiseq_num	noiseq_den	noiseq_theta	noiseq_prob	noiseq_p	limma_adjp_fdr	deseq_adjp_fdr	edger_adjp_fdr	basic_adjp_fdr	noiseq_adjp_fdr	lfc_meta	lfc_var	lfc_varbymed	p_meta	p_var
ENSG00000179344	ENSG00000179344	ENST00000399084	16	5	major histocompatibility complex, class II, DQ beta 1 [Source:HGNC Symbol;Acc:HGNC:4944]	protein_coding	786	6	-	32659467	32668383	HLA-DQB1	HLA-DQB1	ENSG00000179344.5	645.5	-5.676	0.0000	-5.668	0.0000	-7.612	0.0155	4.4060	0	0.0597	5.6120	5.9251	7.993	-5.552	-5.227	0.0001	0.0515	4151.00	0.8485	-6.689	0.0000	6.782	12.458	6.575	36.280	0.0000	3.5810	-4.604	1.291	0.0001	104.486	4.928	4.5901	1.0000	0.0000	0.0155	0.0000	0.0000	0.0515	0.0000	-5.952	1.301e+00	-2.186e-01	3.393e-05	3.454e-09
ENSG00000112139	ENSG00000112139	ENST00000515437	16	5	MAM domain containing glycosylphosphatidylinositol anchor 1 [Source:HGNC Symbol;Acc:HGNC:19267]	protein_coding	388	6	-	37630679	37699306	MDGA1	MDGA1	ENSG00000112139.5	1438.71428571429	-5.037	0.0015	-4.942	0.0398	-2.844	0.2584	0.7507	1	-3.3850	-0.1629	10.6461	14.783	-3.222	-2.249	0.0366	0.4206	149.80	1.1640	-4.327	0.0000	2.708	7.744	1.813	10.640	0.0011	-1.5710	-2.141	-3.687	0.0421	5.239	3.114	0.7443	0.8075	0.1925	0.2584	0.0015	0.0398	0.4206	0.8357	-3.895	9.304e-01	-2.389e-01	1.440e-02	5.749e-04
ENSG00000203972	ENSG00000203972	ENST00000545705	10	1	glycine-N-acyltransferase like 3 [Source:HGNC Symbol;Acc:HGNC:21349]	protein_coding	468	6	+	49499923	49528078	GLYATL3	GLYATL3	ENSG00000203972.1	667.5	-4.718	0.0498	-4.599	0.0363	-2.962	0.1916	0.4092	1	-5.7280	-3.3770	0.5218	6.712	-2.351	-3.493	0.0023	0.1675	27.99	1.5560	-3.032	0.0024	-1.881	2.837	-0.443	10.870	0.0010	-4.5450	-2.444	-3.577	0.0218	2.711	2.042	0.5496	0.6925	0.3075	0.1915	0.0498	0.0363	0.1675	0.9068	-3.901	8.898e-02	-2.281e-02	8.413e-03	1.355e-04
ENSG00000196526	ENSG00000196526	ENST00000358461	10	6	actin filament associated protein 1 [Source:HGNC Symbol;Acc:HGNC:24017]	protein_coding	2193	4	-	7758714	7939926	AFAP1	AFAP1	ENSG00000196526.6	1911	-3.294	0.0252	-3.293	0.0574	-2.538	0.2974	2.2170	1	0.5888	2.6700	2.1638	11.578	-2.081	-2.168	0.0405	0.4335	982.20	0.9856	-3.342	0.0008	6.891	10.185	4.492	9.604	0.0019	1.8040	-1.999	-4.138	0.0565	23.858	5.131	2.1394	0.9761	0.0239	0.2974	0.0252	0.0574	0.4335	0.8357	-2.956	1.507e-01	-5.097e-02	1.976e-02	1.013e-03
ENSG00000175592	ENSG00000175592	ENST00000312562	9	7	FOS like 1, AP-1 transcription factor subunit [Source:HGNC Symbol;Acc:HGNC:13718]	protein_coding	816	11	-	65892049	65900573	FOSL1	FOSL1	ENSG00000175592.7	496	-3.097	0.0000	-3.081	0.0000	-2.221	0.1738	0.7531	1	0.1522	1.8020	3.6253	4.212	-1.650	-2.045	0.0561	0.4818	267.80	0.5692	-5.441	0.0000	5.087	8.184	2.640	30.640	0.0000	1.0720	-2.528	-3.082	0.0181	6.345	3.764	1.1519	0.9081	0.0919	0.1738	0.0000	0.0000	0.4818	0.8357	-2.752	2.007e-01	-7.295e-02	6.027e-03	1.090e-04
ENSG00000122877	ENSG00000122877	ENST00000637191	16	1	early growth response 2 [Source:HGNC Symbol;Acc:HGNC:3239]	protein_coding	418	10	-	62811996	62919900	EGR2	EGR2	ENSG00000122877.1	1140.25	-2.789	0.0117	-2.779	0.0110	-2.011	0.2749	0.6118	1	-1.2330	-0.0550	0.8328	5.118	-1.178	-1.878	0.0731	0.5187	96.53	0.7679	-3.632	0.0003	3.932	6.721	1.188	14.120	0.0002	-0.7511	-2.078	-3.767	0.0480	3.917	2.563	0.7938	0.8090	0.1910	0.2749	0.0117	0.0110	0.5187	0.8357	-2.514	2.515e-01	-1.000e-01	1.616e-02	7.617e-04

Repeat with batch in the model.

t_cf_eosinophil_de_batchvisit <- all_pairwise(t_eosinophils, model_batch = TRUE,
                                              parallel = parallel, filter = TRUE,
                                              methods = methods)

## 
##    tumaco_cure tumaco_failure 
##             17              9 
## 
## 3 2 1 
## 9 9 8

t_cf_eosinophil_de_batchvisit

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: batch in model/limma.

## The primary analysis performed 10 comparisons.

## The logFC agreement among the methods follows:

##                 tmc_flr___
## limma_vs_deseq      0.8678
## limma_vs_edger      0.8696
## limma_vs_basic      0.9676
## limma_vs_noiseq    -0.8444
## deseq_vs_edger      0.9998
## deseq_vs_basic      0.8961
## deseq_vs_noiseq    -0.8351
## edger_vs_basic      0.8977
## edger_vs_noiseq    -0.8396
## basic_vs_noiseq    -0.8713

t_cf_eosinophil_table_batchvisit <- combine_de_tables(
  t_cf_eosinophil_de_batchvisit, keepers = t_cf_contrast,
  excel = glue("{cf_prefix}/Eosinophils/t_eosinophil_cf_table_batchvisit-v{ver}.xlsx"))
t_cf_eosinophil_table_batchvisit

## A set of combined differential expression results.

##                           table deseq_sigup deseq_sigdown edger_sigup
## 1 tumaco_failure_vs_tumaco_cure          99            35         103
##   edger_sigdown limma_sigup limma_sigdown
## 1            24          35            15

## `geom_line()`: Each group consists of only one observation.
## i Do you need to adjust the group aesthetic?

## Plot describing unique/shared genes in a differential expression table.

t_cf_eosinophil_sig_batchvisit <- extract_significant_genes(
  t_cf_eosinophil_table_batchvisit,
  excel = glue("{cf_prefix}/Eosinophils/t_eosinophil_cf_sig_batchvisit-v{ver}.xlsx"))
t_cf_eosinophil_sig_batchvisit

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##         limma_up limma_down edger_up edger_down deseq_up deseq_down basic_up
## outcome       35         15      103         24       99         35        0
##         basic_down
## outcome          0

dim(t_cf_eosinophil_sig_batchvisit$deseq$ups[[1]])

## [1] 99 59

dim(t_cf_eosinophil_sig_batchvisit$deseq$downs[[1]])

## [1] 35 59

knitr::kable(head(t_cf_eosinophil_sig_batchvisit$deseq$ups[[1]]))

	ensembl_gene_id	ensembl_transcript_id	version	transcript_version	description	gene_biotype	cds_length	chromosome_name	strand	start_position	end_position	hgnc_symbol	uniprot_gn_symbol	transcript	mean_cds_len	deseq_logfc	deseq_adjp	edger_logfc	edger_adjp	limma_logfc	limma_adjp	noiseq_logfc	noiseq_adjp	basic_num	basic_den	basic_numvar	basic_denvar	basic_logfc	basic_t	basic_p	basic_adjp	deseq_basemean	deseq_lfcse	deseq_stat	deseq_p	deseq_num	deseq_den	edger_logcpm	edger_lr	edger_p	limma_ave	limma_t	limma_b	limma_p	noiseq_num	noiseq_den	noiseq_theta	noiseq_prob	noiseq_p	limma_adjp_fdr	deseq_adjp_fdr	edger_adjp_fdr	basic_adjp_fdr	noiseq_adjp_fdr	lfc_meta	lfc_var	lfc_varbymed	p_meta	p_var
ENSG00000165949	ENSG00000165949	ENST00000611954	12	4	interferon alpha inducible protein 27 [Source:HGNC Symbol;Acc:HGNC:5397]	protein_coding	180	14	+	94104836	94116698	IFI27	IFI27	ENSG00000165949.4	287.454545454545	5.668	0e+00	5.624	0e+00	4.397	0.0591	-1.2302	1	-0.1291	-3.3160	7.379	1.696	3.186	3.323	0.0077	0.2504	54.63	0.8188	6.923	0	7.505	1.8370	0.4882	47.95	0	-2.6550	3.858	-1.1070	0.0007	2.125	4.985	-1.156	0.9137	0.0863	0.0591	0e+00	0e+00	0.2504	0.8357	5.230	0.000e+00	0.000e+00	2.264e-04	1.537e-07
ENSG00000187569	ENSG00000187569	ENST00000345088	3	3	developmental pluripotency associated 3 [Source:HGNC Symbol;Acc:HGNC:19199]	protein_coding	480	12	+	7711433	7717559	DPPA3	DPPA3	ENSG00000187569.3	480	5.537	4e-04	5.386	1e-04	4.404	0.0351	-0.6819	1	-0.6872	-4.3010	7.263	2.839	3.614	3.662	0.0035	0.1989	23.10	1.1660	4.747	0	5.915	0.3782	-0.6406	25.31	0	-3.3720	4.263	-0.7244	0.0002	2.088	3.349	-1.177	0.9103	0.0897	0.0351	4e-04	1e-04	0.1989	0.8357	5.088	3.906e-02	7.677e-03	7.965e-05	1.843e-08
ENSG00000136235	ENSG00000136235	ENST00000479625	16	1	glycoprotein nmb [Source:HGNC Symbol;Acc:HGNC:4462]	protein_coding	undefined	7	+	23235967	23275108	GPNMB	GPNMB	ENSG00000136235.1	1447.5	5.426	2e-04	5.360	0e+00	2.104	0.5947	-1.2238	1	-1.1190	-3.6950	12.101	2.665	2.576	2.102	0.0621	0.4987	53.03	1.0740	5.053	0	7.515	2.0897	0.4259	29.50	0	-3.2370	1.413	-4.4990	0.1695	2.100	4.905	-1.087	0.9259	0.0741	0.5947	2e-04	0e+00	0.4987	0.8357	4.060	2.227e+00	5.485e-01	5.650e-02	9.577e-03
ENSG00000089127	ENSG00000089127	ENST00000540589	13	2	2’-5’-oligoadenylate synthetase 1 [Source:HGNC Symbol;Acc:HGNC:8086]	protein_coding	68	12	+	112906783	112933222	OAS1	OAS1	ENSG00000089127.2	682.8	4.820	0e+00	4.830	0e+00	3.978	0.1341	-2.1220	1	1.9510	-1.3010	8.237	1.116	3.252	3.284	0.0092	0.2669	184.60	0.7144	6.746	0	8.840	4.0197	2.1430	56.76	0	-0.4596	3.160	-1.9060	0.0040	2.522	10.978	-2.205	0.9849	0.0151	0.1341	0e+00	0e+00	0.2669	0.8357	4.652	5.441e-02	1.170e-02	1.328e-03	5.293e-06
ENSG00000111335	ENSG00000111335	ENST00000551603	12	1	2’-5’-oligoadenylate synthetase 2 [Source:HGNC Symbol;Acc:HGNC:8087]	protein_coding	183	12	+	112978395	113011723	OAS2	OAS2	ENSG00000111335.1	1319.5	4.447	0e+00	4.461	0e+00	3.601	0.2035	-3.0969	0	4.0180	0.8398	12.517	3.025	3.178	2.537	0.0293	0.3953	1028.00	0.8157	5.451	0	11.444	6.9971	4.5830	38.11	0	1.7970	2.779	-2.7210	0.0100	5.227	44.719	-2.937	1.0000	0.0000	0.2035	0e+00	0e+00	0.3953	0.0000	4.110	6.269e-02	1.525e-02	3.340e-03	3.347e-05
ENSG00000137959	ENSG00000137959	ENST00000450498	16	1	interferon induced protein 44 like [Source:HGNC Symbol;Acc:HGNC:17817]	protein_coding	699	1	+	78619922	78646145	IFI44L	IFI44L	ENSG00000137959.1	783.333333333333	4.200	0e+00	4.213	0e+00	3.902	0.0422	-3.6168	0	5.5560	1.7930	6.584	3.318	3.763	3.909	0.0020	0.1645	1932.00	0.7590	5.534	0	12.253	8.0528	5.4890	40.52	0	3.0090	4.118	0.2407	0.0003	7.837	96.140	-3.316	1.0000	0.0000	0.0422	0e+00	0e+00	0.1645	0.0000	4.204	7.374e-02	1.754e-02	1.151e-04	3.976e-08

knitr::kable(head(t_cf_eosinophil_sig_batchvisit$deseq$downs[[1]]))

	ensembl_gene_id	ensembl_transcript_id	version	transcript_version	description	gene_biotype	cds_length	chromosome_name	strand	start_position	end_position	hgnc_symbol	uniprot_gn_symbol	transcript	mean_cds_len	deseq_logfc	deseq_adjp	edger_logfc	edger_adjp	limma_logfc	limma_adjp	noiseq_logfc	noiseq_adjp	basic_num	basic_den	basic_numvar	basic_denvar	basic_logfc	basic_t	basic_p	basic_adjp	deseq_basemean	deseq_lfcse	deseq_stat	deseq_p	deseq_num	deseq_den	edger_logcpm	edger_lr	edger_p	limma_ave	limma_t	limma_b	limma_p	noiseq_num	noiseq_den	noiseq_theta	noiseq_prob	noiseq_p	limma_adjp_fdr	deseq_adjp_fdr	edger_adjp_fdr	basic_adjp_fdr	noiseq_adjp_fdr	lfc_meta	lfc_var	lfc_varbymed	p_meta	p_var
ENSG00000189430	ENSG00000189430	ENST00000338835	13	9	natural cytotoxicity triggering receptor 1 [Source:HGNC Symbol;Acc:HGNC:6731]	protein_coding	864	19	+	54906148	54916140	NCR1	NCR1	ENSG00000189430.9	798.5	-5.820	0.0002	-5.752	0.0019	-3.214	0.2894	0.7731	1	-4.4710	-1.0630	1.835	11.567	-3.408	-3.624	0.0014	0.1472	93.09	1.1550	-5.038	0e+00	1.4989	7.319	1.0570	19.19	0e+00	-2.621	-2.436	-3.3850	0.0220	3.639	2.130	1.1623	0.9112	0.0888	0.2894	0.0002	0.0019	0.1472	0.8357	-5.110	1.583e+00	-3.098e-01	7.344e-03	1.615e-04
ENSG00000179344	ENSG00000179344	ENST00000399084	16	5	major histocompatibility complex, class II, DQ beta 1 [Source:HGNC Symbol;Acc:HGNC:4944]	protein_coding	786	6	-	32659467	32668383	HLA-DQB1	HLA-DQB1	ENSG00000179344.5	645.5	-5.667	0.0000	-5.653	0.0006	-5.936	0.0332	4.4060	0	0.0597	5.6120	5.925	7.993	-5.552	-5.227	0.0001	0.0515	4151.00	0.8879	-6.382	0e+00	8.4098	14.076	6.5750	21.99	0e+00	3.581	-4.357	0.8177	0.0002	104.486	4.928	4.5901	1.0000	0.0000	0.0332	0.0000	0.0006	0.0515	0.0000	-5.464	1.039e-01	-1.902e-02	6.255e-05	1.123e-08
ENSG00000162669	ENSG00000162669	ENST00000427444	16	1	helicase for meiosis 1 [Source:HGNC Symbol;Acc:HGNC:20193]	protein_coding	589	1	-	91260766	91404856	HFM1	HFM1	ENSG00000162669.1	1528.4	-4.617	0.0012	-4.588	0.0054	-2.665	0.1784	1.2318	1	-3.0190	-0.0765	2.110	8.425	-2.942	-3.444	0.0021	0.1672	207.70	1.0200	-4.526	0e+00	4.2354	8.853	2.1480	16.51	0e+00	-1.450	-2.912	-2.4740	0.0073	5.418	2.307	1.1461	0.9059	0.0941	0.1784	0.0012	0.0054	0.1672	0.8357	-3.922	2.638e-01	-6.727e-02	2.452e-03	1.764e-05
ENSG00000167634	ENSG00000167634	ENST00000328092	12	9	NLR family pyrin domain containing 7 [Source:HGNC Symbol;Acc:HGNC:22947]	protein_coding	3030	19	-	54923509	54966312	NLRP7	NLRP7	ENSG00000167634.9	2077.44444444444	-4.061	0.0071	-4.004	0.0317	-1.875	0.4968	0.3273	1	-4.1970	-2.2500	0.258	8.472	-1.947	-2.682	0.0153	0.3229	27.08	1.0170	-3.995	1e-04	0.8671	4.928	-0.7663	12.24	5e-04	-3.363	-1.721	-4.2100	0.0971	2.616	2.085	0.6656	0.7700	0.2300	0.4968	0.0071	0.0317	0.3229	0.8626	-3.179	1.172e+00	-3.687e-01	3.255e-02	3.126e-03
ENSG00000196526	ENSG00000196526	ENST00000358461	10	6	actin filament associated protein 1 [Source:HGNC Symbol;Acc:HGNC:24017]	protein_coding	2193	4	-	7758714	7939926	AFAP1	AFAP1	ENSG00000196526.6	1911	-3.879	0.0038	-3.877	0.0088	-2.357	0.3886	2.2170	1	0.5888	2.6700	2.164	11.578	-2.081	-2.168	0.0405	0.4335	982.20	0.9252	-4.192	0e+00	6.7815	10.660	4.4950	15.42	1e-04	1.804	-2.067	-4.0740	0.0489	23.858	5.131	2.1394	0.9761	0.0239	0.3886	0.0038	0.0088	0.4335	0.8357	-3.320	3.403e-01	-1.025e-01	1.633e-02	7.942e-04
ENSG00000277150	ENSG00000277150	ENST00000622749	1	1	coagulation factor VIII associated 3 [Source:HGNC Symbol;Acc:HGNC:31850]	protein_coding	1116	X	-	155456914	155458672	F8A3	F8A1	ENSG00000277150.1	1116	-3.788	0.0153	-3.704	0.0589	-1.712	0.3848	0.3180	1	-4.5610	-2.3170	1.574	6.415	-2.244	-3.020	0.0059	0.2249	21.72	1.0170	-3.724	2e-04	1.9848	5.772	-1.5770	10.78	1e-03	-3.506	-2.080	-3.8470	0.0476	2.591	2.078	0.5401	0.6847	0.3153	0.3849	0.0154	0.0588	0.2249	0.9109	-2.842	9.067e-01	-3.191e-01	1.628e-02	7.364e-04

Repeat with visit in the condition contrast.

visitcf_factor <- paste0("v", pData(t_eosinophils)[["visitnumber"]],
                         pData(t_eosinophils)[["finaloutcome"]])
t_eosinophil_visitcf <- set_expt_conditions(t_eosinophils, fact = visitcf_factor)

## The numbers of samples by condition are:

## 
##    v1cure v1failure    v2cure v2failure    v3cure v3failure 
##         5         3         6         3         6         3

t_cf_eosinophil_visits_de_sva <- all_pairwise(t_eosinophil_visitcf, model_batch = "svaseq",
                                              parallel = parallel, filter = TRUE,
                                              methods = methods)

## 
##    v1cure v1failure    v2cure v2failure    v3cure v3failure 
##         5         3         6         3         6         3

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

## Setting 373 low elements to zero.

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

t_cf_eosinophil_visits_de_sva

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

t_cf_eosinophil_visits_table_sva <- combine_de_tables(
   t_cf_eosinophil_visits_de_sva, keepers = visitcf_contrasts,
  excel = glue("{cf_prefix}/Eosinophils/t_eosinophil_visitcf_table_sva-v{ver}.xlsx"))

## The keepers has no elements in the coefficients.

## Here are the keepers: v1_failure, v1_cure, v2_failure, v2_cure, v3_failure, v3_cure

## Here are the coefficients: v1failure, v1cure, v2cure, v1cure, v2failure, v1cure, v3cure, v1cure, v3failure, v1cure, v2cure, v1failure, v2failure, v1failure, v3cure, v1failure, v3failure, v1failure, v2failure, v2cure, v3cure, v2cure, v3failure, v2cure, v3cure, v2failure, v3failure, v2failure, v3failure, v3cure

## Error in extract_keepers(extracted, keepers, table_names, all_coefficients, : Unable to find the set of contrasts to keep, fix this and try again.

t_cf_eosinophil_visits_table_sva

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

t_cf_eosinophil_visits_sig_sva <- extract_significant_genes(
  t_cf_eosinophil_visits_table_sva,
  excel = glue("{cf_prefix}/Eosinophils/t_eosinophil_visitcf_sig_sva-v{ver}.xlsx"))

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

t_cf_eosinophil_visits_sig_sva

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

dim(t_cf_eosinophil_visits_sig_sva$deseq$ups[[1]])

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

dim(t_cf_eosinophil_visits_sig_sva$deseq$downs[[1]])

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

knitr::kable(head(t_cf_eosinophil_visits_sig_sva$deseq$ups[[1]]))

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'head': object 't_cf_eosinophil_visits_sig_sva' not found

knitr::kable(head(t_cf_eosinophil_visits_sig_sva$deseq$downs[[1]]))

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'head': object 't_cf_eosinophil_visits_sig_sva' not found

8.0.3.0.1 Repeat test of different models: Eosinophils

We wish to ensure that my model simplification did not do anything incorrect to the data for all three cell types, I already did this for the neutrophils, let us repeat for the eosinophils. I am therefore (mostly) copy/pasting the neutrophil section here.a

## The original pairwise invocation with sva:
#t_cf_eosinophil_de_sva <- all_pairwise(t_eosinophils, model_batch = "svaseq",
#                                       filter = TRUE, parallel=FALSE, methods = methods)
test_eosinophils <- normalize_expt(t_eosinophils, filter = "simple")

## Removing 2652 low-count genes (17300 remaining).

test_eo_design <- pData(test_eosinophils)
test_formula <- as.formula("~ 0 + finaloutcome + visitnumber")
test_model <- model.matrix(test_formula, data = test_eo_design)
null_formula <- as.formula("~ 0 + visitnumber")
null_model <- model.matrix(null_formula, data = test_eo_design)

linear_mtrx <- exprs(test_eosinophils)
l2_mtrx <- log2(linear_mtrx + 1)
chosen_surrogates <- sva::num.sv(dat = l2_mtrx, mod = test_model)

surrogate_result <- sva::svaseq(
  dat = linear_mtrx, n.sv = chosen_surrogates, mod = test_model, mod0 = null_model)

## Number of significant surrogate variables is:  3 
## Iteration (out of 5 ):1  2  3  4  5

model_adjust <- as.matrix(surrogate_result[["sv"]])

colnames(model_adjust) <- c("SV1", "SV2", "SV3")
rownames(model_adjust) <- rownames(pData(test_eosinophils))
longer_model <- as.formula("~ finaloutcome + visitnumber + SV1 + SV2 + SV3")
eo_design_svs <- cbind(test_eo_design, model_adjust)
summarized <- DESeq2::DESeqDataSetFromMatrix(countData = linear_mtrx,
                                             colData = eo_design_svs,
                                             design = longer_model)

## converting counts to integer mode

deseq_run <- DESeq2::DESeq(summarized)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

deseq_table <- as.data.frame(DESeq2::results(object = deseq_run,
                                             contrast = c("finaloutcome", "failure", "cure"),
                                             format = "DataFrame"))

big_table <- t_cf_eosinophil_table_sva[["data"]][["outcome"]]
only_deseq <- big_table[, c("deseq_logfc", "deseq_adjp")]
merged <- merge(deseq_table, only_deseq, by = "row.names")
rownames(merged) <- merged[["Row.names"]]
merged[["Row.names"]] <- NULL

cor_value <- cor.test(merged[["log2FoldChange"]], merged[["deseq_logfc"]])
cor_value

## 
##  Pearson's product-moment correlation
## 
## data:  merged[["log2FoldChange"]] and merged[["deseq_logfc"]]
## t = 228, df = 10530, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9084 0.9149
## sample estimates:
##    cor 
## 0.9117

logfc_plotter <- plot_linear_scatter(merged[, c("log2FoldChange", "deseq_logfc")])
logfc_plot <- logfc_plotter[["scatter"]] +
  xlab("DESeq2 log2FC: Visit explicitly in model") +
  ylab("DESeq2 log2FC: Default pairwise comparison") +
  ggtitle(glue("Comparing results from models: {prettyNum(cor_value[['estimate']])} (pearson)"))
pp(file = "figures/compare_cf_and_visit_in_model_eosinophil_logfc.svg")
logfc_plot
dev.off()

## png 
##   2

logfc_plot

cor_value <- cor.test(merged[["padj"]], merged[["deseq_adjp"]], method = "spearman")

## Warning in cor.test.default(merged[["padj"]], merged[["deseq_adjp"]], method =
## "spearman"): Cannot compute exact p-value with ties

cor_value

## 
##  Spearman's rank correlation rho
## 
## data:  merged[["padj"]] and merged[["deseq_adjp"]]
## S = 3.5e+10, p-value <2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##    rho 
## 0.8214

adjp_plotter <- plot_linear_scatter(merged[, c("padj", "deseq_adjp")])
adjp_plot <- adjp_plotter[["scatter"]] +
  xlab("DESeq2 adjp: Visit explicitly in model") +
  ylab("DESeq2 adjp: Default pairwise comparison") +
  ggtitle(glue("Comparing results from models: {prettyNum(cor_value[['estimate']])} (spearman)"))
pp(file = "images/compare_cf_and_visit_in_model_eosinophil_adjp.svg")
adjp_plot
dev.off()

## png 
##   2

adjp_plot

previous_sig_idx <- big_table[["deseq_adjp"]] <= 0.05 & abs(big_table[["deseq_logfc"]] >= 1.0)
summary(previous_sig_idx)

##    Mode   FALSE    TRUE 
## logical   10416     116

previous_genes <- rownames(big_table)[previous_sig_idx]

new_sig_idx <- abs(deseq_table[["log2FoldChange"]]) >= 1.0 & deseq_table[["padj"]] < 0.05
new_genes <- rownames(deseq_table)[new_sig_idx]
na_idx <- is.na(new_genes)
new_genes <- new_genes[!na_idx]

Vennerable::Venn(list("previous" = previous_genes, "new" = new_genes))

## A Venn object on 2 sets named
## previous,new 
##  00  10  01  11 
##   0  38 193  78

test_new <- simple_gprofiler(new_genes)
test_new

## A set of ontologies produced by gprofiler using 271
## genes against the hsapiens annotations and significance cutoff 0.05.
## There are: 
## 11 MF
## 186 BP
## 4 KEGG
## 6 REAC
## 5 WP
## 7 TF
## 0 MIRNA
## 1 HPA
## 0 CORUM
## 0 HP hits.

test_old <- simple_gprofiler(previous_genes)
test_old

## A set of ontologies produced by gprofiler using 116
## genes against the hsapiens annotations and significance cutoff 0.05.
## There are: 
## 26 MF
## 112 BP
## 3 KEGG
## 7 REAC
## 4 WP
## 72 TF
## 1 MIRNA
## 0 HPA
## 0 CORUM
## 0 HP hits.

new_annotated <- merge(fData(t_eosinophils), deseq_table, by = "row.names")
rownames(new_annotated) <- new_annotated[["Row.names"]]
new_annotated[["Row.names"]] <- NULL
write_xlsx(data = new_annotated, excel = "excel/eosinophil_visit_in_model_sva_cf_new.xlsx")

## write_xlsx() wrote excel/eosinophil_visit_in_model_sva_cf_new.xlsx.
## The cursor is on sheet first, row: 17303 column: 23.

old_annotated <- merge(fData(t_eosinophils), big_table, by = "row.names")
rownames(old_annotated) <- old_annotated[["Row.names"]]
old_annotated[["Row.names"]] <- NULL
write_xlsx(data = old_annotated, excel = "excel/eosinophil_visit_in_model_sva_cf_old.xlsx")

## write_xlsx() wrote excel/eosinophil_visit_in_model_sva_cf_old.xlsx.
## The cursor is on sheet first, row: 10535 column: 76.

8.0.3.1 C/F celltype volcano plots with specific labels

8.0.4 Figure 4B: Volcano plots

Note to self, when I rendered the html, stupid R ran out of temp files and so did not actually print the darn html document, as a result I modified the render function to try to make sure there is a clean directory in which to work; testing now. If it continues to not work, I will need to remove some of the images created in this document.

num_color <- color_choices[["clinic_cf"]][["tumaco_failure"]]
den_color <- color_choices[["clinic_cf"]][["tumaco_cure"]]
wanted_genes <- c("FI44L", "IFI27", "PRR5", "PRR5-ARHGAP8", "RHCE",
                  "FBXO39", "RSAD2", "SMTNL1", "USP18", "AFAP1")

cf_monocyte_table <- t_cf_monocyte_table_sva[["data"]][["outcome"]]
cf_monocyte_volcano <- plot_volcano_condition_de(
  cf_monocyte_table, "outcome", label = wanted_genes,
  fc_col = "deseq_logfc", p_col = "deseq_adjp", line_position = NULL,
  color_high = num_color, color_low = den_color, label_size = 6)
pp(file = "figures/cf_monocyte_volcano_labeled.svg")
cf_monocyte_volcano[["plot"]]
dev.off()

## png 
##   2

cf_monocyte_volcano[["plot"]]

cf_monocyte_volcano_top10 <- plot_volcano_condition_de(
  cf_monocyte_table, "outcome", label = 10,
  fc_col = "deseq_logfc", p_col = "deseq_adjp", line_position = NULL,
  color_high = num_color, color_low = den_color, label_size = 6)
pp(file = glue("images/cf_monocyte_volcano_labeled_top10-v{ver}.svg"))
cf_monocyte_volcano_top10[["plot"]]
dev.off()

## png 
##   2

cf_monocyte_volcano_top10[["plot"]]

cf_eosinophil_table <- t_cf_eosinophil_table_sva[["data"]][["outcome"]]
cf_eosinophil_volcano <- plot_volcano_condition_de(
  cf_eosinophil_table, "outcome", label = wanted_genes,
  fc_col = "deseq_logfc", p_col = "deseq_adjp", line_position = NULL,
  color_high = num_color, color_low = den_color, label_size = 6)
pp(file = "figures/cf_eosinophil_volcano_labeled.svg")
cf_eosinophil_volcano[["plot"]]
dev.off()

## png 
##   2

cf_eosinophil_volcano[["plot"]]

cf_eosinophil_volcano_top10 <- plot_volcano_condition_de(
  cf_eosinophil_table, "outcome", label = 10,
  fc_col = "deseq_logfc", p_col = "deseq_adjp", line_position = NULL,
  color_high = num_color, color_low = den_color, label_size = 6)
pp(file = glue("images/cf_eosinophil_volcano_labeled_top10-v{ver}.svg"))
cf_eosinophil_volcano_top10[["plot"]]
dev.off()

## png 
##   2

cf_eosinophil_volcano_top10[["plot"]]

cf_neutrophil_table <- t_cf_neutrophil_table_sva[["data"]][["outcome"]]
cf_neutrophil_volcano <- plot_volcano_condition_de(
  cf_neutrophil_table, "outcome", label = wanted_genes,
  fc_col = "deseq_logfc", p_col = "deseq_adjp", line_position = NULL,
  color_high = num_color, color_low = den_color, label_size = 6)
pp(file = "figures/cf_neutrophil_volcano_labeled.svg")
cf_neutrophil_volcano[["plot"]]
dev.off()

## png 
##   2

cf_neutrophil_volcano[["plot"]]

cf_neutrophil_volcano_top10 <- plot_volcano_condition_de(
  cf_neutrophil_table, "outcome", label = 10,
  fc_col = "deseq_logfc", p_col = "deseq_adjp", line_position = NULL,
  color_high = num_color, color_low = den_color, label_size = 6)
pp(file = glue("images/cf_neutrophil_volcano_labeled_top10-v{ver}.svg"))
cf_neutrophil_volcano_top10[["plot"]]
dev.off()

## png 
##   2

cf_neutrophil_volcano_top10[["plot"]]

8.0.5 Eosinophil time comparisons

t_cf_eosinophil_v1_de_sva <- all_pairwise(tv1_eosinophils, model_batch = "svaseq",
                                          parallel = parallel, filter = TRUE,
                                          methods = methods)

## 
##    tumaco_cure tumaco_failure 
##              5              3

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

## Setting 57 low elements to zero.

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

t_cf_eosinophil_v1_de_sva

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

## The logFC agreement among the methods follows:

##                 tmc_flr___
## limma_vs_deseq      0.8464
## limma_vs_edger      0.8498
## limma_vs_basic      0.9207
## limma_vs_noiseq    -0.8343
## deseq_vs_edger      0.9996
## deseq_vs_basic      0.8703
## deseq_vs_noiseq    -0.7805
## edger_vs_basic      0.8716
## edger_vs_noiseq    -0.7849
## basic_vs_noiseq    -0.8407

t_cf_eosinophil_v1_table_sva <- combine_de_tables(
  t_cf_eosinophil_v1_de_sva, keepers = t_cf_contrast,
  excel = glue("{cf_prefix}/Eosinophils/t_eosinophil_v1_cf_table_sva-v{ver}.xlsx"))
t_cf_eosinophil_v1_table_sva

## A set of combined differential expression results.

##                           table deseq_sigup deseq_sigdown edger_sigup
## 1 tumaco_failure_vs_tumaco_cure          13            19          11
##   edger_sigdown limma_sigup limma_sigdown
## 1            10           0             0

## `geom_line()`: Each group consists of only one observation.
## i Do you need to adjust the group aesthetic?

## Plot describing unique/shared genes in a differential expression table.

t_cf_eosinophil_v1_sig_sva <- extract_significant_genes(
  t_cf_eosinophil_v1_table_sva,
  excel = glue("{cf_prefix}/Eosinophils/t_eosinophil_v1_cf_sig_sva-v{ver}.xlsx"))
t_cf_eosinophil_v1_table_sva

## A set of combined differential expression results.

##                           table deseq_sigup deseq_sigdown edger_sigup
## 1 tumaco_failure_vs_tumaco_cure          13            19          11
##   edger_sigdown limma_sigup limma_sigdown
## 1            10           0             0

## `geom_line()`: Each group consists of only one observation.
## i Do you need to adjust the group aesthetic?
## Plot describing unique/shared genes in a differential expression table.

dim(t_cf_eosinophil_v1_sig_sva$deseq$ups[[1]])

## [1] 13 59

dim(t_cf_eosinophil_v1_sig_sva$deseq$downs[[1]])

## [1] 19 59

t_cf_eosinophil_v2_de_sva <- all_pairwise(tv2_eosinophils, model_batch = "svaseq",
                                          parallel = parallel, filter = TRUE,
                                          methods = methods)

## 
##    tumaco_cure tumaco_failure 
##              6              3

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

## Setting 90 low elements to zero.

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

t_cf_eosinophil_v2_de_sva

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

## The logFC agreement among the methods follows:

##                 tmc_flr___
## limma_vs_deseq      0.8540
## limma_vs_edger      0.8581
## limma_vs_basic      0.8936
## limma_vs_noiseq    -0.7843
## deseq_vs_edger      0.9996
## deseq_vs_basic      0.8446
## deseq_vs_noiseq    -0.8331
## edger_vs_basic      0.8465
## edger_vs_noiseq    -0.8359
## basic_vs_noiseq    -0.8616

t_cf_eosinophil_v2_table_sva <- combine_de_tables(
  t_cf_eosinophil_v2_de_sva, keepers = t_cf_contrast,
  excel = glue("{cf_prefix}/Eosinophils/t_eosinophil_v2_cf_table_sva-v{ver}.xlsx"))
t_cf_eosinophil_v2_table_sva

## A set of combined differential expression results.

##                           table deseq_sigup deseq_sigdown edger_sigup
## 1 tumaco_failure_vs_tumaco_cure           9             4          10
##   edger_sigdown limma_sigup limma_sigdown
## 1             1           0             0

## `geom_line()`: Each group consists of only one observation.
## i Do you need to adjust the group aesthetic?

## Plot describing unique/shared genes in a differential expression table.

t_cf_eosinophil_v2_sig_sva <- extract_significant_genes(
  t_cf_eosinophil_v2_table_sva,
  excel = glue("{cf_prefix}/Eosinophils/t_eosinophil_v2_cf_sig_sva-v{ver}.xlsx"))
t_cf_eosinophil_v2_sig_sva

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##         limma_up limma_down edger_up edger_down deseq_up deseq_down basic_up
## outcome        0          0       10          1        9          4        0
##         basic_down
## outcome          0

dim(t_cf_eosinophil_v2_sig_sva$deseq$ups[[1]])

## [1]  9 59

dim(t_cf_eosinophil_v2_sig_sva$deseq$downs[[1]])

## [1]  4 59

t_cf_eosinophil_v3_de_sva <- all_pairwise(tv3_eosinophils, model_batch = "svaseq",
                                          parallel = parallel, filter = TRUE,
                                          methods = methods)

## 
##    tumaco_cure tumaco_failure 
##              6              3

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

## Setting 48 low elements to zero.

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

t_cf_eosinophil_v3_de_sva

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

## The logFC agreement among the methods follows:

##                 tmc_flr___
## limma_vs_deseq      0.9137
## limma_vs_edger      0.9138
## limma_vs_basic      0.8620
## limma_vs_noiseq    -0.7838
## deseq_vs_edger      1.0000
## deseq_vs_basic      0.8019
## deseq_vs_noiseq    -0.8252
## edger_vs_basic      0.8021
## edger_vs_noiseq    -0.8258
## basic_vs_noiseq    -0.8696

t_cf_eosinophil_v3_table_sva <- combine_de_tables(
  t_cf_eosinophil_v3_de_sva, keepers = t_cf_contrast,
  excel = glue("{cf_prefix}/Eosinophils/t_eosinophil_v3_cf_table_sva-v{ver}.xlsx"))
t_cf_eosinophil_v3_table_sva

## A set of combined differential expression results.

##                           table deseq_sigup deseq_sigdown edger_sigup
## 1 tumaco_failure_vs_tumaco_cure          68            29          73
##   edger_sigdown limma_sigup limma_sigdown
## 1            10           0             0

## `geom_line()`: Each group consists of only one observation.
## i Do you need to adjust the group aesthetic?

## Plot describing unique/shared genes in a differential expression table.

t_cf_eosinophil_v3_sig_sva <- extract_significant_genes(
  t_cf_eosinophil_v3_table_sva,
  excel = glue("{cf_prefix}/Eosinophils/t_eosinophil_v3_cf_sig_sva-v{ver}.xlsx"))
t_cf_eosinophil_v3_sig_sva

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##         limma_up limma_down edger_up edger_down deseq_up deseq_down basic_up
## outcome        0          0       73         10       68         29        0
##         basic_down
## outcome          0

dim(t_cf_eosinophil_v3_sig_sva$deseq$ups[[1]])

## [1] 68 59

dim(t_cf_eosinophil_v3_sig_sva$deseq$downs[[1]])

## [1] 29 59

8.0.5.1 Eosinophils: Compare sva to batch-in-visit

sva_aucc <- calculate_aucc(t_cf_eosinophil_table_sva[["data"]][[1]],
                           tbl2 = t_cf_eosinophil_table_batchvisit[["data"]][[1]],
                           py = "deseq_adjp", ly = "deseq_logfc")
sva_aucc

## These two tables have an aucc value of: 0.576029928864987 and correlation:

## 
##  Pearson's product-moment correlation
## 
## data:  tbl[[lx]] and tbl[[ly]]
## t = 152, df = 10530, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.823 0.835
## sample estimates:
##    cor 
## 0.8291

shared_ids <- rownames(t_cf_eosinophil_table_sva[["data"]][[1]]) %in%
  rownames(t_cf_eosinophil_table_batchvisit[["data"]][[1]])
first <- t_cf_eosinophil_table_sva[["data"]][[1]][shared_ids, ]
second <- t_cf_eosinophil_table_batchvisit[["data"]][[1]][rownames(first), ]
cor.test(first[["deseq_logfc"]], second[["deseq_logfc"]])

## 
##  Pearson's product-moment correlation
## 
## data:  first[["deseq_logfc"]] and second[["deseq_logfc"]]
## t = 152, df = 10530, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.823 0.835
## sample estimates:
##    cor 
## 0.8291

8.0.5.2 Compare monocyte CF, neutrophil CF, eosinophil CF

t_mono_neut_sva_aucc <- calculate_aucc(t_cf_monocyte_table_sva[["data"]][["outcome"]],
                                       tbl2 = t_cf_neutrophil_table_sva[["data"]][["outcome"]],
                                       py = "deseq_adjp", ly = "deseq_logfc")
t_mono_neut_sva_aucc

## These two tables have an aucc value of: 0.204316386168083 and correlation:

## 
##  Pearson's product-moment correlation
## 
## data:  tbl[[lx]] and tbl[[ly]]
## t = 43, df = 8577, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4028 0.4376
## sample estimates:
##    cor 
## 0.4203

t_mono_eo_sva_aucc <- calculate_aucc(t_cf_monocyte_table_sva[["data"]][["outcome"]],
                                     tbl2 = t_cf_eosinophil_table_sva[["data"]][["outcome"]],
                                     py = "deseq_adjp", ly = "deseq_logfc")
t_mono_eo_sva_aucc

## These two tables have an aucc value of: 0.0963678364630121 and correlation:

## 
##  Pearson's product-moment correlation
## 
## data:  tbl[[lx]] and tbl[[ly]]
## t = 22, df = 9765, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2015 0.2393
## sample estimates:
##    cor 
## 0.2205

t_neut_eo_sva_aucc <- calculate_aucc(t_cf_neutrophil_table_sva[["data"]][["outcome"]],
                                     tbl2 = t_cf_eosinophil_table_sva[["data"]][["outcome"]],
                                     py = "deseq_adjp", ly = "deseq_logfc")
t_neut_eo_sva_aucc

## These two tables have an aucc value of: 0.20148477670576 and correlation:

## 
##  Pearson's product-moment correlation
## 
## data:  tbl[[lx]] and tbl[[ly]]
## t = 42, df = 8571, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3973 0.4323
## sample estimates:
##   cor 
## 0.415

8.0.6 By visit

For these contrasts, we want to see fail_v1 vs. cure_v1, fail_v2 vs. cure_v2 etc. As a result, we will need to juggle the data slightly and add another set of contrasts.

8.0.6.1 Cure/Fail by visits, all cell types

t_visit_cf_all_de_sva <- all_pairwise(t_visitcf, model_batch = "svaseq",
                                      parallel = parallel, filter = TRUE,
                                      methods = methods)

## 
##    v1_cure v1_failure    v2_cure v2_failure    v3_cure v3_failure 
##         30         24         20         15         17         17

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

## Setting 17160 low elements to zero.

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

t_visit_cf_all_de_sva

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

t_visit_cf_all_table_sva <- combine_de_tables(
  t_visit_cf_all_de_sva, keepers = visitcf_contrasts,
  excel = glue("{cf_prefix}/t_all_visitcf_table_sva-v{ver}.xlsx"))
t_visit_cf_all_table_sva

## A set of combined differential expression results.

##                   table deseq_sigup deseq_sigdown edger_sigup edger_sigdown
## 1 v1_failure_vs_v1_cure          26            76          26            58
## 2 v2_failure_vs_v2_cure          51            41          43            28
## 3 v3_failure_vs_v3_cure          77            32          33            25
##   limma_sigup limma_sigdown
## 1           9            17
## 2           3             0
## 3           3             0

## Plot describing unique/shared genes in a differential expression table.

t_visit_cf_all_sig_sva <- extract_significant_genes(
  t_visit_cf_all_table_sva,
  excel = glue("{cf_prefix}/t_all_visitcf_sig_sva-v{ver}.xlsx"))
t_visit_cf_all_sig_sva

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##      limma_up limma_down edger_up edger_down deseq_up deseq_down basic_up
## v1cf        9         17       26         58       26         76        0
## v2cf        3          0       43         28       51         41        0
## v3cf        3          0       33         25       77         32        0
##      basic_down
## v1cf          0
## v2cf          0
## v3cf          0

8.0.6.2 Cure/Fail by visit, Monocytes

visitcf_factor <- paste0("v", pData(t_monocytes)[["visitnumber"]], "_",
                         pData(t_monocytes)[["finaloutcome"]])
t_monocytes_visitcf <- set_expt_conditions(t_monocytes, fact = visitcf_factor)

## The numbers of samples by condition are:

## 
##    v1_cure v1_failure    v2_cure v2_failure    v3_cure v3_failure 
##          8          8          7          6          6          7

t_visit_cf_monocyte_de_sva <- all_pairwise(t_monocytes_visitcf, model_batch = "svaseq",
                                           parallel = parallel, filter = TRUE,
                                           methods = methods)

## 
##    v1_cure v1_failure    v2_cure v2_failure    v3_cure v3_failure 
##          8          8          7          6          6          7

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

## Setting 700 low elements to zero.

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

t_visit_cf_monocyte_de_sva

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

t_visit_cf_monocyte_table_sva <- combine_de_tables(
  t_visit_cf_monocyte_de_sva, keepers = visitcf_contrasts,
  excel = glue("{cf_prefix}/Monocytes/t_monocyte_visitcf_table_sva-v{ver}.xlsx"))
t_visit_cf_monocyte_table_sva

## A set of combined differential expression results.

##                   table deseq_sigup deseq_sigdown edger_sigup edger_sigdown
## 1 v1_failure_vs_v1_cure          15            10          10            13
## 2 v2_failure_vs_v2_cure           0             0           0             0
## 3 v3_failure_vs_v3_cure           0             0           0             0
##   limma_sigup limma_sigdown
## 1           1             1
## 2           0             0
## 3           0             0

## `geom_line()`: Each group consists of only one observation.
## i Do you need to adjust the group aesthetic?

## Plot describing unique/shared genes in a differential expression table.

t_visit_cf_monocyte_sig_sva <- extract_significant_genes(
  t_visit_cf_monocyte_table_sva,
  excel = glue("{cf_prefix}/Monocytes/t_monocyte_visitcf_sig_sva-v{ver}.xlsx"))
t_visit_cf_monocyte_sig_sva

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##      limma_up limma_down edger_up edger_down deseq_up deseq_down basic_up
## v1cf        1          1       10         13       15         10        0
## v2cf        0          0        0          0        0          0        0
## v3cf        0          0        0          0        0          0        0
##      basic_down
## v1cf          0
## v2cf          0
## v3cf          0

t_v1fc_deseq_ma <- t_visit_cf_monocyte_table_sva[["plots"]][["v1cf"]][["deseq_ma_plots"]][["plot"]]
dev <- pp(file = "images/monocyte_cf_de_v1_maplot.png")
t_v1fc_deseq_ma

## NULL

closed <- dev.off()
t_v1fc_deseq_ma

## NULL

t_v2fc_deseq_ma <- t_visit_cf_monocyte_table_sva[["plots"]][["v2cf"]][["deseq_ma_plots"]][["plot"]]
dev <- pp(file = "images/monocyte_cf_de_v2_maplot.png")
t_v2fc_deseq_ma

## NULL

closed <- dev.off()
t_v2fc_deseq_ma

## NULL

t_v3fc_deseq_ma <- t_visit_cf_monocyte_table_sva[["plots"]][["v3cf"]][["deseq_ma_plots"]][["plot"]]
dev <- pp(file = "images/monocyte_cf_de_v3_maplot.png")
t_v3fc_deseq_ma

## NULL

closed <- dev.off()
t_v3fc_deseq_ma

## NULL

One query from Alejandro is to look at the genes shared up/down across visits. I am not entirely certain we have enough samples for this to work, but let us find out.

I am thinking this is a good place to use the AUCC curves I learned about thanks to Julie Cridland.

Note that the following is all monocyte samples, this should therefore potentially be moved up and a version of this with only the Tumaco samples put here?

v1cf <- t_visit_cf_monocyte_table_sva[["data"]][["v1cf"]]
v2cf <- t_visit_cf_monocyte_table_sva[["data"]][["v2cf"]]
v3cf <- t_visit_cf_monocyte_table_sva[["data"]][["v3cf"]]

v1_sig <- c(
  rownames(t_visit_cf_monocyte_sig_sva[["deseq"]][["ups"]][["v1cf"]]),
  rownames(t_visit_cf_monocyte_sig_sva[["deseq"]][["downs"]][["v1cf"]]))
length(v1_sig)

## [1] 25

v2_sig <- c(
  rownames(t_visit_cf_monocyte_sig_sva[["deseq"]][["ups"]][["v2cf"]]),
  rownames(t_visit_cf_monocyte_sig_sva[["deseq"]][["downs"]][["v2cf"]]))
length(v2_sig)

## [1] 0

v3_sig <- c(
  rownames(t_visit_cf_monocyte_sig_sva[["deseq"]][["ups"]][["v2cf"]]),
  rownames(t_visit_cf_monocyte_sig_sva[["deseq"]][["downs"]][["v2cf"]]))
length(v3_sig)

## [1] 0

t_monocyte_visit_aucc_v2v1 <- calculate_aucc(v1cf, tbl2 = v2cf,
                                             py = "deseq_adjp", ly = "deseq_logfc")
dev <- pp(file = "images/monocyte_visit_v2v1_aucc.png")
t_monocyte_visit_aucc_v2v1[["plot"]]
closed <- dev.off()
t_monocyte_visit_aucc_v2v1[["plot"]]

t_monocyte_visit_aucc_v3v1 <- calculate_aucc(v1cf, tbl2 = v3cf,
                                             py = "deseq_adjp", ly = "deseq_logfc")
dev <- pp(file = "images/monocyte_visit_v3v1_aucc.png")
t_monocyte_visit_aucc_v3v1[["plot"]]
closed <- dev.off()
t_monocyte_visit_aucc_v3v1[["plot"]]

8.0.6.3 Cure/Fail by visit, Neutrophils

visitcf_factor <- paste0("v", pData(t_neutrophils)[["visitnumber"]], "_",
                         pData(t_neutrophils)[["finaloutcome"]])
t_neutrophil_visitcf <- set_expt_conditions(t_neutrophils, fact = visitcf_factor)

## The numbers of samples by condition are:

## 
##    v1_cure v1_failure    v2_cure v2_failure    v3_cure v3_failure 
##          8          8          7          6          5          7

t_visit_cf_neutrophil_de_sva <- all_pairwise(t_neutrophil_visitcf, model_batch = "svaseq",
                                             parallel = parallel, filter = TRUE,
                                             methods = methods)

## 
##    v1_cure v1_failure    v2_cure v2_failure    v3_cure v3_failure 
##          8          8          7          6          5          7

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

## Setting 685 low elements to zero.

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

t_visit_cf_neutrophil_de_sva

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

t_visit_cf_neutrophil_table_sva <- combine_de_tables(
  t_visit_cf_neutrophil_de_sva, keepers = visitcf_contrasts,
  excel = glue("{cf_prefix}/Neutrophils/t_neutrophil_visitcf_table_sva-v{ver}.xlsx"))
t_visit_cf_neutrophil_table_sva

## A set of combined differential expression results.

##                   table deseq_sigup deseq_sigdown edger_sigup edger_sigdown
## 1 v1_failure_vs_v1_cure          12             6           6             6
## 2 v2_failure_vs_v2_cure           2             6           2             3
## 3 v3_failure_vs_v3_cure           2             2           0             2
##   limma_sigup limma_sigdown
## 1           1             0
## 2           0             0
## 3           0             0

## Plot describing unique/shared genes in a differential expression table.

t_visit_cf_neutrophil_sig_sva <- extract_significant_genes(
  t_visit_cf_neutrophil_table_sva,
  excel = glue("{cf_prefix}/Neutrophils/t_neutrophil_visitcf_sig_sva-v{ver}.xlsx"))
t_visit_cf_neutrophil_sig_sva

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##      limma_up limma_down edger_up edger_down deseq_up deseq_down basic_up
## v1cf        1          0        6          6       12          6        0
## v2cf        0          0        2          3        2          6        0
## v3cf        0          0        0          2        2          2        0
##      basic_down
## v1cf          0
## v2cf          0
## v3cf          0

8.0.6.4 Cure/Fail by visit, Eosinophils

visitcf_factor <- paste0("v", pData(t_eosinophils)[["visitnumber"]], "_",
                         pData(t_eosinophils)[["finaloutcome"]])
t_eosinophil_visitcf <- set_expt_conditions(t_eosinophils, fact = visitcf_factor)

## The numbers of samples by condition are:

## 
##    v1_cure v1_failure    v2_cure v2_failure    v3_cure v3_failure 
##          5          3          6          3          6          3

t_visit_cf_eosinophil_de_sva <- all_pairwise(t_eosinophil_visitcf, model_batch = "svaseq",
                                             parallel = parallel, filter = TRUE,
                                             methods = methods)

## 
##    v1_cure v1_failure    v2_cure v2_failure    v3_cure v3_failure 
##          5          3          6          3          6          3

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

## Setting 373 low elements to zero.

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

t_visit_cf_eosinophil_de_sva

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

t_visit_cf_eosinophil_table_sva <- combine_de_tables(
  t_visit_cf_eosinophil_de_sva, keepers = visitcf_contrasts,
  excel = glue("{cf_prefix}/Eosinophils/t_eosinophil_visitcf_table_sva-v{ver}.xlsx"))
t_visit_cf_eosinophil_table_sva

## A set of combined differential expression results.

##                   table deseq_sigup deseq_sigdown edger_sigup edger_sigdown
## 1 v1_failure_vs_v1_cure           9            11           2             3
## 2 v2_failure_vs_v2_cure           4             3           5             2
## 3 v3_failure_vs_v3_cure          14             7          17             2
##   limma_sigup limma_sigdown
## 1           0             1
## 2           0             0
## 3           0             0

## Plot describing unique/shared genes in a differential expression table.

t_visit_cf_eosinophil_sig_sva <- extract_significant_genes(
  t_visit_cf_eosinophil_table_sva,
  excel = glue("{cf_prefix}/Eosinophils/t_eosinophil_visitcf_sig_sva-v{ver}.xlsx"))
t_visit_cf_eosinophil_sig_sva

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##      limma_up limma_down edger_up edger_down deseq_up deseq_down basic_up
## v1cf        0          1        2          3        9         11        0
## v2cf        0          0        5          2        4          3        0
## v3cf        0          0       17          2       14          7        0
##      basic_down
## v1cf          0
## v2cf          0
## v3cf          0

8.1 Persistence in visit 3

Having put some SL read mapping information in the sample sheet, Maria Adelaida added a new column using it with the putative persistence state on a per-sample basis. One question which arised from that: what differences are observable between the persistent yes vs. no samples on a per-cell-type basis among the visit 3 samples.

8.1.1 Setting up

First things first, create the datasets.

persistence_expt <- subset_expt(t_clinical, subset = "persistence=='Y'|persistence=='N'") %>%
  subset_expt(subset = 'visitnumber==3') %>%
  set_expt_conditions(fact = 'persistence')

## subset_expt(): There were 123, now there are 83 samples.

## subset_expt(): There were 83, now there are 30 samples.

## The numbers of samples by condition are:

## 
##  N  Y 
##  6 24

## persistence_biopsy <- subset_expt(persistence_expt, subset = "typeofcells=='biopsy'")
persistence_monocyte <- subset_expt(persistence_expt, subset = "typeofcells=='monocytes'")

## subset_expt(): There were 30, now there are 12 samples.

persistence_neutrophil <- subset_expt(persistence_expt, subset = "typeofcells=='neutrophils'")

## subset_expt(): There were 30, now there are 10 samples.

persistence_eosinophil <- subset_expt(persistence_expt, subset = "typeofcells=='eosinophils'")

## subset_expt(): There were 30, now there are 8 samples.

8.1.2 Take a look

See if there are any patterns which look usable.

## All
persistence_norm <- normalize_expt(persistence_expt, transform = "log2", convert = "cpm",
                                   norm = "quant", filter = TRUE)

## Removing 8563 low-count genes (11389 remaining).

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

plot_pca(persistence_norm)[["plot"]]

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

## Removing 8563 low-count genes (11389 remaining).

## Setting 1544 low elements to zero.

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

plot_pca(persistence_nb)[["plot"]]

## Biopsies
##persistence_biopsy_norm <- normalize_expt(persistence_biopsy, transform = "log2", convert = "cpm",
##                                   norm = "quant", filter = TRUE)
##plot_pca(persistence_biopsy_norm)[["plot"]]
## Insufficient data

## Monocytes
persistence_monocyte_norm <- normalize_expt(persistence_monocyte, transform = "log2", convert = "cpm",
                                            norm = "quant", filter = TRUE)

## Removing 9623 low-count genes (10329 remaining).

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

plot_pca(persistence_monocyte_norm)[["plot"]]

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

## Removing 9623 low-count genes (10329 remaining).

## Setting 47 low elements to zero.

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

plot_pca(persistence_monocyte_nb)[["plot"]]

## Neutrophils
persistence_neutrophil_norm <- normalize_expt(persistence_neutrophil, transform = "log2", convert = "cpm",
                                              norm = "quant", filter = TRUE)

## Removing 11558 low-count genes (8394 remaining).

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

plot_pca(persistence_neutrophil_norm)[["plot"]]

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

## Removing 11558 low-count genes (8394 remaining).

## Setting 46 low elements to zero.

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

plot_pca(persistence_neutrophil_nb)[["plot"]]

## Eosinophils
persistence_eosinophil_norm <- normalize_expt(persistence_eosinophil, transform = "log2", convert = "cpm",
                                              norm = "quant", filter = TRUE)

## Removing 9922 low-count genes (10030 remaining).

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

plot_pca(persistence_eosinophil_norm)[["plot"]]

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

## Removing 9922 low-count genes (10030 remaining).

## Setting 25 low elements to zero.

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

plot_pca(persistence_eosinophil_nb)[["plot"]]

8.1.3 persistence DE

persistence_de_sva <- all_pairwise(persistence_expt, filter = TRUE, methods = methods,
                                   parallel = parallel, model_batch = "svaseq")

## 
##  N  Y 
##  6 24

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

## Setting 1544 low elements to zero.

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

persistence_de_sva

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

## The logFC agreement among the methods follows:

##                  Y_vs_N
## limma_vs_deseq   0.8112
## limma_vs_edger   0.8765
## limma_vs_basic   0.8217
## limma_vs_noiseq -0.6444
## deseq_vs_edger   0.9605
## deseq_vs_basic   0.7178
## deseq_vs_noiseq -0.5868
## edger_vs_basic   0.7791
## edger_vs_noiseq -0.6432
## basic_vs_noiseq -0.7960

persistence_table_sva <- combine_de_tables(
  persistence_de_sva,
  excel = glue("{xlsx_prefix}/DE_Persistence/persistence_all_de_sva-v{ver}.xlsx"))
persistence_table_sva

## A set of combined differential expression results.

##    table deseq_sigup deseq_sigdown edger_sigup edger_sigdown limma_sigup
## 1 Y_vs_N          55            44          26            49           7
##   limma_sigdown
## 1            22

## `geom_line()`: Each group consists of only one observation.
## i Do you need to adjust the group aesthetic?

## Plot describing unique/shared genes in a differential expression table.

persistence_monocyte_de_sva <- all_pairwise(persistence_monocyte, filter = TRUE,
                                            parallel = parallel, model_batch = "svaseq",
                                            methods = methods)

## 
##  N  Y 
##  2 10

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

## Setting 47 low elements to zero.

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

persistence_monocyte_de_sva

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

## The logFC agreement among the methods follows:

##                  Y_vs_N
## limma_vs_deseq   0.9260
## limma_vs_edger   0.9270
## limma_vs_basic   0.9858
## limma_vs_noiseq -0.9013
## deseq_vs_edger   1.0000
## deseq_vs_basic   0.9237
## deseq_vs_noiseq -0.9181
## edger_vs_basic   0.9245
## edger_vs_noiseq -0.9193
## basic_vs_noiseq -0.9039

persistence_monocyte_table_sva <- combine_de_tables(
  persistence_monocyte_de_sva,
  excel = glue("{xlsx_prefix}/DE_Persistence/persistence_monocyte_de_sva-v{ver}.xlsx"))
persistence_monocyte_table_sva

## A set of combined differential expression results.

##    table deseq_sigup deseq_sigdown edger_sigup edger_sigdown limma_sigup
## 1 Y_vs_N           1             0           0             1           0
##   limma_sigdown
## 1             0

## Only Y_vs_N_up has information, cannot create an UpSet.

## Plot describing unique/shared genes in a differential expression table.

## NULL

persistence_neutrophil_de_sva <- all_pairwise(persistence_neutrophil, filter = TRUE,
                                              parallel = parallel, model_batch = "svaseq",
                                              methods = methods)

## 
## N Y 
## 3 7

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

## Setting 46 low elements to zero.

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

persistence_neutrophil_de_sva

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

## The logFC agreement among the methods follows:

##                  Y_vs_N
## limma_vs_deseq   0.9407
## limma_vs_edger   0.9408
## limma_vs_basic   0.8808
## limma_vs_noiseq -0.7817
## deseq_vs_edger   0.9985
## deseq_vs_basic   0.8270
## deseq_vs_noiseq -0.7593
## edger_vs_basic   0.8283
## edger_vs_noiseq -0.7602
## basic_vs_noiseq -0.9003

persistence_neutrophil_table_sva <- combine_de_tables(
  persistence_neutrophil_de_sva,
  excel = glue("{xlsx_prefix}/DE_Persistence/persistence_neutrophil_de_sva-v{ver}.xlsx"))
persistence_neutrophil_table_sva

## A set of combined differential expression results.

##    table deseq_sigup deseq_sigdown edger_sigup edger_sigdown limma_sigup
## 1 Y_vs_N          26            49          17            35           0
##   limma_sigdown
## 1             0

## `geom_line()`: Each group consists of only one observation.
## i Do you need to adjust the group aesthetic?

## Plot describing unique/shared genes in a differential expression table.

## There are insufficient samples (1) in the 'N' category.
##persistence_eosinophil_de_sva <- all_pairwise(persistence_eosinophil, filter = TRUE,
##                                              parallel = parallel, model_batch = "svaseq",
##                                              methods = methods)
##persistence_eosinophil_de_sva
##persistence_eosinophil_table_sva <- combine_de_tables(
##  persistence_eosinophil_de_sva,
##  excel = glue("{xlsx_prefix}/DE_Persistence/persistence_eosinophil_de_sva-v{ver}.xlsx"))

8.2 Comparing visits without regard to cure/fail

8.2.1 All cell types

t_visit_all_de_sva <- all_pairwise(t_visit, filter = TRUE, methods = methods,
                                   parallel = parallel, model_batch = "svaseq")

## 
##  3  2  1 
## 34 35 40

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

## Setting 9636 low elements to zero.

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

t_visit_all_de_sva

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

t_visit_all_table_sva <- combine_de_tables(
  t_visit_all_de_sva, keepers = visit_contrasts,
  excel = glue("{xlsx_prefix}/DE_Visits/t_all_visit_table_sva-v{ver}.xlsx"))
t_visit_all_table_sva

## A set of combined differential expression results.

##      table deseq_sigup deseq_sigdown edger_sigup edger_sigdown limma_sigup
## 1 c2_vs_c1          25             9          20            10          19
## 2 c3_vs_c1          20            20          18            16          21
## 3 c3_vs_c2           0             2           0             2           0
##   limma_sigdown
## 1             7
## 2             7
## 3             0

## Plot describing unique/shared genes in a differential expression table.

t_visit_all_sig_sva <- extract_significant_genes(
  t_visit_all_table_sva,
  excel = glue("{xlsx_prefix}/DE_Visits/t_all_visit_sig_sva-v{ver}.xlsx"))
t_visit_all_sig_sva

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##      limma_up limma_down edger_up edger_down deseq_up deseq_down basic_up
## v2v1       19          7       20         10       25          9        0
## v3v1       21          7       18         16       20         20        0
## v3v2        0          0        0          2        0          2        0
##      basic_down
## v2v1          0
## v3v1          0
## v3v2          0

8.2.2 Monocyte samples

t_visit_monocytes <- set_expt_conditions(t_monocytes, fact = "visitnumber")

## The numbers of samples by condition are:

## 
##  3  2  1 
## 13 13 16

t_visit_monocyte_de_sva <- all_pairwise(t_visit_monocytes, filter = TRUE,
                                        parallel = parallel, model_batch = "svaseq",
                                        methods = methods)

## 
##  3  2  1 
## 13 13 16

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

## Setting 658 low elements to zero.

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

t_visit_monocyte_de_sva

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

t_visit_monocyte_table_sva <- combine_de_tables(
  t_visit_monocyte_de_sva, keepers = visit_contrasts,
  excel = glue("{xlsx_prefix}/DE_Visits/Monocytes/t_monocyte_visit_table_sva-v{ver}.xlsx"))
t_visit_monocyte_table_sva

## A set of combined differential expression results.

##      table deseq_sigup deseq_sigdown edger_sigup edger_sigdown limma_sigup
## 1 c2_vs_c1           1             2           1             1           0
## 2 c3_vs_c1           2             1           1             1           0
## 3 c3_vs_c2           0             0           0             0           0
##   limma_sigdown
## 1             0
## 2             0
## 3             0

## `geom_line()`: Each group consists of only one observation.
## i Do you need to adjust the group aesthetic?

## Plot describing unique/shared genes in a differential expression table.

t_visit_monocyte_sig_sva <- extract_significant_genes(
  t_visit_monocyte_table_sva,
  excel = glue("{xlsx_prefix}/DE_Visits/Monocytes/t_monocyte_visit_sig_sva-v{ver}.xlsx"))
t_visit_monocyte_sig_sva

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##      limma_up limma_down edger_up edger_down deseq_up deseq_down basic_up
## v2v1        0          0        1          1        1          2        0
## v3v1        0          0        1          1        2          1        0
## v3v2        0          0        0          0        0          0        0
##      basic_down
## v2v1          0
## v3v1          0
## v3v2          0

8.2.3 Neutrophil samples

t_visit_neutrophils <- set_expt_conditions(t_neutrophils, fact = "visitnumber")

## The numbers of samples by condition are:

## 
##  3  2  1 
## 12 13 16

t_visit_neutrophil_de_sva <- all_pairwise(t_visit_neutrophils, filter = TRUE,
                                          parallel = parallel, model_batch = "svaseq",
                                          methods = methods)

## 
##  3  2  1 
## 12 13 16

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

## Setting 593 low elements to zero.

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

t_visit_neutrophil_de_sva

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

t_visit_neutrophil_table_sva <- combine_de_tables(
  t_visit_neutrophil_de_sva, keepers = visit_contrasts,
  excel = glue("{xlsx_prefix}/DE_Visits/Neutrophils/t_neutrophil_visit_table_sva-v{ver}.xlsx"))
t_visit_neutrophil_table_sva

## A set of combined differential expression results.

##      table deseq_sigup deseq_sigdown edger_sigup edger_sigdown limma_sigup
## 1 c2_vs_c1         111            88         111            88         116
## 2 c3_vs_c1         127            45         122            44          93
## 3 c3_vs_c2           1             0           0             0           0
##   limma_sigdown
## 1            52
## 2            67
## 3             0

## Plot describing unique/shared genes in a differential expression table.

t_visit_neutrophil_sig_sva <- extract_significant_genes(
  t_visit_neutrophil_table_sva,
  excel = glue("{xlsx_prefix}/DE_Visits/Neutrophils/t_neutrophil_visit_sig_sva-v{ver}.xlsx"))
t_visit_neutrophil_sig_sva

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##      limma_up limma_down edger_up edger_down deseq_up deseq_down basic_up
## v2v1      116         52      111         88      111         88       83
## v3v1       93         67      122         44      127         45       52
## v3v2        0          0        0          0        1          0        0
##      basic_down
## v2v1         35
## v3v1         17
## v3v2          0

8.2.4 Eosinophil samples

t_visit_eosinophils <- set_expt_conditions(t_eosinophils, fact="visitnumber")

## The numbers of samples by condition are:

## 
## 3 2 1 
## 9 9 8

t_visit_eosinophil_de <- all_pairwise(t_visit_eosinophils, filter = TRUE,
                                      parallel = parallel, model_batch = "svaseq",
                                      methods = methods)

## 
## 3 2 1 
## 9 9 8

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

## Setting 271 low elements to zero.

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

t_visit_eosinophil_de

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

t_visit_eosinophil_table <- combine_de_tables(
  t_visit_eosinophil_de, keepers = visit_contrasts,
  excel = glue("{xlsx_prefix}/DE_Visits/Eosinophils/t_eosinophil_visit_table_sva-v{ver}.xlsx"))
t_visit_eosinophil_table

## A set of combined differential expression results.

##      table deseq_sigup deseq_sigdown edger_sigup edger_sigdown limma_sigup
## 1 c2_vs_c1           0             0           0             0           0
## 2 c3_vs_c1           0             0           0             0           0
## 3 c3_vs_c2           0             0           0             0           0
##   limma_sigdown
## 1             0
## 2             0
## 3             0

## Only  has information, cannot create an UpSet.

## Plot describing unique/shared genes in a differential expression table.

## NULL

t_visit_eosinophil_sig <- extract_significant_genes(
  t_visit_eosinophil_table,
  excel = glue("{xlsx_prefix}/DE_Visits/Eosinophils/t_eosinophil_visit_sig_sva-v{ver}.xlsx"))
## No significant genes observed.

9 Explore ROC

Alejandro showed some ROC curves for eosinophil data showing sensitivity vs. specificity of a couple genes which were observed in v1 eosinophils vs. all-times eosinophils across cure/fail. I am curious to better understand how this was done and what utility it might have in other contexts.

To that end, I want to try something similar myself. In order to properly perform the analysis with these various tools, I need to reconfigure the data in a pretty specific format:

1. Single df with 1 row per set of observations (sample in this case I think)
2. The outcome column(s) need to be 1 (or more?) metadata factor(s) (cure/fail or a paste0 of relevant queries (eo_v1_cure, eo_v123_cure, etc)
3. The predictor column(s) are the measurements (rpkm of 1 or more genes), 1 column each gene.

If I intend to use this for our tx data, I will likely need a utility function to create the properly formatted input df.

For the purposes of my playing, I will choose three genes from the eosinophil C/F table, one which is significant, one which is not, and an arbitrary.

The input genes will therefore be chosen from the data structure: t_cf_eosinophil_table_sva:

ENSG00000198178, ENSG00000179344, ENSG00000182628

eo_rpkm <- normalize_expt(tv1_eosinophils, convert = "rpkm", column = "cds_length")

## There appear to be 5355 genes without a length.

10 An external dataset

This paper is DOI:10.1126/scitranslmed.aax4204

Variable gene expression and parasite load predict treatment outcome in cutaneous leishmaniasis.

One query from Maria Adelaida is to see how this data fits with ours. I have read this paper a couple of times now and I get confused on a couple of points every time, which I will explain in a moment. The expermental design is key to my confusion and key to what I think is being missed in our interpretation of the results:

1. The PCA is not cure vs. fail but healthy skin vs. CL lesion.

10.1 Only the Scott data

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

## Removing 7327 low-count genes (14154 remaining).

plot_pca(external_norm)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by cure, failure
## Shapes are defined by female, male.

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

## Removing 7327 low-count genes (14154 remaining).

## Setting 171 low elements to zero.

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

plot_pca(external_nb)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by cure, failure
## Shapes are defined by female, male.

external_de <- all_pairwise(external_cf, filter = TRUE, methods = methods,
                            parallel = parallel, model_batch = "svaseq")

## 
##    cure failure 
##      14       7

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

## Setting 171 low elements to zero.

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

external_de

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

## The logFC agreement among the methods follows:

##                 falr_vs_cr
## limma_vs_deseq      0.8487
## limma_vs_edger      0.8497
## limma_vs_basic      0.4180
## limma_vs_noiseq    -0.3286
## deseq_vs_edger      0.9997
## deseq_vs_basic      0.3908
## deseq_vs_noiseq    -0.3756
## edger_vs_basic      0.3914
## edger_vs_noiseq    -0.3768
## basic_vs_noiseq    -0.9338

external_table <- combine_de_tables(external_de, excel = "excel/scott_table.xlsx")
external_table

## A set of combined differential expression results.

##             table deseq_sigup deseq_sigdown edger_sigup edger_sigdown
## 1 failure_vs_cure           0             0           0             0
##   limma_sigup limma_sigdown
## 1           0             0

## Only  has information, cannot create an UpSet.

## Plot describing unique/shared genes in a differential expression table.

## NULL

external_sig <- extract_significant_genes(external_table, excel = "excel/scott_sig.xlsx")
external_sig

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##                 limma_up limma_down edger_up edger_down deseq_up deseq_down
## failure_vs_cure        0          0        0          0        0          0
##                 basic_up basic_down
## failure_vs_cure        0          0

external_top100 <- extract_significant_genes(external_table, n = 100)
external_up <- external_top100[["deseq"]][["ups"]][["failure_vs_cure"]]
external_down <- external_top100[["deseq"]][["downs"]][["failure_vs_cure"]]

11 An explicit comparison of methods.

I think I am getting a significantly different result from Scott, so I am going to do an explicit side-by-side comparison of our results at each step. In order to do this, I am using the capsule they kindly provided with their publication.

I am copy/pasting material from their publication with some modification which I will note as I go.

Here is their block ‘r packages’

Note/Spoiler alert: It actually turns out our results are basically relatively similar, I just didn’t understand what comparisons are actually in paper vs those I have primary interest. In addition, we handled gene IDs differently (gene card vs. EnsemblID) which has a surprisingly big effect.

Oh, I just realized that when I did these analyses, I did them in a completely separate tree and compared the results post-facto. This assumption remains in this document and therefore is unlikely to work properly in the containerized environment I am attempting to create. Given that the primary goal of this section is to show to myself that I compared the two datasets as thoroughly as I could, perhaps I should just disable them for the container and allow the reader to perform the exercise de-novo.

library(tidyverse)
library(ggthemes)
library(reshape2)
library(edgeR)
library(patchwork)
library(vegan)
library(DT)
library(tximport)
library(gplots)
library(FinCal)
library(ggrepel)
library(gt)
library(ggExtra)
library(EnsDb.Hsapiens.v86)
library(stringr)
library(cowplot)
library(ggpubr)

I have a separate tree in which I copied the capsule and data. I performed exactly their steps kallisto quant steps within it and put the output data into the same place within it. I did change the commands slightly because I downloaded the files from SRA and so don’t have them with names like ‘host_CL01’, but instead ‘PRJNA…’. But the samples are in the same order, so I sent the output files to the same final filenames. Here is an example from the first sample:

cd preprocessing
module add kallisto
kallisto index -i Homo_sapiens.GRCh38.cdna.all.Index Homo_sapiens.GRCh38.cdna.all.fa
# Map reads to the indexed reference transcriptome for HOST
# first the healthy subjects (HS)
export LESS = '--buffers 0 -B'
kallisto quant -i Homo_sapiens.GRCh38.cdna.all.Index -o host_HS01 -t 24 -b 60 \
         --single -l 250 -s 30 <(less SRR8668755/*-trimmed.fastq.xz) 2>host_HS01.log 1>&2 &

12 Block ‘sample_info’

I am going to change the path very slightly in the following block simply because I put the capsule in a separate directory and do not want to copy it here. Otherwise it is unmodified. Also, the function gt::tab_header() annoys the crap out of me.

import <- read_tsv("../scott_2019/capsule-6534016/data/studydesign.txt")
import %>% dplyr::filter(disease == "cutaneous") %>%
  dplyr::select(-2) %>%  gt() %>%
  tab_header(title = md("Clinical metadata from patients with cutaneous leishmaniasis (CL)"),
             subtitle = md("`(n=21)`")) %>%  cols_align(align = "center", columns = TRUE)
targets.lesion <- import
targets.onlypatients <- targets.lesion[8:28,] # only CL lesions (n=21)

# Making factors that will be used for pairwise comparisons:
# HS vs. CL lesions as a factor:
disease.lesion <- factor(targets.lesion$disease)
# Cure vs. Failure lesions as a factor:
treatment.lesion <- factor(targets.onlypatients$treatment_outcome)

13 Importing the data and annotations

They did use a slightly different annotation set, Ensembl revision 86. Once again I am modifying the paths slightly to reflect where I put the capsule.

# capturing Ensembl transcript IDs (tx) and gene symbols ("gene_name") from
# EnsDb.Hsapiens.v86 annotation package
Tx <- as.data.frame(transcripts(EnsDb.Hsapiens.v86,
                                columns=c(listColumns(EnsDb.Hsapiens.v86, "tx"),
                                          "gene_name")))

Tx <- dplyr::rename(Tx, target_id = tx_id)
row.names(Tx) <- NULL
Tx <- Tx[,c(6,12)]

# getting file paths for Kallisto outputs
paths.all <- file.path("../scott_2019/capsule-6534016/data/readMapping/human", targets.lesion$sample, "abundance.h5")
paths.patients <- file.path("../scott_2019/capsule-6534016/data/readMapping/human", targets.onlypatients$sample, "abundance.h5")

# importing .h5 Kallisto data and collapsing transcript-level data to genes
Txi.lesion.coding <- tximport(paths.all,
                              type = "kallisto",
                              tx2gene = Tx,
                              txOut = FALSE,
                              ignoreTxVersion = TRUE,
                              countsFromAbundance = "lengthScaledTPM")

# importing againg, but this time just the CL patients
Txi.lesion.coding.onlypatients <- tximport(paths.patients,
                                           type = "kallisto",
                                           tx2gene = Tx,
                                           txOut = FALSE,
                                           ignoreTxVersion = TRUE,
                                           countsFromAbundance = "lengthScaledTPM")

14 Filtering and normalization

The block ‘visualizationDatasets’ follows unchanged. In the next block I will add another plot or perhaps 2

# First make a DGEList from the counts:
Txi.lesion.coding.DGEList <- DGEList(Txi.lesion.coding$counts)
colnames(Txi.lesion.coding.DGEList$counts) <- targets.lesion$sample
colnames(Txi.lesion.coding$counts) <- targets.lesion$sample

Txi.lesion.coding.DGEList.OP <- DGEList(Txi.lesion.coding.onlypatients$counts)
colnames(Txi.lesion.coding.DGEList.OP) <- targets.onlypatients$sample

# Convert to counts per million:
Txi.lesion.coding.DGEList.cpm <- edgeR::cpm(Txi.lesion.coding.DGEList, log = TRUE)
Txi.lesion.coding.DGEList.OP.cpm <- edgeR::cpm(Txi.lesion.coding.DGEList.OP, log = TRUE)

keepers.coding <- rowSums(Txi.lesion.coding.DGEList.cpm>1)>=7
keepers.coding.OP <- rowSums(Txi.lesion.coding.DGEList.OP.cpm>1)>=7

Txi.lesion.coding.DGEList.filtered <- Txi.lesion.coding.DGEList[keepers.coding,]
Txi.lesion.coding.DGEList.OP.filtered <- Txi.lesion.coding.DGEList.OP[keepers.coding.OP,]

# convert back to cpm:
Txi.lesion.coding.DGEList.LogCPM.filtered <- edgeR::cpm(Txi.lesion.coding.DGEList.filtered,
                                                        log=TRUE)
Txi.lesion.coding.DGEList.LogCPM.OP.filtered <- edgeR::cpm(Txi.lesion.coding.DGEList.OP.filtered,
                                                           log=TRUE)

# Normalizing data:
calcNorm1 <- calcNormFactors(Txi.lesion.coding.DGEList.filtered, method = "TMM")
calcNorm2 <- calcNormFactors(Txi.lesion.coding.DGEList.OP.filtered, method = "TMM")

Txi.lesion.coding.DGEList.LogCPM.filtered.norm <- edgeR::cpm(calcNorm1, log=TRUE)
colnames(Txi.lesion.coding.DGEList.LogCPM.filtered.norm) <- targets.lesion$sample
Txi.lesion.coding.DGEList.OP.LogCPM.filtered.norm <- edgeR::cpm(calcNorm2, log=TRUE)
colnames(Txi.lesion.coding.DGEList.OP.LogCPM.filtered.norm) <- targets.onlypatients$sample
# Raw dataset:
V1 <- as.data.frame(Txi.lesion.coding.DGEList.cpm)
colnames(V1) <- targets.lesion$sample
V1 <- melt(V1)
colnames(V1) <- c("sample","expression")

# Filtered dataset:
V1.1 <- as.data.frame(Txi.lesion.coding.DGEList.LogCPM.filtered)
colnames(V1.1) <- targets.lesion$sample
V1.1 <- melt(V1.1)
colnames(V1.1) <- c("sample","expression")

# Filtered-normalized dataset:
V1.1.1 <- as.data.frame(Txi.lesion.coding.DGEList.LogCPM.filtered.norm)
colnames(V1.1.1) <- targets.lesion$sample
V1.1.1 <- melt(V1.1.1)
colnames(V1.1.1) <- c("sample","expression")

# plotting:
ggplot(V1, aes(x=sample, y=expression, fill=sample)) +
  geom_violin(trim = TRUE, show.legend = TRUE) +
  stat_summary(fun.y = "median", geom = "point", shape = 95, size = 10, color = "black") +
  theme_bw() +
  theme(legend.position = "none", axis.title=element_text(size=7),
        axis.title.x=element_blank(), axis.text=element_text(size=5),
        axis.text.x = element_text(angle = 90, hjust = 1),
        plot.title = element_text(size = 7)) +
  ggtitle("Raw dataset") +
  ggplot(V1.1, aes(x=sample, y=expression, fill=sample)) +
  geom_violin(trim = TRUE, show.legend = TRUE) +
  stat_summary(fun.y = "median", geom = "point", shape = 95, size = 10, color = "black") +
  theme_bw() +
  theme(legend.position = "none", axis.title=element_text(size=7),
        axis.title.x=element_blank(), axis.text=element_text(size=5),
        axis.text.x = element_text(angle = 90, hjust = 1),
        plot.title = element_text(size = 7)) +
  ggtitle("Filtered dataset") +
  ggplot(V1.1.1, aes(x=sample, y=expression, fill=sample)) +
  geom_violin(trim = TRUE, show.legend = TRUE) +
  stat_summary(fun.y = "median", geom = "point", shape = 95, size = 10, color = "black") +
  theme_bw() +
  theme(legend.position = "none", axis.title=element_text(size=7),
        axis.title.x=element_blank(), axis.text=element_text(size=5),
        axis.text.x = element_text(angle = 90, hjust = 1),
        plot.title = element_text(size = 7)) +
  ggtitle("Filtered and normalized dataset")

15 The unfiltered data

The following block in their dataset recreated the matrix without filtering and will use that for differential expression. It is a little hard to follow for me because they subset based on the sample numbers (8 to 28, which if I am not mistaken just drops the healthy samples?).

DataNotFiltered_Norm_OP <- calcNormFactors(Txi.lesion.coding.DGEList[,8:28],
                                           method = "TMM")
DataNotFiltered_Norm_log2CPM_OP <- edgeR::cpm(DataNotFiltered_Norm_OP, log=TRUE)
colnames(DataNotFiltered_Norm_log2CPM_OP) <- targets.onlypatients$sample
CPM_normData_notfiltered_OP <- 2^(DataNotFiltered_Norm_log2CPM_OP)
#uncomment the next line to produce raw data that was uploaded to the Gene Expression Omnibus (GEO) for publication.
#write.table(Txi.lesion.coding$counts, file = "Amorim_GEO_raw.txt", sep = "\t", quote = FALSE)

# Including all the individuals (HS and CL patients) for public domain submission:
DataNotFiltered_Norm <- calcNormFactors(Txi.lesion.coding.DGEList, method = "TMM")
DataNotFiltered_Norm_log2CPM <- edgeR::cpm(DataNotFiltered_Norm, log=TRUE)
colnames(DataNotFiltered_Norm_log2CPM) <- targets.lesion$sample
CPM_normData_notfiltered <- 2^(DataNotFiltered_Norm_log2CPM)
#uncomment the next line to produce the normalized data file that was uploaded to the Gene Expression Omnibus (GEO) for publication.
#write.table(DataNotFiltered_Norm_log2CPM, "Amorim_GEO_normalized.txt", sep = "\t", quote = FALSE)

16 The scott exploratory analysis

The following block generated a couple of the figures in the paper and comprise a pretty straightforward PCA. I am going to make a following block containing the same image with the cure/fail visualization using the same method/data.

pca.res <- prcomp(t(Txi.lesion.coding.DGEList.LogCPM.filtered.norm), scale.=F, retx=T)
pc.var <- pca.res$sdev^2
pc.per <- round(pc.var/sum(pc.var)*100, 1)
data.frame <- as.data.frame(pca.res$x)

# Calculate distance between samples by permanova:
allsamples.dist <- vegdist(t(2^Txi.lesion.coding.DGEList.LogCPM.filtered.norm),
                           method = "bray")

vegan <- adonis2(allsamples.dist~targets.lesion$disease,
                 data=targets.lesion,
                 permutations = 999, method="bray")

targets.lesion$disease
ggplot(data.frame, aes(x=PC1, y=PC2, color=factor(targets.lesion$disease))) +
  geom_point(size=5, shape=20) +
  theme_calc() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.text.x = element_text(size = 15, vjust = 0.5),
        axis.text.y = element_text(size = 15), axis.title = element_text(size = 15),
        legend.position="none") +
  scale_color_manual(values = c("#073F80","#EB512C")) +
  annotate("text", x=-50, y=80, label=paste("Permanova Pr(>F) =",
                                            vegan[1,5]), size=3, fontface="bold") +
  xlab(paste("PC1 -",pc.per[1],"%")) +
  ylab(paste("PC2 -",pc.per[2],"%")) +
  xlim(-200,110)

16.0.1 My most similar pca

I just realized that somewhere along the way in creating this container, I messed up this analysis pretty badly:

1. I dropped the 7 control samples.
2. I am comparing cure/fail but these analyses are all control/cutaneous.

When I originally did this on my workstation I had an actual 1:1 comparison and saw that our results were quite similar. I need to bring that back into this in order to show that neither we nor they are crazy people.

Either way, I think the main takeaway is that their dataset does not spend much time looking at cure/fail but instead control/infected for a reason.

Note, the fun aspects of the experiment (time to cure, size of lesion, etc) are not annotated in the metadata provided by SRA, but instead may be found in the capsule kindly provided by the lab. As a result, I copied that file into the sample_sheets/ directory and have added it to the expressionset. There is an important caveat, though: I did not include the non-diseased samples for this comparison; as a result the disease metadata factor is boring (e.g. it is only cutaneous).

external_cf[["accession"]] <- pData(external_cf)[["sample"]]
disease_factor <- pData(external_cf)[["disease"]]
table(disease_factor)

## disease_factor
## cutaneous 
##        21

external_disease <- set_expt_conditions(external_cf, fact = disease_factor)

## The numbers of samples by condition are:

## 
## cutaneous 
##        21

external_l2cpm <- normalize_expt(external_cf, filter = TRUE,
                                convert = "cpm", transform = "log2")

## Removing 7327 low-count genes (14154 remaining).

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

plot_pca(external_l2cpm, plot_labels = "repel")

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by cure, failure
## Shapes are defined by female, male.

Use the following block if you wish to bring together SRA-downloaded data with the experimental design from the Scott paper. It requires running the blocks above in which I loaded the capsule-derived metadata.

test <- pData(external_cf)
test_import <- as.data.frame(import)
test_import[["accession"]] <- pData(external_cf[["accession"]])
test_merged <- merge(test, import, by = "accession")

This is real comparison point to their cure/fail analysis.

16.1 Cure/Fail PCA using the same prcomp result

I am just copy/pasting their code again, but changing the color factor so that cure is purple, failure is red, and na(uninfected) is black.

The following plot should be the first direct comparison point between the two analysis pipelines. Thus, if you look back a few block at my invocation of plot_pca(external_norm), you will see a green/orange plot which is functionally identical if you note:

1. The x and y axes are flipped, which ok whatever it is PCA.
2. I excluded the healthy samples.
3. I dropped to gene level and used hisat.

With those caveats in mind, it is trivial to find the same relationshipes in the samples. E.g. the bottom red/purple individual samples are in the same relative position as my top orange/green pair. the same 4 samples are relative x-axis outliers (my right green, their left purple). The last 6 samples (my orange, their red) are all in the relative orientation.

I think I can further prove the similarity of our inputs via a direct comparison of the datastructures: Txi.lesion.coding.DGEList.LogCPM.filtered.norm (ugh what a name) vs. external_cf. In order to make that comparison, I need to rename my rows to the genecard IDs and the columns.

their_norm_exprs <- Txi.lesion.coding.DGEList.LogCPM.filtered.norm

my_hgnc_ids <- make.names(fData(external_cf)[["hgnc_symbol"]], unique = TRUE)
my_renamed <- set_expt_genenames(external_cf, ids = my_hgnc_ids)
my_norm <- normalize_expt(my_renamed, filter = TRUE, transform = "log2", convert = "cpm")
my_norm_exprs <- as.data.frame(exprs(my_norm))

our_exprs <- merge(their_norm_exprs, my_norm_exprs, by = "row.names")
rownames(our_exprs) <- our_exprs[["Row.names"]]
our_exprs[["Row.names"]] <- NULL
dim(our_exprs)

## I fully expected a correlation heatmap of the combined
## data to show a set of paired samples across the board.
## That is absolutely not true.
correlations <- plot_corheat(our_exprs)
correlations[["scatter"]]
correlations[["plot"]]

color_fact <- factor(targets.lesion$treatment_outcome)
levels(color_fact)
## Added by atb to see cure/fail on the same dataset
ggplot(data.frame, aes(x=PC1, y=PC2, color=color_fact)) +
  geom_point(size=5, shape=20) +
  theme_calc() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.text.x = element_text(size = 15, vjust = 0.5),
        axis.text.y = element_text(size = 15), axis.title = element_text(size = 15),
        legend.position="none") +
  scale_color_manual(values = c("purple", "red","black")) +
  annotate("text", x=-50, y=80, label=paste("Permanova Pr(>F) =",
                                            vegan[1,5]), size=3, fontface="bold") +
  xlab(paste("PC1 -",pc.per[1],"%")) +
  ylab(paste("PC2 -",pc.per[2],"%")) +
  xlim(-200,110)

17 DE comparisons

The following is their comparison of healthy tissue vs. CL lesion and Failure vs. Cure. I am going to follow it with my analagous examination using limma. Note, each of the pairs of variables created in the following block is xxx followed by xxx.treat; the former is healthy vs lesion and the latter is the fail vs cure set.

# Model matrices:
# CL lesions vs. HS:
design.lesion <- model.matrix(~0 + disease.lesion)
colnames(design.lesion) <- levels(disease.lesion)

# Failure vs. Cure:
design.lesion.treatment <- model.matrix(~0 + treatment.lesion)
colnames(design.lesion.treatment) <- levels(treatment.lesion)

myDGEList.lesion.coding <- DGEList(calcNorm1$counts)
myDGEList.OP.NotFil <- DGEList(CPM_normData_notfiltered_OP)

# Model mean-variance trend and fit linear model to data.
# Use VOOM function from Limma package to model the mean-variance relationship
normData.lesion.coding <- voom(myDGEList.lesion.coding, design.lesion)
normData.OP.NotFil <- voom(myDGEList.OP.NotFil, design.lesion.treatment)

colnames(normData.lesion.coding) <- targets.lesion$sample
colnames(normData.OP.NotFil) <- targets.onlypatients$sample

# fit a linear model to your data
fit.lesion.coding <- lmFit(normData.lesion.coding, design.lesion)
fit.lesion.coding.treatment <- lmFit(normData.OP.NotFil, design.lesion.treatment)

# contrast matrix
contrast.matrix.lesion <- makeContrasts(CL.vs.CON = cutaneous - control,
                                        levels=design.lesion)
contrast.matrix.lesion.treat <- makeContrasts(failure.vs.cure = failure - cure,
                                              levels=design.lesion.treatment)

# extract the linear model fit
fits.lesion.coding <- contrasts.fit(fit.lesion.coding,
                                    contrast.matrix.lesion)
fits.lesion.coding.treat <- contrasts.fit(fit.lesion.coding.treatment,
                                          contrast.matrix.lesion.treat)

# get bayesian stats for your linear model fit
ebFit.lesion.coding <- eBayes(fits.lesion.coding)
ebFit.lesion.coding.treat <- eBayes(fits.lesion.coding.treat)

# TopTable ----
allHits.lesion.coding <- topTable(ebFit.lesion.coding,
                                  adjust ="BH", coef=1,
                                  number=34935, sort.by="logFC")
allHits.lesion.coding.treat <- topTable(ebFit.lesion.coding.treat,
                                        adjust ="BH", coef=1,
                                        number=34776, sort.by="logFC")
myTopHits <- rownames_to_column(allHits.lesion.coding, "geneID")
myTopHits.treat <- rownames_to_column(allHits.lesion.coding.treat, "geneID")

# mutate the format of numeric values:
myTopHits <- mutate(myTopHits, log10Pval = round(-log10(adj.P.Val),2),
                    adj.P.Val = round(adj.P.Val, 2),
                    B = round(B, 2),
                    AveExpr = round(AveExpr, 2),
                    t = round(t, 2),
                    logFC = round(logFC, 2),
                    geneID = geneID)

myTopHits.treat <- mutate(myTopHits.treat, log10Pval = round(-log10(adj.P.Val),2),
                          adj.P.Val = round(adj.P.Val, 2),
                          B = round(B, 2),
                          AveExpr = round(AveExpr, 2),
                          t = round(t, 2),
                          logFC = round(logFC, 2),
                          geneID = geneID)
#save(myTopHits, file = "myTopHits")
#save(myTopHits.treat, file = "myTopHits.treat")

18 Perform my analagous limma analysis

my_filt <- normalize_expt(my_renamed, filter = "simple")
limma_cf <- limma_pairwise(my_filt, model_batch = FALSE)

my_table <- limma_cf[["all_tables"]][["failure_vs_cure"]]
their_table <- myTopHits.treat

dim(my_table)
dim(myTopHits.treat)
our_table <- merge(my_table, myTopHits.treat, by.x = "row.names", by.y = "geneID")
dim(our_table)
comparison <- plot_linear_scatter(our_table[, c("logFC.x", "logFC.y")])
comparison$scatter
comparison$correlation
comparison$lm_model

Ok, so there is a constituitive difference in our results, and it is significant. What does that mean for the set of genes observed?

With that said, in my most recent manual run of this, the results are quite good, I got a 0.75 correlation; I bet the primary outliers (on the axes) are just genes for which we got different gene<->tx mappings due to me using hisat and their usage of kallisto.

I guess I can test this hypothesis by just swapping in their counts into my data structure.

test_counts <- as.data.frame(myDGEList.lesion.coding[["counts"]])
test_counts[["host_HS01"]] <- NULL
test_counts[["host_HS02"]] <- NULL
test_counts[["host_HS03"]] <- NULL
test_counts[["host_HS04"]] <- NULL
test_counts[["host_HS05"]] <- NULL
test_counts[["host_HS06"]] <- NULL
test_counts[["host_HS07"]] <- NULL

dim(test_counts)
dim(exprs(my_test))
## Oh, that surprises me, the kallisto data has ~ 6k fewer genes?

19 See if there are shared DE genes

!!NOTE!! I am using a non-adjusted p-value filter here because I want to use the same filter they used for the volcano plot.

my_filter <- abs(my_table[["logFC"]]) > 1.0 & my_table[["P.Value"]] <= 0.05
sum(my_filter)
their_filter <- abs(their_table[["logFC"]]) > 1.0 & their_table[["P.Value"]] <= 0.05
sum(their_filter)

my_shared <- rownames(my_table)[my_filter] %in% their_table[their_filter, "geneID"]
sum(my_shared)

shared <- rownames(my_table)[my_filter]
shared[my_shared]

both <- list(
  "us" = rownames(my_table)[my_filter],
  "them" = their_table[their_filter, "geneID"])
tt <- UpSetR::fromList(both)
UpSetR::upset(tt)

19.1 Compare the two datasets directly

only_tmrc3 <- subset_expt(tmrc3_external, subset = "condition=='Colombia'") %>%
  set_expt_conditions(fact = "finaloutcome")

## subset_expt(): There were 39, now there are 18 samples.

## The numbers of samples by condition are:

## 
##    cure failure 
##      13       5

only_tmrc3_de <- all_pairwise(only_tmrc3, model_batch = "svaseq",
                              parallel = parallel, filter = TRUE,
                              methods = methods)

## 
##    cure failure 
##      13       5

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

## Setting 225 low elements to zero.

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

only_tmrc3_de

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

## The logFC agreement among the methods follows:

##                 falr_vs_cr
## limma_vs_deseq      0.7154
## limma_vs_edger      0.8311
## limma_vs_basic      0.9061
## limma_vs_noiseq    -0.8065
## deseq_vs_edger      0.9247
## deseq_vs_basic      0.7781
## deseq_vs_noiseq    -0.7249
## edger_vs_basic      0.8963
## edger_vs_noiseq    -0.8438
## basic_vs_noiseq    -0.8787

only_tmrc3_table <- combine_de_tables(only_tmrc3_de)
only_tmrc3_table

## A set of combined differential expression results.

##             table deseq_sigup deseq_sigdown edger_sigup edger_sigdown
## 1 failure_vs_cure          27            26          28            15
##   limma_sigup limma_sigdown
## 1           1             0

## `geom_line()`: Each group consists of only one observation.
## i Do you need to adjust the group aesthetic?

## Plot describing unique/shared genes in a differential expression table.

only_tmrc3_top100 <- extract_significant_genes(only_tmrc3_table, n = 100)
only_tmrc3_up <- only_tmrc3_top100[["deseq"]][["ups"]][["failure_vs_cure"]]
only_tmrc3_down <- only_tmrc3_top100[["deseq"]][["downs"]][["failure_vs_cure"]]

tmrc3_external_de <- all_pairwise(tmrc3_external, model_batch = "svaseq",
                                  parallel = parallel, filter = "simple",
                                  methods = methods)

## 
##   Brazil Colombia 
##       21       18

## Potentially check over the experimental design, there appear to be missing values.

## Warning in plot_pca(pre_batch, plot_labels = FALSE, ...): There are NA values
## in the component data.  This can lead to weird plotting errors.

## Removing 4594 low-count genes (14577 remaining).

## Setting 3653 low elements to zero.

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

## Potentially check over the experimental design, there appear to be missing values.

## Warning in plot_pca(post_batch, plot_labels = FALSE, ...): There are NA values
## in the component data.  This can lead to weird plotting errors.

tmrc3_external_table <- combine_de_tables(tmrc3_external_de, excel = "excel/tmrc3_scott_biopsies.xlsx")

## Error : colNames must be a unique vector (case sensitive)

tmrc3_external_sig <- extract_significant_genes(tmrc3_external_table, excel = "excel/tmrc3_scott_biopsies_sig.xlsx")

tmrc3_external_cf <- set_expt_conditions(tmrc3_external, fact = "finaloutcome") %>%
  set_expt_batches(fact = "lab")

## Warning in `[<-.factor`(`*tmp*`, null_ids, value = "null"): invalid factor
## level, NA generated

## The numbers of samples by condition are:

## 
##    cure failure 
##      13       5

## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': 'names' attribute [3] must be the same length as the vector [2]

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

## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'tmrc3_external_cf' not found

plot_pca(tmrc3_external_cf_norm)

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

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

## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'tmrc3_external_cf' not found

plot_pca(tmrc3_external_cf_nb)

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

tmrc3_external_cf_de <- all_pairwise(tmrc3_external_cf, model_batch = "svaseq",
                                     parallel = parallel, filter = TRUE,
                                     methods = methods)

## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': object 'tmrc3_external_cf' not found

tmrc3_external_cf_de

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

tmrc3_external_cf_table <- combine_de_tables(tmrc3_external_cf_de, excel = "excel/tmrc3_scott_cf_table.xlsx")

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

tmrc3_external_cf_table

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

tmrc3_external_cf_sig <- extract_significant_genes(tmrc3_external_cf_table, excel = "excel/tmrc3_scott_cf_sig.xlsx")

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

tmrc3_external_cf_sig

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

tmrc3_external_species <- set_expt_conditions(tmrc3_external, fact = "ParasiteSpecies") %>%
  set_expt_colors(color_choices[["parasite"]])

## The numbers of samples by condition are:

## 
## lvbraziliensis   lvpanamensis  notapplicable 
##             22             14              3

## Warning in set_expt_colors(., color_choices[["parasite"]]): Colors for the
## following categories are not being used: lvguyanensis.

20 Compare the l2FC values

Let us look at the top/bottom 100 genes of these two datasets and see if they have any similarities.

Note to self, set up s4 dispatch on compare_de_tables!

compared <- compare_de_tables(only_tmrc3_table, external_table, first_table = 1, second_table = 1)
compared$scatter

compared$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[[xcol]] and df[[ycol]]
## t = 14, df = 13240, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1033 0.1368
## sample estimates:
##    cor 
## 0.1201

21 Compare visits by celltype and C/F

I assume this request came out of the review process, but I am not quite sure where to put it. If I understand it correctly, the goal is to look across visits for combinations of cure and fail (not fail/cure, but v2/v1) and across cell types.

Thus, in order to do this, I will need to combine those three parameters or set up a more complex model to handle this.

t_cellvisitcf <- set_expt_conditions(t_clinical_nobiop, fact = "cell_visit_cf")

## The numbers of samples by condition are:

## 
##    eosinophils_1_cure eosinophils_1_failure    eosinophils_2_cure 
##                     5                     3                     6 
## eosinophils_2_failure    eosinophils_3_cure eosinophils_3_failure 
##                     3                     6                     3 
##      monocytes_1_cure   monocytes_1_failure      monocytes_2_cure 
##                     8                     8                     7 
##   monocytes_2_failure      monocytes_3_cure   monocytes_3_failure 
##                     6                     6                     7 
##    neutrophils_1_cure neutrophils_1_failure    neutrophils_2_cure 
##                     8                     8                     7 
## neutrophils_2_failure    neutrophils_3_cure neutrophils_3_failure 
##                     6                     5                     7

t_cellvisitcf_de <- all_pairwise(t_cellvisitcf, keepers = visittype_contrasts,
                                 model_batch = "svaseq", filter = TRUE, parallel = parallel,
                                 methods = methods)

## 
##    eosinophils_1_cure eosinophils_1_failure    eosinophils_2_cure 
##                     5                     3                     6 
## eosinophils_2_failure    eosinophils_3_cure eosinophils_3_failure 
##                     3                     6                     3 
##      monocytes_1_cure   monocytes_1_failure      monocytes_2_cure 
##                     8                     8                     7 
##   monocytes_2_failure      monocytes_3_cure   monocytes_3_failure 
##                     6                     6                     7 
##    neutrophils_1_cure neutrophils_1_failure    neutrophils_2_cure 
##                     8                     8                     7 
## neutrophils_2_failure    neutrophils_3_cure neutrophils_3_failure 
##                     6                     5                     7

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

## Setting 5938 low elements to zero.

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

t_cellvisitcf_de

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

t_cellvisitcf_mono_table <- combine_de_tables(
  t_cellvisitcf_de, keepers = visittype_contrasts_mono,
  excel = glue("{xlsx_prefix}/DE_Visits/Cure_Fail/monocyte_visit_cf_combined_table_sva-v{ver}.xlsx"))
t_cellvisitcf_mono_table

## A set of combined differential expression results.

##                                        table deseq_sigup deseq_sigdown
## 1       monocytes_2_cure_vs_monocytes_1_cure           0             2
## 2 monocytes_2_failure_vs_monocytes_1_failure           2             2
## 3       monocytes_3_cure_vs_monocytes_1_cure           1             3
## 4 monocytes_3_failure_vs_monocytes_1_failure           1             3
##   edger_sigup edger_sigdown limma_sigup limma_sigdown
## 1           0             0           0             0
## 2           1             1           1             0
## 3           0             0           0             0
## 4           0             0           0             0

## Plot describing unique/shared genes in a differential expression table.

t_cellvisitcf_mono_sig <- extract_significant_genes(
  t_cellvisitcf_mono_table,
  excel = glue("{xlsx_prefix}/DE_Visits/Cure_Fail/monocyte_visit_cf_combined_sig_sva-v{ver}.xlsx"))
t_cellvisitcf_mono_sig

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##                   limma_up limma_down edger_up edger_down deseq_up deseq_down
## v2v1_mono_cure           0          0        0          0        0          2
## v2v1_mono_failure        1          0        1          1        2          2
## v3v1_mono_cure           0          0        0          0        1          3
## v3v1_mono_failure        0          0        0          0        1          3
##                   basic_up basic_down
## v2v1_mono_cure           0          0
## v2v1_mono_failure        0          0
## v3v1_mono_cure           0          0
## v3v1_mono_failure        0          0

t_cellvisitcf_neut_table <- combine_de_tables(
  t_cellvisitcf_de, keepers = visittype_contrasts_ne,
  excel = glue("{xlsx_prefix}/DE_Visits/Cure_Fail/neutrophil_visit_cf_combined_table_sva-v{ver}.xlsx"))
t_cellvisitcf_neut_table

## A set of combined differential expression results.

##                                            table deseq_sigup deseq_sigdown
## 1       neutrophils_2_cure_vs_neutrophils_1_cure          85           132
## 2 neutrophils_2_failure_vs_neutrophils_1_failure         127           150
## 3       neutrophils_3_cure_vs_neutrophils_1_cure         105           195
## 4 neutrophils_3_failure_vs_neutrophils_1_failure          87            24
##   edger_sigup edger_sigdown limma_sigup limma_sigdown
## 1          75           118          90            31
## 2         116           175         105            42
## 3         110           157         114            39
## 4          77            20          56            36

## Plot describing unique/shared genes in a differential expression table.

t_cellvisitcf_neut_sig <- extract_significant_genes(
  t_cellvisitcf_neut_table,
  excel = glue("{xlsx_prefix}/DE_Visits/Cure_Fail/neutrophil_visit_cf_combined_sig_sva-v{ver}.xlsx"))
t_cellvisitcf_neut_sig

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##                 limma_up limma_down edger_up edger_down deseq_up deseq_down
## v2v1_ne_cure          90         31       75        118       85        132
## v2v1_ne_failure      105         42      116        175      127        150
## v3v1_ne_cure         114         39      110        157      105        195
## v3v1_ne_failure       56         36       77         20       87         24
##                 basic_up basic_down
## v2v1_ne_cure           2          0
## v2v1_ne_failure        4          0
## v3v1_ne_cure           0          0
## v3v1_ne_failure        0          0

t_cellvisitcf_eo_table <- combine_de_tables(
  t_cellvisitcf_de, keepers = visittype_contrasts_eo,
  excel = glue("{xlsx_prefix}/DE_Visits/Cure_Fail/eosinophil_visit_cf_combined_table_sva-v{ver}.xlsx"))
t_cellvisitcf_eo_table

## A set of combined differential expression results.

##                                            table deseq_sigup deseq_sigdown
## 1       eosinophils_2_cure_vs_eosinophils_1_cure           5             1
## 2 eosinophils_2_failure_vs_eosinophils_1_failure           1             5
## 3       eosinophils_3_cure_vs_eosinophils_1_cure           9             1
## 4 eosinophils_3_failure_vs_eosinophils_1_failure           0             8
##   edger_sigup edger_sigdown limma_sigup limma_sigdown
## 1           0             0           0             0
## 2           0             0           0             0
## 3           0             0           0             1
## 4           0             0           0             0

## Plot describing unique/shared genes in a differential expression table.

t_cellvisitcf_eo_sig <- extract_significant_genes(
  t_cellvisitcf_eo_table,
  excel = glue("{xlsx_prefix}/DE_Visits/Cure_Fail/eosinophil_visit_cf_combined_sig_sva-v{ver}.xlsx"))
t_cellvisitcf_eo_sig

## A set of genes deemed significant according to limma, edger, deseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##                 limma_up limma_down edger_up edger_down deseq_up deseq_down
## v2v1_eo_cure           0          0        0          0        5          1
## v2v1_eo_failure        0          0        0          0        1          5
## v3v1_eo_cure           0          1        0          0        9          1
## v3v1_eo_failure        0          0        0          0        0          8
##                 basic_up basic_down
## v2v1_eo_cure           0          0
## v2v1_eo_failure        0          0
## v3v1_eo_cure           0          0
## v3v1_eo_failure        0          0

tmp <- loadme(filename = savefile)

Bibliography

Chung, Matthew, Vincent M. Bruno, David A. Rasko, Christina A. Cuomo, José F. Muñoz, Jonathan Livny, Amol C. Shetty, Anup Mahurkar, and Julie C. Dunning Hotopp. 2021. “Best Practices on the Differential Expression Analysis of Multi-Species RNA-seq.” Genome Biology 22 (April): 121. https://doi.org/10.1186/s13059-021-02337-8.

Hoffman, Gabriel E, and Panos Roussos. 2020. “Dream: Powerful Differential Expression Analysis for Repeated Measures Designs.” Bioinformatics 37 (2): 192–201. https://doi.org/10.1093/bioinformatics/btaa687.

Leng, Ning, John A. Dawson, James A. Thomson, Victor Ruotti, Anna I. Rissman, Bart M. G. Smits, Jill D. Haag, Michael N. Gould, Ron M. Stewart, and Christina Kendziorski. 2013. “EBSeq: An Empirical Bayes Hierarchical Model for Inference in RNA-seq Experiments.” Bioinformatics 29 (8): 1035–43. https://doi.org/10.1093/bioinformatics/btt087.

Love, Michael I., Wolfgang Huber, and Simon Anders. 2014. “Moderated Estimation of Fold Change and Dispersion for RNA-Seq Data with DESeq2.” bioRxiv. https://doi.org/10.1101/002832.

McCarthy, Davis J., Yunshun Chen, and Gordon K. Smyth. 2012. “Differential Expression Analysis of Multifactor RNA-Seq Experiments with Respect to Biological Variation.” Nucleic Acids Research 40 (10): 4288–97. https://doi.org/10.1093/nar/gks042.

Molania, Ramyar, Momeneh Foroutan, Johann A. Gagnon-Bartsch, Luke C. Gandolfo, Aryan Jain, Abhishek Sinha, Gavriel Olshansky, Alexander Dobrovic, Anthony T. Papenfuss, and Terence P. Speed. 2023. “Removing Unwanted Variation from Large-Scale RNA Sequencing Data with PRPS.” Nature Biotechnology 41 (1): 82–95. https://doi.org/10.1038/s41587-022-01440-w.

Risso, Davide, John Ngai, Terence P. Speed, and Sandrine Dudoit. 2014. “Normalization of RNA-seq Data Using Factor Analysis of Control Genes or Samples.” Nature Biotechnology 32 (9): 896–902. https://doi.org/10.1038/nbt.2931.

Ritchie, Matthew E., Belinda Phipson, Di Wu, Yifang Hu, Charity W. Law, Wei Shi, and Gordon K. Smyth. 2015. “Limma Powers Differential Expression Analyses for RNA-sequencing and Microarray Studies.” Nucleic Acids Research 43 (7): e47–47. https://doi.org/10.1093/nar/gkv007.

Tarazona, Sonia, Pedro Furió-Tarí, David Turrà, Antonio Di Pietro, María José Nueda, Alberto Ferrer, and Ana Conesa. 2015. “Data Quality Aware Analysis of Differential Expression in RNA-seq with NOISeq R/Bioc Package.” Nucleic Acids Research 43 (21): e140. https://doi.org/10.1093/nar/gkv711.