TMRC3 202412: Differential Expression analyses

1 Changelog

• 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 will have to trust me when I say those blocks all work?
• 202309: Moved all GSEA analyses to 04lrt_gsea_gsva.Rmd
• 202309 next day: Moving GSEA 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.

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.

2.3 GSEA

The GSEA analyses will follow each DE analysis during this document.

Most (all?) of the GSEA analyses used in this paper were done via gProfiler rather than goseq/clusterProfiler/topGO/GOstats. Primarily because it is so easy to invoke gprofiler.

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("tumaco_failure", "tumaco_cure"),
  "cure" = c("tumaco_cure", "cali_cure"))
## In other cases, we have cure/fail for both places.
clinic_cf_contrasts <- list(
  "cali" = c("cali_failure", "cali_cure"),
  "tumaco" = c("tumaco_failure", "tumaco_cure"),
  "cure" = c("tumaco_cure", "cali_cure"),
  "fail" = c("tumaco_failure", "cali_failure"))
cf_contrast <- list(
  "outcome" = c("tumaco_failure", "tumaco_cure"))
t_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"))

3 Compare samples by clinic

3.1 DE: Compare clinics, all samples

Perform a svaseq-guided comparison of the two clinics. Ideally this will give some clue about just how strong the clinic-based batch effect really is and what its causes are.

tc_clinic_type <- tc_valid %>%
  set_expt_conditions(fact = "clinic") %>%
  set_expt_batches(fact = "typeofcells")

## The numbers of samples by condition are:

## 
##   cali tumaco 
##     61    123

## The number of samples by batch are:

## 
##      biopsy eosinophils   monocytes neutrophils 
##          18          41          63          62

table(pData(tc_clinic_type)[["condition"]])

## 
##   cali tumaco 
##     61    123

tc_all_clinic_de_sva <- all_pairwise(tc_clinic_type, model_batch = "svaseq",
                                     filter = TRUE, parallel = parallel, methods = methods)

## 
##   cali tumaco 
##     61    123

## Error in checkForRemoteErrors(val): one node produced an error: c("Error in all_adjusters(input, estimate_type = model_batch, surrogates = surrogates) : \n  If an expt is not passed, then design _must_ be.\n", "noiseq")

tc_all_clinic_de_sva

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

tc_all_clinic_de_sva[["deseq"]][["contrasts_performed"]]

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

tc_all_clinic_table_sva <- combine_de_tables(
  tc_all_clinic_de_sva, keepers = clinic_contrasts,
  excel = glue("{clinic_prefix}/tc_all_clinic_table_sva-v{ver}.xlsx"))

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

tc_all_clinic_table_sva

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

tc_all_clinic_sig_sva <- extract_significant_genes(
  tc_all_clinic_table_sva,
  excel = glue("{clinic_prefix}/compare_clinics/tc_clinic_type_sig_sva-v{ver}.xlsx"))

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

tc_all_clinic_sig_sva

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

3.1.1 GSEA: comparing the clinics

increased_tumaco_categories_up <- simple_gprofiler(
  tc_all_clinic_sig_sva[["deseq"]][["ups"]][["clinics"]],
  excel = glue("{gsea_prefix}/tumaco_cateogies_up-v{ver}.xlsx"))

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

increased_tumaco_categories_up

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

increased_tumaco_categories_up[["pvalue_plots"]][["BP"]]

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

increased_cali_categories <- simple_gprofiler(
  tc_all_clinic_sig_sva[["deseq"]][["downs"]][["clinics"]],
  excel = glue("{gsea_prefix}/cali_cateogies_up-v{ver}.xlsx"))

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

increased_cali_categories

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

increased_cali_categories[["pvalue_plots"]][["BP"]]

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

3.1.2 Visualize clinic differences

Let us take a quick look at the results of the comparison of Tumaco/Cali

Note: I keep re-introducing an error which causes these (volcano and MA) plots to be reversed with respect to the logFC values. Pay careful attention to these and make sure that they agree with the numbers of genes observed in the contrast.

## Check that up is up
summary(tc_all_clinic_table_sva[["data"]][["clinics"]][["deseq_logfc"]])

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

## I think we can assume that most genes are down when considering Tumaco/Cali.
sum(tc_all_clinic_table_sva$data$clinics$deseq_logfc < -1.0 &
      tc_all_clinic_table_sva$data$clinics$deseq_adjp < 0.05)

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

tc_all_clinic_table_sva[["plots"]][["clinics"]][["deseq_vol_plots"]]

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

## Ok, so it says 1794 up, but that is clearly the down side...  Something is definitely messed up.
## The points are on the correct sides of the plot, but the categories of up/down are reversed.
## Theresa noted that she colors differently, and I think better: left side gets called
## 'increased in denominator', right side gets called 'increased in numerator';
## these two groups are colored according to their condition colors, and everything else is gray.
## I am checking out Theresa's helper_functions.R to get a sense of how she handles this, I think
## I can use a variant of her idea pretty easily:
##  1.  Add a column 'Significance', which is a factor, and contains either 'Not enriched',
##      'Enriched in x', or 'Enriched in y' according to the logfc/adjp.
##  2.  use the significance column for the geom_point color/fill in the volcano plot.
## My change to this idea would be to extract the colors from the input expressionset.

There appear to be many more genes which are increased in the Tumaco samples with respect to the Cali samples.

3.2 DE: Compare clinics, eosinophil samples

The remaining cell types all have pretty strong clinic-based variance; but I am not certain if it is consistent across cell types.

table(pData(tc_eosinophils)[["condition"]])

## 
##      cali_cure    tumaco_cure tumaco_failure 
##             15             17              9

tc_eosinophils_clinic_de_nobatch <- all_pairwise(tc_eosinophils, parallel = parallel,
                                                 model_batch = FALSE, filter = TRUE,
                                                 methods = methods)

## 
##      cali_cure    tumaco_cure tumaco_failure 
##             15             17              9

## Error in all_pairwise(tc_eosinophils, parallel = parallel, model_batch = FALSE, : object 'sv_model' not found

tc_eosinophils_clinic_de_nobatch

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

tc_eosinophils_clinic_de_nobatch[["deseq"]][["contrasts_performed"]]

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

