TMRC3 202305: Differential Expression analyses, Tumaco only

1 Changelog

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

2.1 Naming conventions

I am going to try to standardize how I name the various data structures created in this document. Most of the large data created are either sets of differential expression analyses, their combined results, or the set of results deemed ‘significant’.

Hopefully by now they all follow these guidelines:

{clinic(s)}sample-subset}{primary-question(s)}{datatype}{batch-method}

• {clinic}: This is either tc or t for Tumaco and Cali, or just Tumaco.
• {sample-subset}: Things like ‘all’ or ‘monocytes’.
• {primary-question}: Shorthand name for the primary contrasts performed, thus ‘clinics’ would suggest a comparison of Tumaco vs. Cali. ‘visits’ would compare v2/v1, etc.
• {datatype}: de, table, sig
• {batch-type}: nobatch, batch{factor}, sva. {factor} in this instance should be a column from the metadata.

With this in mind, ‘tc_biopsies_clinic_de_sva’ should be the Tumaco+Cali biopsy data after performing the differential expression analyses comparing the clinics using sva.

I suspect there remain some exceptions and/or errors.

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

clinic_contrasts <- list(
  "clinics" = c("Cali", "Tumaco"))
## In some cases we have no Cali failure samples, so there remain only 2
## contrasts that are likely of interest
tc_cf_contrasts <- list(
  "tumaco" = c("Tumacofailure", "Tumacocure"),
  "cure" = c("Tumacocure", "Calicure"))
## In other cases, we have cure/fail for both places.
clinic_cf_contrasts <- list(
  "cali" = c("Califailure", "Calicure"),
  "tumaco" = c("Tumacofailure", "Tumacocure"),
  "cure" = c("Tumacocure", "Calicure"),
  "fail" = c("Tumacofailure", "Califailure"))
cf_contrast <- list(
  "outcome" = c("Tumacofailure", "Tumacocure"))
t_cf_contrast <- list(
  "outcome" = c("failure", "cure"))
visitcf_contrasts <- list(
  "v1cf" = c("v1failure", "v1cure"),
  "v2cf" = c("v2failure", "v2cure"),
  "v3cf" = c("v3failure", "v3cure"))
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"))

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.

xlsx_prefix <- "analyses/4_tumaco/DE_Cure_vs_Fail"

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", filter = TRUE)

## 
##    cure failure 
##      67      56

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

## Setting 17282 low elements to zero.

## transform_counts: Found 17282 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.

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

## The primary analysis performed 6 comparisons.

## The logFC agreement among the methods follows:

##                failure_vs_cure
## limma_vs_deseq          0.8064
## limma_vs_edger          0.8179
## limma_vs_basic          0.8702
## deseq_vs_edger          0.9844
## deseq_vs_basic          0.8243
## edger_vs_basic          0.8286