tc_eosinophils_clinic_table_nobatch <- combine_de_tables(
  tc_eosinophils_clinic_de_nobatch, keepers = tc_cf_contrasts,
  excel = glue("{clinic_cf_prefix}/Eosinophils/tc_eosinophils_clinic_table_nobatch-v{ver}.xlsx"))

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

tc_eosinophils_clinic_table_nobatch

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

tc_eosinophils_clinic_sig_nobatch <- extract_significant_genes(
  tc_eosinophils_clinic_table_nobatch,
  excel = glue("{clinic_cf_prefix}/Eosinophils/tc_eosinophils_clinic_sig_nobatch-v{ver}.xlsx"))

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

tc_eosinophils_clinic_sig_nobatch

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

tc_eosinophils_clinic_de_sva <- all_pairwise(tc_eosinophils, model_batch = "svaseq",
                                             parallel = parallel, filter = TRUE, methods = methods)

## 
##      cali_cure    tumaco_cure tumaco_failure 
##             15             17              9

## Error in checkForRemoteErrors(val): one node produced an error: c("Error in all_adjusters(input, estimate_type = model_batch, surrogates = surrogates) : \n  If an expt is not passed, then design _must_ be.\n", "noiseq")

tc_eosinophils_clinic_de_sva

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

tc_eosinophils_clinic_de_sva[["deseq"]][["contrasts_performed"]]

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

tc_eosinophils_clinic_table_sva <- combine_de_tables(
  tc_eosinophils_clinic_de_sva, keepers = tc_cf_contrasts,
  excel = glue("{clinic_cf_prefix}/Eosinophils/tc_eosinophils_clinic_table_sva-v{ver}.xlsx"))

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

tc_eosinophils_clinic_table_sva

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

tc_eosinophils_clinic_sig_sva <- extract_significant_genes(
  tc_eosinophils_clinic_table_sva,
  excel = glue("{clinic_cf_prefix}/Eosinophils/tc_eosinophils_clinic_sig_sva-v{ver}.xlsx"))

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

tc_eosinophils_clinic_sig_sva

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

3.3 DE: Compare clinics, biopsy samples

Interestingly to me, the biopsy samples appear to have the least location-based variance. But we can perform an explicit DE and see how well that hypothesis holds up.

Note that these data include cure and fail samples for

table(pData(tc_biopsies)[["condition"]])

## 
##      cali_cure    tumaco_cure tumaco_failure 
##              4              9              5

tc_biopsies_clinic_de_sva <- all_pairwise(tc_biopsies, parallel = parallel,
                                          model_batch = "svaseq", filter = TRUE,
                                          methods = methods)

## 
##      cali_cure    tumaco_cure tumaco_failure 
##              4              9              5

## Error in checkForRemoteErrors(val): one node produced an error: c("Error in all_adjusters(input, estimate_type = model_batch, surrogates = surrogates) : \n  If an expt is not passed, then design _must_ be.\n", "noiseq")

tc_biopsies_clinic_de_sva

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

tc_biopsies_clinic_de_sva[["deseq"]][["contrasts_performed"]]

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

tc_biopsies_clinic_table_sva <- combine_de_tables(
  tc_biopsies_clinic_de_sva, keepers = tc_cf_contrasts,
  excel = glue("{clinic_cf_prefix}/Biopsies/tc_biopsies_clinic_table_sva-v{ver}.xlsx"))

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

tc_biopsies_clinic_table_sva

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

tc_biopsies_clinic_sig_sva <- extract_significant_genes(
  tc_biopsies_clinic_table_sva,
  excel = glue("{clinic_cf_prefix}/Biopsies/tc_biopsies_clinic_sig_sva-v{ver}.xlsx"))

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

tc_biopsies_clinic_sig_sva

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

3.4 DE: Compare clinics, monocyte samples

At least for the moment, I am only looking at the differences between no-batch vs. sva across clinics for the monocyte samples. This was chosen mostly arbitrarily.

3.4.1 DE: Compare clinics, monocytes without batch estimation

Our baseline is the comparison of the monocytes samples without batch in the model or surrogate estimation. In theory at least, this should correspond to the PCA plot above when no batch estimation was performed.

table(pData(tc_monocytes)[["condition"]])

## 
##      cali_cure   cali_failure    tumaco_cure tumaco_failure 
##             18              3             21             21

tc_monocytes_de_nobatch <- all_pairwise(tc_monocytes, model_batch = FALSE,
                                        parallel = parallel, filter = TRUE,
                                        methods = methods)

## 
##      cali_cure   cali_failure    tumaco_cure tumaco_failure 
##             18              3             21             21

## Error in all_pairwise(tc_monocytes, model_batch = FALSE, parallel = parallel, : object 'sv_model' not found

tc_monocytes_de_nobatch

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

tc_monocytes_table_nobatch <- combine_de_tables(
  tc_monocytes_de_nobatch, keepers = clinic_cf_contrasts,
  excel = glue("{clinic_cf_prefix}/Monocytes/tc_monocytes_clinic_table_nobatch-v{ver}.xlsx"))

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

tc_monocytes_table_nobatch

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

tc_monocytes_sig_nobatch <- extract_significant_genes(
  tc_monocytes_table_nobatch,
  excel = glue("{clinic_cf_prefix}/Monocytes/tc_monocytes_clinic_sig_nobatch-v{ver}.xlsx"))

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

tc_monocytes_sig_nobatch

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

3.4.2 DE: Compare clinics, monocytes with svaseq

In contrast, the following comparison should give a view of the data corresponding to the svaseq PCA plot above. In the best case scenario, we should therefore be able to see some significane differences between the Tumaco cure and fail samples.

tc_monocytes_de_sva <- all_pairwise(tc_monocytes, model_batch = "svaseq",
                                    parallel = parallel, filter = TRUE,
                                    methods = methods)

## 
##      cali_cure   cali_failure    tumaco_cure tumaco_failure 
##             18              3             21             21

## Error in checkForRemoteErrors(val): one node produced an error: c("Error in all_adjusters(input, estimate_type = model_batch, surrogates = surrogates) : \n  If an expt is not passed, then design _must_ be.\n", "noiseq")

tc_monocytes_de_sva

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

tc_monocytes_table_sva <- combine_de_tables(
  tc_monocytes_de_sva, keepers = clinic_cf_contrasts,
  excel = glue("{clinic_cf_prefix}/Monocytes/tc_monocytes_clinic_table_sva-v{ver}.xlsx"))

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

tc_monocytes_table_sva

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

tc_monocytes_sig_sva <- extract_significant_genes(
  tc_monocytes_table_sva,
  excel = glue("{clinic_cf_prefix}/Monocytes/tc_monocytes_clinic_sig_sva-v{ver}.xlsx"))

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

tc_monocytes_sig_sva

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

3.4.3 DE Compare: How similar are the no-batch vs. SVA results?

The following block shows that these two results are exceedingly different, sugesting that the Cali cure/fail and Tumaco cure/fail cannot easily be considered in the same analysis. I did some playing around with my calculate_aucc function in this block and found that it is in some important way broken, at least if one expands the top-n genes to more than 20% of the number of genes in the data.

cali_table <- tc_monocytes_table_nobatch[["data"]][["cali"]]

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

table <- tc_monocytes_table_nobatch[["data"]][["tumaco"]]

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

cali_aucc <- calculate_aucc(cali_table, table, px = "deseq_adjp", py = "deseq_adjp",
                            lx = "deseq_logfc", ly = "deseq_logfc")

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

cali_aucc

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

cali_table_sva <- tc_monocytes_table_sva[["data"]][["cali"]]

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

tumaco_table_sva <- tc_monocytes_table_sva[["data"]][["tumaco"]]

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

cali_aucc_sva <- calculate_aucc(cali_table_sva, tumaco_table_sva, px = "deseq_adjp",
                                py = "deseq_adjp", lx = "deseq_logfc", ly = "deseq_logfc")

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

cali_aucc_sva

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

3.5 DE: Compare clinics, neutrophil samples

tc_neutrophils_de_nobatch <- all_pairwise(tc_neutrophils, parallel = parallel,
                                          model_batch = FALSE, filter = TRUE,
                                          methods = methods)

## 
##      cali_cure   cali_failure    tumaco_cure tumaco_failure 
##             18              3             20             21

## Error in all_pairwise(tc_neutrophils, parallel = parallel, model_batch = FALSE, : object 'sv_model' not found

tc_neutrophils_de_nobatch

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

tc_neutrophils_table_nobatch <- combine_de_tables(
  tc_neutrophils_de_nobatch, keepers = clinic_cf_contrasts,
  excel = glue("{clinic_cf_prefix}/Neutrophils/tc_neutrophils_table_nobatch-v{ver}.xlsx"))

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

tc_neutrophils_table_nobatch

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

tc_neutrophils_sig_nobatch <- extract_significant_genes(
  tc_neutrophils_table_nobatch,
  excel = glue("{clinic_cf_prefix}/Neutrophils/tc_neutrophils_sig_nobatch-v{ver}.xlsx"))

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

tc_neutrophils_sig_nobatch

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

tc_neutrophils_de_sva <- all_pairwise(tc_neutrophils, parallel = parallel,
                                      model_batch = "svaseq", filter = TRUE,
                                      methods = methods)

## 
##      cali_cure   cali_failure    tumaco_cure tumaco_failure 
##             18              3             20             21

## Error in checkForRemoteErrors(val): one node produced an error: c("Error in all_adjusters(input, estimate_type = model_batch, surrogates = surrogates) : \n  If an expt is not passed, then design _must_ be.\n", "noiseq")

tc_neutrophils_de_sva

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

tc_neutrophils_table_sva <- combine_de_tables(
  tc_neutrophils_de_sva, keepers = clinic_cf_contrasts,
  excel = glue("{clinic_cf_prefix}/Neutrophils/tc_neutrophils_table_sva-v{ver}.xlsx"))

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

tc_neutrophils_table_sva

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

tc_neutrophils_sig_sva <- extract_significant_genes(
  tc_neutrophils_table_sva,
  excel = glue("{clinic_cf_prefix}/Neutrophils/tc_neutrophils_sig_sva-v{ver}.xlsx"))

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

tc_neutrophils_sig_sva

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

3.6 GSEA: Extract clinic-specific genes

Given the above comparisons, we can extract some gene sets which resulted from those DE analyses and eventually perform some ontology/KEGG/reactome/etc searches. This reminds me, I want to make my extract_significant_ functions to return gene-set data structures and my various ontology searches to take them as inputs. This should help avoid potential errors when extracting up/down genes.

clinic_sigenes_up <- rownames(tc_all_clinic_sig_sva[["deseq"]][["ups"]][["clinics"]])

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

clinic_sigenes_down <- rownames(tc_all_clinic_sig_sva[["deseq"]][["downs"]][["clinics"]])

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

clinic_sigenes <- c(clinic_sigenes_up, clinic_sigenes_down)

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

tc_eosinophils_sigenes_up <- rownames(tc_eosinophils_clinic_sig_sva[["deseq"]][["ups"]][["cure"]])

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

tc_eosinophils_sigenes_down <- rownames(tc_eosinophils_clinic_sig_sva[["deseq"]][["downs"]][["cure"]])

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

tc_monocytes_sigenes_up <- rownames(tc_monocytes_sig_sva[["deseq"]][["ups"]][["cure"]])

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

tc_monocytes_sigenes_down <- rownames(tc_monocytes_sig_sva[["deseq"]][["downs"]][["cure"]])

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

tc_neutrophils_sigenes_up <- rownames(tc_neutrophils_sig_sva[["deseq"]][["ups"]][["cure"]])

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

tc_neutrophils_sigenes_down <- rownames(tc_neutrophils_sig_sva[["deseq"]][["downs"]][["cure"]])

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

tc_eosinophils_sigenes <- c(tc_eosinophils_sigenes_up,
                            tc_eosinophils_sigenes_down)

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

tc_monocytes_sigenes <- c(tc_monocytes_sigenes_up,
                          tc_monocytes_sigenes_down)

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

tc_neutrophils_sigenes <- c(tc_neutrophils_sigenes_up,
                            tc_neutrophils_sigenes_down)

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