t_cf_clinical_table_sva <- combine_de_tables(
  t_cf_clinical_de_sva, keepers = t_cf_contrast,
#  rda = glue("rda/t_clinical_cf_table_sva-v{ver}.rda"),
  excel = glue("{xlsx_prefix}/All_Samples/t_clinical_cf_tables_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/All_Samples/t_clinical_cf_tables_sva-v202305.xlsx before writing the tables.

## Adding venn plots for outcome.

t_cf_clinical_table_sva

## A set of combined differential expression results.

##             table deseq_sigup deseq_sigdown edger_sigup edger_sigdown limma_sigup
## 1 failure_vs_cure          93           183         103           157          50
##   limma_sigdown
## 1            38

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

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("{xlsx_prefix}/All_Samples/t_clinical_cf_sig_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/All_Samples/t_clinical_cf_sig_sva-v202305.xlsx before writing the tables.

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 basic_down
## outcome       50         38      103        157       93        183       16          4

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

## [1] 93 50

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

## [1] 183  50

3.2 gProfiler search of all samples

The following gProfiler searches use the all_gprofiler() function instead of simple_gprofiler(). As a result, the results are separated by {contrast}_{direction}. Thus ‘outcome_down’.

The same plots are available as the previous gProfiler searches, but in many of the following runs, I used the dotplot() function to get a slightly different view of the results.

t_cf_clinical_gp <- all_gprofiler(t_cf_clinical_sig_sva)
## Wikipathways of the up c/f genes
enrichplot::dotplot(t_cf_clinical_gp[["outcome_up"]][["WP_enrich"]])

## Transcription factor database of the up c/f genes
enrichplot::dotplot(t_cf_clinical_gp[["outcome_up"]][["TF_enrich"]])

## Reactome of the up c/f genes
enrichplot::dotplot(t_cf_clinical_gp[["outcome_up"]][["REAC_enrich"]])

## GO of the down c/f genes
enrichplot::dotplot(t_cf_clinical_gp[["outcome_down"]][["GO_enrich"]])

t_cf_clinical_gp[["outcome_up"]][["pvalue_plots"]][["BP"]]

## Reactome of the down c/f genes
enrichplot::dotplot(t_cf_clinical_gp[["outcome_down"]][["GO_enrich"]])

enrichplot::dotplot(t_cf_clinical_gp[["outcome_up"]][["GO_enrich"]])

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", filter = TRUE)

## 
## first later 
##    40    69

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

## Setting 9616 low elements to zero.

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

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

## Deleting the file excel/tv1_vs_later_tables-v202305.xlsx before writing the tables.

## Adding venn plots for later_vs_first.

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

## Deleting the file excel/tv1_vs_later_sig-v202305.xlsx before writing the tables.

tv1later_gp <- all_gprofiler(tv1_vs_later_sig)
tv1later_gp[[1]]$pvalue_plots$BP

tv1later_gp[[2]]$pvalue_plots$BP

5 Sex comparison

Can we observe consistent difference in the fe/male samples?

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

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

t_sex

## An expressionSet containing experiment with 19923
## genes and 123 samples. There are 156 metadata columns and
## 14 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", filter = TRUE)

## 
## female   male 
##     22    101

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

## Setting 17259 low elements to zero.

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

t_sex_de

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

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

## The primary analysis performed 6 comparisons.

## The logFC agreement among the methods follows:

##                male_vs_female
## limma_vs_deseq         0.8603
## limma_vs_edger         0.8672
## limma_vs_basic         0.9482
## deseq_vs_edger         0.9909
## deseq_vs_basic         0.8712
## edger_vs_basic         0.8758

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

## Deleting the file excel/t_sex_table-v202305.xlsx before writing the tables.

## Adding venn plots for male_vs_female.

t_sex_table

## A set of combined differential expression results.

##            table deseq_sigup deseq_sigdown edger_sigup edger_sigdown limma_sigup
## 1 male_vs_female         128            96         116            95          53
##   limma_sigdown
## 1            74

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

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

## Deleting the file excel/t_sex_sig-v202305.xlsx before writing the tables.

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 basic_up
## male_vs_female       53         74      116         95      128         96       15
##                basic_down
## male_vs_female         10

t_sex_gp <- all_gprofiler(t_sex_sig)
t_sex_gp

## Running gProfiler on every set of significant genes found:

##                     GO KEGG REAC WP TF MIRNA HPA CORUM HP
## male_vs_female_up   63    0    4  1 46     0   0     0  2
## male_vs_female_down 22    0    0  0  0     0   0     0  0

What if we limit the question to only the people who cured?

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", filter = TRUE)

## 
## female   male 
##     13     54

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

## Setting 8959 low elements to zero.

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

t_sex_cure_de

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

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

## The primary analysis performed 6 comparisons.

## The logFC agreement among the methods follows:

##                male_vs_female
## limma_vs_deseq         0.7789
## limma_vs_edger         0.8383
## limma_vs_basic         0.9281
## deseq_vs_edger         0.9276
## deseq_vs_basic         0.7981
## edger_vs_basic         0.8480

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

## Deleting the file excel/t_sex_cure_table-v202305.xlsx before writing the tables.

## Adding venn plots for male_vs_female.

t_sex_cure_table

## A set of combined differential expression results.

##            table deseq_sigup deseq_sigdown edger_sigup edger_sigdown limma_sigup
## 1 male_vs_female         172           129         161           143          63
##   limma_sigdown
## 1           107

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

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

## Deleting the file excel/t_sex_cure_sig-v202305.xlsx before writing the tables.

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 basic_up
## male_vs_female       63        107      161        143      172        129       12
##                basic_down
## male_vs_female          5

t_sex_cure_gp <- all_gprofiler(t_sex_cure_sig)
t_sex_cure_gp

## Running gProfiler on every set of significant genes found:

##                     GO KEGG REAC WP TF MIRNA HPA CORUM HP
## male_vs_female_up   63    2    7  0 35     0   0     0  9
## male_vs_female_down  0    0    0  0  0     0   0     0  0

t_sex_cure_gp[[1]][["pvalue_plots"]][["BP"]]

6 Ethnicity comparisons

The set of ethnicities observed in Tumaco is a bit different than that in Cali. Can we see differences among those groups? Note that this is confounded with cure/fail.

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

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

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

## Setting 15817 low elements to zero.

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

t_ethnicity_de

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

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

## The primary analysis performed 6 comparisons.

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

## Deleting the file excel/t_ethnicity_table-v202305.xlsx before writing the tables.

## Adding venn plots for indigena_vs_afrocol.

## Adding venn plots for mestiza_vs_afrocol.

## Adding venn plots for mestiza_vs_indigena.

t_ethnicity_table

## A set of combined differential expression results.

##                 table deseq_sigup deseq_sigdown edger_sigup edger_sigdown limma_sigup
## 1 indigena_vs_afrocol         162           236         186           216         164
## 2  mestiza_vs_afrocol          56            92          51            96          41
## 3 mestiza_vs_indigena          83            97          67           108          58
##   limma_sigdown
## 1           146
## 2            53
## 3            56

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

## Deleting the file excel/t_ethnicity_sig-v202305.xlsx before writing the tables.

t_ethnicity_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 basic_up
## indigena_vs_afrocol      164        146      186        216      162        236       16
## mestiza_vs_afrocol        41         53       51         96       56         92        2
## mestiza_vs_indigena       58         56       67        108       83         97        2
##                     basic_down
## indigena_vs_afrocol         17
## mestiza_vs_afrocol           9
## mestiza_vs_indigena          2

t_ethnicity_gp <- all_gprofiler(t_ethnicity_sig)
t_ethnicity_gp

## Running gProfiler on every set of significant genes found:

##                          GO KEGG REAC WP TF MIRNA HPA CORUM HP
## indigena_vs_afrocol_up   57    1    2  0  0     0   1     0  0
## indigena_vs_afrocol_down 26    0    0  0  0     0   0     0  0
## mestiza_vs_afrocol_up     6    0    0  0  0     0   0     0  0
## mestiza_vs_afrocol_down  21    0    4  3  1     0   0     0  0
## mestiza_vs_indigena_up   20    0    2  0  7     0   0     0  0
## mestiza_vs_indigena_down 10    0    1  0  5     5   0     0  0

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", filter = TRUE)

## 
##    cure failure 
##      30      24

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

## Setting 7615 low elements to zero.

## transform_counts: Found 7615 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.

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

## The primary analysis performed 6 comparisons.

## The logFC agreement among the methods follows:

##                failure_vs_cure
## limma_vs_deseq          0.7393
## limma_vs_edger          0.7826
## limma_vs_basic          0.6896
## deseq_vs_edger          0.9534
## deseq_vs_basic          0.6956
## edger_vs_basic          0.7235

t_cf_clinical_v1_table_sva <- combine_de_tables(
  t_cf_clinical_v1_de_sva, keepers = t_cf_contrast,
#  rda = glue("rda/t_clinical_v1_cf_table_sva-v{ver}.rda"),
  excel = glue("{xlsx_prefix}/All_Samples/t_clinical_v1_cf_tables_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/All_Samples/t_clinical_v1_cf_tables_sva-v202305.xlsx before writing the tables.

## Adding venn plots for outcome.

t_cf_clinical_v1_table_sva

## A set of combined differential expression results.

##             table deseq_sigup deseq_sigdown edger_sigup edger_sigdown limma_sigup
## 1 failure_vs_cure          28            74          28            54           3
##   limma_sigdown
## 1             3

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

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

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/All_Samples/t_clinical_v1_cf_sig_sva-v202305.xlsx before writing the tables.

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

## [1] 28 50

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

## [1] 74 50

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", filter = TRUE)

## 
##    cure failure 
##      20      15

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

## Setting 2848 low elements to zero.

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

t_cf_clinical_v2_table_sva <- combine_de_tables(
  t_cf_clinical_v2_de_sva, keepers = t_cf_contrast,
#  rda = glue("rda/t_clinical_v2_cf_table_sva-v{ver}.rda"),
  excel = glue("{xlsx_prefix}/All_Samples/t_clinical_v2_cf_tables_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/All_Samples/t_clinical_v2_cf_tables_sva-v202305.xlsx before writing the tables.

## Adding venn plots for outcome.

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

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/All_Samples/t_clinical_v2_cf_sig_sva-v202305.xlsx before writing the tables.

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

## [1] 51 50

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

## [1] 15 50

6.0.1.3 Cure/Fail, Tumaco Visit 3

t_cf_clinical_v3_de_sva <- all_pairwise(tv3_samples, model_batch = "svaseq", filter = TRUE)

## 
##    cure failure 
##      17      17

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

## Setting 1878 low elements to zero.

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

t_cf_clinical_v3_table_sva <- combine_de_tables(
  t_cf_clinical_v3_de_sva, keepers = t_cf_contrast,
#  rda = glue("rda/t_clinical_v3_cf_table_sva-v{ver}.rda"),
  excel = glue("{xlsx_prefix}/All_Samples/t_clinical_v3_cf_tables_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/All_Samples/t_clinical_v3_cf_tables_sva-v202305.xlsx before writing the tables.

## Adding venn plots for outcome.

t_cf_clinical_v3_sig_sva <- extract_significant_genes(
  t_cf_clinical_v3_table_sva,
  excel = glue("{xlsx_prefix}/All_Samples/t_clinical_v3_cf_sig_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/All_Samples/t_clinical_v3_cf_sig_sva-v202305.xlsx before writing the tables.

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

## [1] 120  50

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

## [1] 62 50

6.0.1.4 Visit 1 gProfiler searches

It looks like there are very few groups in the visit 1 significant genes.

t_cf_clinical_v1_sig_sva_gp <- all_gprofiler(t_cf_clinical_v1_sig_sva)

## Wikipathways of the up c/f genes
enrichplot::dotplot(t_cf_clinical_v1_sig_sva_gp[["outcome_up"]][["GO_enrich"]])

enrichplot::dotplot(t_cf_clinical_v1_sig_sva_gp[["outcome_down"]][["GO_enrich"]])

6.0.1.5 Visit 2 gProfiler searches

Up: 74 GO, 4 KEGG, 6 reactome, 4 WP, 56 TF, 1 miRNA, 0 HP/HPA/CORUM. Down: 19 GO, 1 KEGG, 1 HP, 2 HPA, 0 reactome/wp/tf/corum

t_cf_clinical_v2_sig_sva_gp <- all_gprofiler(t_cf_clinical_v2_sig_sva)

## Wikipathways of the up c/f genes
enrichplot::dotplot(t_cf_clinical_v2_sig_sva_gp[["outcome_up"]][["GO_enrich"]])

enrichplot::dotplot(t_cf_clinical_v2_sig_sva_gp[["outcome_up"]][["REAC_enrich"]])

enrichplot::dotplot(t_cf_clinical_v2_sig_sva_gp[["outcome_up"]][["TF_enrich"]])

enrichplot::dotplot(t_cf_clinical_v2_sig_sva_gp[["outcome_down"]][["GO_enrich"]])

6.0.1.6 Visit 3 gProfiler searches

Up: 120 genes; 141 GO, 1 KEGG, 5 Reactome, 2 WP, 30 TF, 1 miRNA, 0 HPA/CORUM/HP Down: 62 genes; 30 GO, 2 KEGG, 1 Reactome, 0 WP/TF/miRNA/HPA/CORUM/HP,

t_cf_clinical_v3_sig_sva_gp <- all_gprofiler(t_cf_clinical_v3_sig_sva)

## Wikipathways of the up c/f genes
enrichplot::dotplot(t_cf_clinical_v3_sig_sva_gp[["outcome_up"]][["GO_enrich"]])

enrichplot::dotplot(t_cf_clinical_v3_sig_sva_gp[["outcome_up"]][["REAC_enrich"]])

enrichplot::dotplot(t_cf_clinical_v3_sig_sva_gp[["outcome_up"]][["TF_enrich"]])

enrichplot::dotplot(t_cf_clinical_v3_sig_sva_gp[["outcome_down"]][["GO_enrich"]])

6.0.2 Repeat no biopsies

The biopsy samples are problematic for a few reasons, so let us repeat without them.

t_cf_clinical_nobiop_de_sva <- all_pairwise(t_clinical_nobiop,
                                            model_batch = "svaseq", filter = TRUE)

## 
##    cure failure 
##      58      51

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

## Setting 9578 low elements to zero.

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

t_cf_clinical_nobiop_table_sva <- combine_de_tables(
  t_cf_clinical_nobiop_de_sva, keepers = t_cf_contrast,
#  rda = glue("rda/t_clinical_nobiop_cf_table_sva-v{ver}.rda"),
  excel = glue("{xlsx_prefix}/No_Biopsies/t_clinical_nobiop_cf_tables_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/No_Biopsies/t_clinical_nobiop_cf_tables_sva-v202305.xlsx before writing the tables.

## Adding venn plots for outcome.

t_cf_clinical_nobiop_sig_sva <- extract_significant_genes(
  t_cf_clinical_nobiop_table_sva,
  excel = glue("{xlsx_prefix}/No_Biopsies/t_clinical_nobiop_cf_sig_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/No_Biopsies/t_clinical_nobiop_cf_sig_sva-v202305.xlsx before writing the tables.

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

## [1] 137  50

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

## [1] 73 50

6.0.2.1 gProfiler: Clinical no biopsies

Up: 137 genes; 88 GO, 0 KEGG, 6 Reactome, 1 WP, 46 TF, 1 miRNA, 0 others Down: 73 genes; 78 GO, 1 KEGG, 1 Reactome, 9 TF, 0 others

t_cf_clinical_nobiop_sig_sva_gp <- all_gprofiler(t_cf_clinical_nobiop_sig_sva)

enrichplot::dotplot(t_cf_clinical_nobiop_sig_sva_gp[["outcome_up"]][["GO_enrich"]])

enrichplot::dotplot(t_cf_clinical_nobiop_sig_sva_gp[["outcome_up"]][["TF_enrich"]])

enrichplot::dotplot(t_cf_clinical_nobiop_sig_sva_gp[["outcome_down"]][["GO_enrich"]])

6.0.3 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.3.1 Cure/Fail, Biopsies

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

## 
##    Tumaco_cure Tumaco_failure 
##              9              5

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

## Setting 145 low elements to zero.

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

t_cf_biopsy_table_sva <- combine_de_tables(
  t_cf_biopsy_de_sva, keepers = cf_contrast,
#  rda = glue("rda/t_biopsy_cf_table_sva-v{ver}.rda"),
  excel = glue("{xlsx_prefix}/Biopsies/t_biopsy_cf_tables_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Biopsies/t_biopsy_cf_tables_sva-v202305.xlsx before writing the tables.

## Adding venn plots for outcome.

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

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Biopsies/t_cf_biopsy_sig_sva-v202305.xlsx before writing the tables.

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

## [1] 17 50

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

## [1] 11 50

6.0.3.2 gProfiler: Biopsies

Up: 17 genes; 74 GO, 3 KEGG, 1 Reactome, 3 WP, 1 TF, 0 others Down: 11 genes; 2 GO, 0 others

t_cf_biopsy_sig_sva_gp <- all_gprofiler(t_cf_biopsy_sig_sva)

enrichplot::dotplot(t_cf_biopsy_sig_sva_gp[["outcome_up"]][["GO_enrich"]])

enrichplot::dotplot(t_cf_biopsy_sig_sva_gp[["outcome_up"]][["WP_enrich"]])

enrichplot::dotplot(t_cf_biopsy_sig_sva_gp[["outcome_down"]][["GO_enrich"]])

6.0.3.3 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",
                                     filter = TRUE)

## 
##    Tumaco_cure Tumaco_failure 
##             21             21

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

## Setting 730 low elements to zero.

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

t_cf_monocyte_tables_sva <- combine_de_tables(
  t_cf_monocyte_de_sva, keepers = cf_contrast,
#  rda = glue("rda/t_monocyte_cf_table_sva-v{ver}.rda"),
  excel = glue("{xlsx_prefix}/Monocytes/t_monocyte_cf_tables_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Monocytes/t_monocyte_cf_tables_sva-v202305.xlsx before writing the tables.

## Adding venn plots for outcome.

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

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Monocytes/t_monocyte_cf_sig_sva-v202305.xlsx before writing the tables.

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

## [1] 60 50

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

## [1] 53 50

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

## 
##    Tumaco_cure Tumaco_failure 
##             21             21 
## 
##  3  2  1 
## 13 13 16

t_cf_monocyte_tables_batchvisit <- combine_de_tables(
  t_cf_monocyte_de_batchvisit, keepers = cf_contrast,
#  rda = glue("rda/t_monocyte_cf_table_batchvisit-v{ver}.rda"),
  excel = glue("{xlsx_prefix}/Monocytes/t_monocyte_cf_tables_batchvisit-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Monocytes/t_monocyte_cf_tables_batchvisit-v202305.xlsx before writing the tables.
## Adding venn plots for outcome.

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

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Monocytes/t_monocyte_cf_sig_batchvisit-v202305.xlsx before writing the tables.

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

## [1] 43 50

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

## [1] 93 50

6.0.3.4 gProfiler: Monocytes

Now that I am looking back over these results, I am not compeltely certain why I only did the gprofiler search for the sva data…

Up: 60 genes; 12 GO, 1 KEGG, 1 WP, 4 TF, 0 others Down: 53 genes; 26 GO, 1 KEGG, 1 Reactome, 2 TF, 0 others

t_cf_monocyte_sig_sva_gp <- all_gprofiler(t_cf_monocyte_sig_sva)

enrichplot::dotplot(t_cf_monocyte_sig_sva_gp[["outcome_up"]][["GO_enrich"]])

enrichplot::dotplot(t_cf_monocyte_sig_sva_gp[["outcome_up"]][["TF_enrich"]])

enrichplot::dotplot(t_cf_monocyte_sig_sva_gp[["outcome_down"]][["GO_enrich"]])

t_cf_monocyte_sig_batch_gp <- all_gprofiler(t_cf_monocyte_sig_batchvisit)
enrichplot::dotplot(t_cf_monocyte_sig_batch_gp[["outcome_up"]][["GO_enrich"]])

enrichplot::dotplot(t_cf_monocyte_sig_batch_gp[["outcome_up"]][["HP_enrich"]])

6.0.4 Individual visits, Monocytes

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

6.0.4.1 Visit 1

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

## 
##    Tumaco_cure Tumaco_failure 
##              8              8

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

## Setting 187 low elements to zero.

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

t_cf_monocyte_v1_tables_sva <- combine_de_tables(
  t_cf_monocyte_v1_de_sva, keepers = cf_contrast,
#  rda = glue("rda/t_monocyte_v1_cf_table_sva-v{ver}.rda"),
  excel = glue("{xlsx_prefix}/Monocytes/t_monocyte_v1_cf_tables_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Monocytes/t_monocyte_v1_cf_tables_sva-v202305.xlsx before writing the tables.

## Adding venn plots for outcome.

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

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Monocytes/t_monocyte_v1_cf_sig_sva-v202305.xlsx before writing the tables.

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

## [1] 14 50

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

## [1] 52 50

6.0.4.2 Visit 2

t_cf_monocyte_v2_de_sva <- all_pairwise(tv2_monocytes, model_batch = "svaseq", filter = TRUE)

## 
##    Tumaco_cure Tumaco_failure 
##              7              6

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

## Setting 115 low elements to zero.

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

t_cf_monocyte_v2_tables_sva <- combine_de_tables(
  t_cf_monocyte_v2_de_sva, keepers = cf_contrast,
#  rda = glue("rda/t_monocyte_v2_cf_table_sva-v{ver}.rda"),
  excel = glue("{xlsx_prefix}/Monocytes/t_monocyte_v2_cf_tables_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Monocytes/t_monocyte_v2_cf_tables_sva-v202305.xlsx before writing the tables.

## Adding venn plots for outcome.

t_cf_monocyte_v2_sig_sva <- extract_significant_genes(
  t_cf_monocyte_v2_tables_sva,
  excel = glue("{xlsx_prefix}/Monocytes/t_monocyte_v2_cf_sig_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Monocytes/t_monocyte_v2_cf_sig_sva-v202305.xlsx before writing the tables.

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

## [1]  0 50

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

## [1]  1 50

6.0.4.3 Visit 3

t_cf_monocyte_v3_de_sva <- all_pairwise(tv3_monocytes, model_batch = "svaseq", filter = TRUE)

## 
##    Tumaco_cure Tumaco_failure 
##              6              7

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

## Setting 55 low elements to zero.

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

t_cf_monocyte_v3_tables_sva <- combine_de_tables(
  t_cf_monocyte_v3_de_sva, keepers = cf_contrast,
#  rda = glue("rda/t_monocyte_v3_cf_table_sva-v{ver}.rda"),
  excel = glue("{xlsx_prefix}/Monocytes/t_monocyte_v3_cf_tables_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Monocytes/t_monocyte_v3_cf_tables_sva-v202305.xlsx before writing the tables.

## Adding venn plots for outcome.

t_cf_monocyte_v3_sig_sva <- extract_significant_genes(
  t_cf_monocyte_v3_tables_sva,
  excel = glue("{xlsx_prefix}/Monocytes/t_monocyte_v3_cf_sig_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Monocytes/t_monocyte_v3_cf_sig_sva-v202305.xlsx before writing the tables.

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

## [1]  0 50

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

## [1]  4 50

6.0.4.4 Monocytes: Compare sva to batch-in-model

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

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

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

shared_ids <- rownames(t_cf_monocyte_tables_sva[["data"]][[1]]) %in%
  rownames(t_cf_monocyte_tables_batchvisit[["data"]][[1]])
first <- t_cf_monocyte_tables_sva[["data"]][[1]][shared_ids, ]
second <- t_cf_monocyte_tables_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 = 180, df = 10857, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8611 0.8705
## sample estimates:
##    cor 
## 0.8659

6.0.4.4.1 gProfiler: Monocytes by visit, V1

V1: Up: 14 genes; No categories V1: Down: 52 genes; 20 GO, 5 TF

t_cf_monocyte_v1_sig_sva_gp <- all_gprofiler(t_cf_monocyte_v1_sig_sva)

enrichplot::dotplot(t_cf_monocyte_v1_sig_sva_gp[["outcome_down"]][["GO_enrich"]])

6.0.4.4.2 gProfiler: Monocytes by visit, V2

V2: Up: 1 gene V2: Down: 0 genes.

6.0.4.4.3 gProfiler: Monocytes by visit, V3

V3: Up: 4 genes. V3: Down: 0 genes.

6.0.5 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", filter = TRUE)

## 
##    Tumaco_cure Tumaco_failure 
##             20             21

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

## Setting 750 low elements to zero.

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

t_cf_neutrophil_tables_sva <- combine_de_tables(
  t_cf_neutrophil_de_sva, keepers = cf_contrast,
#  rda = glue("rda/t_neutrophil_cf_table_sva-v{ver}.rda"),
  excel = glue("{xlsx_prefix}/Neutrophils/t_neutrophil_cf_tables_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Neutrophils/t_neutrophil_cf_tables_sva-v202305.xlsx before writing the tables.

## Adding venn plots for outcome.

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

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Neutrophils/t_neutrophil_cf_sig_sva-v202305.xlsx before writing the tables.

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

## [1] 84 50

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

## [1] 29 50

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

## 
##    Tumaco_cure Tumaco_failure 
##             20             21 
## 
##  3  2  1 
## 12 13 16

t_cf_neutrophil_tables_batchvisit <- combine_de_tables(
  t_cf_neutrophil_de_batchvisit, keepers = cf_contrast,
#  rda = glue("rda/t_neutrophil_cf_table_batchvisit-v{ver}.rda"),
  excel = glue("{xlsx_prefix}/Neutrophils/t_neutrophil_cf_tables_batchvisit-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Neutrophils/t_neutrophil_cf_tables_batchvisit-v202305.xlsx before writing the tables.
## Adding venn plots for outcome.

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

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Neutrophils/t_neutrophil_cf_sig_batchvisit-v202305.xlsx before writing the tables.

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

## [1] 92 50

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

## [1] 47 50

6.0.5.1 gProfiler: Neutrophils

Up: 84 genes; 5 GO, 2 Reactome, 3 TF, no others. Down: 29 genes: 12 GO, 1 Reactome, 1 TF, 1 miRNA, 11 HP, 0 others

t_cf_neutrophil_sig_sva_gp <- all_gprofiler(t_cf_neutrophil_sig_sva)

enrichplot::dotplot(t_cf_neutrophil_sig_sva_gp[["outcome_up"]][["GO_enrich"]])

enrichplot::dotplot(t_cf_neutrophil_sig_sva_gp[["outcome_up"]][["TF_enrich"]])

enrichplot::dotplot(t_cf_neutrophil_sig_sva_gp[["outcome_down"]][["GO_enrich"]])

enrichplot::dotplot(t_cf_neutrophil_sig_sva_gp[["outcome_down"]][["HP_enrich"]])

6.0.5.2 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",
                                              filter = TRUE)

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

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

## Setting 686 low elements to zero.

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

t_cf_neutrophil_visits_tables_sva <- combine_de_tables(
  t_cf_neutrophil_visits_de_sva, keepers = visitcf_contrasts,
#  rda = glue("rda/t_neutrophil_visitcf_table_sva-v{ver}.rda"),
  excel = glue("{xlsx_prefix}/Neutrophils/t_neutrophil_visitcf_tables_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Neutrophils/t_neutrophil_visitcf_tables_sva-v202305.xlsx before writing the tables.

## Adding venn plots for v1cf.

## Adding venn plots for v2cf.

## Adding venn plots for v3cf.

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

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Neutrophils/t_neutrophil_visitcf_sig_sva-v202305.xlsx before writing the tables.

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

## [1] 12 50

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

## [1]  6 50

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

## 
##    Tumaco_cure Tumaco_failure 
##              8              8

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

## Setting 145 low elements to zero.

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

t_cf_neutrophil_v1_tables_sva <- combine_de_tables(
  t_cf_neutrophil_v1_de_sva, keepers = cf_contrast,
#  rda = glue("rda/t_neutrophil_v1_cf_table_sva-v{ver}.rda"),
  excel = glue("{xlsx_prefix}/Neutrophils/t_neutrophil_v1_cf_tables_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Neutrophils/t_neutrophil_v1_cf_tables_sva-v202305.xlsx before writing the tables.

## Adding venn plots for outcome.

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

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Neutrophils/t_neutrophil_v1_cf_sig_sva-v202305.xlsx before writing the tables.

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

## [1]  5 50

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

## [1]  8 50

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

## 
##    Tumaco_cure Tumaco_failure 
##              7              6

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

## Setting 78 low elements to zero.

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

t_cf_neutrophil_v2_tables_sva <- combine_de_tables(
  t_cf_neutrophil_v2_de_sva,
  keepers = cf_contrast,
#  rda = glue("rda/t_neutrophil_v2_cf_table_sva-v{ver}.rda"),
  excel = glue("{xlsx_prefix}/Neutrophils/t_neutrophil_v2_cf_tables_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Neutrophils/t_neutrophil_v2_cf_tables_sva-v202305.xlsx before writing the tables.

## Adding venn plots for outcome.

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

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Neutrophils/t_neutrophil_v2_cf_sig_sva-v202305.xlsx before writing the tables.

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

## [1]  9 50

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

## [1]  3 50

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

## 
##    Tumaco_cure Tumaco_failure 
##              5              7

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

## Setting 83 low elements to zero.

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

t_cf_neutrophil_v3_tables_sva <- combine_de_tables(
  t_cf_neutrophil_v3_de_sva, keepers = cf_contrast,
#  rda = glue("rda/t_neutrophil_v3_cf_table_sva-v{ver}.rda"),
  excel = glue("{xlsx_prefix}/Neutrophils/t_neutrophil_v3_cf_tables_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Neutrophils/t_neutrophil_v3_cf_tables_sva-v202305.xlsx before writing the tables.

## Adding venn plots for outcome.

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

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Neutrophils/t_neutrophil_v3_cf_sig_sva-v202305.xlsx before writing the tables.

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

## [1]  5 50

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

## [1]  4 50

6.0.5.2.1 gProfiler: Neutrophils by visit, V1

V1: Up: 5 genes V1: Down: 8 genes; 14 GO.

t_cf_neutrophil_v1_sig_sva_gp <- all_gprofiler(t_cf_neutrophil_v1_sig_sva)

enrichplot::dotplot(t_cf_neutrophil_v1_sig_sva_gp[["outcome_down"]][["GO_enrich"]])

6.0.5.2.2 gProfiler: Neutrophils by visit, V2

Up: 5 genes; 3 GO, 10 TF. Down: 1 gene.

6.0.5.3 Neutrophils: Compare sva to batch-in-model

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

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

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

shared_ids <- rownames(t_cf_neutrophil_tables_sva[["data"]][[1]]) %in%
  rownames(t_cf_neutrophil_tables_batchvisit[["data"]][[1]])
first <- t_cf_neutrophil_tables_sva[["data"]][[1]][shared_ids, ]
second <- t_cf_neutrophil_tables_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 = 192, df = 9097, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8915 0.8996
## sample estimates:
##    cor 
## 0.8956

6.0.6 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", filter = TRUE)

## 
##    Tumaco_cure Tumaco_failure 
##             17              9

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

## Setting 325 low elements to zero.

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

t_cf_eosinophil_tables_sva <- combine_de_tables(
  t_cf_eosinophil_de_sva, keepers = cf_contrast,
#  rda = glue("rda/t_eosinophil_cf_table_sva-v{ver}.rda"),
  excel = glue("{xlsx_prefix}/Eosinophils/t_eosinophil_cf_tables_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Eosinophils/t_eosinophil_cf_tables_sva-v202305.xlsx before writing the tables.

## Adding venn plots for outcome.

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

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Eosinophils/t_eosinophil_cf_sig_sva-v202305.xlsx before writing the tables.

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

## [1] 116  50

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

## [1] 74 50

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

## 
##    Tumaco_cure Tumaco_failure 
##             17              9 
## 
## 3 2 1 
## 9 9 8

t_cf_eosinophil_tables_batchvisit <- combine_de_tables(
  t_cf_eosinophil_de_batchvisit, keepers = cf_contrast,
#  rda = glue("rda/t_eosinophil_cf_table_batchvisit-v{ver}.rda"),
  excel = glue("{xlsx_prefix}/Eosinophils/t_eosinophil_cf_tables_batchvisit-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Eosinophils/t_eosinophil_cf_tables_batchvisit-v202305.xlsx before writing the tables.
## Adding venn plots for outcome.

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

## Deleting the file excel/t_eosinophil_cf_sig_batchvisit-v202305.xlsx before writing the tables.

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

## [1] 99 50

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

## [1] 35 50

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",
                                              filter = TRUE)

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

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

## Setting 374 low elements to zero.

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

t_cf_eosinophil_visits_tables_sva <- combine_de_tables(
  t_cf_eosinophil_visits_de_sva, keepers = visitcf_contrasts,
#  rda = glue("rda/t_eosinophil_visitcf_table_sva-v{ver}.rda"),
  excel = glue("{xlsx_prefix}/Eosinophils/t_eosinophil_visitcf_tables_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Eosinophils/t_eosinophil_visitcf_tables_sva-v202305.xlsx before writing the tables.

## Adding venn plots for v1cf.

## Adding venn plots for v2cf.

## Adding venn plots for v3cf.

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

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Eosinophils/t_eosinophil_visitcf_sig_sva-v202305.xlsx before writing the tables.

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

## [1]  9 50

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

## [1] 11 50

6.0.6.1 C/F celltype volcano plots with specific labels

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_tables_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 = glue("images/cf_monocyte_volcano_labeled-v{ver}.svg"))
cf_monocyte_volcano$plot
dev.off()

## png 
##   2

cf_monocyte_volcano$plot

cf_eosinophil_table <- t_cf_eosinophil_tables_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 = glue("images/cf_eosinophil_volcano_labeled-v{ver}.svg"))
cf_eosinophil_volcano$plot
dev.off()

## png 
##   2

cf_eosinophil_volcano$plot

cf_neutrophil_table <- t_cf_neutrophil_tables_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 = glue("images/cf_neutrophil_volcano_labeled-v{ver}.svg"))
cf_neutrophil_volcano$plot
dev.off()

## png 
##   2

cf_neutrophil_volcano$plot

6.0.6.2 gProfiler: Eosinophils

Up: 116 genes; 123 GO, 2 KEGG, 7 Reactome, 5 WP, 69 TF, 1 miRNA, 0 others Down: 74 genes; 5 GO, 1 Reactome, 4 TF, 0 others

t_cf_eosinophil_sig_sva_gp <- all_gprofiler(t_cf_eosinophil_sig_sva)

enrichplot::dotplot(t_cf_eosinophil_sig_sva_gp[["outcome_up"]][["GO_enrich"]])

enrichplot::dotplot(t_cf_eosinophil_sig_sva_gp[["outcome_up"]][["REAC_enrich"]])

enrichplot::dotplot(t_cf_eosinophil_sig_sva_gp[["outcome_up"]][["WP_enrich"]])

enrichplot::dotplot(t_cf_eosinophil_sig_sva_gp[["outcome_up"]][["TF_enrich"]])

enrichplot::dotplot(t_cf_eosinophil_sig_sva_gp[["outcome_down"]][["GO_enrich"]])

enrichplot::dotplot(t_cf_eosinophil_sig_sva_gp[["outcome_down"]][["TF_enrich"]])

6.0.7 Eosinophil time comparisons

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

## 
##    Tumaco_cure Tumaco_failure 
##              5              3

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

## Setting 57 low elements to zero.

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

t_cf_eosinophil_v1_tables_sva <- combine_de_tables(
  t_cf_eosinophil_v1_de_sva, keepers = cf_contrast,
#  rda = glue("rda/t_eosinophil_v1_cf_table_sva-v{ver}.rda"),
  excel = glue("{xlsx_prefix}/Eosinophils/t_eosinophil_v1_cf_tables_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Eosinophils/t_eosinophil_v1_cf_tables_sva-v202305.xlsx before writing the tables.

## Adding venn plots for outcome.

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

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Eosinophils/t_eosinophil_v1_cf_sig_sva-v202305.xlsx before writing the tables.

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

## [1] 13 50

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

## [1] 19 50

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

## 
##    Tumaco_cure Tumaco_failure 
##              6              3

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

## Setting 90 low elements to zero.

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

t_cf_eosinophil_v2_tables_sva <- combine_de_tables(
  t_cf_eosinophil_v2_de_sva, keepers = cf_contrast,
#  rda = glue("rda/t_eosinophil_v2_cf_table_sva-v{ver}.rda"),
  excel = glue("{xlsx_prefix}/Eosinophils/t_eosinophil_v2_cf_tables_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Eosinophils/t_eosinophil_v2_cf_tables_sva-v202305.xlsx before writing the tables.

## Adding venn plots for outcome.

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

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Eosinophils/t_eosinophil_v2_cf_sig_sva-v202305.xlsx before writing the tables.

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

## [1]  9 50

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

## [1]  4 50

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

## 
##    Tumaco_cure Tumaco_failure 
##              6              3

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

## Setting 48 low elements to zero.

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

t_cf_eosinophil_v3_tables_sva <- combine_de_tables(
  t_cf_eosinophil_v3_de_sva, keepers = cf_contrast,
#  rda = glue("rda/t_eosinophil_v3_cf_table_sva-v{ver}.rda"),
  excel = glue("{xlsx_prefix}/Eosinophils/t_eosinophil_v3_cf_tables_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Eosinophils/t_eosinophil_v3_cf_tables_sva-v202305.xlsx before writing the tables.

## Adding venn plots for outcome.

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

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Eosinophils/t_eosinophil_v3_cf_sig_sva-v202305.xlsx before writing the tables.

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

## [1] 68 50

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

## [1] 29 50

6.0.7.1 gProfiler: Eosinophils V1

Up: 13 genes, no hits. Down: 19 genes; 11 GO, 1 Reactome, 1 TF

t_cf_eosinophil_v1_sig_sva_gp <- all_gprofiler(t_cf_eosinophil_v1_sig_sva)

enrichplot::dotplot(t_cf_eosinophil_sig_sva_gp[["outcome_down"]][["GO_enrich"]])

enrichplot::dotplot(t_cf_eosinophil_sig_sva_gp[["outcome_down"]][["TF_enrich"]])

6.0.7.2 gProfiler: Eosinophils V2

Up: 9 genes; 23 GO, 2 KEGG, 2 Reactome, 4 WP Down: 4 genes; no hits

t_cf_eosinophil_v2_sig_sva_gp <- all_gprofiler(t_cf_eosinophil_v2_sig_sva)

enrichplot::dotplot(t_cf_eosinophil_sig_sva_gp[["outcome_up"]][["GO_enrich"]])

enrichplot::dotplot(t_cf_eosinophil_sig_sva_gp[["outcome_up"]][["WP_enrich"]])

6.0.7.3 gProfiler: Eosinophils V3

Up: 68 genes; 95 GO, 2 KEGG, 12 Reactome, 3 WP, 63 TF, 1 miRNA Down: 29 genes; 3 GO, 1 WP, 1 TF, 3 miRNA

t_cf_eosinophil_v3_sig_sva_gp <- all_gprofiler(t_cf_eosinophil_v3_sig_sva)

enrichplot::dotplot(t_cf_eosinophil_sig_sva_gp[["outcome_up"]][["GO_enrich"]])

enrichplot::dotplot(t_cf_eosinophil_sig_sva_gp[["outcome_up"]][["WP_enrich"]])

6.0.7.4 Eosinophils: Compare sva to batch-in-visit

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

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

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

shared_ids <- rownames(t_cf_eosinophil_tables_sva[["data"]][[1]]) %in%
  rownames(t_cf_eosinophil_tables_batchvisit[["data"]][[1]])
first <- t_cf_eosinophil_tables_sva[["data"]][[1]][shared_ids, ]
second <- t_cf_eosinophil_tables_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 = 10528, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8232 0.8352
## sample estimates:
##    cor 
## 0.8293

6.0.7.5 Compare monocyte CF, neutrophil CF, eosinophil CF

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

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

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

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

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

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

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

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

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

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

6.0.8.1 Cure/Fail by visits, all cell types

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

## 
##    v1cure v1failure    v2cure v2failure    v3cure v3failure 
##        30        24        20        15        17        17

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

## Setting 17117 low elements to zero.

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

t_visit_cf_all_tables_sva <- combine_de_tables(
  t_visit_cf_all_de_sva, keepers = visitcf_contrasts,
#  rda = glue("rda/t_all_visitcf_table_sva-v{ver}.rda"),
  excel = glue("analyses/4_tumaco/DE_Cure_vs_Fail/t_all_visitcf_tables_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/t_all_visitcf_tables_sva-v202305.xlsx before writing the tables.

## Adding venn plots for v1cf.

## Adding venn plots for v2cf.

## Adding venn plots for v3cf.

t_visit_cf_all_sig_sva <- extract_significant_genes(
  t_visit_cf_all_tables_sva,
  excel = glue("analyses/4_tumaco/DE_Cure_vs_Fail/t_all_visitcf_sig_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/t_all_visitcf_sig_sva-v202305.xlsx before writing the tables.

t_visit_cf_all_gp <- all_gprofiler(t_visit_cf_all_sig_sva)

6.0.8.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",
                                           filter = TRUE)

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

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

## Setting 688 low elements to zero.

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

t_visit_cf_monocyte_tables_sva <- combine_de_tables(
  t_visit_cf_monocyte_de_sva, keepers = visitcf_contrasts,
#  rda = glue("rda/t_monocyte_visitcf_table_sva-v{ver}.rda"),
  excel = glue("{xlsx_prefix}/Monocytes/t_monocyte_visitcf_tables_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Monocytes/t_monocyte_visitcf_tables_sva-v202305.xlsx before writing the tables.

## Adding venn plots for v1cf.

## Adding venn plots for v2cf.

## Adding venn plots for v3cf.

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

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Monocytes/t_monocyte_visitcf_sig_sva-v202305.xlsx before writing the tables.

t_v1fc_deseq_ma <- t_visit_cf_monocyte_tables_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_tables_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_tables_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_tables_sva[["data"]][["v1cf"]]
v2cf <- t_visit_cf_monocyte_tables_sva[["data"]][["v2cf"]]
v3cf <- t_visit_cf_monocyte_tables_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"]]

6.0.8.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",
                                             filter = TRUE)

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

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

## Setting 686 low elements to zero.

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

t_visit_cf_neutrophil_tables_sva <- combine_de_tables(
  t_visit_cf_neutrophil_de_sva, keepers = visitcf_contrasts,
#  rda = glue("rda/t_neutrophil_visitcf_table_sva-v{ver}.rda"),
  excel = glue("{xlsx_prefix}/Neutrophils/t_neutrophil_visitcf_tables_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Neutrophils/t_neutrophil_visitcf_tables_sva-v202305.xlsx before writing the tables.

## Adding venn plots for v1cf.

## Adding venn plots for v2cf.

## Adding venn plots for v3cf.

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

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Neutrophils/t_neutrophil_visitcf_sig_sva-v202305.xlsx before writing the tables.

6.0.8.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",
                                             filter = TRUE)

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

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

## Setting 374 low elements to zero.

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

t_visit_cf_eosinophil_tables_sva <- combine_de_tables(
  t_visit_cf_eosinophil_de_sva, keepers = visitcf_contrasts,
#  rda = glue("rda/t_eosinophil_visitcf_table_sva-v{ver}.rda"),
  excel = glue("{xlsx_prefix}/Eosinophils/t_eosinophil_visitcf_tables_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Eosinophils/t_eosinophil_visitcf_tables_sva-v202305.xlsx before writing the tables.

## Adding venn plots for v1cf.

## Adding venn plots for v2cf.

## Adding venn plots for v3cf.

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

## Deleting the file analyses/4_tumaco/DE_Cure_vs_Fail/Eosinophils/t_eosinophil_visitcf_sig_sva-v202305.xlsx before writing the tables.

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

6.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 97 samples.

## subset_expt(): There were 97, 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.

6.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 8537 low-count genes (11386 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 8537 low-count genes (11386 remaining).

## Setting 1538 low elements to zero.

## transform_counts: Found 1538 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 9597 low-count genes (10326 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 9597 low-count genes (10326 remaining).

## Setting 46 low elements to zero.

## transform_counts: Found 46 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 11531 low-count genes (8392 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 11531 low-count genes (8392 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 9895 low-count genes (10028 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 9895 low-count genes (10028 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

6.1.3 persistence DE

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

## 
##  N  Y 
##  6 24

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

## Setting 1538 low elements to zero.

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

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

## Deleting the file analyses/4_tumaco/DE_Persistence/persistence_all_de_sva-v202305.xlsx before writing the tables.

## Adding venn plots for Y_vs_N.

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

## 
##  N  Y 
##  2 10

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

## Setting 46 low elements to zero.

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

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

## Deleting the file analyses/4_tumaco/DE_Persistence/persistence_monocyte_de_sva-v202305.xlsx before writing the tables.

## Adding venn plots for Y_vs_N.

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

## 
## N Y 
## 3 7

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

## Setting 46 low elements to zero.

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

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

## Deleting the file analyses/4_tumaco/DE_Persistence/persistence_neutrophil_de_sva-v202305.xlsx before writing the tables.

## Adding venn plots for Y_vs_N.

persistence_eosinophil_de_sva <- all_pairwise(persistence_eosinophil, filter = TRUE, model_batch = "svaseq")

## 
## N Y 
## 1 7

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

## Setting 25 low elements to zero.

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

persistence_eosinophil_table_sva <- combine_de_tables(
  persistence_eosinophil_de_sva,
  excel = glue("analyses/4_tumaco/DE_Persistence/persistence_eosinophil_de_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Persistence/persistence_eosinophil_de_sva-v202305.xlsx before writing the tables.

## Adding venn plots for Y_vs_N.

6.2 Comparing visits without regard to cure/fail

6.2.1 All cell types

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

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

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

## Setting 9614 low elements to zero.

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

t_visit_all_table_sva <- combine_de_tables(
  t_visit_all_de_sva, keepers = visit_contrasts,
#  rda = glue("rda/t_all_visit_table_sva-v{ver}.rda"),
  excel = glue("analyses/4_tumaco/DE_Visits/t_all_visit_tables_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Visits/t_all_visit_tables_sva-v202305.xlsx before writing the tables.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Adding venn plots for v2v1.

## Adding venn plots for v3v1.

## Adding venn plots for v3v2.

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

## Deleting the file analyses/4_tumaco/DE_Visits/t_all_visit_sig_sva-v202305.xlsx before writing the tables.

6.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, model_batch = "svaseq")

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

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

## Setting 648 low elements to zero.

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

t_visit_monocyte_table_sva <- combine_de_tables(
  t_visit_monocyte_de_sva, keepers = visit_contrasts,
#  rda = glue("rda/t_monocyte_visit_table_sva-v{ver}.rda"),
  excel = glue("analyses/4_tumaco/DE_Visits/Monocytes/t_monocyte_visit_tables_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Visits/Monocytes/t_monocyte_visit_tables_sva-v202305.xlsx before writing the tables.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Adding venn plots for v2v1.

## Adding venn plots for v3v1.

## Adding venn plots for v3v2.

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

## Deleting the file analyses/4_tumaco/DE_Visits/Monocytes/t_monocyte_visit_sig_sva-v202305.xlsx before writing the tables.

6.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, model_batch = "svaseq")

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

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

## Setting 589 low elements to zero.

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

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

## Deleting the file analyses/4_tumaco/DE_Visits/Neutrophils/t_neutrophil_visit_table_sva-v202305.xlsx before writing the tables.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Adding venn plots for v2v1.

## Adding venn plots for v3v1.

## Adding venn plots for v3v2.

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

## Deleting the file analyses/4_tumaco/DE_Visits/Neutrophils/t_neutrophil_visit_sig_sva-v202305.xlsx before writing the tables.

6.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, model_batch = "svaseq")

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

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

## Setting 272 low elements to zero.

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

t_visit_eosinophil_table <- combine_de_tables(
  t_visit_eosinophil_de, keepers = visit_contrasts,
#  rda = glue("rda/t_eosinophil_visit_table_sva-v{ver}.rda"),
  excel = glue("analyses/4_tumaco/DE_Visits/Eosinophils/t_eosinophil_visit_table_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Visits/Eosinophils/t_eosinophil_visit_table_sva-v202305.xlsx before writing the tables.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.

## Adding venn plots for v2v1.

## Adding venn plots for v3v1.

## Adding venn plots for v3v2.

t_visit_eosinophil_sig <- extract_significant_genes(
  t_visit_eosinophil_table,
  excel = glue("analyses/4_tumaco/DE_Visits/Eosinophils/t_eosinophil_visit_sig_sva-v{ver}.xlsx"))

## Deleting the file analyses/4_tumaco/DE_Visits/Eosinophils/t_eosinophil_visit_sig_sva-v202305.xlsx before writing the tables.

7 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_tables_sva:

ENSG00000198178, ENSG00000179344, ENSG00000182628

eo_rpkm <- normalize_expt(tv1_eosinophils, convert = "rpkm", column = "cds_length")

## There appear to be 5391 genes without a length.

test <- all_pairwise(tmrc3_external, model_batch = "svaseq", filter = "simple")

## 
##   Brazil Colombia 
##       21       18

## Removing 4594 low-count genes (14576 remaining).

## Setting 3707 low elements to zero.

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

test_table <- combine_de_tables(test, excel = "excel/tmrc3_scott_biopsies.xlsx")

## Adding venn plots for Colombia_vs_Brazil.

test_sig <- extract_significant_genes(test_table, excel = "excel/tmrc3_scott_biopsies_sig.xlsx")

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

## Skipping this because it is taking too long.
##if (!isTRUE(get0("skip_load"))) {
##  pander::pander(sessionInfo())
##  message(paste0("This is hpgltools commit: ", get_git_commit()))
##  message(paste0("Saving to ", savefile))
##  tmp <- sm(saveme(filename=savefile))
##}

tmp <- loadme(filename = savefile)