3.7 GSEA: gProfiler of genes deemed up/down when comparing Cali and Tumaco

I was curious to try to understand why the two clinics appear to be so different vis a vis their PCA/DE; so I thought that gProfiler might help boil those results down to something more digestible.

3.7.1 GSEA: Compare clinics, all samples

Note that in the following block I used the function simple_gprofiler(), but later in this document I will use all_gprofiler(). The first invocation limits the search to a single table, while the second will iterate over every result in a pairwise differential expression analysis.

In this instance, we are looking at the vector of gene IDs deemed significantly different between the two clinics in either the up or down direction.

One other thing worth noting, the new version of gProfiler provides some fun interactive plots. I will add an example here.

tc_eosinophil_gprofiler <- simple_gprofiler(
  tc_eosinophils_sigenes_up,
  excel = glue("{gsea_prefix}/eosinophil_clinics_tumaco_up-v{ver}.xlsx"))

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

tc_eosinophil_gprofiler

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

clinic_gp <- simple_gprofiler(
  clinic_sigenes,
  excel = glue("{gsea_prefix}/both_clinics_cali_up-v{ver}.xlsx"))

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

clinic_gp$pvalue_plots$REAC

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

clinic_gp$pvalue_plots$BP

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

clinic_gp$pvalue_plots$TF

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

clinic_gp$interactive_plots$GO

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

3.7.2 GSEA: Compare clinics, Eosinophil samples

In the following block, I am looking at the gProfiler over represented groups observed across clinics in only the Eosinophils. First I do so for all genes(up or down), followed by only the up and down groups. Each of the following will include only the Reactome and GO:BP plots. These searches did not have too many other hits, excepting the transcription factor database.

tc_eosinophils_gp <- simple_gprofiler(
  tc_eosinophils_sigenes,
  excel = glue("{gsea_prefix}/eosinophil_clinics-v{ver}.xlsx"))

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

tc_eosinophils_gp

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

tc_eosinophils_gp$pvalue_plots$REAC

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

tc_eosinophils_gp$pvalue_plots$BP

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

tc_eosinophils_up_gp <- simple_gprofiler(
  tc_eosinophils_sigenes_up,
  excel = glue("{gsea_prefix}/eosinophil_clinics_tumaco_up-v{ver}.xlsx"))

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

tc_eosinophils_up_gp

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

tc_eosinophils_up_gp$pvalue_plots$REAC

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

tc_eosinophils_down_gp <- simple_gprofiler(
  tc_eosinophils_sigenes_down,
  excel = glue("{gsea_prefix}/eosinophil_clinics_cali_up-v{ver}.xlsx"))

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

tc_eosinophils_down_gp

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

tc_eosinophils_down_gp$pvalue_plots$REAC

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

3.7.3 GSEA: Compare clinics, Monocyte samples

In the following block I repeated the above query, but this time looking at the monocyte samples.

tc_monocytes_up_gp <- simple_gprofiler(
  tc_monocytes_sigenes,
  excel = glue("{gsea_prefix}/monocyte_clinics-v{ver}.xlsx"))

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

tc_monocytes_up_gp

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

tc_monocytes_up_gp$pvalue_plots$REAC

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

tc_monocytes_up_gp$pvalue_plots$BP

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

tc_monocytes_down_gp <- simple_gprofiler(
  tc_monocytes_sigenes_down,
  excel = glue("{gsea_prefix}/monocyte_clinics_cali_up-v{ver}.xlsx"))

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

tc_monocytes_down_gp$pvalue_plots$REAC

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

tc_monocytes_down_gp$pvalue_plots$BP

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

3.7.3.1 GSEA: Compare clinics, Neutrophil samples

Ibid. This time looking at the Neutrophils. Thus the first two images should be a superset of the second and third pairs of images; assuming that the genes in the up/down list do not cause the groups to no longer be significant. Interestingly, the reactome search did not return any hits for the increased search.

tc_neutrophils_gp <- simple_gprofiler(
  tc_neutrophils_sigenes,
  excel = glue("{gsea_prefix}/neutrophil_clinics-v{ver}.xlsx"))

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

## tc_neutrophils_gp$pvalue_plots$REAC ## no hits
tc_neutrophils_gp$pvalue_plots$BP

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

tc_neutrophils_gp$pvalue_plots$TF

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

tc_neutrophils_up_gp <- simple_gprofiler(
  tc_neutrophils_sigenes_up,
  excel = glue("{gsea_prefix}/neutrophil_clinics_tumaco_up-v{ver}.xlsx"))

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

## tc_neutrophils_up_gp$pvalue_plots$REAC ## No hits
tc_neutrophils_up_gp$pvalue_plots$BP

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

tc_neutrophils_down_gp <- simple_gprofiler(
  tc_neutrophils_sigenes_down,
  excel = glue("{gsea_prefix}/neutrophil_clinics_cali_up-v{ver}.xlsx"))

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

tc_neutrophils_down_gp$pvalue_plots$REAC

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

tc_neutrophils_down_gp$pvalue_plots$BP

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

4 Compare DE: How similar are Tumaco C/F vs. Cali C/F

The following expands the cross-clinic query above to also test the neutrophils. Once again, I think it will pretty strongly support the hypothesis that the two clinics are not compatible.

We are concerned that the clinic-based batch effect may make our results essentially useless. One way to test this concern is to compare the set of genes observed different between the Cali Cure/Fail vs. the Tumaco Cure/Fail.

cali_table_nobatch <- tc_neutrophils_table_nobatch[["data"]][["cali"]]

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

tumaco_table_nobatch <- tc_neutrophils_table_nobatch[["data"]][["tumaco"]]

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

cali_merged_nobatch <- merge(cali_table_nobatch, tumaco_table_nobatch, by="row.names")

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

cor.test(cali_merged_nobatch[, "deseq_logfc.x"], cali_merged_nobatch[, "deseq_logfc.y"])

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

cali_aucc_nobatch <- calculate_aucc(cali_table_nobatch, tumaco_table_nobatch, px = "deseq_adjp",
                                    py = "deseq_adjp", lx = "deseq_logfc", ly = "deseq_logfc")

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

cali_aucc_nobatch$plot

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

5 Tumaco and Cali, cure vs. fail

In all of the above, we are looking to understand the differences between the two location. Let us now step back and perform the original question: fail/cure without regard to location.

I performed this query with a few different parameters, notably with(out) sva and again using each cell type, including biopsies. The main reasion I am keeping these comparisons is in the relatively weak hope that there will be sufficient signal in the full dataset that it might be able to overcome the apparently ridiculous batch effect from the two clinics.

5.1 All cell types together, with(out) SVA

table(pData(tc_valid)[["condition"]])

## 
##    cure failure 
##     122      62

tc_all_cf_de_sva <- all_pairwise(tc_valid, filter = TRUE, methods = methods,
                                 parallel = parallel, model_batch = "svaseq")

## 
##    cure failure 
##     122      62

## Error in checkForRemoteErrors(val): one node produced an error: c("Error in all_adjusters(input, estimate_type = model_batch, surrogates = surrogates) : \n  If an expt is not passed, then design _must_ be.\n", "noiseq")

tc_all_cf_table_sva <- combine_de_tables(
  tc_all_cf_de_sva, keepers = t_cf_contrast,
  excel = glue("{cf_prefix}/All_Samples/tc_valid_cf_table_sva-v{ver}.xlsx"))

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

tc_all_cf_sig_sva <- extract_significant_genes(
  tc_all_cf_table_sva,
  excel = glue("{cf_prefix}/All_Samples/tc_valid_cf_sig_sva-v{ver}.xlsx"))

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

tc_all_cf_de_batch <- all_pairwise(tc_valid, filter = TRUE, methods = methods,
                                   parallel = parallel, model_batch = TRUE)

## 
##    cure failure 
##     122      62 
## 
##  1  2  3 
## 83 50 51

## Error in all_pairwise(tc_valid, filter = TRUE, methods = methods, parallel = parallel, : object 'sv_model' not found

tc_all_cf_table_batch <- combine_de_tables(
  tc_all_cf_de_batch, keepers = t_cf_contrast,
  excel = glue("{cf_prefix}/All_Samples/tc_valid_cf_table_batch-v{ver}.xlsx"))

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

tc_all_cf_sig_batch <- extract_significant_genes(
  tc_all_cf_table_batch,
  excel = glue("{cf_prefix}/All_Samples/tc_valid_cf_sig_batch-v{ver}.xlsx"))

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

5.2 All cell types except biopsies

I am not sure if this is the best choice, but I call the set of all samples excluding biopsies ‘clinical’.

table(pData(tc_clinical_nobiop)[["condition"]])

## 
##    cure failure 
##     109      57

tc_clinical_cf_de_sva <- all_pairwise(tc_clinical_nobiop, filter = TRUE,
                                      parallel = parallel, model_batch = "svaseq",
                                      methods = methods)

## 
##    cure failure 
##     109      57

## Error in checkForRemoteErrors(val): one node produced an error: c("Error in all_adjusters(input, estimate_type = model_batch, surrogates = surrogates) : \n  If an expt is not passed, then design _must_ be.\n", "noiseq")

tc_clinical_cf_de_sva

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

tc_clinical_cf_table_sva <- combine_de_tables(
  tc_clinical_cf_de_sva, keepers = t_cf_contrast,
  excel = glue("{cf_prefix}/Clinical_Samples/tc_clinical_cf_table_sva-v{ver}.xlsx"))

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

tc_clinical_cf_table_sva

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

tc_clinical_cf_sig_sva <- extract_significant_genes(
  tc_clinical_cf_table_sva, according_to = "deseq",
  excel = glue("{cf_prefix}/Clinical_Samples/tc_clinical_cf_sig_sva-v{ver}.xlsx"))

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

tc_clinical_cf_sig_sva

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

tc_clinical_cf_de_batch <- all_pairwise(tc_clinical_nobiop, filter = TRUE,
                                        parallel = parallel, model_batch = TRUE,
                                        methods = methods)

## 
##    cure failure 
##     109      57 
## 
## eosinophils   monocytes neutrophils 
##          41          63          62

## Error in all_pairwise(tc_clinical_nobiop, filter = TRUE, parallel = parallel, : object 'sv_model' not found

tc_clinical_cf_de_batch

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

tc_clinical_cf_table_batch <- combine_de_tables(
  tc_clinical_cf_de_batch, keepers = t_cf_contrast,
  excel = glue("{cf_prefix}/Clinical_Samples/tc_clinical_cf_table_batch-v{ver}.xlsx"))

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

tc_clinical_cf_table_batch

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

tc_clinical_cf_sig_batch <- extract_significant_genes(
  tc_clinical_cf_table_batch, according_to = "deseq",
  excel = glue("{cf_prefix}/Clinical_Samples/tc_clinical_cf_sig_batch-v{ver}.xlsx"))

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

tc_clinical_cf_sig_batch

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

5.2.1 A portion of Supplemental Figure 11.

num_color <- color_choices[["cf"]][["cure"]]
den_color <- color_choices[["cf"]][["failure"]]
tc_clinical_cf_table <- tc_clinical_cf_table_sva[["data"]][["outcome"]]

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

tc_clinical_cf_volcano_top10 <- plot_volcano_condition_de(
  tc_clinical_cf_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)

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

pp(file = "figures/s11c_tc_clinical_cf_volcano_labeled_top10.svg")
tc_clinical_cf_volcano_top10[["plot"]]

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

dev.off()

## png 
##   2

tc_clinical_cf_volcano_top10[["plot"]]

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

5.3 Biopsies, with(out) SVA

In the following block, we repeat the same question, but using only the biopsy samples from both clinics.

tc_biopsies_cf <- set_expt_conditions(tc_biopsies, fact = "finaloutcome")

## The numbers of samples by condition are:

## 
##    cure failure 
##      13       5

tc_biopsies_cf_de_sva <- all_pairwise(tc_biopsies_cf, filter = TRUE, methods = methods,
                                      parallel = parallel, model_batch = "svaseq")

## 
##    cure failure 
##      13       5

## Error in checkForRemoteErrors(val): one node produced an error: c("Error in all_adjusters(input, estimate_type = model_batch, surrogates = surrogates) : \n  If an expt is not passed, then design _must_ be.\n", "noiseq")

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

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

tc_biopsies_cf_sig_sva <- extract_significant_genes(
  tc_biopsies_cf_table_sva,
  excel = glue("{cf_prefix}/All_Samples/tc_biopsies_cf_sig_sva-v{ver}.xlsx"))

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

tc_biopsies_cf_de_batch <- all_pairwise(tc_biopsies_cf, filter = TRUE, methods = methods,
                                        parallel = parallel, model_batch = TRUE)

## 
##    cure failure 
##      13       5 
## 
##  1 
## 18

## Error in all_pairwise(tc_biopsies_cf, filter = TRUE, methods = methods, : object 'sv_model' not found

tc_biopsies_cf_table_batch <- combine_de_tables(
  tc_biopsies_cf_de_batch, keepers = t_cf_contrast,
  excel = glue("{cf_prefix}/All_Samples/tc_biopsies_cf_table_batch-v{ver}.xlsx"))

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

tc_biopsies_cf_sig_batch <- extract_significant_genes(
  tc_biopsies_cf_table_batch,
  excel = glue("{cf_prefix}/All_Samples/tc_biopsies_cf_sig_batch-v{ver}.xlsx"))

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

5.4 Eosinophils, with(out) SVA

In the following block, we repeat the same question, but using only the Eosinophil samples from both clinics.

tc_eosinophils_cf <- set_expt_conditions(tc_eosinophils, fact = "finaloutcome")

## The numbers of samples by condition are:

## 
##    cure failure 
##      32       9

tc_eosinophils_cf_de_sva <- all_pairwise(tc_eosinophils_cf, filter = TRUE, methods = methods,
                                         parallel = parallel, model_batch = "svaseq")

## 
##    cure failure 
##      32       9

## Error in checkForRemoteErrors(val): one node produced an error: c("Error in all_adjusters(input, estimate_type = model_batch, surrogates = surrogates) : \n  If an expt is not passed, then design _must_ be.\n", "noiseq")

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

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

tc_eosinophils_cf_sig_sva <- extract_significant_genes(
  tc_eosinophils_cf_table_sva,
  excel = glue("{cf_prefix}/All_Samples/tc_eosinophils_cf_sig_sva-v{ver}.xlsx"))

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

tc_eosinophils_cf_de_batch <- all_pairwise(tc_eosinophils_cf, filter = TRUE,
                                           parallel = parallel, model_batch = TRUE,
                                           methods = methods)

## 
##    cure failure 
##      32       9 
## 
##  3  2  1 
## 13 14 14

## Error in all_pairwise(tc_eosinophils_cf, filter = TRUE, parallel = parallel, : object 'sv_model' not found

tc_eosinophils_cf_table_batch <- combine_de_tables(
  tc_eosinophils_cf_de_batch, keepers = t_cf_contrast,
  excel = glue("{cf_prefix}/All_Samples/tc_eosinophils_cf_table_batch-v{ver}.xlsx"))

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

tc_eosinophils_cf_sig_batch <- extract_significant_genes(
  tc_eosinophils_cf_table_batch,
  excel = glue("{cf_prefix}/All_Samples/tc_eosinophils_cf_sig_batch-v{ver}.xlsx"))

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

5.5 Monocytes, with(out) SVA

Repeat yet again, this time with the monocyte samples. The idea is to see if there is a cell type which is particularly good (or bad) at discriminating the two clinics.

tc_monocytes_cf <- set_expt_conditions(tc_monocytes, fact = "finaloutcome")

## The numbers of samples by condition are:

## 
##    cure failure 
##      39      24

tc_monocytes_cf_de_sva <- all_pairwise(tc_monocytes_cf, filter = TRUE, methods = methods,
                                       parallel = parallel, model_batch = "svaseq")

## 
##    cure failure 
##      39      24

## Error in checkForRemoteErrors(val): one node produced an error: c("Error in all_adjusters(input, estimate_type = model_batch, surrogates = surrogates) : \n  If an expt is not passed, then design _must_ be.\n", "noiseq")

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

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

tc_monocytes_cf_sig_sva <- extract_significant_genes(
  tc_monocytes_cf_table_sva,
  excel = glue("{cf_prefix}/All_Samples/tc_monocytes_cf_sig_sva-v{ver}.xlsx"))

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

tc_monocytes_cf_de_batch <- all_pairwise(tc_monocytes_cf, filter = TRUE, methods = methods,
                                         parallel = parallel, model_batch = TRUE)

## 
##    cure failure 
##      39      24 
## 
##  3  2  1 
## 19 18 26

## Error in all_pairwise(tc_monocytes_cf, filter = TRUE, methods = methods, : object 'sv_model' not found

tc_monocytes_cf_table_batch <- combine_de_tables(
  tc_monocytes_cf_de_batch, keepers = t_cf_contrast,
  excel = glue("{cf_prefix}/All_Samples/tc_monocytes_cf_table_batch-v{ver}.xlsx"))

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

tc_monocytes_cf_sig_batch <- extract_significant_genes(
  tc_monocytes_cf_table_batch,
  excel = glue("{cf_prefix}/All_Samples/tc_monocytes_cf_sig_batch-v{ver}.xlsx"))

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

5.6 Neutrophils, with(out) SVA

Last try, this time using the Neutrophil samples.

tc_neutrophils_cf <- set_expt_conditions(tc_neutrophils, fact = "finaloutcome")

## The numbers of samples by condition are:

## 
##    cure failure 
##      38      24

tc_neutrophils_cf_de_sva <- all_pairwise(tc_neutrophils_cf, parallel = parallel,
                                         filter = TRUE, model_batch = "svaseq",
                                         methods = methods)

## 
##    cure failure 
##      38      24

## Error in checkForRemoteErrors(val): one node produced an error: c("Error in all_adjusters(input, estimate_type = model_batch, surrogates = surrogates) : \n  If an expt is not passed, then design _must_ be.\n", "noiseq")

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

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

tc_neutrophils_cf_sig_sva <- extract_significant_genes(
  tc_neutrophils_cf_table_sva,
  excel = glue("{cf_prefix}/All_Samples/tc_neutrophils_cf_sig_sva-v{ver}.xlsx"))

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

tc_neutrophils_cf_de_batch <- all_pairwise(tc_neutrophils_cf, filter = TRUE,
                                           parallel = parallel, model_batch = TRUE,
                                           methods = methods)

## 
##    cure failure 
##      38      24 
## 
##  3  2  1 
## 19 18 25

## Error in all_pairwise(tc_neutrophils_cf, filter = TRUE, parallel = parallel, : object 'sv_model' not found

tc_neutrophils_cf_table_batch <- combine_de_tables(
  tc_neutrophils_cf_de_batch, keepers = t_cf_contrast,
  excel = glue("{cf_prefix}/All_Samples/tc_neutrophils_cf_table_batch-v{ver}.xlsx"))

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

tc_neutrophils_cf_sig_batch <- extract_significant_genes(
  tc_neutrophils_cf_table_batch,
  excel = glue("{cf_prefix}/All_Samples/tc_neutrophils_cf_sig_batch-v{ver}.xlsx"))

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

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

v1_vs_later <- all_pairwise(tc_v1vs, model_batch = "svaseq", methods = methods,
                            parallel = parallel, filter = TRUE)

## 
## first later 
##    65   101

## Error in checkForRemoteErrors(val): one node produced an error: c("Error in all_adjusters(input, estimate_type = model_batch, surrogates = surrogates) : \n  If an expt is not passed, then design _must_ be.\n", "noiseq")

v1_vs_later_table <- combine_de_tables(
  v1_vs_later, keepers = visit_v1later,
  excel = glue("{visit_prefix}/v1_vs_later_tables-v{ver}.xlsx"))

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

v1_vs_later_sig <- extract_significant_genes(
  v1_vs_later_table,
  excel = glue("{visit_prefix}/v1_vs_later_sig-v{ver}.xlsx"))

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

6.0.0.1 GSEA: V1 vs other visits.

v1later_gp <- all_gprofiler(v1_vs_later_sig)

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

v1later_gp[[1]]$pvalue_plots$REAC

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

v1later_gp[[2]]$pvalue_plots$REAC

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

7 Sex comparison

tc_sex_de <- all_pairwise(tc_sex, model_batch = "svaseq", methods = methods,
                          parallel = parallel, filter = TRUE)

## 
## female   male 
##     28    156

## Error in checkForRemoteErrors(val): one node produced an error: c("Error in all_adjusters(input, estimate_type = model_batch, surrogates = surrogates) : \n  If an expt is not passed, then design _must_ be.\n", "noiseq")

tc_sex_table <- combine_de_tables(
  tc_sex_de, excel = glue("{sex_prefix}/tc_sex_table-v{ver}.xlsx"))

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

tc_sex_sig <- extract_significant_genes(
  tc_sex_table, excel = glue("{sex_prefix}/tc_sex_sig-v{ver}.xlsx"))

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

tc_sex_gp <- all_gprofiler(tc_sex_sig)

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

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

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

tc_sex_cure_de <- all_pairwise(tc_sex_cure, model_batch = "svaseq",
                               parallel = parallel, filter = TRUE,
                               methods = methods)

## 
## female   male 
##     19    103

## Error in checkForRemoteErrors(val): one node produced an error: c("Error in all_adjusters(input, estimate_type = model_batch, surrogates = surrogates) : \n  If an expt is not passed, then design _must_ be.\n", "noiseq")

tc_sex_cure_de

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

tc_sex_cure_table <- combine_de_tables(
  tc_sex_cure_de, excel = glue("{sex_prefix}/tc_sex_cure_table-v{ver}.xlsx"))

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

tc_sex_cure_table

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

tc_sex_cure_sig <- extract_significant_genes(
  tc_sex_cure_table, excel = glue("{sex_prefix}/tc_sex_cure_sig-v{ver}.xlsx"))

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

tc_sex_cure_sig

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

7.0.0.1 GSEA: Sex comparisons both clinics

tc_sex_cure_gp <- all_gprofiler(tc_sex_cure_sig)

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

tc_sex_cure_gp

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

tc_sex_cure_gp[[1]][["pvalue_plots"]][["BP"]]

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

tc_sex_cure_gp[[2]][["pvalue_plots"]][["BP"]]

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

8 Ethnicity comparisons

tc_ethnicity_de <- all_pairwise(tc_etnia_expt, model_batch = "svaseq",
                                parallel = parallel, filter = TRUE,
                                methods = methods)

## 
##  afrocol indigena  mestiza 
##       91       46       47

## Error in checkForRemoteErrors(val): one node produced an error: c("Error in all_adjusters(input, estimate_type = model_batch, surrogates = surrogates) : \n  If an expt is not passed, then design _must_ be.\n", "noiseq")

tc_ethnicity_de

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

tc_ethnicity_table <- combine_de_tables(
  tc_ethnicity_de, keepers = ethnicity_contrasts,
  excel = glue("{eth_prefix}/tc_ethnicity_table-v{ver}.xlsx"))

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

tc_ethnicity_table

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

tc_ethnicity_table[["plots"]][["mestizo_indigenous"]][["deseq_ma_plots"]]

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

tc_ethnicity_table[["plots"]][["mestizo_afrocol"]][["deseq_ma_plots"]]

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

tc_ethnicity_table[["plots"]][["indigenous_afrocol"]][["deseq_ma_plots"]]

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

tc_ethnicity_sig <- extract_significant_genes(
  tc_ethnicity_table, excel = glue("{eth_prefix}/tc_ethnicity_sig-v{ver}.xlsx"))

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

ethnicity_cure <- subset_expt(tc_etnia_expt, subset = "finaloutcome=='cure'")

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

ethnicity_cure_de <- all_pairwise(ethnicity_cure, model_batch = "svaseq",
                                  parallel = parallel, filter = TRUE,
                                  methods = methods)

## 
##  afrocol indigena  mestiza 
##       39       36       47

## Error in checkForRemoteErrors(val): one node produced an error: c("Error in all_adjusters(input, estimate_type = model_batch, surrogates = surrogates) : \n  If an expt is not passed, then design _must_ be.\n", "noiseq")

ethnicity_cure_table <- combine_de_tables(
  ethnicity_cure_de, keepers = ethnicity_contrasts,
  excel = glue("{eth_prefix}/ethnicity_cure_table-v{ver}.xlsx"))

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

ethnicity_cure_table

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

ethnicity_cure_table[["plots"]][["mestizo_indigenous"]][["deseq_ma_plots"]]

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

ethnicity_cure_table[["plots"]][["mestizo_afrocol"]][["deseq_ma_plots"]]

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

ethnicity_cure_table[["plots"]][["indigenous_afrocol"]][["deseq_ma_plots"]]

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

ethnicity_cure_sig <- extract_significant_genes(
  ethnicity_cure_table, excel = glue("{eth_prefix}/ethnicity_cure_sig-v{ver}.xlsx"))

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

ethnicity_cure_sig

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

8.0.0.1 GSEA: Ethnicity differences

Performed once with both clinics and again with only Tumaco.

tc_ethnicity_gp <- all_gprofiler(tc_ethnicity_sig)

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

pander::pander(sessionInfo())

R version 4.4.1 (2024-06-14)

Platform: x86_64-conda-linux-gnu

locale: C

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

other attached packages: ruv(v.0.9.7.1), DOSE(v.3.28.2), forcats(v.1.0.0), dplyr(v.1.1.4), hpgltools(v.1.0), Matrix(v.1.6-5), glue(v.1.7.0), SummarizedExperiment(v.1.32.0), GenomicRanges(v.1.54.1), GenomeInfoDb(v.1.38.6), IRanges(v.2.36.0), S4Vectors(v.0.40.2), MatrixGenerics(v.1.14.0), matrixStats(v.1.2.0), Biobase(v.2.62.0) and BiocGenerics(v.0.48.1)

loaded via a namespace (and not attached): RColorBrewer(v.1.1-3), jsonlite(v.1.8.8), magrittr(v.2.0.3), rmarkdown(v.2.25), fs(v.1.6.3), zlibbioc(v.1.48.0), vctrs(v.0.6.5), memoise(v.2.0.1), RCurl(v.1.98-1.14), htmltools(v.0.5.7), S4Arrays(v.1.2.0), broom(v.1.0.5), SparseArray(v.1.2.4), sass(v.0.4.8), bslib(v.0.6.1), desc(v.1.4.3), htmlwidgets(v.1.6.4), plyr(v.1.8.9), testthat(v.3.2.1), plotly(v.4.10.4), cachem(v.1.0.8), mime(v.0.12), lifecycle(v.1.0.4), iterators(v.1.0.14), pkgconfig(v.2.0.3), R6(v.2.5.1), fastmap(v.1.1.1), GenomeInfoDbData(v.1.2.11), shiny(v.1.8.0), digest(v.0.6.34), colorspace(v.2.1-0), AnnotationDbi(v.1.64.1), DESeq2(v.1.42.0), rprojroot(v.2.0.4), pkgload(v.1.3.4), RSQLite(v.2.3.5), fansi(v.1.0.6), httr(v.1.4.7), abind(v.1.4-5), mgcv(v.1.9-1), compiler(v.4.4.1), pander(v.0.6.5), doParallel(v.1.0.17), bit64(v.4.0.5), backports(v.1.4.1), BiocParallel(v.1.36.0), DBI(v.1.2.2), DelayedArray(v.0.28.0), corpcor(v.1.6.10), HDO.db(v.0.99.1), tools(v.4.4.1), zip(v.2.3.1), httpuv(v.1.6.14), nlme(v.3.1-164), GOSemSim(v.2.28.1), promises(v.1.2.1), grid(v.4.4.1), reshape2(v.1.4.4), fgsea(v.1.28.0), generics(v.0.1.3), sva(v.3.50.0), gtable(v.0.3.4), preprocessCore(v.1.64.0), tidyr(v.1.3.1), data.table(v.1.15.0), utf8(v.1.2.4), XVector(v.0.42.0), foreach(v.1.5.2), pillar(v.1.9.0), stringr(v.1.5.1), yulab.utils(v.0.1.7), limma(v.3.58.1), genefilter(v.1.84.0), later(v.1.3.2), splines(v.4.4.1), lattice(v.0.22-5), survival(v.3.5-8), bit(v.4.0.5), annotate(v.1.80.0), tidyselect(v.1.2.0), GO.db(v.3.18.0), locfit(v.1.5-9.8), Biostrings(v.2.70.2), knitr(v.1.45), gridExtra(v.2.3), edgeR(v.4.0.16), xfun(v.0.42), statmod(v.1.5.0), brio(v.1.1.4), stringi(v.1.8.3), lazyeval(v.0.2.2), yaml(v.2.3.8), evaluate(v.0.23), codetools(v.0.2-19), tibble(v.3.2.1), qvalue(v.2.34.0), BiocManager(v.1.30.25), graph(v.1.80.0), cli(v.3.6.2), xtable(v.1.8-4), munsell(v.0.5.0), jquerylib(v.0.1.4), Rcpp(v.1.0.12), png(v.0.1-8), XML(v.3.99-0.16.1), parallel(v.4.4.1), ellipsis(v.0.3.2), ggplot2(v.3.5.0), blob(v.1.2.4), bitops(v.1.0-7), viridisLite(v.0.4.2), GSEABase(v.1.64.0), scales(v.1.3.0), openxlsx(v.4.2.5.2), purrr(v.1.0.2), crayon(v.1.5.2), rlang(v.1.1.3), cowplot(v.1.1.3), fastmatch(v.1.1-4) and KEGGREST(v.1.42.0)

message("This is hpgltools commit: ", get_git_commit())

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

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

## > git reset e7016a27043e88fcf9b6f1dffe682dcc7a6477c5

## This is hpgltools commit: Mon Dec 2 16:30:15 2024 -0500: e7016a27043e88fcf9b6f1dffe682dcc7a6477c5

message("Saving to ", savefile)

## Saving to 03differential_expression_both.rda.xz

# tmp <- sm(saveme(filename = savefile))

tmp <- loadme(filename = savefile)