TMRC2 202509: Macrophage Differential Expression.

1 Changelog

• 202401-202405: Cleanups, formatting, ensuring that everything works in the container.
• 202310: Cleaning up to make everything pass within a containerized environment.
• 202310: Received a set of colors and contrasts of interest for a barplot of significance.
• 20230410: Making some changes to improve the differential expression plots as well as prepare for some different pathway/GSEA/GSVA analyses on the data.

2 Introduction

Having established that the TMRC2 macrophage data looks robust and illustrative of a couple of interesting questions, let us perform a couple of differential analyses of it.

Also note that as of 202212, we received a new set of samples which now include some which are a completely different cell type, U937. As their ATCC page states, they are malignant cells taken from the pleural effusion of a 37 year old white male with histiocytic lymphoma and which exhibit the morphology of monocytes. Thus, this document now includes some comparisons of the cell types as well as the various macrophage donors (given that there are now more donors too).

2.1 Human data

I am moving the dataset manipulations here so that I can look at them all together before running the various DE analyses.

2.2 Create sets focused on drug, celltype, strain, and combinations

Let us start by playing with the metadata a little and create sets with the condition set to:

• Drug treatment
• Cell type (macrophage or U937)
• Donor
• Infection Strain
• Some useful combinations thereof

In addition, keep mental track of which datasets are comprised of all samples vs. those which are only macrophage vs. those which are only U937. (Thus, the usage of all_human vs. hs_macr vs. u937 as prefixes for the data structures.)

Ideally, these recreations of the data should perhaps be in the datastructures worksheet.

all_human <- sanitize_metadata(hs_macrophage, columns = "drug") %>%
  set_conditions(fact = "drug", colors = color_choices[["drug"]]) %>%
  set_batches(fact = "typeofcells")

## Recasting the data.frame to DataFrame.

## The numbers of samples by condition are:

## 
## antimony     none 
##       34       34

## The number of samples by batch are:

## 
## Macrophages        U937 
##          54          14

## The following 3 lines were copy/pasted to datastructures and should be removed soon.
no_strain_idx <- colData(all_human)[["strainid"]] == "none"
##colData(all_human)[["strainid"]] <- paste0("s", colData(all_human)[["strainid"]],
##                                         "_", colData(all_human)[["macrophagezymodeme"]])
colData(all_human)[no_strain_idx, "strainid"] <- "none"
table(colData(all_human)[["strainid"]])

## 
## 10763 10772 10977 11026 11075 11126 12251 12309 12355 12367  2169  7158  none 
##     2     8     2     2     2     8     7     8     2     7     8     2    10

all_human_types <- set_conditions(all_human, fact = "typeofcells") %>%
  set_batches(fact = "drug")

## The numbers of samples by condition are:

## 
## Macrophages        U937 
##          54          14

## The number of samples by batch are:

## 
## antimony     none 
##       34       34

type_zymo_fact <- paste0(colData(all_human_types)[["condition"]], "_",
                         colData(all_human_types)[["macrophagezymodeme"]])
type_zymo <- set_conditions(all_human_types, fact = type_zymo_fact)

## The numbers of samples by condition are:

## 
## Macrophages_none  Macrophages_z22  Macrophages_z23        U937_none 
##                8               23               23                2 
##         U937_z22         U937_z23 
##                6                6

type_drug_fact <- paste0(colData(all_human_types)[["condition"]], "_",
                         colData(all_human_types)[["drug"]])
type_drug <- set_conditions(all_human_types, fact = type_drug_fact)

## The numbers of samples by condition are:

## 
## Macrophages_antimony     Macrophages_none        U937_antimony 
##                   27                   27                    7 
##            U937_none 
##                    7

strain_fact <- colData(all_human_types)[["strainid"]]
table(strain_fact)

## strain_fact
## 10763 10772 10977 11026 11075 11126 12251 12309 12355 12367  2169  7158  none 
##     2     8     2     2     2     8     7     8     2     7     8     2    10

new_conditions <- paste0(colData(hs_macrophage)[["macrophagetreatment"]], "_",
                         colData(hs_macrophage)[["macrophagezymodeme"]])
## Note the sanitize() call is redundant with the addition of sanitize() in the
## datastructures file, but I don't want to wait to rerun that.
hs_macr <- set_conditions(hs_macrophage, fact = new_conditions) %>%
  sanitize_metadata(column = "drug") %>%
  subset_se(subset = "typeofcells!='U937'") %>%
  set_se_colors(color_choices[["treatment_zymo"]])

## The numbers of samples by condition are:

## 
##    inf_sb_z22    inf_sb_z23       inf_z22       inf_z23    uninf_none 
##            15            14            14            15             5 
## uninf_sb_none 
##             5

## Recasting the data.frame to DataFrame.

2.2.1 Separate Macrophage samples

Once again, we should reconsider where the following block is placed, but these datastructures are likely to be used in many of the following analyses.

hs_macr_drug_se <- set_conditions(hs_macr, fact = "drug", colors = color_choices[["drug"]])

## The numbers of samples by condition are:

## 
## antimony     none 
##       27       27

hs_macr_strain_se <- set_conditions(hs_macr, fact = "macrophagezymodeme",
                                      colors = color_choices[["zymo"]]) %>%
  subset_se(subset = "macrophagezymodeme != 'none'")

## The numbers of samples by condition are:

## 
## none  z22  z23 
##    8   23   23

table(colData(hs_macr)[["strainid"]])

## 
## 10763 10772 10977 11026 11075 11126 12251 12309 12355 12367  2169  7158  none 
##     2     6     2     2     2     6     5     6     2     5     6     2     8

2.2.2 Refactor U937 samples

The U937 samples were separated in the datastructures file, but we want to use the combination of drug/zymodeme with them pretty much exclusively.

new_conditions <- paste0(colData(hs_u937)[["macrophagetreatment"]], "_",
                         colData(hs_u937)[["macrophagezymodeme"]])
u937_se <- set_conditions(hs_u937, fact = new_conditions,
                            colors = color_choices[["treatment_zymo"]])

## The numbers of samples by condition are:

## 
##    inf_sb_z22    inf_sb_z23       inf_z22       inf_z23    uninf_none 
##             3             3             3             3             1 
## uninf_sb_none 
##             1

2.3 Contrasts used in this document

Given the various ways we have chopped up this dataset, there are a few general types of contrasts we will perform, which will then be combined into greater complexity:

• drug treatment: Antimonal treated or not.
• strains used: Uninfected, z2.3, and z2.2.
• cellltypes: U937 or macrophage.
• donors: The person from whom the macrophages were taken.

In the end, our actual goal is to consider the variable effects of drug+strain and see if we can discern patterns which lead to better or worse drug treatment outcome.

There is a set of contrasts in which we are primarily interested in this data, these follow. I created one ratio of ratios contrast which I think has the potential to ask our biggest question.

## Each of the following lists has the name of the contrast as the key
## followed by a two element vector comprised of the numerator and
## denominator as the value.  In the case of this first contrast, that
## is comprised of a string which manually defines a series of more
## complex contrasts than the usual/simple pairwise.
tmrc2_human_extra <- "z23drugnodrug_vs_z22drugnodrug = (conditioninf_sb_z23 - conditioninf_z23) - (conditioninf_sb_z22 - conditioninf_z22), z23z22drug_vs_z23z22nodrug = (conditioninf_sb_z23 - conditioninf_sb_z22) - (conditioninf_z23 - conditioninf_z22)"
tmrc2_human_keepers <- list(
  "z23nosb_vs_uninf" = c("inf_z23", "uninf_none"),
  "z22nosb_vs_uninf" = c("inf_z22", "uninf_none"),
  "z23nosb_vs_z22nosb" = c("inf_z23", "inf_z22"),
  "z23sb_vs_z22sb" = c("inf_sb_z23", "inf_sb_z22"),
  "z23sb_vs_z23nosb" = c("inf_sb_z23", "inf_z23"),
  "z22sb_vs_z22nosb" = c("inf_sb_z22", "inf_z22"),
  "z23sb_vs_sb" = c("inf_sb_z23", "uninf_sb_none"),
  "z22sb_vs_sb" = c("inf_sb_z22", "uninf_sb_none"),
  "z23sb_vs_uninf" = c("inf_sb_z23", "uninf_none"),
  "z22sb_vs_uninf" = c("inf_sb_z22", "uninf_none"),
  "sb_vs_uninf" = c("uninf_sb_none", "uninf_none"),
  "extra_z2322" = c("z23drugnodrug", "z22drugnodrug"),
  "extra_drugnodrug" = c("z23z22drug", "z23z22nodrug"))
single_tmrc2_keeper <- list(
  "z22sb_vs_sb" = c("inf_sb_z22", "uninf_sb_none"))
tmrc2_drug_keepers <- list(
  "drug" = c("antimony", "none"))
tmrc2_type_keepers <- list(
  "type" = c("U937", "Macrophages"))
tmrc2_strain_keepers <- list(
  "strain" = c("z23", "z22"))
type_zymo_extra <- "zymos_vs_types = (conditionU937_z23 - conditionU937_z22) - (conditionMacrophages_z23 - conditionMacrophages_z22)"
tmrc2_typezymo_keepers <- list(
  "u937_macr" = c("Macrophages_none", "U937_none"),
  "zymo_macr" = c("Macrophages_z23", "Macrophages_z22"),
  "zymo_u937" = c("U937_z23", "U937_z22"),
  "z23_types" = c("U937_z23", "Macrophages_z23"),
  "z22_types" = c("U937_z22", "Macrophages_z22"),
  "zymos_types" = c("zymos_vs_types"))
tmrc2_typedrug_keepers <- list(
  "type_nodrug" = c("U937_none", "Macrophages_none"),
  "type_drug" = c("U937_antimony", "Macrophages_antimony"),
  "macr_drugs" = c("Macrophages_antimony", "Macrophages_none"),
  "u937_drugs" = c("U937_antimony", "U937_none"))
u937_keepers <- list(
  "z23nosb_vs_uninf" = c("inf_z23", "uninf_none"),
  "z22nosb_vs_uninf" = c("inf_z22", "uninf_none"),
  "z23nosb_vs_z22nosb" = c("inf_z23", "inf_z22"),
  "z23sb_vs_z22sb" = c("inf_sb_z23", "inf_sb_z22"),
  "z23sb_vs_z23nosb" = c("inf_sb_z23", "inf_z23"),
  "z22sb_vs_z22nosb" = c("inf_sb_z22", "inf_z22"),
  "z23sb_vs_sb" = c("inf_sb_z23", "uninf_sb_none"),
  "z22sb_vs_sb" = c("inf_sb_z22", "uninf_sb_none"),
  "z23sb_vs_uninf" = c("inf_sb_z23", "uninf_none"),
  "z22sb_vs_uninf" = c("inf_sb_z22", "uninf_none"),
  "sb_vs_uninf" = c("uninf_sb_none", "uninf_none"))
## If some cases, when the set of significant genes was chosen, an
## additional filter was added to exclude genes with expression values
## less than 'high_expression' according to the
## 'high_expression_column' in the table.
high_expression <- 128
high_expression_column <- "deseq_basemean"

combined_to_tsv <- function(combined, celltype = "all") {
  keepers <- combined[["keepers"]]
  for (k in seq_len(length(keepers))) {
    kname <- names(keepers)[k]
    numerator <- keepers[[k]][1]
    denominator <- keepers[[k]][2]
    filename <- glue("analyses/macrophage_de/tsv_tables/tmrc2_{celltype}_{kname}_n{numerator}_d{denominator}-v{ver}.xlsx")
    kdata <- combined[["data"]][[kname]]
    if (is.null(kdata[["basic_num"]])) {
      next
    }
    wanted <- c("hgnc_symbol", "deseq_logfc", "deseq_adjp",
                "deseq_basemean", "deseq_num", "deseq_den")
    wanted_data <- kdata[, wanted]
    colnames(wanted_data) <- c("hgnc_symbol", "deseq_logfc", "deseq_adjp",
                               "deseq_mean", "deseq_numerator", "deseq_denominator")
    write_xlsx(data = wanted_data, excel = filename)
  }
}

write_all_gp <- function(all_gp, suffix = NULL) {
  all_written <- list()
  for (g in seq_len(length(all_gp))) {
    name <- names(all_gp)[g]
    datum <- all_gp[[name]]
    filename <- glue("analyses/macrophage_de/gprofiler/{name}_gprofiler-v{ver}.xlsx")
    if (!is.null(suffix)) {
      filename <- glue("analyses/macrophage_de/gprofiler/{name}_gprofiler{suffix}-v{ver}.xlsx")
    }
    written <- sm(write_gprofiler_data(datum, excel = filename))
    all_written[[g]] <- written
  }
  return(all_written)
}

2.4 Primary queries

There is a series of initial questions which make some sense to me, but these do not necessarily match the set of questions which are most pressing. I am hoping to pull both of these sets of queries in one.

Before extracting these groups of queries, let us invoke the all_pairwise() function and get all of the likely contrasts along with one or more extras that might prove useful (the ‘extra’ argument).

The structure of these blocks will all basically be identical:

• Perform a set of pairwise contrasts of all the conditions against each other. Optionally use sva.
• Given that result, dump it in its entirety to an xlsx file in the analyses/ directory.
• Given those combined tables, extract from them the set deemed ‘significant’ by whatever criteria we want to try. (Usually |lfc| >= 1.0, adjusted p <= 0.05; but potentially also expression >= x and sometimes a set of less stringent values (|lfc| >= 0.6))
• Given one or more gene sets deemed ‘significant’ pass them to gProfiler2 and see what pops out.

2.4.1 Combined U937 and Macrophages: Compare drug effects

When we have the u937 cells in the same dataset as the macrophages, that provides an interesting opportunity to see if we can observe drug-dependant effects which are shared across both cell types.

Note to self: given the changes to hpgltools I may need to specify the statistical model string when I am using svaseq for some/many/all of these comparisons.

drug_de <- all_pairwise(all_human, filter = TRUE, model_svs = "svaseq",
                        model_fstring = "~ 0 + condition")

## antimony     none 
##       34       34

## Removing 9198 low-count genes (12283 remaining).

## Basic step 0/3: Normalizing data.

## Basic step 0/3: Converting data.

## I think this is failing? SummarizedExperiment

## Basic step 0/3: Transforming data.

## Setting 85798 entries to zero.

## This received a matrix of SVs.

## converting counts to integer mode

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## Warning in createContrastL(objFlt$formula, objFlt$data, L): Contrasts with only
## a single non-zero term are already evaluated by default.

## conditions
## antimony     none 
##       34       34

## conditions
## antimony     none 
##       34       34

## conditions
## antimony     none 
##       34       34

drug_de

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

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

## The primary analysis performed 1 comparisons.

## The logFC agreement among the methods follows:

##                 nn_vs_ntmn
## basic_vs_deseq      0.8651
## basic_vs_dream      0.8738
## basic_vs_ebseq      0.8367
## basic_vs_edger      0.8671
## basic_vs_limma      0.8771
## basic_vs_noiseq     0.8726
## deseq_vs_dream      0.9719
## deseq_vs_ebseq      0.9344
## deseq_vs_edger      0.9988
## deseq_vs_limma      0.9665
## deseq_vs_noiseq     0.9738
## dream_vs_ebseq      0.9770
## dream_vs_edger      0.9741
## dream_vs_limma      0.9953
## dream_vs_noiseq     0.9616
## ebseq_vs_edger      0.9383
## ebseq_vs_limma      0.9728
## ebseq_vs_noiseq     0.9486
## edger_vs_limma      0.9688
## edger_vs_noiseq     0.9745
## limma_vs_noiseq     0.9584

drug_table <- combine_de_tables(
  drug_de, keepers = tmrc2_drug_keepers,
  excel = glue("analyses/macrophage_de/de_tables/macrophage_drug_comparison-v{ver}.xlsx"))

## Looking for subscript invalid names, start of extract_keepers.

## Error in combine_extracted_plots(entry_name, combined, wanted_denominator, : object 'color_low' not found

drug_table

## Error: object 'drug_table' not found

combined_to_tsv(drug_table, celltype = "all")

## Error: object 'drug_table' not found

drug_sig <- extract_significant_genes(
  drug_table,
  excel = glue("analyses/macrophage_de/sig_tables/macrophage_drug_sig-v{ver}.xlsx"))

## Error: object 'drug_table' not found

drug_sig

## Error: object 'drug_sig' not found

drug_highsig <- extract_significant_genes(
  drug_table, min_mean_exprs = high_expression, exprs_column = high_expression_column,
  excel = glue("analyses/macrophage_de/sig_tables/macrophage_drug_highsig-v{ver}.xlsx"))

## Error: object 'drug_table' not found

drug_highsig

## Error: object 'drug_highsig' not found

drug_lesssig <- extract_significant_genes(
  drug_table, lfc = 0.6,
  excel = glue("analyses/macrophage_de/sig_tables/macrophage_drug_lesssig-v{ver}.xlsx"))

## Error: object 'drug_table' not found

drug_lesssig

## Error: object 'drug_lesssig' not found

2.4.1.1 gProfiler2 results of the significant drug genes

all_drug_gp <- all_gprofiler(drug_sig, enrich_id_column = "hgnc_symbol")

## Error: object 'drug_sig' not found

all_drug_gp

## Error: object 'all_drug_gp' not found

written <- write_all_gp(all_drug_gp)

## Error: object 'all_drug_gp' not found

all_drug_lesssig <- all_gprofiler(drug_lesssig, enrich_id_column = "hgnc_symbol")

## Error: object 'drug_lesssig' not found

written <- write_all_gp(all_drug_lesssig, suffix = "_lfc0.6_")

## Error: object 'all_drug_lesssig' not found

2.4.2 Combined U937 and Macrophages: compare cell types

There are a couple of ways one might want to directly compare the two cell types.

• Given that the variance between the two celltypes is so huge, just compare all samples.
• One might want to compare them with the interaction effects of drug/zymodeme.

type_de <- all_pairwise(all_human_types, filter = TRUE, model_fstring = "~ 0 + condition",
                        model_svs = "svaseq")

## Macrophages        U937 
##          54          14

## Removing 9198 low-count genes (12283 remaining).

## Basic step 0/3: Normalizing data.

## Basic step 0/3: Converting data.

## I think this is failing? SummarizedExperiment

## Basic step 0/3: Transforming data.

## Setting 85798 entries to zero.

## This received a matrix of SVs.

## converting counts to integer mode

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## Warning in createContrastL(objFlt$formula, objFlt$data, L): Contrasts with only
## a single non-zero term are already evaluated by default.

## conditions
## Macrophages        U937 
##          54          14

## conditions
## Macrophages        U937 
##          54          14

## conditions
## Macrophages        U937 
##          54          14

type_de

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

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

## The primary analysis performed 1 comparisons.

## The logFC agreement among the methods follows:

##                 U937_vs_Mc
## basic_vs_deseq      0.8587
## basic_vs_dream      0.8853
## basic_vs_ebseq      0.8125
## basic_vs_edger      0.8576
## basic_vs_limma      0.8982
## basic_vs_noiseq     0.9114
## deseq_vs_dream      0.9938
## deseq_vs_ebseq      0.9805
## deseq_vs_edger      0.9974
## deseq_vs_limma      0.9849
## deseq_vs_noiseq     0.9835
## dream_vs_ebseq      0.9682
## dream_vs_edger      0.9947
## dream_vs_limma      0.9932
## dream_vs_noiseq     0.9910
## ebseq_vs_edger      0.9836
## ebseq_vs_limma      0.9490
## ebseq_vs_noiseq     0.9652
## edger_vs_limma      0.9859
## edger_vs_noiseq     0.9829
## limma_vs_noiseq     0.9817

type_table <- combine_de_tables(
  type_de, keepers = tmrc2_type_keepers,
  excel = glue("analyses/macrophage_de/de_tables/macrophage_type_comparison-v{ver}.xlsx"))

## Looking for subscript invalid names, start of extract_keepers.

## Error in combine_extracted_plots(entry_name, combined, wanted_denominator, : object 'color_low' not found

type_table

## Error: object 'type_table' not found

combined_to_tsv(type_table, celltype = "all")

## Error: object 'type_table' not found

type_sig <- extract_significant_genes(
  type_table,
  excel = glue("analyses/macrophage_de/sig_tables/macrophage_type_sig-v{ver}.xlsx"))

## Error: object 'type_table' not found

type_sig

## Error: object 'type_sig' not found

type_highsig <- extract_significant_genes(
  type_table, min_mean_exprs = high_expression, exprs_column = high_expression_column,
  excel = glue("analyses/macrophage_de/sig_tables/macrophage_type_highsig-v{ver}.xlsx"))

## Error: object 'type_table' not found

type_highsig

## Error: object 'type_highsig' not found

type_lesssig <- extract_significant_genes(
  type_table, lfc = 0.6,
  excel = glue("analyses/macrophage_de/sig_tables/macrophage_type_lesssig-v{ver}.xlsx"))

## Error: object 'type_table' not found

type_sig

## Error: object 'type_sig' not found

2.4.2.1 Combined factors of interest: celltype+zymodeme

Given the above explicit comparison of all samples comprising the two cell types, now let us look at the drug treatment+zymodeme status with all samples, macrophages and U937.

type_zymo_de <- all_pairwise(type_zymo, filter = TRUE, model_svs = "svaseq",
                             model_fstring = "~ 0 + condition",
                             extra_contrasts = type_zymo_extra)

## Macrophages_none  Macrophages_z22  Macrophages_z23        U937_none 
##                8               23               23                2 
##         U937_z22         U937_z23 
##                6                6

## Removing 9198 low-count genes (12283 remaining).

## Basic step 0/3: Normalizing data.

## Basic step 0/3: Converting data.

## I think this is failing? SummarizedExperiment

## Basic step 0/3: Transforming data.

## Setting 85798 entries to zero.

## This received a matrix of SVs.

## converting counts to integer mode

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## The contrast zymos is not in the results.

## If this is not an extra contrast, then this is an error.

## Warning in createContrastL(objFlt$formula, objFlt$data, L): Contrasts with only
## a single non-zero term are already evaluated by default.

## conditions
## Macrophages_none  Macrophages_z22  Macrophages_z23        U937_none 
##                8               23               23                2 
##         U937_z22         U937_z23 
##                6                6

## conditions
## Macrophages_none  Macrophages_z22  Macrophages_z23        U937_none 
##                8               23               23                2 
##         U937_z22         U937_z23 
##                6                6

## conditions
## Macrophages_none  Macrophages_z22  Macrophages_z23        U937_none 
##                8               23               23                2 
##         U937_z22         U937_z23 
##                6                6

type_zymo_de

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

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

## The primary analysis performed 15 comparisons.

Strangely, as of 20250903, the following line throws an error in the container: “Error: subscript contains invalid names.”

But when I run it manually I get no error. I assume this means that I fell behind while maintaining hpgltools HEAD with the container?

In addition, the functions actually did return successfully.

type_zymo_table <- combine_de_tables(
  type_zymo_de, keepers = tmrc2_typezymo_keepers,
  excel = glue("analyses/macrophage_de/de_tables/macrophage_type_zymo_comparison-v{ver}.xlsx"))

## Looking for subscript invalid names, start of extract_keepers.

## Error in combine_extracted_plots(entry_name, combined, wanted_denominator, : object 'color_low' not found

type_zymo_table

## Error: object 'type_zymo_table' not found

combined_to_tsv(type_zymo_table, celltype = "all")

## Error: object 'type_zymo_table' not found

type_zymo_sig <- extract_significant_genes(
  type_zymo_table,
  excel = glue("analyses/macrophage_de/sig_tables/macrophage_type_zymo_sig-v{ver}.xlsx"))

## Error: object 'type_zymo_table' not found

type_zymo_sig

## Error: object 'type_zymo_sig' not found

type_zymo_highsig <- extract_significant_genes(
  type_zymo_table, min_mean_exprs = high_expression, exprs_column = high_expression_column,
  excel = glue("analyses/macrophage_de/sig_tables/macrophage_type_zymo_highsig-v{ver}.xlsx"))

## Error: object 'type_zymo_table' not found

type_zymo_lesssig <- extract_significant_genes(
  type_zymo_table, lfc = 0.6,
  excel = glue("analyses/macrophage_de/sig_tables/macrophage_type_zymo_lesssig-v{ver}.xlsx"))

## Error: object 'type_zymo_table' not found

type_zymo_lesssig

## Error: object 'type_zymo_lesssig' not found

2.4.2.2 Combined factors of interest: celltype+drug

The ‘type_drug’ datastructure is the same as above, but the condition is created from the concatenation of the cell type and drug treatment.

type_drug_de <- all_pairwise(type_drug, filter = TRUE, model_svs = "svaseq",
                             model_fstring = "~ 0 + condition")

## Macrophages_antimony     Macrophages_none        U937_antimony 
##                   27                   27                    7 
##            U937_none 
##                    7

## Removing 9198 low-count genes (12283 remaining).

## Basic step 0/3: Normalizing data.

## Basic step 0/3: Converting data.

## I think this is failing? SummarizedExperiment

## Basic step 0/3: Transforming data.

## Setting 85798 entries to zero.

## This received a matrix of SVs.

## converting counts to integer mode

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## Warning in createContrastL(objFlt$formula, objFlt$data, L): Contrasts with only
## a single non-zero term are already evaluated by default.

## conditions
## Macrophages_antimony     Macrophages_none        U937_antimony 
##                   27                   27                    7 
##            U937_none 
##                    7

## conditions
## Macrophages_antimony     Macrophages_none        U937_antimony 
##                   27                   27                    7 
##            U937_none 
##                    7

## conditions
## Macrophages_antimony     Macrophages_none        U937_antimony 
##                   27                   27                    7 
##            U937_none 
##                    7

type_drug_de

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

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

## The primary analysis performed 6 comparisons.

type_drug_table <- combine_de_tables(
  type_drug_de, keepers = tmrc2_typedrug_keepers,
  excel = glue("analyses/macrophage_de/de_tables/macrophage_type_drug_comparison-v{ver}.xlsx"))

## Looking for subscript invalid names, start of extract_keepers.

## Error in combine_extracted_plots(entry_name, combined, wanted_denominator, : object 'color_low' not found

type_drug_table

## Error: object 'type_drug_table' not found

#combined_to_tsv(type_drug_table, celltype = "all")

type_drug_sig <- extract_significant_genes(
  type_drug_table,
  excel = glue("analyses/macrophage_de/sig_tables/macrophage_type_drug_sig-v{ver}.xlsx"))

## Error: object 'type_drug_table' not found

type_drug_sig

## Error: object 'type_drug_sig' not found

type_drug_highsig <- extract_significant_genes(
  type_drug_table, min_mean_exprs = high_expression, exprs_column = high_expression_column,
  excel = glue("analyses/macrophage_de/sig_tables/macrophage_type_drug_highsig-v{ver}.xlsx"))

## Error: object 'type_drug_table' not found

type_drug_highsig

## Error: object 'type_drug_highsig' not found

type_drug_lesssig <- extract_significant_genes(
  type_drug_table, lfc = 0.6,
  excel = glue("analyses/macrophage_de/sig_tables/macrophage_type_drug_lesssig-v{ver}.xlsx"))

## Error: object 'type_drug_table' not found

type_drug_lesssig

## Error: object 'type_drug_lesssig' not found

3 Individual cell types

At this point, I think it is fair to say that the two cell types are sufficiently different that they do not really belong together in a single analysis.

3.1 drug or strain effects, single cell type

One of the queries Najib asked which I think I misinterpreted was to look at drug and/or strain effects. My interpretation is somewhere below and was not what he was looking for. Instead, he was looking to see all(macrophage) drug/nodrug and all(macrophage) z23/z22 and compare them to each other. It may be that this is still a wrong interpretation, if so the most likely comparison is either:

• (z23drug/z22drug) / (z23nodrug/z22nodrug), or perhaps
• (z23drug/z23nodrug) / (z22drug/z22nodrug),

I am not sure those confuse me, and at least one of them is below

3.2 Macrophages

In these blocks we will explicitly query only one factor at a time, drug and strain. The eventual goal is to look for effects of drug treatment and/or strain treatment which are shared?

3.2.1 Macrophage Drug only

Thus we will start with the pure drug query. In this block we will look only at the drug/nodrug effect.

hs_macr_drug_de <- all_pairwise(hs_macr_drug_se, filter = TRUE, model_svs = "svaseq",
                                model_fstring = "~ 0 + condition")

## antimony     none 
##       27       27

## Removing 9725 low-count genes (11756 remaining).

## Basic step 0/3: Normalizing data.

## Basic step 0/3: Converting data.

## I think this is failing? SummarizedExperiment

## Basic step 0/3: Transforming data.

## Setting 40036 entries to zero.

## This received a matrix of SVs.

## converting counts to integer mode

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## Warning in createContrastL(objFlt$formula, objFlt$data, L): Contrasts with only
## a single non-zero term are already evaluated by default.

## conditions
## antimony     none 
##       27       27

## conditions
## antimony     none 
##       27       27

## conditions
## antimony     none 
##       27       27

hs_macr_drug_de

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

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

## The primary analysis performed 1 comparisons.

## The logFC agreement among the methods follows:

##                 nn_vs_ntmn
## basic_vs_deseq      0.8930
## basic_vs_dream      0.8901
## basic_vs_ebseq      0.8651
## basic_vs_edger      0.8942
## basic_vs_limma      0.8948
## basic_vs_noiseq     0.9048
## deseq_vs_dream      0.9951
## deseq_vs_ebseq      0.9647
## deseq_vs_edger      0.9997
## deseq_vs_limma      0.9911
## deseq_vs_noiseq     0.9837
## dream_vs_ebseq      0.9725
## dream_vs_edger      0.9952
## dream_vs_limma      0.9961
## dream_vs_noiseq     0.9880
## ebseq_vs_edger      0.9643
## ebseq_vs_limma      0.9694
## ebseq_vs_noiseq     0.9891
## edger_vs_limma      0.9911
## edger_vs_noiseq     0.9837
## limma_vs_noiseq     0.9852

hs_macr_drug_table <- combine_de_tables(
  hs_macr_drug_de, keepers = tmrc2_drug_keepers,
  excel = glue("analyses/macrophage_de/de_tables/macrophage_onlydrug_table-v{ver}.xlsx"))

## Looking for subscript invalid names, start of extract_keepers.

## Error in combine_extracted_plots(entry_name, combined, wanted_denominator, : object 'color_low' not found

hs_macr_drug_table

## Error: object 'hs_macr_drug_table' not found

#combined_to_tsv(hs_macr_drug_table, celltype = "macrophage")

hs_macr_drug_sig <- extract_significant_genes(
  hs_macr_drug_table,
  excel = glue("analyses/macrophage_de/sig_tables/macrophageonly_drug_sig-v{ver}.xlsx"))

## Error: object 'hs_macr_drug_table' not found

hs_macr_drug_sig

## Error: object 'hs_macr_drug_sig' not found

hs_macr_drug_highsig <- extract_significant_genes(
  hs_macr_drug_table, min_mean_exprs = high_expression, exprs_column = high_expression_column,
  excel = glue("analyses/macrophage_de/sig_tables/macrophageonly_drug_highsig-v{ver}.xlsx"))

## Error: object 'hs_macr_drug_table' not found

hs_macr_drug_highsig

## Error: object 'hs_macr_drug_highsig' not found

## Creating the following to see how it affects gProfiler.
hs_macr_drug_lesssig <- extract_significant_genes(
  hs_macr_drug_table, lfc = 0.6,
  excel = glue("analyses/macrophage_de/sig_tables/macrophageonly_drug_sig_lfc0.6-v{ver}.xlsx"))

## Error: object 'hs_macr_drug_table' not found

3.2.2 Macrophage Strain only

In a similar fashion, let us look for effects which are observed when we consider only the strain used during infection.

hs_macr_strain_de <- all_pairwise(hs_macr_strain_se, filter = TRUE, model_svs = "svaseq",
                                  model_fstring = "~ 0 + condition")

## z22 z23 
##  23  23

## Removing 9761 low-count genes (11720 remaining).

## Basic step 0/3: Normalizing data.

## Basic step 0/3: Converting data.

## I think this is failing? SummarizedExperiment

## Basic step 0/3: Transforming data.

## Setting 32467 entries to zero.

## This received a matrix of SVs.

## converting counts to integer mode

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## Warning in createContrastL(objFlt$formula, objFlt$data, L): Contrasts with only
## a single non-zero term are already evaluated by default.

## conditions
## z22 z23 
##  23  23

## conditions
## z22 z23 
##  23  23

## conditions
## z22 z23 
##  23  23

hs_macr_strain_de

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

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

## The primary analysis performed 1 comparisons.

## The logFC agreement among the methods follows:

##                 z23_vs_z22
## basic_vs_deseq      0.8273
## basic_vs_dream      0.8408
## basic_vs_ebseq      0.8021
## basic_vs_edger      0.8313
## basic_vs_limma      0.8629
## basic_vs_noiseq     0.1082
## deseq_vs_dream      0.9850
## deseq_vs_ebseq      0.9721
## deseq_vs_edger      0.9991
## deseq_vs_limma      0.9637
## deseq_vs_noiseq     0.2183
## dream_vs_ebseq      0.9734
## dream_vs_edger      0.9876
## dream_vs_limma      0.9782
## dream_vs_noiseq     0.2520
## ebseq_vs_edger      0.9726
## ebseq_vs_limma      0.9614
## ebseq_vs_noiseq     0.4180
## edger_vs_limma      0.9668
## edger_vs_noiseq     0.2127
## limma_vs_noiseq     0.2754

hs_macr_strain_table <- combine_de_tables(
  hs_macr_strain_de, keepers = tmrc2_strain_keepers,
  excel = glue("analyses/macrophage_de/de_tables/macrophage_onlystrain_table-v{ver}.xlsx"))

## Looking for subscript invalid names, start of extract_keepers.

## Error in combine_extracted_plots(entry_name, combined, wanted_denominator, : object 'color_low' not found

hs_macr_strain_table

## Error: object 'hs_macr_strain_table' not found

combined_to_tsv(hs_macr_strain_table, celltype = "macrophage")

## Error: object 'hs_macr_strain_table' not found

hs_macr_strain_sig <- extract_significant_genes(
  hs_macr_strain_table,
  excel = glue("analyses/macrophage_de/sig_tables/macrophageonly_onlystrain_sig-v{ver}.xlsx"))

## Error: object 'hs_macr_strain_table' not found

hs_macr_strain_sig

## Error: object 'hs_macr_strain_sig' not found

hs_macr_strain_highsig <- extract_significant_genes(
  hs_macr_strain_table, min_mean_exprs = high_expression, exprs_column = high_expression_column,
  excel = glue("analyses/macrophage_de/sig_tables/macrophageonly_onlystrain_highsig-v{ver}.xlsx"))

## Error: object 'hs_macr_strain_table' not found

hs_macr_strain_highsig

## Error: object 'hs_macr_strain_highsig' not found

hs_macr_strain_lesssig <- extract_significant_genes(
  hs_macr_strain_table, lfc = 0.6,
  excel = glue("analyses/macrophage_de/sig_tables/macrophageonly_onlystrain_lesssig-v{ver}.xlsx"))

## Error: object 'hs_macr_strain_table' not found

hs_macr_strain_lesssig

## Error: object 'hs_macr_strain_lesssig' not found

3.2.3 Compare Drug and Strain Effects

Now let us consider the above two comparisons together. First, I will plot the logFC values of them against each other (drug on x-axis and strain on the y-axis). Then we can extract the significant genes in a few combined categories of interest. I assume these will focus exclusively on the categories which include the introduction of the drug.

drug_strain_comp_df <- merge(hs_macr_drug_table[["data"]][["drug"]],
                             hs_macr_strain_table[["data"]][["strain"]],
                             by = "row.names")

## Error: object 'hs_macr_drug_table' not found

drug_strain_comp_plot <- plot_linear_scatter(
  drug_strain_comp_df[, c("deseq_logfc.x", "deseq_logfc.y")])

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'as.data.frame': object 'drug_strain_comp_df' not found

## Contrasts: antimony/none, z23/z22; x-axis: drug, y-axis: strain
## top left: higher no drug, z23; top right: higher drug z23
## bottom left: higher no drug, z22; bottom right: higher drug z22
drug_strain_comp_plot[["scatter"]]

## Error: object 'drug_strain_comp_plot' not found

As I noted in the comments above, some quadrants of the scatter plot are likely to be of greater interest to us than others (the right side). Because I get confused sometimes, the following block will explicitly name the categories of likely interest, then ask which genes are shared among them, and finally use UpSetR to extract the various gene intersection/union categories.

higher_drug <- hs_macr_drug_sig[["deseq"]][["downs"]][[1]]

## Error: object 'hs_macr_drug_sig' not found

higher_nodrug <- hs_macr_drug_sig[["deseq"]][["ups"]][[1]]

## Error: object 'hs_macr_drug_sig' not found

higher_z23 <- hs_macr_strain_sig[["deseq"]][["ups"]][[1]]

## Error: object 'hs_macr_strain_sig' not found

higher_z22 <- hs_macr_strain_sig[["deseq"]][["downs"]][[1]]

## Error: object 'hs_macr_strain_sig' not found

sum(rownames(higher_drug) %in% rownames(higher_z23))

## Error: object 'higher_drug' not found

sum(rownames(higher_drug) %in% rownames(higher_z22))

## Error: object 'higher_drug' not found

sum(rownames(higher_nodrug) %in% rownames(higher_z23))

## Error: object 'higher_nodrug' not found

sum(rownames(higher_nodrug) %in% rownames(higher_z22))

## Error: object 'higher_nodrug' not found

drug_z23_lst <- list("drug" = rownames(higher_drug),
                     "z23" = rownames(higher_z23))

## Error: object 'higher_drug' not found

upset_input <- UpSetR::fromList(drug_z23_lst)

## Error: object 'drug_z23_lst' not found

higher_drug_z23 <- upset(upset_input, text.scale = 2)

## Error: object 'upset_input' not found

higher_drug_z23

## Error: object 'higher_drug_z23' not found

drug_z23_shared_genes <- overlap_groups(drug_z23_lst)

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

shared_genes_drug_z23 <- overlap_geneids(drug_z23_shared_genes, "drug:z23")

## Error: object 'drug_z23_shared_genes' not found

shared_genes_drug_z23 <- attr(drug_z23_shared_genes, "elements")[drug_z23_shared_genes[["drug:z23"]]]

## Error: object 'drug_z23_shared_genes' not found

drug_z22_lst <- list("drug" = rownames(higher_drug),
                     "z22" = rownames(higher_z22))

## Error: object 'higher_drug' not found

higher_drug_z22 <- upset(UpSetR::fromList(drug_z22_lst), text.scale = 2)

## Error: object 'drug_z22_lst' not found

higher_drug_z22

## Error: object 'higher_drug_z22' not found

drug_z22_shared_genes <- overlap_groups(drug_z22_lst)

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

shared_genes_drug_z22 <- overlap_geneids(drug_z22_shared_genes, "drug:z22")

## Error: object 'drug_z22_shared_genes' not found

shared_genes_drug_z22 <- attr(drug_z22_shared_genes, "elements")[drug_z22_shared_genes[["drug:z22"]]]

## Error: object 'drug_z22_shared_genes' not found

3.2.4 Perform gProfiler on drug/strain effect shared genes

Now that we have some populations of genes which are shared across the drug/strain effects, let us pass them to some GSEA analyses and see what pops out.

wanted <- drug_z23_shared_genes[["drug:z23"]]

## Error: object 'drug_z23_shared_genes' not found

shared_genes_drug_z23 <- attr(drug_z23_shared_genes, "elements")[wanted]

## Error: object 'drug_z23_shared_genes' not found

shared_drug_z23_gp <- simple_gprofiler(shared_genes_drug_z23)

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

shared_drug_z23_gp[["pvalue_plots"]][["MF"]]

## Error: object 'shared_drug_z23_gp' not found

shared_drug_z23_gp[["pvalue_plots"]][["BP"]]

## Error: object 'shared_drug_z23_gp' not found

shared_drug_z23_gp[["pvalue_plots"]][["REAC"]]

## Error: object 'shared_drug_z23_gp' not found

wanted <- drug_z22_shared_genes[["drug:z22"]]

## Error: object 'drug_z22_shared_genes' not found

shared_genes_drug_z22 <- attr(drug_z22_shared_genes, "elements")[wanted]

## Error: object 'drug_z22_shared_genes' not found

shared_drug_z22_gp <- simple_gprofiler(shared_genes_drug_z22)

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

shared_drug_z22_gp[["pvalue_plots"]][["BP"]]

## Error: object 'shared_drug_z22_gp' not found

4 Our main question of interest

The data structure hs_macr contains our primary macrophages, which are, as shown above, the data we can really sink our teeth into.

Note, we expect some errors when running the combine_de_tables() because not all methods I use are comfortable using the ratio or ratios contrasts we added in the ‘extras’ argument. As a result, when we combine them into the larger output tables, those peculiar contrasts fail. This does not stop it from writing the rest of the results, however.

hs_macr_de_noextra <- all_pairwise(hs_macr, model_svs = "svaseq",
                                   model_fstring = "~ 0 + condition", filter = TRUE)

##    inf_sb_z22    inf_sb_z23       inf_z22       inf_z23    uninf_none 
##            12            11            11            12             4 
## uninf_sb_none 
##             4

## Removing 9725 low-count genes (11756 remaining).

## Basic step 0/3: Normalizing data.

## Basic step 0/3: Converting data.

## I think this is failing? SummarizedExperiment

## Basic step 0/3: Transforming data.

## Setting 40036 entries to zero.

## This received a matrix of SVs.

## converting counts to integer mode

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## Warning in createContrastL(objFlt$formula, objFlt$data, L): Contrasts with only
## a single non-zero term are already evaluated by default.

## conditions
##    inf_sb_z22    inf_sb_z23       inf_z22       inf_z23    uninf_none 
##            12            11            11            12             4 
## uninf_sb_none 
##             4

## conditions
##    inf_sb_z22    inf_sb_z23       inf_z22       inf_z23    uninf_none 
##            12            11            11            12             4 
## uninf_sb_none 
##             4

## conditions
##    inf_sb_z22    inf_sb_z23       inf_z22       inf_z23    uninf_none 
##            12            11            11            12             4 
## uninf_sb_none 
##             4

hs_macr_de <- all_pairwise(hs_macr, model_svs = "svaseq", model_fstring = "~ 0 + condition",
                           filter = TRUE, extra_contrasts = tmrc2_human_extra)

##    inf_sb_z22    inf_sb_z23       inf_z22       inf_z23    uninf_none 
##            12            11            11            12             4 
## uninf_sb_none 
##             4

## Removing 9725 low-count genes (11756 remaining).

## Basic step 0/3: Normalizing data.

## Basic step 0/3: Converting data.

## I think this is failing? SummarizedExperiment

## Basic step 0/3: Transforming data.

## Setting 40036 entries to zero.

## This received a matrix of SVs.

## converting counts to integer mode

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## The contrast z23drugnodrug is not in the results.

## If this is not an extra contrast, then this is an error.

## The contrast z23z22drug is not in the results.

## If this is not an extra contrast, then this is an error.

## Warning in createContrastL(objFlt$formula, objFlt$data, L): Contrasts with only
## a single non-zero term are already evaluated by default.

## conditions
##    inf_sb_z22    inf_sb_z23       inf_z22       inf_z23    uninf_none 
##            12            11            11            12             4 
## uninf_sb_none 
##             4

## conditions
##    inf_sb_z22    inf_sb_z23       inf_z22       inf_z23    uninf_none 
##            12            11            11            12             4 
## uninf_sb_none 
##             4

## conditions
##    inf_sb_z22    inf_sb_z23       inf_z22       inf_z23    uninf_none 
##            12            11            11            12             4 
## uninf_sb_none 
##             4

hs_macr_de

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

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

## The primary analysis performed 15 comparisons.

Write out the results.

hs_single_table <- combine_de_tables(
  hs_macr_de_noextra, keepers = single_tmrc2_keeper,
  excel = glue("analyses/macrophage_de/de_tables/hs_macr_drug_zymo_z22sb_sb-v{ver}.xlsx"))

## Looking for subscript invalid names, start of extract_keepers.

## Error in combine_extracted_plots(entry_name, combined, wanted_denominator, : object 'color_low' not found

hs_single_table

## Error: object 'hs_single_table' not found

hs_macr_table <- combine_de_tables(
  hs_macr_de, keepers = tmrc2_human_keepers,
  excel = glue("analyses/macrophage_de/de_tables/hs_macr_drug_zymo_table_macr_only-v{ver}.xlsx"))

## Looking for subscript invalid names, start of extract_keepers.

## Error in combine_extracted_plots(entry_name, combined, wanted_denominator, : object 'color_low' not found

hs_macr_table

## Error: object 'hs_macr_table' not found

combined_to_tsv(hs_macr_table, "macrophage")

## Error: object 'hs_macr_table' not found

hs_macr_sig <- extract_significant_genes(
  hs_macr_table,
  excel = glue("analyses/macrophage_de/sig_tables/hs_macr_drug_zymo_sig-v{ver}.xlsx"))

## Error: object 'hs_macr_table' not found

hs_macr_sig

## Error: object 'hs_macr_sig' not found

hs_macr_highsig <- extract_significant_genes(
  hs_macr_table, min_mean_exprs = high_expression, exprs_column = high_expression_column,
  excel = glue("analyses/macrophage_de/sig_tables/hs_macr_drug_zymo_highsig-v{ver}.xlsx"))

## Error: object 'hs_macr_table' not found

hs_macr_highsig

## Error: object 'hs_macr_highsig' not found

hs_macr_lesssig <- extract_significant_genes(
  hs_macr_table, lfc = 0.6,
  excel = glue("analyses/macrophage_de/sig_tables/hs_macr_drug_zymo_sig_lfc0.6-v{ver}.xlsx"))

## Error: object 'hs_macr_table' not found

hs_macr_lesssig

## Error: object 'hs_macr_lesssig' not found

4.1 gene group upset

4.1.1 2.3 vs 2.2 up and down vs. uninfected

This is my version of the Venn diagram which includes the text:

“Differentially expressed genes in macrophages infected with subpopulations 2.2 or 2.3. Volcano plots contrast of: A. Venn diagram for upregulated and downregulated genes by infection with 2.3 and 2.2 strains. B. infected cells with 2.3 strains and uninfected cells; C. infected cells with 2.2 strains and uninfected cells; D. infected cells with 2.3 strains and infected cells with 2.2 strains”

The following upset plot is currently Figure 2E.

nodrug_upset <- upsetr_combined_de(hs_macr_table,
                                   desired_contrasts = c("z22nosb_vs_uninf", "z23nosb_vs_uninf"))

## Error: object 'hs_macr_table' not found

pp(file = "images/nodrug_upset.svg")
nodrug_upset[["plot"]]

## Error: object 'nodrug_upset' not found

dev.off()

## png 
##   2

nodrug_upset

## Error: object 'nodrug_upset' not found

4.1.1.1 A point of interest while Olga visits Umd

Najib and Olga asked about pulling the 9 gene IDs which are in the peculiar situation of increased expression in z2.2/uninf and decreased in z2.3/uninf. In the previous upset plot, these are visible in the 6th bar. I can access these via the attr() function, which I should admit I can never remember how to use, so I am going to use the code under the ‘Compare(no)Sb z2.3/z2.2 treatment’ heading to remember how to extract these genes.

all_groups <- nodrug_upset[["groups"]]

## Error: object 'nodrug_upset' not found

wanted_group <- "z23nosb_vs_uninf_down:z22nosb_vs_uninf_up"
gene_idx <- all_groups[[wanted_group]]

## Error: object 'all_groups' not found

wanted_genes <- attr(all_groups, "elements")[gene_idx]

## Error: object 'all_groups' not found

wanted_genes

## Error: object 'wanted_genes' not found

gene_symbol_idx <- rownames(rowData(hs_macr)) %in% as.character(wanted_genes)

## Error: object 'wanted_genes' not found

rowData(hs_macr)[gene_symbol_idx, "hgnc_symbol"]

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

• ABCB5: ATB Binding Cassette Subfamily B Member #5, wide range of functions in this diverse paralogous family. Associated with skin diseases (melanoma and Epidermolysis Bullosa; participate in ATP-dependent transmembrane transport).
• RFX4: Regulatory Factor X #4: transcription factor.
• CA14: Carbonic anhydrase #14: Zync metalloenzyme catalyzes reversible hydration of CO2. This gene looks pretty neat, but not really relevant to anything we are likely to care about.
• EGR1: Early Growth Response Protein #1: Another Tx factor (zinc-finger) – important for cell survival/proliferation/cell death. Presumably important for healing?
• MCF2L: MCF.2 Cell Line Derived Transforming Sequence Like? guanine nucleotide exchange factor interacting with GTP-bound Rac1. Apparently associated with ostroarthritis; potentially relevant to regulation of RHOA and CDC42 signalling.
• DNASE1L3: Deoxyribonuclease I family member: not inhibited by actin, breaks down DNA during apoptosis. Important during necrosis.
• FOS: Proto-Oncogene, AP-1 Transcription Factor: leucine zipper dimerizes with JUN family proteins, forming tx factor complex AP-1. Important for cell proliferation, differentiation, and transformation.
• IFITM10: Interferon-Induced Transmembrane Protein #10
• PKD1L3: Polycystin 1 Like #3, Transient Receptor Potential Channel Interacting: 11 transmembrane domain protein which might help create cation channels.

As some comparison points, the Venn in the current figure has:

• 387 up z2.3
• 259 up z2.2
• 83 shared up z2.3 and z2.2
• 247 down z2.3
• 3 down z2.2
• 3 shared down z2.3 and z2.2

4.1.2 2.2 and 2.3 with SbV vs 2.2 and 2.3 without SbV

This is my version of the Venn with the text:

“Differentially expressed genes in macrophages infected with subpopulations 2.2 or 2.3, in presence of SbV. Volcano plots contrast of: A. infected cells with 2.3 strains + SbV and infected cells with 2.3 strains; B. infected cells with 2.2 strains + SbV and infected cells with 2.2 strains; C. infected cells with 2.3 strains + SbV and infected cells with 2.2 strains + SbV. D. Venn diagram for upregulated and downregulated genes by infection with 2.3+SbV and 2.2+SbV strains.”

A query from Olga (20240801): Please include in the upset in figure 3 the contrast of uninfected cells + SbV vs uninfected without SbV.

## I keep mis-interpreting this text, it is z2.3/z2.3SbV and z2.2/z2.2SbV
drugnodrug_upset <- upsetr_combined_de(hs_macr_table,
                                       desired_contrasts = c("z23sb_vs_z23nosb", "z22sb_vs_z22nosb"))

## Error: object 'hs_macr_table' not found

pp(file = "images/drugnodrug_upset.pdf")
drugnodrug_upset[["plot"]]

## Error: object 'drugnodrug_upset' not found

dev.off()

## png 
##   2

drugnodrug_upset

## Error: object 'drugnodrug_upset' not found

drugnodrug_uninf_contrasts <- c("z23sb_vs_z23nosb", "z22sb_vs_z22nosb", "sb_vs_uninf")
drugnodrug_upset_with_uninf <- upsetr_combined_de(hs_macr_table,
                                       desired_contrasts = drugnodrug_uninf_contrasts)

## Error: object 'hs_macr_table' not found

pp(file = "figures/drugnodrug_with_uninf_upset.svg")
drugnodrug_upset_with_uninf[["plot"]]

## Error: object 'drugnodrug_upset_with_uninf' not found

dev.off()

## png 
##   2

drugnodrug_upset_with_uninf

## Error: object 'drugnodrug_upset_with_uninf' not found

For some comparison points, the venn image has:

• 222 up z2.3 SbV
• 134 up z2.2 SbV
• 182 down z2.3 SbV
• 396 down z2.2 SbV
• 605 shared down z2.2 and z2.3 SbV
• 34 shared down z2.2 SbV and up z2.3 SbV
• 363 shared up z2.2 SbV and z2.3 SbV

4.1.3 Compare z2.2SbV vs SbV and z2.3SbV and SbV

drug_upset <- upsetr_combined_de(hs_macr_table,
                                 desired_contrasts = c("z22sb_vs_sb", "z23sb_vs_sb"))

## Error: object 'hs_macr_table' not found

pp(file = "images/drug_upset.pdf")
drug_upset[["plot"]]

## Error: object 'drug_upset' not found

dev.off()

## png 
##   2

drug_upset

## Error: object 'drug_upset' not found

4.2 Significance barplot of interest

Olga kindly sent a set of particularly interesting contrasts and colors for a significance barplot, they include the following:

• z2.3 vs. uninfected.
• z2.2 vs. uninfected.
• z2.3 vs z2.2
• z2.3Sbv vs z2.3
• z2.2Sbv vs z2.2
• z2.3Sbv vs z2.2Sbv
• Sbv vs uninfected.

The existing set of ‘keepers’ exvised to these is taken from the extant set of ‘tmrc2_human_keepers’ and is as follows:

barplot_keepers <- list(
  ## z2.3 vs uninfected
  "z23nosb_vs_uninf" = c("inf_z23", "uninf_none"),
  ## z2.2 vs uninfected
  "z22nosb_vs_uninf" = c("inf_z22", "uninf_none"),
  ## z2.3 vs z2.2
  "z23nosb_vs_z22nosb" = c("inf_z23", "inf_z22"),
  ## z2.3Sbv vs z2.3
  "z23sb_vs_z23nosb" = c("inf_sb_z23", "inf_z23"),
  ## z2.2Sbv vs z2.2
  "z22sb_vs_z22nosb" = c("inf_sb_z22", "inf_z22"),
  ## z2.3Sbv vs z2.2Sbv
  "z23sb_vs_z22sb" = c("inf_sb_z23", "inf_sb_z22"),
  ## Sbv vs uninfected.
  "sb_vs_uninf" = c("uninf_sb_none", "uninf_none"))
barplot_combined <- combine_de_tables(
  hs_macr_de, keepers = barplot_keepers,
  excel = glue("analyses/macrophage_de/de_tables/hs_macr_drug_zymo_7contrasts-v{ver}.xlsx"))

## Looking for subscript invalid names, start of extract_keepers.

## Error in combine_extracted_plots(entry_name, combined, wanted_denominator, : object 'color_low' not found

Now let us use the colors suggested by Olga to make a barplot of these…

color_list <-  c( "#de8bf9", "#ad07e3","#410257", "#ffa0a0", "#f94040", "#a00000")
barplot_sig <- extract_significant_genes(
  barplot_combined, color_list = color_list, according_to = "deseq",
  excel = glue("analyses/macrophage_de/sig_tables/hs_macr_drug_zymo_7contrasts_sig-v{ver}.xlsx"))

## Error: object 'barplot_combined' not found

barplot_sig

## Error: object 'barplot_sig' not found

5 PROPER

In our last meeting there were some questions about the statistical power of different future experimental designs. One thing I can do is to use PROPER to estimate the power of an extant dataset and infer from that the likely power of other designs.

In order to use proper, one must feed it one or more DE tables.

power_estimate <- simple_proper(hs_single_table)

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

power_estimate[[1]][["power_plot"]]

## Error: object 'power_estimate' not found

power_estimate[[1]][["powertd_plot"]]

## Error: object 'power_estimate' not found

power_estimate[[1]][["powerfd_plot"]]

## Error: object 'power_estimate' not found

6 Our main questions in U937

Let us do the same comparisons in the U937 samples, though I will not do the extra contrasts, primarily because I think the dataset is less likely to support them.

u937_de <- all_pairwise(u937_se, model_svs = "svaseq",
                        filter = TRUE, model_fstring = "~ 0 + condition")

##    inf_sb_z22    inf_sb_z23       inf_z22       inf_z23    uninf_none 
##             3             3             3             3             1 
## uninf_sb_none 
##             1

## Removing 10730 low-count genes (10751 remaining).

## Basic step 0/3: Normalizing data.

## Basic step 0/3: Converting data.

## I think this is failing? SummarizedExperiment

## Basic step 0/3: Transforming data.

## Setting 938 entries to zero.

## This received a matrix of SVs.

## converting counts to integer mode

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## Warning in createContrastL(objFlt$formula, objFlt$data, L): Contrasts with only
## a single non-zero term are already evaluated by default.

## conditions
##    inf_sb_z22    inf_sb_z23       inf_z22       inf_z23    uninf_none 
##             3             3             3             3             1 
## uninf_sb_none 
##             1

## conditions
##    inf_sb_z22    inf_sb_z23       inf_z22       inf_z23    uninf_none 
##             3             3             3             3             1 
## uninf_sb_none 
##             1

## conditions
##    inf_sb_z22    inf_sb_z23       inf_z22       inf_z23    uninf_none 
##             3             3             3             3             1 
## uninf_sb_none 
##             1 
## Error in NOISeq::noiseqbio(norm_input, k = k, norm = norm, factor = condition_column,  : 
##   ERROR: To run NOISeqBIO at least two replicates per condition are needed.
##          Please, run NOISeq if there are not enough replicates in your experiment.

u937_de

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

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

## The primary analysis performed 15 comparisons.

u937_table <- combine_de_tables(
  u937_de, keepers = u937_keepers,
  excel = glue("analyses/macrophage_de/de_tables/u937_drug_zymo_table-v{ver}.xlsx"))

## Looking for subscript invalid names, start of extract_keepers.

## Error in combine_extracted_plots(entry_name, combined, wanted_denominator, : object 'color_low' not found

u937_table

## Error: object 'u937_table' not found

combined_to_tsv(u937_table, celltype = "u937")

## Error: object 'u937_table' not found

u937_sig <- extract_significant_genes(
  u937_table,
  excel = glue("analyses/macrophage_de/sig_tables/u937_drug_zymo_sig-v{ver}.xlsx"))

## Error: object 'u937_table' not found

u937_sig

## Error: object 'u937_sig' not found

u937_highsig <- extract_significant_genes(
  u937_table, min_mean_exprs = high_expression, exprs_column = high_expression_column,
  excel = glue("analyses/macrophage_de/sig_tables/u937_drug_zymo_highsig-v{ver}.xlsx"))

## Error: object 'u937_table' not found

u937_highsig

## Error: object 'u937_highsig' not found

u937_lesssig <- extract_significant_genes(
  u937_table, lfc = 0.6,
  excel = glue("analyses/macrophage_de/sig_tables/u937_drug_zymo_lesssig-v{ver}.xlsx"))

## Error: object 'u937_table' not found

u937_lesssig

## Error: object 'u937_lesssig' not found

7 Compare (no)Sb z2.3/z2.2 treatments among macrophages

In the following block, I will jump back to the macrophage samples and look for genes which are shared/unique when comparing z2.3/z2.2 for the drug treated samples and the untreated samples.

upset_plots_hs_macr <- upsetr_sig(
  hs_macr_sig, both = TRUE,
  contrasts = c("z23sb_vs_z22sb", "z23nosb_vs_z22nosb"))

## Error: object 'hs_macr_sig' not found

upset_plots_hs_macr[["both"]]

## Error: object 'upset_plots_hs_macr' not found

groups <- upset_plots_hs_macr[["both_groups"]]

## Error: object 'upset_plots_hs_macr' not found

shared_genes <- attr(groups, "elements")[groups[[2]]] %>%
  gsub(pattern = "^gene:", replacement = "")

## Error: object 'groups' not found

length(shared_genes)

## Error: object 'shared_genes' not found

shared_gp <- simple_gprofiler(shared_genes)

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

shared_gp[["pvalue_plots"]][["MF"]]

## Error: object 'shared_gp' not found

shared_gp[["pvalue_plots"]][["BP"]]

## Error: object 'shared_gp' not found

shared_gp[["pvalue_plots"]][["REAC"]]

## Error: object 'shared_gp' not found

drug_genes <- attr(groups, "elements")[groups[["z23sb_vs_z22sb"]]] %>%
  gsub(pattern = "^gene:", replacement = "")

## Error: object 'groups' not found

drugonly_gp <- simple_gprofiler(drug_genes)

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

drugonly_gp[["pvalue_plots"]][["BP"]]

## Error: object 'drugonly_gp' not found

I want to try something, directly include the u937 data in this. Thus, in the following block I will repeat but compare all samples and the U937 using the same logic.

both_sig <- hs_macr_sig

## Error: object 'hs_macr_sig' not found

names(both_sig[["deseq"]][["ups"]]) <- paste0("macr_", names(both_sig[["deseq"]][["ups"]]))

## Error: object 'both_sig' not found

names(both_sig[["deseq"]][["downs"]]) <- paste0("macr_", names(both_sig[["deseq"]][["downs"]]))

## Error: object 'both_sig' not found

u937_deseq <- u937_sig[["deseq"]]

## Error: object 'u937_sig' not found

names(u937_deseq[["ups"]]) <- paste0("u937_", names(u937_deseq[["ups"]]))

## Error: object 'u937_deseq' not found

names(u937_deseq[["downs"]]) <- paste0("u937_", names(u937_deseq[["downs"]]))

## Error: object 'u937_deseq' not found

both_sig[["deseq"]][["ups"]] <- c(both_sig[["deseq"]][["ups"]], u937_deseq[["ups"]])

## Error: object 'both_sig' not found

both_sig[["deseq"]][["downs"]] <- c(both_sig[["deseq"]][["ups"]], u937_deseq[["downs"]])

## Error: object 'both_sig' not found

summary(both_sig[["deseq"]][["ups"]])

## Error: object 'both_sig' not found

upset_plots_both <- upsetr_sig(
  both_sig, both = TRUE,
  contrasts = c("macr_z23sb_vs_z22sb", "macr_z23nosb_vs_z22nosb",
                "u937_z23sb_vs_z22sb", "u937_z23nosb_vs_z22nosb"))

## Error: object 'both_sig' not found

upset_plots_both[["both"]]

## Error: object 'upset_plots_both' not found

7.1 Compare DE results from macrophages and U937 samples

Looking a bit more closely at these, I think the u937 data is too sparse to effectively compare.

macr_u937_comparison <- compare_de_results(hs_macr_table, u937_table)

## Error: object 'hs_macr_table' not found

macr_u937_comparison[["lfc_heat"]]

## Error: object 'macr_u937_comparison' not found

macr_u937_venns <- compare_significant_contrasts(hs_macr_sig, second_sig_tables = u937_sig,
                                                 contrasts = "z23sb_vs_z23nosb")

## Error: object 'hs_macr_sig' not found

macr_u937_venns[["up_plot"]]

## Error: object 'macr_u937_venns' not found

macr_u937_venns[["down_plot"]]

## Error: object 'macr_u937_venns' not found

macr_u937_venns_v2 <- compare_significant_contrasts(
  hs_macr_sig, second_sig_tables = u937_sig, contrasts = "z22sb_vs_z22nosb")

## Error: object 'hs_macr_sig' not found

macr_u937_venns_v2[["up_plot"]]

## Error: object 'macr_u937_venns_v2' not found

macr_u937_venns_v2[["down_plot"]]

## Error: object 'macr_u937_venns_v2' not found

macr_u937_venns_v3 <- compare_significant_contrasts(
  hs_macr_sig, second_sig_tables = u937_sig, contrasts = "sb_vs_uninf")

## Error: object 'hs_macr_sig' not found

macr_u937_venns_v3[["up_plot"]]

## Error: object 'macr_u937_venns_v3' not found

macr_u937_venns_v3[["down_plot"]]

## Error: object 'macr_u937_venns_v3' not found

7.2 Compare macrophage/u937 with respect to z2.3/z2.2

comparison_df <- merge(hs_macr_table[["data"]][["z23sb_vs_z22sb"]],
                       u937_table[["data"]][["z23sb_vs_z22sb"]],
                       by = "row.names")

## Error: object 'hs_macr_table' not found

macru937_z23z22_plot <- plot_linear_scatter(comparison_df[, c("deseq_logfc.x", "deseq_logfc.y")])

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'as.data.frame': object 'comparison_df' not found

macru937_z23z22_plot[["scatter"]]

## Error: object 'macru937_z23z22_plot' not found

comparison_df <- merge(hs_macr_table[["data"]][["z23nosb_vs_z22nosb"]],
                       u937_table[["data"]][["z23nosb_vs_z22nosb"]],
                       by = "row.names")

## Error: object 'hs_macr_table' not found

macru937_z23z22_plot <- plot_linear_scatter(comparison_df[, c("deseq_logfc.x", "deseq_logfc.y")])

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'as.data.frame': object 'comparison_df' not found

macru937_z23z22_plot[["scatter"]]

## Error: object 'macru937_z23z22_plot' not found

8 Add donor to the contrasts, no sva

In the following block, I will change the sample condition to include the donor.

no_power_fact <- paste0(colData(hs_macr)[["donor"]], "_",
                        colData(hs_macr)[["condition"]])
table(colData(hs_macr)[["donor"]])

## 
## d01 d02 d09 d81 
##  13  14  13  14

table(no_power_fact)

## no_power_fact
##    d01_inf_sb_z22    d01_inf_sb_z23       d01_inf_z22       d01_inf_z23 
##                 3                 3                 2                 3 
##    d01_uninf_none d01_uninf_sb_none    d02_inf_sb_z22    d02_inf_sb_z23 
##                 1                 1                 3                 3 
##       d02_inf_z22       d02_inf_z23    d02_uninf_none d02_uninf_sb_none 
##                 3                 3                 1                 1 
##    d09_inf_sb_z22    d09_inf_sb_z23       d09_inf_z22       d09_inf_z23 
##                 3                 2                 3                 3 
##    d09_uninf_none d09_uninf_sb_none    d81_inf_sb_z22    d81_inf_sb_z23 
##                 1                 1                 3                 3 
##       d81_inf_z22       d81_inf_z23    d81_uninf_none d81_uninf_sb_none 
##                 3                 3                 1                 1

hs_nopower <- set_conditions(hs_macr, fact = no_power_fact)

## The numbers of samples by condition are:

## 
##    d01_inf_sb_z22    d01_inf_sb_z23       d01_inf_z22       d01_inf_z23 
##                 3                 3                 2                 3 
##    d01_uninf_none d01_uninf_sb_none    d02_inf_sb_z22    d02_inf_sb_z23 
##                 1                 1                 3                 3 
##       d02_inf_z22       d02_inf_z23    d02_uninf_none d02_uninf_sb_none 
##                 3                 3                 1                 1 
##    d09_inf_sb_z22    d09_inf_sb_z23       d09_inf_z22       d09_inf_z23 
##                 3                 2                 3                 3 
##    d09_uninf_none d09_uninf_sb_none    d81_inf_sb_z22    d81_inf_sb_z23 
##                 1                 1                 3                 3 
##       d81_inf_z22       d81_inf_z23    d81_uninf_none d81_uninf_sb_none 
##                 3                 3                 1                 1

hs_nopower <- subset_se(hs_nopower, subset = "macrophagezymodeme!='none'")
hs_nopower_nosva_de <- all_pairwise(hs_nopower, model_svs = FALSE, filter = TRUE, model_fstring = "~ 0 + condition")

## d01_inf_sb_z22 d01_inf_sb_z23    d01_inf_z22    d01_inf_z23 d02_inf_sb_z22 
##              3              3              2              3              3 
## d02_inf_sb_z23    d02_inf_z22    d02_inf_z23 d09_inf_sb_z22 d09_inf_sb_z23 
##              3              3              3              3              2 
##    d09_inf_z22    d09_inf_z23 d81_inf_sb_z22 d81_inf_sb_z23    d81_inf_z22 
##              3              3              3              3              3 
##    d81_inf_z23 
##              3

## Removing 9761 low-count genes (11720 remaining).

## Basic step 0/3: Normalizing data.

## Basic step 0/3: Converting data.

## I think this is failing? SummarizedExperiment

## Basic step 0/3: Transforming data.

## Setting 32467 entries to zero.

## converting counts to integer mode

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## Warning in createContrastL(objFlt$formula, objFlt$data, L): Contrasts with only
## a single non-zero term are already evaluated by default.

## conditions
## d01_inf_sb_z22 d01_inf_sb_z23    d01_inf_z22    d01_inf_z23 d02_inf_sb_z22 
##              3              3              2              3              3 
## d02_inf_sb_z23    d02_inf_z22    d02_inf_z23 d09_inf_sb_z22 d09_inf_sb_z23 
##              3              3              3              3              2 
##    d09_inf_z22    d09_inf_z23 d81_inf_sb_z22 d81_inf_sb_z23    d81_inf_z22 
##              3              3              3              3              3 
##    d81_inf_z23 
##              3

## conditions
## d01_inf_sb_z22 d01_inf_sb_z23    d01_inf_z22    d01_inf_z23 d02_inf_sb_z22 
##              3              3              2              3              3 
## d02_inf_sb_z23    d02_inf_z22    d02_inf_z23 d09_inf_sb_z22 d09_inf_sb_z23 
##              3              3              3              3              2 
##    d09_inf_z22    d09_inf_z23 d81_inf_sb_z22 d81_inf_sb_z23    d81_inf_z22 
##              3              3              3              3              3 
##    d81_inf_z23 
##              3

## conditions
## d01_inf_sb_z22 d01_inf_sb_z23    d01_inf_z22    d01_inf_z23 d02_inf_sb_z22 
##              3              3              2              3              3 
## d02_inf_sb_z23    d02_inf_z22    d02_inf_z23 d09_inf_sb_z22 d09_inf_sb_z23 
##              3              3              3              3              2 
##    d09_inf_z22    d09_inf_z23 d81_inf_sb_z22 d81_inf_sb_z23    d81_inf_z22 
##              3              3              3              3              3 
##    d81_inf_z23 
##              3

nopower_keepers <- list(
  "d01_zymo" = c("d01infz23", "d01infz22"),
  "d01_sbzymo" = c("d01infsbz23", "d01infsbz22"),
  "d02_zymo" = c("d02infz23", "d02infz22"),
  "d02_sbzymo" = c("d02infsbz23", "d02infsbz22"),
  "d09_zymo" = c("d09infz23", "d09infz22"),
  "d09_sbzymo" = c("d09infsbz23", "d09infsbz22"),
  "d81_zymo" = c("d81infz23", "d81infz22"),
  "d81_sbzymo" = c("d81infsbz23", "d81infsbz22"))
hs_nopower_nosva_table <- combine_de_tables(
  hs_nopower_nosva_de, keepers = nopower_keepers,
  excel = glue("analyses/macrophage_de/de_tables/hs_nopower_table-v{ver}.xlsx"))

## Looking for subscript invalid names, start of extract_keepers.

## The keepers has no elements in the coefficients.

## Here are the keepers: d01infz23, d01infz22, d01infsbz23, d01infsbz22, d02infz23, d02infz22, d02infsbz23, d02infsbz22, d09infz23, d09infz22, d09infsbz23, d09infsbz22, d81infz23, d81infz22, d81infsbz23, d81infsbz22

## Here are the coefficients: d81_inf_z23, d81_inf_z22, d81_inf_z23, d81_inf_sb_z23, d81_inf_z22, d81_inf_sb_z23, d81_inf_z23, d81_inf_sb_z22, d81_inf_z22, d81_inf_sb_z22, d81_inf_sb_z23, d81_inf_sb_z22, d81_inf_z23, d09_inf_z23, d81_inf_z22, d09_inf_z23, d81_inf_sb_z23, d09_inf_z23, d81_inf_sb_z22, d09_inf_z23, d81_inf_z23, d09_inf_z22, d81_inf_z22, d09_inf_z22, d81_inf_sb_z23, d09_inf_z22, d81_inf_sb_z22, d09_inf_z22, d09_inf_z23, d09_inf_z22, d81_inf_z23, d09_inf_sb_z23, d81_inf_z22, d09_inf_sb_z23, d81_inf_sb_z23, d09_inf_sb_z23, d81_inf_sb_z22, d09_inf_sb_z23, d09_inf_z23, d09_inf_sb_z23, d09_inf_z22, d09_inf_sb_z23, d81_inf_z23, d09_inf_sb_z22, d81_inf_z22, d09_inf_sb_z22, d81_inf_sb_z23, d09_inf_sb_z22, d81_inf_sb_z22, d09_inf_sb_z22, d09_inf_z23, d09_inf_sb_z22, d09_inf_z22, d09_inf_sb_z22, d09_inf_sb_z23, d09_inf_sb_z22, d81_inf_z23, d02_inf_z23, d81_inf_z22, d02_inf_z23, d81_inf_sb_z23, d02_inf_z23, d81_inf_sb_z22, d02_inf_z23, d09_inf_z23, d02_inf_z23, d09_inf_z22, d02_inf_z23, d09_inf_sb_z23, d02_inf_z23, d09_inf_sb_z22, d02_inf_z23, d81_inf_z23, d02_inf_z22, d81_inf_z22, d02_inf_z22, d81_inf_sb_z23, d02_inf_z22, d81_inf_sb_z22, d02_inf_z22, d09_inf_z23, d02_inf_z22, d09_inf_z22, d02_inf_z22, d09_inf_sb_z23, d02_inf_z22, d09_inf_sb_z22, d02_inf_z22, d02_inf_z23, d02_inf_z22, d81_inf_z23, d02_inf_sb_z23, d81_inf_z22, d02_inf_sb_z23, d81_inf_sb_z23, d02_inf_sb_z23, d81_inf_sb_z22, d02_inf_sb_z23, d09_inf_z23, d02_inf_sb_z23, d09_inf_z22, d02_inf_sb_z23, d09_inf_sb_z23, d02_inf_sb_z23, d09_inf_sb_z22, d02_inf_sb_z23, d02_inf_z23, d02_inf_sb_z23, d02_inf_z22, d02_inf_sb_z23, d81_inf_z23, d02_inf_sb_z22, d81_inf_z22, d02_inf_sb_z22, d81_inf_sb_z23, d02_inf_sb_z22, d81_inf_sb_z22, d02_inf_sb_z22, d09_inf_z23, d02_inf_sb_z22, d09_inf_z22, d02_inf_sb_z22, d09_inf_sb_z23, d02_inf_sb_z22, d09_inf_sb_z22, d02_inf_sb_z22, d02_inf_z23, d02_inf_sb_z22, d02_inf_z22, d02_inf_sb_z22, d02_inf_sb_z23, d02_inf_sb_z22, d81_inf_z23, d01_inf_z23, d81_inf_z22, d01_inf_z23, d81_inf_sb_z23, d01_inf_z23, d81_inf_sb_z22, d01_inf_z23, d09_inf_z23, d01_inf_z23, d09_inf_z22, d01_inf_z23, d09_inf_sb_z23, d01_inf_z23, d09_inf_sb_z22, d01_inf_z23, d02_inf_z23, d01_inf_z23, d02_inf_z22, d01_inf_z23, d02_inf_sb_z23, d01_inf_z23, d02_inf_sb_z22, d01_inf_z23, d81_inf_z23, d01_inf_z22, d81_inf_z22, d01_inf_z22, d81_inf_sb_z23, d01_inf_z22, d81_inf_sb_z22, d01_inf_z22, d09_inf_z23, d01_inf_z22, d09_inf_z22, d01_inf_z22, d09_inf_sb_z23, d01_inf_z22, d09_inf_sb_z22, d01_inf_z22, d02_inf_z23, d01_inf_z22, d02_inf_z22, d01_inf_z22, d02_inf_sb_z23, d01_inf_z22, d02_inf_sb_z22, d01_inf_z22, d01_inf_z23, d01_inf_z22, d81_inf_z23, d01_inf_sb_z23, d81_inf_z22, d01_inf_sb_z23, d81_inf_sb_z23, d01_inf_sb_z23, d81_inf_sb_z22, d01_inf_sb_z23, d09_inf_z23, d01_inf_sb_z23, d09_inf_z22, d01_inf_sb_z23, d09_inf_sb_z23, d01_inf_sb_z23, d09_inf_sb_z22, d01_inf_sb_z23, d02_inf_z23, d01_inf_sb_z23, d02_inf_z22, d01_inf_sb_z23, d02_inf_sb_z23, d01_inf_sb_z23, d02_inf_sb_z22, d01_inf_sb_z23, d01_inf_z23, d01_inf_sb_z23, d01_inf_z22, d01_inf_sb_z23, d81_inf_z23, d01_inf_sb_z22, d81_inf_z22, d01_inf_sb_z22, d81_inf_sb_z23, d01_inf_sb_z22, d81_inf_sb_z22, d01_inf_sb_z22, d09_inf_z23, d01_inf_sb_z22, d09_inf_z22, d01_inf_sb_z22, d09_inf_sb_z23, d01_inf_sb_z22, d09_inf_sb_z22, d01_inf_sb_z22, d02_inf_z23, d01_inf_sb_z22, d02_inf_z22, d01_inf_sb_z22, d02_inf_sb_z23, d01_inf_sb_z22, d02_inf_sb_z22, d01_inf_sb_z22, d01_inf_z23, d01_inf_sb_z22, d01_inf_z22, d01_inf_sb_z22, d01_inf_sb_z23, d01_inf_sb_z22

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

## extra_contrasts = extra)
hs_nopower_nosva_sig <- extract_significant_genes(
  hs_nopower_nosva_table,
  excel = glue("analyses/macrophage_de/sig_tables/hs_nopower_nosva_sig-v{ver}.xlsx"))

## Error: object 'hs_nopower_nosva_table' not found

d01d02_zymo_nosva_comp <- merge(hs_nopower_nosva_table[["data"]][["d01_zymo"]],
                                hs_nopower_nosva_table[["data"]][["d02_zymo"]],
                                by = "row.names")

## Error: object 'hs_nopower_nosva_table' not found

d0102_zymo_nosva_plot <- plot_linear_scatter(d01d02_zymo_nosva_comp[, c("deseq_logfc.x", "deseq_logfc.y")])

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'as.data.frame': object 'd01d02_zymo_nosva_comp' not found

d0102_zymo_nosva_plot[["scatter"]]

## Error: object 'd0102_zymo_nosva_plot' not found

d0102_zymo_nosva_plot[["correlation"]]

## Error: object 'd0102_zymo_nosva_plot' not found

d0102_zymo_nosva_plot[["lm_rsq"]]

## Error: object 'd0102_zymo_nosva_plot' not found

d09d81_zymo_nosva_comp <- merge(hs_nopower_nosva_table[["data"]][["d09_zymo"]],
                                hs_nopower_nosva_table[["data"]][["d81_zymo"]],
                                by = "row.names")

## Error: object 'hs_nopower_nosva_table' not found

d0981_zymo_nosva_plot <- plot_linear_scatter(d09d81_zymo_nosva_comp[, c("deseq_logfc.x", "deseq_logfc.y")])

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'as.data.frame': object 'd09d81_zymo_nosva_comp' not found

d0981_zymo_nosva_plot[["scatter"]]

## Error: object 'd0981_zymo_nosva_plot' not found

d0981_zymo_nosva_plot[["correlation"]]

## Error: object 'd0981_zymo_nosva_plot' not found

d0981_zymo_nosva_plot[["lm_rsq"]]

## Error: object 'd0981_zymo_nosva_plot' not found

d01d81_zymo_nosva_comp <- merge(hs_nopower_nosva_table[["data"]][["d01_zymo"]],
                                hs_nopower_nosva_table[["data"]][["d81_zymo"]],
                                by = "row.names")

## Error: object 'hs_nopower_nosva_table' not found

d0181_zymo_nosva_plot <- plot_linear_scatter(d01d81_zymo_nosva_comp[, c("deseq_logfc.x", "deseq_logfc.y")])

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'as.data.frame': object 'd01d81_zymo_nosva_comp' not found

d0181_zymo_nosva_plot[["scatter"]]

## Error: object 'd0181_zymo_nosva_plot' not found

d0181_zymo_nosva_plot[["correlation"]]

## Error: object 'd0181_zymo_nosva_plot' not found

d0181_zymo_nosva_plot[["lm_rsq"]]

## Error: object 'd0181_zymo_nosva_plot' not found

upset_plots_nosva <- upsetr_sig(hs_nopower_nosva_sig, both = TRUE,
                                contrasts = c("d01_zymo", "d02_zymo", "d09_zymo", "d81_zymo"))

## Error: object 'hs_nopower_nosva_sig' not found

upset_plots_nosva[["up"]]

## Error: object 'upset_plots_nosva' not found

upset_plots_nosva[["down"]]

## Error: object 'upset_plots_nosva' not found

upset_plots_nosva[["both"]]

## Error: object 'upset_plots_nosva' not found

## The 7th element in the both groups list is the set shared among all donors.
## I don't feel like writing out x:y:z:a
groups <- upset_plots_nosva[["both_groups"]]

## Error: object 'upset_plots_nosva' not found

shared_genes <- attr(groups, "elements")[groups[[7]]] %>%
  gsub(pattern = "^gene:", replacement = "")

## Error: object 'groups' not found

shared_gp <- simple_gprofiler(shared_genes)

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

shared_gp[["pvalue_plots"]][["MF"]]

## Error: object 'shared_gp' not found

shared_gp[["pvalue_plots"]][["BP"]]

## Error: object 'shared_gp' not found

shared_gp[["pvalue_plots"]][["REAC"]]

## Error: object 'shared_gp' not found

shared_gp[["pvalue_plots"]][["WP"]]

## Error: object 'shared_gp' not found

9 Add donor to the contrasts, sva

Same deal as the last block, but this time add SVA into the mix!

hs_nopower_sva_de <- all_pairwise(hs_nopower, model_svs = "svaseq",
                                  model_fstring = "~ 0 + condition", filter = TRUE)

## d01_inf_sb_z22 d01_inf_sb_z23    d01_inf_z22    d01_inf_z23 d02_inf_sb_z22 
##              3              3              2              3              3 
## d02_inf_sb_z23    d02_inf_z22    d02_inf_z23 d09_inf_sb_z22 d09_inf_sb_z23 
##              3              3              3              3              2 
##    d09_inf_z22    d09_inf_z23 d81_inf_sb_z22 d81_inf_sb_z23    d81_inf_z22 
##              3              3              3              3              3 
##    d81_inf_z23 
##              3

## Removing 9761 low-count genes (11720 remaining).

## Basic step 0/3: Normalizing data.

## Basic step 0/3: Converting data.

## I think this is failing? SummarizedExperiment

## Basic step 0/3: Transforming data.

## Setting 32467 entries to zero.

## This received a matrix of SVs.

## converting counts to integer mode

##   the design formula contains one or more numeric variables with integer values,
##   specifying a model with increasing fold change for higher values.
##   did you mean for this to be a factor? if so, first convert
##   this variable to a factor using the factor() function
##   the design formula contains one or more numeric variables with integer values,
##   specifying a model with increasing fold change for higher values.
##   did you mean for this to be a factor? if so, first convert
##   this variable to a factor using the factor() function

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## Warning in createContrastL(objFlt$formula, objFlt$data, L): Contrasts with only
## a single non-zero term are already evaluated by default.

## conditions
## d01_inf_sb_z22 d01_inf_sb_z23    d01_inf_z22    d01_inf_z23 d02_inf_sb_z22 
##              3              3              2              3              3 
## d02_inf_sb_z23    d02_inf_z22    d02_inf_z23 d09_inf_sb_z22 d09_inf_sb_z23 
##              3              3              3              3              2 
##    d09_inf_z22    d09_inf_z23 d81_inf_sb_z22 d81_inf_sb_z23    d81_inf_z22 
##              3              3              3              3              3 
##    d81_inf_z23 
##              3

## conditions
## d01_inf_sb_z22 d01_inf_sb_z23    d01_inf_z22    d01_inf_z23 d02_inf_sb_z22 
##              3              3              2              3              3 
## d02_inf_sb_z23    d02_inf_z22    d02_inf_z23 d09_inf_sb_z22 d09_inf_sb_z23 
##              3              3              3              3              2 
##    d09_inf_z22    d09_inf_z23 d81_inf_sb_z22 d81_inf_sb_z23    d81_inf_z22 
##              3              3              3              3              3 
##    d81_inf_z23 
##              3

## conditions
## d01_inf_sb_z22 d01_inf_sb_z23    d01_inf_z22    d01_inf_z23 d02_inf_sb_z22 
##              3              3              2              3              3 
## d02_inf_sb_z23    d02_inf_z22    d02_inf_z23 d09_inf_sb_z22 d09_inf_sb_z23 
##              3              3              3              3              2 
##    d09_inf_z22    d09_inf_z23 d81_inf_sb_z22 d81_inf_sb_z23    d81_inf_z22 
##              3              3              3              3              3 
##    d81_inf_z23 
##              3

nopower_keepers <- list(
  "d01_zymo" = c("d01infz23", "d01infz22"),
  "d01_sbzymo" = c("d01infsbz23", "d01infsbz22"),
  "d02_zymo" = c("d02infz23", "d02infz22"),
  "d02_sbzymo" = c("d02infsbz23", "d02infsbz22"),
  "d09_zymo" = c("d09infz23", "d09infz22"),
  "d09_sbzymo" = c("d09infsbz23", "d09infsbz22"),
  "d81_zymo" = c("d81infz23", "d81infz22"),
  "d81_sbzymo" = c("d81infsbz23", "d81infsbz22"))
hs_nopower_sva_table <- combine_de_tables(
  hs_nopower_sva_de, keepers = nopower_keepers,
  excel = glue("analyses/macrophage_de/de_tables/hs_nopower_table-v{ver}.xlsx"))

## Looking for subscript invalid names, start of extract_keepers.

## The keepers has no elements in the coefficients.

## Here are the keepers: d01infz23, d01infz22, d01infsbz23, d01infsbz22, d02infz23, d02infz22, d02infsbz23, d02infsbz22, d09infz23, d09infz22, d09infsbz23, d09infsbz22, d81infz23, d81infz22, d81infsbz23, d81infsbz22

## Here are the coefficients: d81_inf_z23, d81_inf_z22, d81_inf_z23, d81_inf_sb_z23, d81_inf_z22, d81_inf_sb_z23, d81_inf_z23, d81_inf_sb_z22, d81_inf_z22, d81_inf_sb_z22, d81_inf_sb_z23, d81_inf_sb_z22, d81_inf_z23, d09_inf_z23, d81_inf_z22, d09_inf_z23, d81_inf_sb_z23, d09_inf_z23, d81_inf_sb_z22, d09_inf_z23, d81_inf_z23, d09_inf_z22, d81_inf_z22, d09_inf_z22, d81_inf_sb_z23, d09_inf_z22, d81_inf_sb_z22, d09_inf_z22, d09_inf_z23, d09_inf_z22, d81_inf_z23, d09_inf_sb_z23, d81_inf_z22, d09_inf_sb_z23, d81_inf_sb_z23, d09_inf_sb_z23, d81_inf_sb_z22, d09_inf_sb_z23, d09_inf_z23, d09_inf_sb_z23, d09_inf_z22, d09_inf_sb_z23, d81_inf_z23, d09_inf_sb_z22, d81_inf_z22, d09_inf_sb_z22, d81_inf_sb_z23, d09_inf_sb_z22, d81_inf_sb_z22, d09_inf_sb_z22, d09_inf_z23, d09_inf_sb_z22, d09_inf_z22, d09_inf_sb_z22, d09_inf_sb_z23, d09_inf_sb_z22, d81_inf_z23, d02_inf_z23, d81_inf_z22, d02_inf_z23, d81_inf_sb_z23, d02_inf_z23, d81_inf_sb_z22, d02_inf_z23, d09_inf_z23, d02_inf_z23, d09_inf_z22, d02_inf_z23, d09_inf_sb_z23, d02_inf_z23, d09_inf_sb_z22, d02_inf_z23, d81_inf_z23, d02_inf_z22, d81_inf_z22, d02_inf_z22, d81_inf_sb_z23, d02_inf_z22, d81_inf_sb_z22, d02_inf_z22, d09_inf_z23, d02_inf_z22, d09_inf_z22, d02_inf_z22, d09_inf_sb_z23, d02_inf_z22, d09_inf_sb_z22, d02_inf_z22, d02_inf_z23, d02_inf_z22, d81_inf_z23, d02_inf_sb_z23, d81_inf_z22, d02_inf_sb_z23, d81_inf_sb_z23, d02_inf_sb_z23, d81_inf_sb_z22, d02_inf_sb_z23, d09_inf_z23, d02_inf_sb_z23, d09_inf_z22, d02_inf_sb_z23, d09_inf_sb_z23, d02_inf_sb_z23, d09_inf_sb_z22, d02_inf_sb_z23, d02_inf_z23, d02_inf_sb_z23, d02_inf_z22, d02_inf_sb_z23, d81_inf_z23, d02_inf_sb_z22, d81_inf_z22, d02_inf_sb_z22, d81_inf_sb_z23, d02_inf_sb_z22, d81_inf_sb_z22, d02_inf_sb_z22, d09_inf_z23, d02_inf_sb_z22, d09_inf_z22, d02_inf_sb_z22, d09_inf_sb_z23, d02_inf_sb_z22, d09_inf_sb_z22, d02_inf_sb_z22, d02_inf_z23, d02_inf_sb_z22, d02_inf_z22, d02_inf_sb_z22, d02_inf_sb_z23, d02_inf_sb_z22, d81_inf_z23, d01_inf_z23, d81_inf_z22, d01_inf_z23, d81_inf_sb_z23, d01_inf_z23, d81_inf_sb_z22, d01_inf_z23, d09_inf_z23, d01_inf_z23, d09_inf_z22, d01_inf_z23, d09_inf_sb_z23, d01_inf_z23, d09_inf_sb_z22, d01_inf_z23, d02_inf_z23, d01_inf_z23, d02_inf_z22, d01_inf_z23, d02_inf_sb_z23, d01_inf_z23, d02_inf_sb_z22, d01_inf_z23, d81_inf_z23, d01_inf_z22, d81_inf_z22, d01_inf_z22, d81_inf_sb_z23, d01_inf_z22, d81_inf_sb_z22, d01_inf_z22, d09_inf_z23, d01_inf_z22, d09_inf_z22, d01_inf_z22, d09_inf_sb_z23, d01_inf_z22, d09_inf_sb_z22, d01_inf_z22, d02_inf_z23, d01_inf_z22, d02_inf_z22, d01_inf_z22, d02_inf_sb_z23, d01_inf_z22, d02_inf_sb_z22, d01_inf_z22, d01_inf_z23, d01_inf_z22, d81_inf_z23, d01_inf_sb_z23, d81_inf_z22, d01_inf_sb_z23, d81_inf_sb_z23, d01_inf_sb_z23, d81_inf_sb_z22, d01_inf_sb_z23, d09_inf_z23, d01_inf_sb_z23, d09_inf_z22, d01_inf_sb_z23, d09_inf_sb_z23, d01_inf_sb_z23, d09_inf_sb_z22, d01_inf_sb_z23, d02_inf_z23, d01_inf_sb_z23, d02_inf_z22, d01_inf_sb_z23, d02_inf_sb_z23, d01_inf_sb_z23, d02_inf_sb_z22, d01_inf_sb_z23, d01_inf_z23, d01_inf_sb_z23, d01_inf_z22, d01_inf_sb_z23, d81_inf_z23, d01_inf_sb_z22, d81_inf_z22, d01_inf_sb_z22, d81_inf_sb_z23, d01_inf_sb_z22, d81_inf_sb_z22, d01_inf_sb_z22, d09_inf_z23, d01_inf_sb_z22, d09_inf_z22, d01_inf_sb_z22, d09_inf_sb_z23, d01_inf_sb_z22, d09_inf_sb_z22, d01_inf_sb_z22, d02_inf_z23, d01_inf_sb_z22, d02_inf_z22, d01_inf_sb_z22, d02_inf_sb_z23, d01_inf_sb_z22, d02_inf_sb_z22, d01_inf_sb_z22, d01_inf_z23, d01_inf_sb_z22, d01_inf_z22, d01_inf_sb_z22, d01_inf_sb_z23, d01_inf_sb_z22

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

## extra_contrasts = extra)
hs_nopower_sva_sig <- extract_significant_genes(
  hs_nopower_sva_table,
  excel = glue("analyses/macrophage_de/sig_tables/hs_nopower_sva_sig-v{ver}.xlsx"))

## Error: object 'hs_nopower_sva_table' not found

d01d02_zymo_sva_comp <- merge(hs_nopower_sva_table[["data"]][["d01_zymo"]],
                              hs_nopower_sva_table[["data"]][["d02_zymo"]],
                              by = "row.names")

## Error: object 'hs_nopower_sva_table' not found

d0102_zymo_sva_plot <- plot_linear_scatter(d01d02_zymo_sva_comp[, c("deseq_logfc.x", "deseq_logfc.y")])

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'as.data.frame': object 'd01d02_zymo_sva_comp' not found

d0102_zymo_sva_plot[["scatter"]]

## Error: object 'd0102_zymo_sva_plot' not found

d0102_zymo_sva_plot[["correlation"]]

## Error: object 'd0102_zymo_sva_plot' not found

d0102_zymo_sva_plot[["lm_rsq"]]

## Error: object 'd0102_zymo_sva_plot' not found

d09d81_zymo_sva_comp <- merge(hs_nopower_sva_table[["data"]][["d09_zymo"]],
                              hs_nopower_sva_table[["data"]][["d81_zymo"]],
                              by = "row.names")

## Error: object 'hs_nopower_sva_table' not found

d0981_zymo_sva_plot <- plot_linear_scatter(d09d81_zymo_sva_comp[, c("deseq_logfc.x", "deseq_logfc.y")])

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'as.data.frame': object 'd09d81_zymo_sva_comp' not found

d0981_zymo_sva_plot[["scatter"]]

## Error: object 'd0981_zymo_sva_plot' not found

d0981_zymo_sva_plot[["correlation"]]

## Error: object 'd0981_zymo_sva_plot' not found

d0981_zymo_sva_plot[["lm_rsq"]]

## Error: object 'd0981_zymo_sva_plot' not found

d01d81_zymo_sva_comp <- merge(hs_nopower_sva_table[["data"]][["d01_zymo"]],
                              hs_nopower_sva_table[["data"]][["d81_zymo"]],
                              by = "row.names")

## Error: object 'hs_nopower_sva_table' not found

d0181_zymo_sva_plot <- plot_linear_scatter(d01d81_zymo_sva_comp[, c("deseq_logfc.x", "deseq_logfc.y")])

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'as.data.frame': object 'd01d81_zymo_sva_comp' not found

d0181_zymo_sva_plot[["scatter"]]

## Error: object 'd0181_zymo_sva_plot' not found

d0181_zymo_sva_plot[["correlation"]]

## Error: object 'd0181_zymo_sva_plot' not found

d0181_zymo_sva_plot[["lm_rsq"]]

## Error: object 'd0181_zymo_sva_plot' not found

upset_plots_sva <- upsetr_sig(hs_nopower_sva_sig, both = TRUE,
                              contrasts = c("d01_zymo", "d02_zymo", "d09_zymo", "d81_zymo"))

## Error: object 'hs_nopower_sva_sig' not found

upset_plots_sva[["up"]]

## Error: object 'upset_plots_sva' not found

upset_plots_sva[["down"]]

## Error: object 'upset_plots_sva' not found

upset_plots_sva[["both"]]

## Error: object 'upset_plots_sva' not found

## The 7th element in the both groups list is the set shared among all donors.
## I don't feel like writing out x:y:z:a
groups <- upset_plots_sva[["both_groups"]]

## Error: object 'upset_plots_sva' not found

shared_genes <- attr(groups, "elements")[groups[[7]]] %>%
  gsub(pattern = "^gene:", replacement = "")

## Error: object 'groups' not found

shared_gp <- simple_gprofiler(shared_genes)

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

shared_gp[["pvalue_plots"]][["MF"]]

## Error: object 'shared_gp' not found

shared_gp[["pvalue_plots"]][["BP"]]

## Error: object 'shared_gp' not found

shared_gp[["pvalue_plots"]][["REAC"]]

## Error: object 'shared_gp' not found

shared_gp[["pvalue_plots"]][["WP"]]

## Error: object 'shared_gp' not found

10 Donor comparison

Now compare the donors to each other directly.

hs_donors <- set_conditions(hs_macr, fact = "donor")

## The numbers of samples by condition are:

## 
## d01 d02 d09 d81 
##  13  14  13  14

donor_de <- all_pairwise(hs_donors, model_svs = "svaseq",
                         model_fstring = "~ 0 + condition", filter = TRUE)

## d01 d02 d09 d81 
##  13  14  13  14

## Removing 9725 low-count genes (11756 remaining).

## Basic step 0/3: Normalizing data.

## Basic step 0/3: Converting data.

## I think this is failing? SummarizedExperiment

## Basic step 0/3: Transforming data.

## Setting 40036 entries to zero.

## This received a matrix of SVs.

## converting counts to integer mode

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## Warning in createContrastL(objFlt$formula, objFlt$data, L): Contrasts with only
## a single non-zero term are already evaluated by default.

## conditions
## d01 d02 d09 d81 
##  13  14  13  14

## conditions
## d01 d02 d09 d81 
##  13  14  13  14

## conditions
## d01 d02 d09 d81 
##  13  14  13  14

donor_de

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

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

## The primary analysis performed 6 comparisons.

donor_table <- combine_de_tables(
  donor_de,
  excel = glue("analyses/macrophage_de/de_tables/donor_tables-v{ver}.xlsx"))

## Looking for subscript invalid names, start of extract_keepers.

## Error in combine_extracted_plots(entry_name, combined, wanted_denominator, : object 'color_low' not found

donor_table

## Error: object 'donor_table' not found

donor_sig <- extract_significant_genes(
  donor_table,
  excel = glue("analyses/macrophage_de/sig_tables/donor_sig-v{ver}.xlsx"))

## Error: object 'donor_table' not found

donor_sig

## Error: object 'donor_sig' not found

11 Primary query contrasts

The final contrast in this list is interesting because it depends on the extra contrasts applied to the all_pairwise() above. In my way of thinking, the primary comparisons to consider are either cross-drug or cross-strain, but not both. However I think in at least a few instances Olga is interested in strain+drug / uninfected+nodrug.

11.1 Write contrast results

Now let us write out the xlsx file containing the above contrasts. The file with the suffix _table-version will therefore contain all genes and the file with the suffix _sig-version will contain only those deemed significant via our default criteria of DESeq2 |logFC| >= 1.0 and adjusted p-value <= 0.05.

11.2 Over representation searches

I decided to make one initially small, but I think quickly big change to the organization of this document: I am moving the GSEA searches up to immediately after the DE. I will then move the plots of the gprofiler results to immediately after the various volcano plots so that it is easier to interpret them.

I am reasonably certain this is the place to check that z23no drug / uninfected has the expected set of genes and that there is or is not a reactome result.

Reproducibility note: Given that this is entirely dependent on an online service, I must assume that the results will change over time; in addition their web servers undergo maintenance regularly, which may result in systematic failure of these analyses. I like gProfiler quite a lot for this type of stuff, but this is an important caveat.

Conversely, the clusterProfiler results later depend on a consistent orgdb annotation set (or reactome or whatever); those versions are fixed by the container installation.

all_gp <- all_gprofiler(hs_macr_sig, enrich_id_column = "hgnc_symbol")

## Error: object 'hs_macr_sig' not found

for (g in seq_len(length(all_gp))) {
  name <- names(all_gp)[g]
  datum <- all_gp[[name]]
  filename <- glue("analyses/macrophage_de/gprofiler/{name}_gprofiler-v{ver}.xlsx")
  written <- sm(write_gprofiler_data(datum, excel = filename))
}

## Error: object 'all_gp' not found

lesssig_all_gp <- all_gprofiler(hs_macr_lesssig, enrich_id_column = "hgnc_symbol")

## Error: object 'hs_macr_lesssig' not found

for (g in seq_len(length(lesssig_all_gp))) {
  name <- names(lesssig_all_gp)[g]
  datum <- lesssig_all_gp[[name]]
  filename <- glue("analyses/macrophage_de/gprofiler/{name}_gprofiler_lesssig-v{ver}.xlsx")
  written <- sm(write_gprofiler_data(datum, excel = filename))
}

## Error: object 'lesssig_all_gp' not found

11.3 Explicit GSEA search vis clusterProfiler

all_cp <- all_cprofiler(hs_macr_sig, hs_macr_table)

## Error: object 'hs_macr_sig' not found

11.4 Specific desires in Reactome results

In previous analyses (I think by Dr. Colmenares), a specific Tryptophan biosynthesis pathway was observed. Partciularly in the 2.3/uninfected comparison. I think my gprofiler analysis is too stringent and therefore not observing this. Olga asked if I could look at that and see if there are trivial settings I can change to highlight this pathway. The two most likely things I can change are the stringencies of the DE analysis and/or gProfiler.

test_z23_uninf_up <- hs_macr_sig[["deseq"]][["ups"]][["z23nosb_vs_uninf"]]

## Error: object 'hs_macr_sig' not found

nrow(test_z23_uninf_up)

## Error: object 'test_z23_uninf_up' not found

test_z23_uninf_down <- hs_macr_sig[["deseq"]][["downs"]][["z23nosb_vs_uninf"]]

## Error: object 'hs_macr_sig' not found

nrow(test_z23_uninf_down)

## Error: object 'test_z23_uninf_down' not found

test_gp_up <- simple_gprofiler(test_z23_uninf_up, enrich_id_column = "hgnc_symbol",
                               threshold = 1.0)

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

test_gp_up

## Error: object 'test_gp_up' not found

written_up <- write_gprofiler_data(test_gp_up, excel = "excel/z23_uninf_gp_up_all.xlsx")

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'as.data.frame': object 'test_gp_up' not found

test_gp_down <- simple_gprofiler(test_z23_uninf_down, enrich_id_column = "hgnc_symbol",
                                 threshold = 1.0)

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

test_gp_down

## Error: object 'test_gp_down' not found

written_down <- write_gprofiler_data(test_gp_down, excel = "excel/z23_uninf_gp_down_all.xlsx")

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'as.data.frame': object 'test_gp_down' not found

11.5 Plot contrasts of interest

One suggestion I received recently was to set the axes for these volcano plots to be static rather than let ggplot choose its own. I am assuming this is only relevant for pairs of contrasts, but that might not be true.

11.6 Individual zymodemes vs. uninfected

The following blocks will be a lot of repetition. In each case I am yanking out the volcano plot for a specific contrast and showing the original followed by a version with different colors/labelling.

11.6.1 Infected with z2.3 no Antimonial vs. Uninfected

plot_colors <- get_se_colors(hs_macr_table[["input"]][["input"]])

## Error in get_se_colors(hs_macr_table[["input"]][["input"]]): could not find function "get_se_colors"

## The original plot from my xlsx file
hs_macr_table[["plots"]][["z23nosb_vs_uninf"]][["deseq_vol_plots"]]

## Error: object 'hs_macr_table' not found

z23nosb_vs_uninf_volcano <- plot_volcano_condition_de(
  input = hs_macr_table[["data"]][["z23nosb_vs_uninf"]],
  fc_col = "deseq_logfc", p_col = "deseq_adjp",
  label = 10, label_column = "hgnc_symbol",
  color_low = plot_colors[["uninfnone"]], color_high = plot_colors[["infz23"]])

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

labeled <- z23nosb_vs_uninf_volcano[["plot"]] +
  scale_x_continuous(limits = c(-6, 21), breaks = c(-6, -4, -2, 0, 2, 4, 6, 8, 10, 20)) +
  ggbreak::scale_x_break(c(10, 19), scales = 0.2, space = 0.02)

## Error: object 'z23nosb_vs_uninf_volcano' not found

pp(file = "figures/fig2a_labeled_with_break.svg")
labeled

## Error: object 'labeled' not found

dev.off()

## png 
##   2

labeled

## Error: object 'labeled' not found

plotly::ggplotly(z23nosb_vs_uninf_volcano[["plot"]])

## Error: object 'z23nosb_vs_uninf_volcano' not found

The following provides some of the over-representation plots from gProfiler2.

all_gp[["z23nosb_vs_uninf_up"]][["pvalue_plots"]][["REAC"]]

## Error: object 'all_gp' not found

## Reactome, zymodeme2.3 without drug vs. uninfected without drug, up.
all_gp[["z23nosb_vs_uninf_up"]][["pvalue_plots"]][["KEGG"]]

## Error: object 'all_gp' not found

## KEGG, zymodeme2.3 without drug vs. uninfected without drug, up.
##all_gp[["z23nosb_vs_uninf_up"]][["pvalue_plots"]][["MF"]]
## MF, zymodeme2.3 without drug vs. uninfected without drug, up.
all_gp[["z23nosb_vs_uninf_up"]][["pvalue_plots"]][["TF"]]

## Error: object 'all_gp' not found

## TF, zymodeme2.3 without drug vs. uninfected without drug, up.
all_gp[["z23nosb_vs_uninf_up"]][["pvalue_plots"]][["WP"]]

## Error: object 'all_gp' not found

## WikiPathways, zymodeme2.3 without drug vs. uninfected without drug, up.
all_gp[["z23nosb_vs_uninf_up"]][["interactive_plots"]][["WP"]]

## Error: object 'all_gp' not found

message("Olga received a query about the following result, I think it is null.")

## Olga received a query about the following result, I think it is null.

all_gp[["z23nosb_vs_uninf_down"]][["pvalue_plots"]][["REAC"]]

## Error: object 'all_gp' not found

message("Is the previous plot null?")

## Is the previous plot null?

## Reactome, zymodeme2.3 without drug vs. uninfected without drug, down.
all_gp[["z23nosb_vs_uninf_down"]][["pvalue_plots"]][["MF"]]

## Error: object 'all_gp' not found

## MF, zymodeme2.3 without drug vs. uninfected without drug, down.
all_gp[["z23nosb_vs_uninf_down"]][["pvalue_plots"]][["TF"]]

## Error: object 'all_gp' not found

## TF, zymodeme2.3 without drug vs. uninfected without drug, down.

We have some other categorical enrichment plots available via enrichplot, let us try a few out for contrasts of interest and see if any of them prove helpful.

First, as a reminder, here are the contrasts which are available to examine, in each case there is an _up and _down enrichment object in the data. Thus in the following list I am going to arbitrarily print out some invocations which extract putatively interesting bits of data.

• z23nosb_vs_uninf: all_gp[[“z23nosb_vs_uninf_up”]][[“BP_enrich”]]
• z22nosb_vs_uninf.
• z23nosb_vs_z22nosb.
• z23sb_vs_z22sb.
• z23sb_vs_z23nosb.
• z22sb_vs_z22nosb.
• z23sb_vs_sb.
• z22sb_vs_sb.
• z23sb_vs_uninf.
• z22sb_vs_uninf.
• sb_vs_uninf.
• extra_z2322.
• extra_drugnodrug.

z23nosb_uninf_up_go <- all_gp[["z23nosb_vs_uninf_up"]][["BP_enrich"]]

## Error: object 'all_gp' not found

z23nosb_uninf_up_go_pair <- pairwise_termsim(z23nosb_uninf_up_go)

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

dotplot(z23nosb_uninf_up_go)

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

emapplot(z23nosb_uninf_up_go_pair)

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

##ssplot(z23nosb_uninf_up_go_pair)
treeplot(z23nosb_uninf_up_go_pair)

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

upsetplot(z23nosb_uninf_up_go)

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

cnetplot(z23nosb_uninf_up_go)

## Error: object 'z23nosb_uninf_up_go' not found

11.6.2 Repeat, but using a less strict set of ‘significant genes’

I am not entirely certain if the Reactome results Olga showed me included both up and down genes? I am going to assume for the moment that it was just up/down, but if that proves intractable I will go back to the manuscript and read more carefully (e.g. I just remembered where the picture came from!)

11.6.2.1 Add a little topgo

In the process of exploring the various parameters used with gProfiler2, I found myself thinking that it would be nice to have some topgo results to compare against. The following block is the result of that thought.

test_genes_up <- hs_macr_lesssig[["deseq"]][["ups"]][["z23nosb_vs_uninf"]]

## Error: object 'hs_macr_lesssig' not found

test_query_up <- simple_gprofiler(test_genes_up, threshold = 0.1)

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

test_query_up[["pvalue_plots"]][["REAC"]]

## Error: object 'test_query_up' not found

pdf(file = "images/test_query_biological_process_z23_vs_uninf_up.pdf", height = 12, width = 9)
test_query_up[["pvalue_plots"]][["BP"]]

## Error: object 'test_query_up' not found

dev.off()

## png 
##   2

enrichplot::dotplot(test_query_up[["BP_enrich"]])

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

test_genes_down <- hs_macr_lesssig[["deseq"]][["downs"]][["z23nosb_vs_uninf"]]

## Error: object 'hs_macr_lesssig' not found

test_query_down <- simple_gprofiler(test_genes_down)

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

test_query_down[["pvalue_plots"]][["REAC"]]

## Error: object 'test_query_down' not found

## I keep getting all sorts of annoying biomart errors.
hs_go <- try(load_biomart_go(archive = FALSE, overwrite = TRUE))

## Using mart: ENSEMBL_MART_ENSEMBL from host: useast.ensembl.org.

## Successfully connected to the hsapiens_gene_ensembl database.

## Error in httr2::req_perform(req) : Failed to perform HTTP request.
## Caused by error in `curl::curl_fetch_memory()` at httr2/R/req-perform.R:201:5:
## ! Timeout was reached [useast.ensembl.org]:
## Operation timed out after 300000 milliseconds with 10690645 bytes received

## Unable to download annotation data.

if ("try-error" %in% class(hs_go)) {
  hs_go <- load_biomart_go(archive = TRUE, month = "04", year = "2020", overwrite = TRUE)
}
test_topgo_up <- simple_topgo(test_genes_up, go_db = hs_go[["go"]], parallel = FALSE)

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

written_topgo <- write_topgo_data(
  test_topgo_up,
  excel = glue("analyses/macrophage_de/ontology_topgo/topgo_z23_uninf_less_strict.xlsx"))

## Error: object 'test_topgo_up' not found

11.6.3 Infected with z2.2 no Antimonial vs. Uninfected

Here is where things will get most repetitive. In each instance I am creating a couple of volcano plots followed by printing some of the gProfiler2 results (when I get the itch).

The following should be a slightly improved version of our extant figure 2B.

## The original plot
hs_macr_table[["plots"]][["z22nosb_vs_uninf"]][["deseq_vol_plots"]]

## Error: object 'hs_macr_table' not found

z22nosb_vs_uninf_volcano <- plot_volcano_condition_de(
  hs_macr_table[["data"]][["z22nosb_vs_uninf"]], "z22nosb_vs_uninf",
  fc_col = "deseq_logfc", p_col = "deseq_adjp",
  label = 10, label_column = "hgnc_symbol",
  color_low = plot_colors[["uninfnone"]], color_high = plot_colors[["infz22"]])

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

labeled <- z22nosb_vs_uninf_volcano[["plot"]] +
  scale_x_continuous(limits = c(-2, 21), breaks = c(-2, 0, 2, 4, 6, 8, 10, 21, 22)) +
  ggbreak::scale_x_break(c(11, 20), scales = 0.2, space = 0.02)

## Error: object 'z22nosb_vs_uninf_volcano' not found

pp(file = "figures/fig2b_labeled_with_break.svg")
labeled

## Error: object 'labeled' not found

dev.off()

## png 
##   2

labeled

## Error: object 'labeled' not found

plotly::ggplotly(z22nosb_vs_uninf_volcano[["plot"]])

## Error: object 'z22nosb_vs_uninf_volcano' not found

Add some pvalue barplots from gProfiler for this contrast.

all_gp[["z22nosb_vs_uninf_up"]][["pvalue_plots"]][["REAC"]]

## Error: object 'all_gp' not found

## Reactome, zymodeme2.2 without drug vs. uninfected without drug, up.
all_gp[["z22nosb_vs_uninf_up"]][["pvalue_plots"]][["MF"]]

## Error: object 'all_gp' not found

## MF, zymodeme2.2 without drug vs. uninfected without drug, up.
all_gp[["z22nosb_vs_uninf_up"]][["pvalue_plots"]][["TF"]]

## Error: object 'all_gp' not found

## TF, zymodeme2.2 without drug vs. uninfected without drug, up.
all_gp[["z22nosb_vs_uninf_up"]][["pvalue_plots"]][["WP"]]

## Error: object 'all_gp' not found

## WikiPathways, zymodeme2.2 without drug vs. uninfected without drug, up.

all_gp[["z22nosb_vs_uninf_down"]][["pvalue_plots"]][["REAC"]]

## Error: object 'all_gp' not found

## Reactome, zymodeme2.2 without drug vs. uninfected without drug, down.
all_gp[["z22nosb_vs_uninf_down"]][["pvalue_plots"]][["MF"]]

## Error: object 'all_gp' not found

## MF, zymodeme2.2 without drug vs. uninfected without drug, down.
all_gp[["z22nosb_vs_uninf_down"]][["pvalue_plots"]][["TF"]]

## Error: object 'all_gp' not found

## TF, zymodeme2.3 without drug vs. uninfected without drug, down.

11.6.4 Infected with z2.3 treated vs. Uninfected treated

I do not think this plot is used at this time.

## The original plot
hs_macr_table[["plots"]][["z23sb_vs_sb"]][["deseq_vol_plots"]]

## Error: object 'hs_macr_table' not found

z23sb_vs_uninfsb_volcano <- plot_volcano_condition_de(
  hs_macr_table[["data"]][["z23sb_vs_sb"]], "z23sb_vs_sb",
  fc_col = "deseq_logfc", p_col = "deseq_adjp",
  label = 10, label_column = "hgnc_symbol",
  color_low = plot_colors[["infsbz23"]], color_high = plot_colors[["uninfsbnone"]])

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

z23sb_vs_uninfsb_volcano[["plot"]]

## Error: object 'z23sb_vs_uninfsb_volcano' not found

plotly::ggplotly(z23sb_vs_uninfsb_volcano[["plot"]])

## Error: object 'z23sb_vs_uninfsb_volcano' not found

11.6.5 Infected with z2.3 untreated vs. z2.2 untreated

This is figure 2C at this time.

## The original plot
hs_macr_table[["plots"]][["z23nosb_vs_z22nosb"]][["deseq_vol_plots"]]

## Error: object 'hs_macr_table' not found

z23nosb_vs_z22nosb_volcano <- plot_volcano_condition_de(
  hs_macr_table[["data"]][["z23nosb_vs_z22nosb"]], "z23nosb_vs_z22nosb",
  fc_col = "deseq_logfc", p_col = "deseq_adjp",
  label = 10, label_column = "hgnc_symbol",
  color_low = plot_colors[["infz23"]], color_high = plot_colors[["infz22"]])

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

labeled <- z23nosb_vs_z22nosb_volcano[["plot"]] +
  scale_x_continuous(breaks = c(-10, -8, -6, -4, -2, 0, 2, 4, 6))

## Error: object 'z23nosb_vs_z22nosb_volcano' not found

pp(file = "figures/fig2c_labeled.svg")
labeled

## Error: object 'labeled' not found

dev.off()

## png 
##   2

labeled

## Error: object 'labeled' not found

11.6.6 Infected with z2.3 treated vs. z2.2 treated

This is currently figure 3C.

FIXME: The axis label isn’t quite right for the ggbreak.

## The original plot
hs_macr_table[["plots"]][["z23sb_vs_z22sb"]][["deseq_vol_plots"]]

## Error: object 'hs_macr_table' not found

z23sb_vs_z22sb_volcano <- plot_volcano_condition_de(
  hs_macr_table[["data"]][["z23sb_vs_z22sb"]], "z23sb_vs_z22sb",
  fc_col = "deseq_logfc", p_col = "deseq_adjp",
  label = 10, label_column = "hgnc_symbol",
  color_high = plot_colors[["infsbz23"]], color_low = plot_colors[["infsbz22"]])

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

labeled <- z23sb_vs_z22sb_volcano[["plot"]] +
  scale_x_continuous(breaks = c(-23, -6, -4, -2, 0, 2, 4, 6)) +
  ggbreak::scale_x_break(c(-5, -22.5), scales = 10, space = 0.02)

## Error: object 'z23sb_vs_z22sb_volcano' not found

pp(file = "figures/fig3c_labeled_breaks.svg")
labeled

## Error: object 'labeled' not found

dev.off()

## png 
##   2

labeled

## Error: object 'labeled' not found

11.6.7 Infected with z2.3 SB treated vs. z2.3 untreated

I think this is currently figure 3A.

FIXME: The axis label for the ggbreak isn’t quite right.

## The original plot
hs_macr_table[["plots"]][["z23sb_vs_z23nosb"]][["deseq_vol_plots"]]

## Error: object 'hs_macr_table' not found

z23sb_vs_z23nosb_volcano <- plot_volcano_condition_de(
  hs_macr_table[["data"]][["z23sb_vs_z23nosb"]], "z23sb_vs_z23nosb",
  fc_col = "deseq_logfc", p_col = "deseq_adjp",
  label = 10, label_column = "hgnc_symbol",
  color_high = plot_colors[["infsbz23"]], color_low = plot_colors[["infz23"]])

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

labeled <- z23sb_vs_z23nosb_volcano[["plot"]] +
  scale_x_continuous(limits = c(-19, 6),
                     breaks = c(-20, -18, -16, -14, -12, -10, -6, -4, -2, 0, 2, 4, 6)) +
  ggbreak::scale_x_break(c(-17, -8), scales = 17, space = 0.02)

## Error: object 'z23sb_vs_z23nosb_volcano' not found

pp(file = "figures/fig3a_labeled_with_break.svg")
labeled

## Error: object 'labeled' not found

dev.off()

## png 
##   2

labeled

## Error: object 'labeled' not found

11.6.8 Infected with z2.3 SB treated vs. z2.3 untreated

## The original plot
hs_macr_table[["plots"]][["z22sb_vs_z22nosb"]][["deseq_vol_plots"]]

## Error: object 'hs_macr_table' not found

z22sb_vs_z22nosb_volcano <- plot_volcano_condition_de(
  hs_macr_table[["data"]][["z22sb_vs_z22nosb"]], "z22sb_vs_z22nosb",
  fc_col = "deseq_logfc", p_col = "deseq_adjp",
  label = 10, label_column = "hgnc_symbol",
  color_high = plot_colors[["infsbz22"]], color_low = plot_colors[["infz22"]])

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

labeled <- z22sb_vs_z22nosb_volcano[["plot"]] +
  scale_x_continuous(breaks = c(-6, -4, -2, 0, 2, 4, 6))

## Error: object 'z22sb_vs_z22nosb_volcano' not found

pp(file = "figures/fig3b_labeled.svg")
labeled

## Error: object 'labeled' not found

dev.off()

## png 
##   2

labeled

## Error: object 'labeled' not found

11.6.9 Infected with z2.3 SB treated vs. uninfected treated

x_limits <- c(-6, 6)
## The original plot
hs_macr_table[["plots"]][["z23sb_vs_sb"]][["deseq_vol_plots"]]

## Error: object 'hs_macr_table' not found

z23sb_vs_sb_volcano <- plot_volcano_condition_de(
  hs_macr_table[["data"]][["z23sb_vs_sb"]], "z23sb_vs_sb",
  fc_col = "deseq_logfc", p_col = "deseq_adjp",
  label = 10, label_column = "hgnc_symbol", invert = TRUE,
  color_low = plot_colors[["infsbz23"]], color_high = plot_colors[["uninfsbnone"]])

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

z23sb_vs_sb_volcano[["plot"]]

## Error: object 'z23sb_vs_sb_volcano' not found

11.6.10 Infected with z2.2 SB treated vs. uninfected treated

## The original plot
hs_macr_table[["plots"]][["z22sb_vs_sb"]][["deseq_vol_plots"]]

## Error: object 'hs_macr_table' not found

z22sb_vs_sb_volcano <- plot_volcano_condition_de(
  hs_macr_table[["data"]][["z22sb_vs_sb"]], "z22sb_vs_sb",
  fc_col = "deseq_logfc", p_col = "deseq_adjp",
  label = 10, label_column = "hgnc_symbol", invert = TRUE,
  color_low = plot_colors[["infsbz22"]], color_high = plot_colors[["uninfsbnone"]])

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

z22sb_vs_sb_volcano[["plot"]]

## Error: object 'z22sb_vs_sb_volcano' not found

11.6.11 Uninfected+SbV vs. Uninfected-SbV

This is currently figure 3D.

FIXME: This needs the BOLA2B ggbreak.

## The original plot
hs_macr_table[["plots"]][["sb_vs_uninf"]][["deseq_vol_plots"]]

## Error: object 'hs_macr_table' not found

sb_vs_uninf_volcano <- plot_volcano_condition_de(
  hs_macr_table[["data"]][["sb_vs_uninf"]], "sb_vs_uninf",
  fc_col = "deseq_logfc", p_col = "deseq_adjp",
  label = 10, label_column = "hgnc_symbol",
  color_high = plot_colors[["uninfsbnone"]], color_low = plot_colors[["uninfnone"]])

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

labeled <- sb_vs_uninf_volcano[["plot"]] +
  scale_x_continuous(breaks = c(-23, -6, -4, -2, 0, 2, 4, 6)) +
  ggbreak::scale_x_break(c(-5, -22.5), scales = 10, space = 0.02)

## Error: object 'sb_vs_uninf_volcano' not found

pp(file = "figures/fig3d_labeled_breaks.svg")
labeled

## Error: object 'labeled' not found

dev.off()

## png 
##   2

labeled

## Error: object 'labeled' not found

11.7 Double-check that gene counts match my perceptions

Check that my perception of the number of significant up/down genes matches what the table/venn says. In the following block I am performing some venn/upset analyses to see if the numbers of genes match what we have in the current version of the manuscript (plus or minus a gene) and thus if my interpretation of the figure/legend text matches what I think it means.

shared <- Vennerable::Venn(list(
  "drug" = rownames(hs_macr_sig[["deseq"]][["ups"]][["z23sb_vs_uninf"]]),
  "nodrug" = rownames(hs_macr_sig[["deseq"]][["ups"]][["z23nosb_vs_uninf"]])))

## Error: object 'hs_macr_sig' not found

pp(file = "images/z23_vs_uninf_venn_up.png")
Vennerable::plot(shared)

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

dev.off()

## png 
##   2

Vennerable::plot(shared)

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

## I see 910 z23sb/uninf and 670 no z23nosb/uninf genes in the venn diagram.
length(shared@IntersectionSets[["10"]]) + length(shared@IntersectionSets[["11"]])

## Error: object 'shared' not found

dim(hs_macr_sig[["deseq"]][["ups"]][["z23sb_vs_uninf"]])

## Error: object 'hs_macr_sig' not found

shared <- Vennerable::Venn(list(
  "drug" = rownames(hs_macr_sig[["deseq"]][["ups"]][["z22sb_vs_uninf"]]),
  "nodrug" = rownames(hs_macr_sig[["deseq"]][["ups"]][["z22nosb_vs_uninf"]])))

## Error: object 'hs_macr_sig' not found

pp(file = "images/z22_vs_uninf_venn_up.png")
Vennerable::plot(shared)

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

dev.off()

## png 
##   2

Vennerable::plot(shared)

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

length(shared@IntersectionSets[["10"]]) + length(shared@IntersectionSets[["11"]])

## Error: object 'shared' not found

dim(hs_macr_sig[["deseq"]][["ups"]][["z22sb_vs_uninf"]])

## Error: object 'hs_macr_sig' not found

Note to self: There is an error in my volcano plot code which takes effect when the numerator and denominator of the all_pairwise contrasts are different than those in combine_de_tables. It is putting the ups/downs on the correct sides of the plot, but calling the down genes ‘up’ and vice-versa. The reason for this is that I did a check for this happening, but used the wrong argument to handle it.

A likely bit of text for these volcano plots:

The set of genes differentially expressed between the zymodeme 2.3 and uninfected samples without druge treatment was quantified with DESeq2 and included surrogate estimates from SVA. Given the criteria of significance of a abs(logFC) >= 1.0 and false discovery rate adjusted p-value <= 0.05, 670 genes were observed as significantly increased between the infected and uninfected samples and 386 were observed as decreased. The most increased genes from the uninfected samples include some which are potentially indicative of a strong innate immune response and the inflammatory response.

In contrast, when the set of genes differentially expressed between the zymodeme 2.2 and uninfected samples was visualized, only 7 genes were observed as decreased and 435 increased. The inflammatory response was significantly less apparent in this set, but instead included genes related to transporter activity and oxidoreductases.

11.8 Direct zymodeme comparisons

An orthogonal comparison to that performed above is to directly compare the zymodeme 2.3 and 2.2 samples with and without antimonial treatment.

11.8.1 Z2.3 / z2.2 without drug

z23nosb_vs_z22nosb_volcano <- plot_volcano_de(
  table = hs_macr_table[["data"]][["z23nosb_vs_z22nosb"]],
  fc_col = "deseq_logfc", p_col = "deseq_adjp",
  shapes_by_state = FALSE, color_by = "fc",  label = 10, label_column = "hgnc_symbol")

## Error in plot_volcano_de(table = hs_macr_table[["data"]][["z23nosb_vs_z22nosb"]], : could not find function "plot_volcano_de"

plotly::ggplotly(z23nosb_vs_z22nosb_volcano[["plot"]])

## Error: object 'z23nosb_vs_z22nosb_volcano' not found

z23sb_vs_z22sb_volcano <- plot_volcano_de(
  table = hs_macr_table[["data"]][["z23sb_vs_z22sb"]],
  fc_col = "deseq_logfc", p_col = "deseq_adjp",
  shapes_by_state = FALSE, color_by = "fc",  label = 10, label_column = "hgnc_symbol")

## Error in plot_volcano_de(table = hs_macr_table[["data"]][["z23sb_vs_z22sb"]], : could not find function "plot_volcano_de"

plotly::ggplotly(z23sb_vs_z22sb_volcano[["plot"]])

## Error: object 'z23sb_vs_z22sb_volcano' not found

z23nosb_vs_z22nosb_volcano[["plot"]] +
  xlim(-10, 10) +
  ylim(0, 60)

## Error: object 'z23nosb_vs_z22nosb_volcano' not found

pp(file = "images/z23nosb_vs_z22nosb_reactome_up.svg",
   image = all_gp[["z23nosb_vs_z22nosb_up"]][["pvalue_plots"]][["REAC"]],
   height = 12, width = 9)

## Error: object 'all_gp' not found

all_gp[["z23nosb_vs_z22nosb_up"]][["pvalue_plots"]][["REAC"]]

## Error: object 'all_gp' not found

## Reactome, zymodeme2.3 without drug vs. uninfected without drug, up.
all_gp[["z23nosb_vs_z22nosb_up"]][["pvalue_plots"]][["KEGG"]]

## Error: object 'all_gp' not found

## KEGG, zymodeme2.3 without drug vs. uninfected without drug, up.
all_gp[["z23nosb_vs_z22nosb_up"]][["pvalue_plots"]][["MF"]]

## Error: object 'all_gp' not found

## MF, zymodeme2.3 without drug vs. uninfected without drug, up.
all_gp[["z23nosb_vs_z22nosb_up"]][["pvalue_plots"]][["TF"]]

## Error: object 'all_gp' not found

## TF, zymodeme2.3 without drug vs. uninfected without drug, up.
all_gp[["z23nosb_vs_z22nosb_up"]][["pvalue_plots"]][["WP"]]

## Error: object 'all_gp' not found

## WikiPathways, zymodeme2.3 without drug vs. uninfected without drug, up.
all_gp[["z23nosb_vs_z22nosb_up"]][["interactive_plots"]][["WP"]]

## Error: object 'all_gp' not found

pp(file = "images/z23nosb_vs_z22nosb_reactome_down.svg",
   image = all_gp[["z23nosb_vs_z22nosb_down"]][["pvalue_plots"]][["REAC"]],
   height = 12, width = 9)

## Error: object 'all_gp' not found

all_gp[["z23nosb_vs_z22nosb_down"]][["pvalue_plots"]][["REAC"]]

## Error: object 'all_gp' not found

## Reactome, zymodeme2.3 without drug vs. uninfected without drug, down.
all_gp[["z23nosb_vs_z22nosb_down"]][["pvalue_plots"]][["MF"]]

## Error: object 'all_gp' not found

## MF, zymodeme2.3 without drug vs. uninfected without drug, down.
all_gp[["z23nosb_vs_z22nosb_down"]][["pvalue_plots"]][["TF"]]

## Error: object 'all_gp' not found

## TF, zymodeme2.3 without drug vs. uninfected without drug, down.

11.8.2 z2.3 / z2.2 with drug

z23sb_vs_z22sb_volcano[["plot"]] +
  xlim(-10, 10) +
  ylim(0, 60)

## Error: object 'z23sb_vs_z22sb_volcano' not found

pp(
  file = "images/z23sb_vs_z22sb_reactome_up.png",
  image = all_gp[["z23sb_vs_z22sb_up"]][["pvalue_plots"]][["REAC"]],
  height = 12, width = 9)

## Error: object 'all_gp' not found

all_gp[["z23sb_vs_z22sb_up"]][["pvalue_plots"]][["REAC"]]

## Error: object 'all_gp' not found

## Reactome, zymodeme2.3 without drug vs. uninfected without drug, up.
all_gp[["z23sb_vs_z22sb_up"]][["pvalue_plots"]][["KEGG"]]

## Error: object 'all_gp' not found

## KEGG, zymodeme2.3 without drug vs. uninfected without drug, up.
all_gp[["z23sb_vs_z22sb_up"]][["pvalue_plots"]][["MF"]]

## Error: object 'all_gp' not found

## MF, zymodeme2.3 without drug vs. uninfected without drug, up.
all_gp[["z23sb_vs_z22sb_up"]][["pvalue_plots"]][["TF"]]

## Error: object 'all_gp' not found

## TF, zymodeme2.3 without drug vs. uninfected without drug, up.
all_gp[["z23sb_vs_z22sb_up"]][["pvalue_plots"]][["WP"]]

## Error: object 'all_gp' not found

## WikiPathways, zymodeme2.3 without drug vs. uninfected without drug, up.
all_gp[["z23sb_vs_z22sb_up"]][["interactive_plots"]][["WP"]]

## Error: object 'all_gp' not found

all_gp[["z23sb_vs_z22sb_down"]][["pvalue_plots"]][["REAC"]]

## Error: object 'all_gp' not found

## Reactome, zymodeme2.3 without drug vs. uninfected without drug, down.
all_gp[["z23sb_vs_z22sb_down"]][["pvalue_plots"]][["MF"]]

## Error: object 'all_gp' not found

## MF, zymodeme2.3 without drug vs. uninfected without drug, down.
all_gp[["z23sb_vs_z22sb_down"]][["pvalue_plots"]][["TF"]]

## Error: object 'all_gp' not found

## TF, zymodeme2.3 without drug vs. uninfected without drug, down.

11.8.3 Venn to see shared/unique genes

Once again I wish to pull out the significant genes and see how my numbers match against the text.

shared <- Vennerable::Venn(list(
  "drug" = rownames(hs_macr_sig[["deseq"]][["ups"]][["z23sb_vs_z22sb"]]),
  "nodrug" = rownames(hs_macr_sig[["deseq"]][["ups"]][["z23nosb_vs_z22nosb"]])))

## Error: object 'hs_macr_sig' not found

pp(file = "images/drug_nodrug_venn_up.png")
Vennerable::plot(shared)

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

dev.off()

## png 
##   2

Vennerable::plot(shared)

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

shared <- Vennerable::Venn(
  list("drug" = rownames(hs_macr_sig[["deseq"]][["downs"]][["z23sb_vs_z22sb"]]),
       "nodrug" = rownames(hs_macr_sig[["deseq"]][["downs"]][["z23nosb_vs_z22nosb"]])))

## Error: object 'hs_macr_sig' not found

pp(file = "images/drug_nodrug_venn_down.png")
Vennerable::plot(shared)

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

dev.off()

## png 
##   2

Vennerable::plot(shared)

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

A slightly different way of looking at the differences between the two zymodeme infections is to directly compare the infected samples with and without drug. Thus, when a volcano plot showing the comparison of the zymodeme 2.3 vs. 2.2 samples was plotted, 484 genes were observed as increased and 422 decreased; these groups include many of the same inflammatory (up) and membrane (down) genes.

Similar patterns were observed when the antimonial was included. Thus, when a Venn diagram of the two sets of increased genes was plotted, a significant number of the genes was observed as increased (313) and decreased (244) in both the untreated and antimonial treated samples.

11.9 Drug effects on each zymodeme infection

Another likely question is to directly compare the treated vs untreated samples for each zymodeme infection in order to visualize the effects of antimonial.

11.9.1 z2.3 with and without drug

z23sb_vs_z23nosb_volcano <- plot_volcano_de(
  table = hs_macr_table[["data"]][["z23sb_vs_z23nosb"]],
  fc_col = "deseq_logfc", p_col = "deseq_adjp",
  shapes_by_state = FALSE, color_by = "fc",  label = 10, label_column = "hgnc_symbol")

## Error in plot_volcano_de(table = hs_macr_table[["data"]][["z23sb_vs_z23nosb"]], : could not find function "plot_volcano_de"

plotly::ggplotly(z23sb_vs_z23nosb_volcano[["plot"]])

## Error: object 'z23sb_vs_z23nosb_volcano' not found

z22sb_vs_z22nosb_volcano <- plot_volcano_de(
  table = hs_macr_table[["data"]][["z22sb_vs_z22nosb"]],
  fc_col = "deseq_logfc", p_col = "deseq_adjp",
  shapes_by_state = FALSE, color_by = "fc",  label = 10, label_column = "hgnc_symbol")

## Error in plot_volcano_de(table = hs_macr_table[["data"]][["z22sb_vs_z22nosb"]], : could not find function "plot_volcano_de"

plotly::ggplotly(z22sb_vs_z22nosb_volcano[["plot"]])

## Error: object 'z22sb_vs_z22nosb_volcano' not found

z23sb_vs_z23nosb_volcano[["plot"]] +
  xlim(-8, 8) +
  ylim(0, 210)

## Error: object 'z23sb_vs_z23nosb_volcano' not found

pp(file = "images/z23sb_vs_z23nosb_reactome_up.png",
   image = all_gp[["z23sb_vs_z23nosb_up"]][["pvalue_plots"]][["REAC"]],
   height = 12, width = 9)

## Error: object 'all_gp' not found

all_gp[["z23sb_vs_z23nosb_up"]][["pvalue_plots"]][["REAC"]]

## Error: object 'all_gp' not found

## Reactome, zymodeme2.3 without drug vs. uninfected without drug, up.
all_gp[["z23sb_vs_z23nosb_up"]][["pvalue_plots"]][["KEGG"]]

## Error: object 'all_gp' not found

## KEGG, zymodeme2.3 without drug vs. uninfected without drug, up.
all_gp[["z23sb_vs_z23nosb_up"]][["pvalue_plots"]][["MF"]]

## Error: object 'all_gp' not found

## MF, zymodeme2.3 without drug vs. uninfected without drug, up.
all_gp[["z23sb_vs_z23nosb_up"]][["pvalue_plots"]][["TF"]]

## Error: object 'all_gp' not found

## TF, zymodeme2.3 without drug vs. uninfected without drug, up.
all_gp[["z23sb_vs_z23nosb_up"]][["pvalue_plots"]][["WP"]]

## Error: object 'all_gp' not found

## WikiPathways, zymodeme2.3 without drug vs. uninfected without drug, up.
all_gp[["z23sb_vs_z23nosb_up"]][["interactive_plots"]][["WP"]]

## Error: object 'all_gp' not found

all_gp[["z23sb_vs_z23nosb_down"]][["pvalue_plots"]][["REAC"]]

## Error: object 'all_gp' not found

## Reactome, zymodeme2.3 without drug vs. uninfected without drug, down.
all_gp[["z23sb_vs_z23nosb_down"]][["pvalue_plots"]][["MF"]]

## Error: object 'all_gp' not found

## MF, zymodeme2.3 without drug vs. uninfected without drug, down.
all_gp[["z23sb_vs_z23nosb_down"]][["pvalue_plots"]][["TF"]]

## Error: object 'all_gp' not found

## TF, zymodeme2.3 without drug vs. uninfected without drug, down.

11.9.2 z2.2 with and without drug

z22sb_vs_z22nosb_volcano[["plot"]] +
  xlim(-8, 8) +
  ylim(0, 210)

## Error: object 'z22sb_vs_z22nosb_volcano' not found

pp(file = "images/z22sb_vs_z22nosb_reactome_up.png",
   image = all_gp[["z22sb_vs_z22nosb_up"]][["pvalue_plots"]][["REAC"]],
   height = 12, width = 9)

## Error: object 'all_gp' not found

all_gp[["z22sb_vs_z22nosb_up"]][["pvalue_plots"]][["REAC"]]

## Error: object 'all_gp' not found

## Reactome, zymodeme2.3 without drug vs. uninfected without drug, up.
all_gp[["z22sb_vs_z22nosb_up"]][["pvalue_plots"]][["KEGG"]]

## Error: object 'all_gp' not found

## KEGG, zymodeme2.3 without drug vs. uninfected without drug, up.
all_gp[["z22sb_vs_z22nosb_up"]][["pvalue_plots"]][["MF"]]

## Error: object 'all_gp' not found

## MF, zymodeme2.3 without drug vs. uninfected without drug, up.
all_gp[["z22sb_vs_z22nosb_up"]][["pvalue_plots"]][["TF"]]

## Error: object 'all_gp' not found

## TF, zymodeme2.3 without drug vs. uninfected without drug, up.
all_gp[["z22sb_vs_z22nosb_up"]][["pvalue_plots"]][["WP"]]

## Error: object 'all_gp' not found

## WikiPathways, zymodeme2.3 without drug vs. uninfected without drug, up.
all_gp[["z22sb_vs_z22nosb_up"]][["interactive_plots"]][["WP"]]

## Error: object 'all_gp' not found

all_gp[["z22sb_vs_z22nosb_down"]][["pvalue_plots"]][["REAC"]]

## Error: object 'all_gp' not found

## Reactome, zymodeme2.3 without drug vs. uninfected without drug, down.
all_gp[["z22sb_vs_z22nosb_down"]][["pvalue_plots"]][["MF"]]

## Error: object 'all_gp' not found

## MF, zymodeme2.3 without drug vs. uninfected without drug, down.
all_gp[["z22sb_vs_z22nosb_down"]][["pvalue_plots"]][["TF"]]

## Error: object 'all_gp' not found

## TF, zymodeme2.3 without drug vs. uninfected without drug, down.

11.9.3 Shared and unique genes after/before drug

shared <- Vennerable::Venn(list(
  "z23" = rownames(hs_macr_sig[["deseq"]][["ups"]][["z23sb_vs_z23nosb"]]),
  "z22" = rownames(hs_macr_sig[["deseq"]][["ups"]][["z22sb_vs_z22nosb"]])))

## Error: object 'hs_macr_sig' not found

pp(file = "images/z23_z22_drug_venn_up.png")
Vennerable::plot(shared)

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

dev.off()

## png 
##   2

Vennerable::plot(shared)

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

shared <- Vennerable::Venn(list(
  "z23" = rownames(hs_macr_sig[["deseq"]][["downs"]][["z23sb_vs_z23nosb"]]),
  "z22" = rownames(hs_macr_sig[["deseq"]][["downs"]][["z22sb_vs_z22nosb"]])))

## Error: object 'hs_macr_sig' not found

pp(file = "images/z23_z22_drug_venn_down.png")
Vennerable::plot(shared)

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

dev.off()

## png 
##   2

Vennerable::plot(shared)

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

Note: I am settig the x and y-axis boundaries by allowing the plotter to pick its own axis the first time, writing down the ranges I observe, and then setting them to the largest of the pair. It is therefore possible that I missed one or more genes which lies outside that range.

The previous plotted contrasts sought to show changes between the two strains z2.3 and z2.2. Conversely, the previous volcano plots seek to directly compare each strain before/after drug treatment.

12 LRT of the Human Macrophage

A slightly different tack to examine the data is to perform a likelihood ratio test in order to look for trends which are shared among genes when examining different conditions in the data.

tmrc2_lrt_strain_drug <- deseq_lrt(hs_macr, interactor_column = "drug",
                                   interest_column = "macrophagezymodeme",
                                   factors = c("drug", "macrophagezymodeme"))

## converting counts to integer mode

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 38 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

## rlog() may take a long time with 50 or more samples,
## vst() is a much faster transformation

## Working with 858 genes.

## Working with 855 genes after filtering: minc > 3

## Joining with `by = join_by(merge)`
## Joining with `by = join_by(merge)`

tmrc2_lrt_strain_drug[["cluster_data"]][["plot"]]

13 Parasite

Let us consider for a moment differences among the parasite transcriptomes for the samples which were not drug treated.

One thing I did in the initial implementation of this document was to repeat the variable ‘up_genes’ for each comparison; I think this time I will make a different variable for each comparison so I can play with them a little further.

comparison <- "z23_vs_z22"
lp_macrophage_de <- all_pairwise(lp_macrophage_nosb, model_svs = "svaseq",
                                 model_fstring = "~ 0 + condition", filter = TRUE)

## z2.2 z2.3 
##   14   15

## Removing 119 low-count genes (8591 remaining).

## Basic step 0/3: Normalizing data.

## Basic step 0/3: Converting data.

## I think this is failing? SummarizedExperiment

## Basic step 0/3: Transforming data.

## Setting 2387 entries to zero.

## This received a matrix of SVs.

## converting counts to integer mode

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## Warning in createContrastL(objFlt$formula, objFlt$data, L): Contrasts with only
## a single non-zero term are already evaluated by default.

## conditions
## z22 z23 
##  14  15

## conditions
## z22 z23 
##  14  15

## conditions
## z22 z23 
##  14  15

tmrc2_parasite_keepers <- list(
  "z23_vs_z22" = c("z23", "z22"))
lp_macrophage_table <- combine_de_tables(
  lp_macrophage_de, keepers = tmrc2_parasite_keepers,
  excel = glue("analyses/macrophage_de/de_tables/parasite_infection_de-v{ver}.xlsx"))

## Looking for subscript invalid names, start of extract_keepers.

## Error in combine_extracted_plots(entry_name, combined, wanted_denominator, : object 'color_low' not found

lp_macrophage_sig <- extract_significant_genes(
  lp_macrophage_table,
  excel = glue("analyses/macrophage_de/sig_tables/parasite_sig-v{ver}.xlsx"))

## Error: object 'lp_macrophage_table' not found

lp_macrophage_table[["plots"]][[comparison]][["deseq_vol_plots"]]

## Error: object 'lp_macrophage_table' not found

lp_macrophage_table[["plots"]][[comparison]][["deseq_ma_plots"]]

## Error: object 'lp_macrophage_table' not found

up_genes_z23z22 <- lp_macrophage_sig[["deseq"]][["ups"]][[comparison]]

## Error: object 'lp_macrophage_sig' not found

dim(up_genes_z23z22)

## Error: object 'up_genes_z23z22' not found

down_genes_z23z22 <- lp_macrophage_sig[["deseq"]][["downs"]][[comparison]]

## Error: object 'lp_macrophage_sig' not found

dim(down_genes_z23z22)

## Error: object 'down_genes_z23z22' not found

lp_z23sb_vs_z22sb_volcano <- plot_volcano_de(
  table = lp_macrophage_table[["data"]][["z23_vs_z22"]],
  fc_col = "deseq_logfc", p_col = "deseq_adjp",
  shapes_by_state = FALSE, color_by = "fc",  label = 10, label_column = "hgnc_symbol")

## Error in plot_volcano_de(table = lp_macrophage_table[["data"]][["z23_vs_z22"]], : could not find function "plot_volcano_de"

plotly::ggplotly(lp_z23sb_vs_z22sb_volcano[["plot"]])

## Error: object 'lp_z23sb_vs_z22sb_volcano' not found

lp_z23sb_vs_z22sb_volcano[["plot"]]

## Error: object 'lp_z23sb_vs_z22sb_volcano' not found

13.1 Ontology search

An important note, I recently added a minimum crossreference argument (defaulting to 40 genes), which causes lots of comparisons to fail for poorly annotated genomes (like panamensis.) Thus, I am relaxing that constraint for these searches.

lp_lengths <- all_lp_annot[, c("gid", "annot_cds_length")]
colnames(lp_lengths)  <- c("ID", "length")

up_goseq <- simple_goseq(up_genes_z23z22, go_db = lp_go,
                         length_db = lp_lengths, min_xref = 15)

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

## View categories over represented in the 2.3 samples
up_goseq[["pvalue_plots"]][["bpp_plot_over"]]

## Error: object 'up_goseq' not found

down_goseq <- simple_goseq(down_genes_z23z22, go_db = lp_go,
                           length_db = lp_lengths, min_xref = 15)

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

## View categories over represented in the 2.2 samples
down_goseq[["pvalue_plots"]][["bpp_plot_over"]]

## Error: object 'down_goseq' not found

created <- dir.create(glue("analyses/macrophage_de/goseq_parasite"))
written_goseq <- write_goseq_data(
  up_goseq,
  excel = glue("analyses/macrophage_de/goseq_parasite/lp_macrophage_increased_z2.3_goseq-v{ver}.xlsx"))

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

written_goseq <- write_goseq_data(
  down_goseq,
  excel = glue("analyses/macrophage_de/goseq_parasite/lp_macrophage_increased_z2.2_goseq-v{ver}.xlsx"))

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

14 GSVA

Note: The following block assumes one is able to download a fresh copy of msigDB, which I am not sure is possible within the constraints of a container (I mean it is trivial to do, but I am not sure if it is ok due to licensing). However, Broad provides a data package of a msigdb release. As a result, the following block will be repeated using that.

hs_infected <- subset_se(hs_macrophage, subset = "macrophagetreatment!='uninf'") %>%
  subset_se(subset = "macrophagetreatment!='uninf_sb'")
hs_gsva_c2 <- simple_gsva(hs_infected)
hs_gsva_c2_meta <- get_msigdb_metadata(hs_gsva_c2, msig_xml = "reference/msigdb_v7.2.xml")
hs_gsva_c2_sig <- get_sig_gsva_categories(
  hs_gsva_c2_meta,
  excel = glue("analyses/macrophage_de/gsva/hs_macrophage_gsva_c2_sig.xlsx"))
hs_gsva_c2_sig[["raw_plot"]]

hs_gsva_c7 <- simple_gsva(hs_infected, signature_category = "c7")
hs_gsva_c7_meta <- get_msigdb_metadata(hs_gsva_c7, msig_xml = "reference/msigdb_v7.2.xml")
hs_gsva_c7_sig <- get_sig_gsva_categories(
  hs_gsva_c7,
  excel = glue("analyses/macrophage_de/gsva/hs_macrophage_gsva_c7_sig.xlsx"))
hs_gsva_c7_sig[["raw_plot"]]

14.1 Repeat using the GSVAdata package.

hs_infected <- subset_se(hs_macrophage, subset = "macrophagetreatment!='uninf'") %>%
  subset_se(subset = "macrophagetreatment!='uninf_sb'")
hs_gsva_c2 <- simple_gsva(hs_infected)

## Error: unable to find an inherited method for function 'annotation' for signature 'object = "SummarizedExperiment"'

##hs_gsva_c2_meta <- get_msigdb_metadata(hs_gsva_c2, msig_xml="reference/msigdb_v7.2.xml")
hs_gsva_c2_sig <- get_sig_gsva_categories(
  hs_gsva_c2,
  excel = glue("analyses/macrophage_de/gsva/hs_macrophage_gsva_c2_sig.xlsx"))

## Error: object 'hs_gsva_c2' not found

hs_gsva_c2_sig[["raw_plot"]]

## Error: object 'hs_gsva_c2_sig' not found

hs_gsva_c7 <- simple_gsva(hs_infected, signature_category = "c7")

## Error: unable to find an inherited method for function 'annotation' for signature 'object = "SummarizedExperiment"'

##hs_gsva_c7_meta <- get_msigdb_metadata(hs_gsva_c7, msig_xml="reference/msigdb_v7.2.xml")
hs_gsva_c7_sig <- get_sig_gsva_categories(
  hs_gsva_c7,
  excel = glue("analyses/macrophage_de/gsva/hs_macrophage_gsva_c7_sig.xlsx"))

## Error: object 'hs_gsva_c7' not found

hs_gsva_c7_sig[["raw_plot"]]

## Error: object 'hs_gsva_c7_sig' not found

15 Try out a new tool

Two reasons: Najib loves him some PCA, this uses wikipathways, which is something I think is neat.

Ok, I spent some time looking through the code and I have some problems with some of the design decisions.

Most importantly, it requires a data.frame() which has the following format:

1. No rownames, instead column #1 is the sample ID.
2. Columns 2-m are the categorical/survival/etc metrics.
3. Columns m-n are 1 gene-per-column with log2 values.

But when I think about it I think I get the idea, they want to be able to do modelling stuff more easily with response factors.

library(pathwayPCA)
library(rWikiPathways)

## 
## Attaching package: 'rWikiPathways'

## The following object is masked from 'package:edgeR':
## 
##     getCounts

downloaded <- downloadPathwayArchive(organism = "Homo sapiens", format = "gmt")
data_path <- system.file("extdata", package = "pathwayPCA")
wikipathways <- read_gmt(paste0(data_path, "/wikipathways_human_symbol.gmt"),
                         description = TRUE)

se <- subset_se(hs_macrophage, subset = "macrophagetreatment!='uninf'") %>%
  subset_se(subset = "macrophagetreatment!='uninf_sb'")
se <- set_conditions(se, fact = "macrophagezymodeme")

## The numbers of samples by condition are:

## 
## none  z22  z23 
##    0   29   29

symbol_column <- "hgnc_symbol"
symbol_vector <- rowData(se)[[symbol_column]]
names(symbol_vector) <- rownames(rowData(se))
symbol_df <- as.data.frame(symbol_vector)

assay_df <- merge(symbol_df, as.data.frame(assay(se)), by = "row.names")
assay_df[["Row.names"]] <- NULL
rownames(assay_df) <- make.names(assay_df[["symbol_vector"]], unique = TRUE)
assay_df[["symbol_vector"]] <- NULL
assay_df <- as.data.frame(t(assay_df))
assay_df[["SampleID"]] <- rownames(assay_df)
assay_df <- dplyr::select(assay_df, "SampleID", everything())

factor_df <- as.data.frame(colData(se))
factor_df[["SampleID"]] <- rownames(factor_df)
factor_df <- dplyr::select(factor_df, "SampleID", everything())
factor_df <- factor_df[, c("SampleID", factors)]

## Error: object 'factors' not found

tt <- CreateOmics(
  assayData_df = assay_df,
  pathwayCollection_ls = wikipathways,
  response = factor_df,
  respType = "categorical",
  minPathSize = 5)

## 3190 genes have variance < epsilon and will be removed. These gene(s) include:

##   [1] "TNMD"         "CYP51A1"      "KRIT1"        "MAD1L1"       "ARF5"        
##   [6] "REXO5"        "FBXL3"        "REX1BD"       "KRT33A"       "TAC1"        
##  [11] "LGALS14"      "SLC13A2"      "TRAPPC6A"     "SELE"         "TFAP2B"      
##  [16] "SS18L2"       "IDS"          "SLC7A14"      "CLDN11"       "MDH1"        
##  [21] "COX15"        "MATR3"        "ISL1"         "INSRR"        "EFCAB1"      
##  [26] "TMSB10"       "OTC"          "HOXC8"        "XK"           "NOP16"       
##  [31] "TNFRSF17"     "GUCA1A"       "NNAT"         "NRIP2"        "MCOLN3"      
##  [36] "SERPINB3"     "MRPS24"       "SEZ6"         "AHRR"         "BORCS8.MEF2B"
##  [41] "KDM4A"        "THUMPD1"      "IFT80"        "ERLEC1"       "PAGE1"       
##  [46] "FRMPD1"       "LNX1"         "IPCEF1"       "ZNF37A"       "TUBA3D"      
##  [51] "SPAG5"        "EXOSC5"       "TIGAR"        "TP53INP2"     "LXN"         
##  [56] "AFM"          "CFHR2"        "UBA5"         "JMJD4"        "PCDHA6"      
##  [61] "PCDHGA2"      "C1QTNF3"      "RNF13"        "ZNF671"       "RRN3"        
##  [66] "CHERP"        "DIMT1"        "NME8"         "PIGS"         "DEFB127"     
##  [71] "FXYD3"        "CMTM1"        "FLT3LG"       "RBM27"        "ANGPT2"      
##  [76] "RNF31"        "SEMA4G"       "NUBP2"        "KCNK16"       "MAGEB2"      
##  [81] "MTAP"         "SERPIND1"     "DDT"          "SEC14L2"      "GGT1"        
##  [86] "PRODH"        "SOX10"        "TIMP3"        "PSMA3"        "SNW1"        
##  [91] "SERPINA4"     "PCK2"         "PRORP"        "TM9SF1"       "RAB5IF"      
##  [96] "CST9L"        "CST4"         "SPINT3"       "EPPIN"        "RBFA"        
## [101] "CEP76"        "H2BW2"        "MCTS2P"       "SRPX"         "F9"          
## [106] "PPP1R2C"      "BRS3"         "TIMP1"        "GLA"          "ACP5"        
## [111] "DHRS12"       "ZNF821"       "CMC2"         "ZNF174"       "CORO7.PAM16" 
## [116] "SALL1"        "AQP9"         "OIP5"         "ARHGEF10"     "CGB2"        
## [121] "CGB3"         "PPP1R37"      "RNASEH2A"     "OAZ1"         "C19orf44"    
## [126] "MED26"        "ZNF419"       "LGALS13"      "CEACAM5"      "BABAM1"      
## [131] "ATP1A3"       "ZNRF4"        "TMEM205"      "WDR83OS"      "PIK3R2"      
## [136] "PDE4C"        "DDX49"        "ERF"          "RASA4"        "TFPI2"       
## [141] "DUS4L"        "NPVF"         "WNT2"         "HOXA5"        "HOXA6"       
## [146] "CHN2"         "MINDY4"       "PSMA2"        "OGN"          "ASPN"        
## [151] "ECM2"         "EXOSC3"       "VSIR"         "RNF43"        "ASPA"        
## [156] "HOXB6"        "SLC16A6"      "RANGRF"       "VTN"          "FOXN1"       
## [161] "UNC119"       "ALDOC"        "ODAM"         "SMR3A"        "CHIC2"       
## [166] "IL2"          "CPZ"          "DBX1"         "SNX15"        "APOC3"       
## [171] "SCGB2A2"      "PTPMT1"       "CALCA"        "MYF6"         "MYF5"        
## [176] "PRR4"         "AKAP3"        "GSG1"         "OGFOD2"       "ART4"        
## [181] "MGP"          "FZD10"        "LPCAT3"       "GYS2"         "MAK"         
## [186] "ASF1A"        "IL17A"        "TSPO2"        "CCNC"         "HDGFL1"      
## [191] "OR12D3"       "MRPL2"        "TMCO6"        "PDE8B"        "IL5"         
## [196] "SMC4"         "NPHP3"        "BCHE"         "ABHD14B"      "ABHD14A.ACY1"
## [201] "TP53I3"       "INO80B"       "REG1A"        "KCNJ13"       "NEU2"        
## [206] "HSPE1"        "ABCB6"        "PNO1"         "ATP6V1B1"     "ANGPTL1"     
## [211] "NCF2"         "PRAMEF1"      "AGMAT"        "TNNI3K"       "TSNAX"       
## [216] "CRYGD"        "ZC2HC1B"      "CCND2"        "FGF23"        "TRIM32"      
## [221] "TGFB3"        "ZNF410"       "GPR75"        "IFIT3"        "NKX2.3"      
## [226] "IFIT2"        "HOXB8"        "HOXB5"        "CRHR1"        "HOXB1"       
## [231] "MLANA"        "IFNA6"        "GRIA2"        "LRP11"        "MAGEB4"      
## [236] "SLC25A2"      "GPR31"        "TNFSF11"      "TRIM6"        "TAS2R10"     
## [241] "IAPP"         "SLITRK3"      "CLCC1"        "GPSM2"        "OBP2A"       
## [246] "TBX22"        "PRM2"         "RWDD3"        "MRM2"         "SLC25A51"    
## [251] "DCAF10"       "WDR83"        "ACTRT1"       "TSFM"         "ORMDL2"      
## [256] "CDK2"         "KBTBD4"       "COL10A1"      "SERPINA7"     "H2BW1"       
## [261] "ESX1"         "B9D2"         "MC3R"         "GCNT7"        "ANKRD60"     
## [266] "C20orf85"     "TP53TG5"      "MAGEA10"      "FASTKD3"      "CRISP2"      
## [271] "AARS2"        "RPS10"        "APOBEC2"      "GCM2"         "TRIM51"      
## [276] "SCGB2A1"      "GPR18"        "IRF1"         "AMELX"        "NPBWR2"      
## [281] "BHLHE23"      "FOSB"         "DEFB126"      "FOXA2"        "NKX2.4"      
## [286] "NKX2.2"       "CSTL1"        "FLRT3"        "MGME1"        "TMEM74B"     
## [291] "CITED1"       "MMP24"        "TMEM115"      "KIRREL2"      "X.5"         
## [296] "STATH"        "HTN1"         "TIMM17B"      "EVI2A"        "OMG"         
## [301] "AVPR2"        "OMD"          "TAS2R3"       "TAS2R4"       "OR7A10"      
## [306] "SLC35E1"      "OR7C2"        "GNG13"        "EMC6"         "OR1E2"       
## [311] "FGL2"         "GNAZ"         "ADORA2A"      "TAS2R16"      "ATP6V1F"     
## [316] "LRRC4"        "LRRC17"       "FEZF1"        "MRPS12"       "HOXD1"       
## [321] "HAT1"         "HOXD9"        "HOXD13"       "ELL3"         "CALML4"      
## [326] "THAP10"       "ACKR4"        "SOX15"        "KLK8"         "NEDD8"       
## [331] "VCY1B"        "VCY"          "CDY2B"        "INS.IGF2"     "KCNA5"       
## [336] "ANGPTL8"      "ACE2"         "GDF1"         "MRPL34"       "LSM7"        
## [341] "ACSBG2"       "BMP15"        "ARPC1B"       "OR11H1"       "CALY"        
## [346] "PPAN"         "HSD17B3"      "PRRG1"        "BPIFA3"       "DEFB118"     
## [351] "GJA9"         "CDX4"         "NAPSA"        "PDLIM4"       "TMEM204"     
## [356] "KRT33B"       "FSHB"         "USP29"        "NR0B2"        "ACTR10"      
## [361] "ABHD12B"      "RTBDN"        "TRIM22"       "TIMM10B"      "SCLY"        
## [366] "FTHL17"       "NIP7"         "VPS4A"        "SCP2D1"       "SSTR4"       
## [371] "APCS"         "TOE1"         "PPP1R3D"      "BHMT2"        "ZBED3"       
## [376] "ANGPTL3"      "STOML3"       "IRS4"         "ERG28"        "GSC"         
## [381] "CAMK1"        "GSTM1"        "TSHB"         "GSTM3"        "FKBP11"      
## [386] "GRP"          "PRH2"         "SOX3"         "BIVM"         "ERCC5"       
## [391] "UGT2A3"       "CSN2"         "LACRT"        "GLS2"         "FAM186B"     
## [396] "BLOC1S1"      "ZC3H10"       "SLC26A10"     "MIP"          "CHST5"       
## [401] "HTR2B"        "TMBIM1"       "RCBTB2"       "KDELR2"       "CIDEB"       
## [406] "NKX2.8"       "NKX2.1"       "SRSF1"        "CHAD"         "GH2"         
## [411] "LIMD2"        "CFC1"         "OR13C9"       "ANGPTL2"      "TLR4"        
## [416] "HINT2"        "YIPF3"        "CLPS"         "FXYD6"        "SLTM"        
## [421] "LHCGR"        "SLC3A1"       "LBX1"         "CUZD1"        "RBP4"        
## [426] "GPR87"        "MSTN"         "CARF"         "IL21"         "GSTCD"       
## [431] "ZCRB1"        "NDUFA9"       "KERA"         "SYCP3"        "CCDC65"      
## [436] "NPFF"         "LPAR6"        "RNF113B"      "SSTR1"        "TSSK4"       
## [441] "RAB15"        "SERF2"        "FGF7"         "CELF6"        "IGSF6"       
## [446] "CHST4"        "ZSCAN32"      "CTRL"         "KIF2B"        "ADCYAP1"     
## [451] "MEP1B"        "SAT2"         "ZNF750"       "ELAC1"        "SLC39A3"     
## [456] "CCDC97"       "TMEM91"       "ZNF593"       "EVA1B"        "DMRTB1"      
## [461] "BARHL2"       "PIGM"         "CRNN"         "JTB"          "CREB3L4"     
## [466] "PYCR2"        "CHAC2"        "ANKRD53"      "TEX261"       "LIPT1"       
## [471] "LYG1"         "GPR17"        "PHOSPHO2"     "RPL32"        "LRTM1"       
## [476] "ZNF660"       "EIF2A"        "AMT"          "ECE2"         "SLC26A1"     
## [481] "CABS1"        "BHMT"         "KCNMB1"       "LRRTM2"       "SLC17A4"     
## [486] "H2BC1"        "HIGD2A"       "TCTE1"        "CLVS2"        "TAAR2"       
## [491] "TAAR6"        "TAAR8"        "WTAP"         "RBAK"         "FERD3L"      
## [496] "TMEM140"      "CLTRN"        "LANCL3"       "SYTL5"        "AKAP4"

## 1103 gene name(s) are invalid. Invalid name(s) include:

##   [1] "NME1.NME2"       "RTEL1.TNFRSF6B"  "STON1.GTF2A1L"   "X.1"            
##   [5] "PTGES3L.AARSD1"  "NKX3.2"          "X.2"             "TMEM189.UBE2V1" 
##   [9] "H1.3"            "X.3"             "H1.1"            "X.4"            
##  [13] "CHURC1.FNTB"     "X.6"             "H3.3B"           "ZNF670.ZNF695"  
##  [17] "X.7"             "ERVK3.1"         "X.8"             "X.9"            
##  [21] "NKX6.2"          "X.10"            "H3.3A"           "NKX6.1"         
##  [25] "NKX6.3"          "NKX3.1"          "X.11"            "X.12"           
##  [29] "H3.4"            "H1.4"            "JMJD7.PLA2G4B"   "X.14"           
##  [33] "KRTAP4.4"        "RAB4B.EGLN2"     "X.15"            "X.16"           
##  [37] "H1.8"            "HLA.DQB1"        "X.17"            "NKX2.6"         
##  [41] "KRTAP9.7"        "KRTAP11.1"       "NKX2.5"          "KRTAP8.1"       
##  [45] "X.19"            "KRTAP19.1"       "H1.5"            "KRTAP6.1"       
##  [49] "H1.10"           "KRTAP5.5"        "KRTAP17.1"       "KRTAP21.2"      
##  [53] "H1.7"            "X.20"            "H1.6"            "H1.2"           
##  [57] "X.21"            "H3.5"            "X.22"            "H1.0"           
##  [61] "HLA.DRB1"        "KRTAP5.3"        "HLA.DQA1"        "X.23"           
##  [65] "H4.16"           "X.25"            "KRTAP4.1"        "HLA.DRB5"       
##  [69] "MT.ND6"          "MT.CO2"          "MT.CYB"          "MT.ND2"         
##  [73] "MT.ND5"          "MT.CO1"          "MT.ND3"          "MT.ND4"         
##  [77] "MT.ND1"          "MT.ATP6"         "MT.CO3"          "X.26"           
##  [81] "HLA.DOA"         "HLA.DMA"         "HLA.DRA"         "X.28"           
##  [85] "HLA.C"           "KRTAP5.11"       "KRTAP5.10"       "HLA.E"          
##  [89] "HLA.G"           "HLA.F"           "X.29"            "CLLU1.AS1"      
##  [93] "X.30"            "TRBV20OR9.2"     "KRTAP5.6"        "KRTAP5.2"       
##  [97] "KRTAP5.1"        "KRTAP19.8"       "HLA.A"           "X.31"           
## [101] "IGKV4.1"         "IGKV6.21"        "IGKV3D.20"       "IGLV10.54"      
## [105] "IGLV5.52"        "IGLV1.51"        "IGLV1.50"        "IGLV1.47"       
## [109] "IGLV7.46"        "IGLV1.44"        "IGLV7.43"        "IGLV1.40"       
## [113] "IGLV3.25"        "IGLV2.23"        "IGLV3.21"        "IGLV3.19"       
## [117] "IGLV3.16"        "IGLV2.14"        "IGLV2.11"        "IGLV3.9"        
## [121] "IGLV3.1"         "TRBV7.3"         "TRBV5.3"         "TRBV10.1"       
## [125] "TRBV6.5"         "TRBV6.6"         "TRBV7.6"         "TRBV5.1"        
## [129] "TRBV20.1"        "TRBV24.1"        "TRBJ2.1"         "TRBJ2.2P"       
## [133] "TRBJ2.3"         "TRBJ2.6"         "TRBJ2.7"         "TRAV12.3"       
## [137] "TRAV8.7"         "IGHD3.10"        "IGHV6.1"         "IGHV1.2"        
## [141] "IGHV1.3"         "IGHV2.5"         "IGHV3.7"         "IGHV3.11"       
## [145] "IGHV3.13"        "IGHV3.15"        "IGHV1.18"        "IGHV3.20"       
## [149] "IGHV3.21"        "IGHV3.23"        "IGHV1.24"        "IGHV2.26"       
## [153] "IGHV4.28"        "IGHV3.33"        "IGHV4.34"        "IGHV4.39"       
## [157] "IGHV3.49"        "IGHV5.51"        "IGHV3.66"        "IGHV3.73"       
## [161] "KRTAP16.1"       "KRTAP3.3"        "KRTAP3.2"        "MT.ND4L"        
## [165] "X.32"            "ZNF625.ZNF20"    "ERV3.1"          "X.33"           
## [169] "RPL17.C18orf32"  "KRTAP1.5"        "ZNF816.ZNF321P"  "IGHV3.64"       
## [173] "HLA.DPB1"        "IGHV4.59"        "IGHV3.74"        "APOC4.APOC2"    
## [177] "X.36"            "X.38"            "ERVMER34.1"      "X.39"           
## [181] "MT.ATP8"         "IGKV3D.7"        "TRBV5.4"         "X.40"           
## [185] "IGKV1OR2.108"    "HLA.DPA1"        "IGHV3.43"        "HLA.DQB2"       
## [189] "TRBV29.1"        "X.41"            "IGKV3OR2.268"    "HLA.B"          
## [193] "HNRNPUL2.BSCL2"  "NKX1.1"          "X.44"            "HLA.DQA2"       
## [197] "IGKV2D.30"       "IGKV1D.8"        "IGKV1.6"         "X.47"           
## [201] "IGKV3.20"        "IGKV1D.33"       "IGKV1.17"        "IGKV1.8"        
## [205] "IGKV1.16"        "HLA.DOB"         "KRTAP5.8"        "IGKV2.24"       
## [209] "IGKV3.11"        "X.48"            "KRTAP5.4"        "IGKV1.9"        
## [213] "X.50"            "IGKV1.33"        "IGKV1.39"        "IGKV2D.28"      
## [217] "HLA.DMB"         "IGKV1D.17"       "ERVW.1"          "PPAN.P2RY11"    
## [221] "IGKV2.30"        "IGKV2D.29"       "IGKV1.12"        "IGKV1.5"        
## [225] "X.51"            "X.52"            "DNAJC25.GNG10"   "KRTAP5.7"       
## [229] "IGKV3.15"        "KRTAP4.2"        "IGKV1.27"        "TRIM39.RPP21"   
## [233] "X.54"            "PRR5.ARHGAP8"    "STIMATE.MUSTN1"  "RBM14.RBM4"     
## [237] "LY75.CD302"      "X.55"            "X.56"            "TNFSF12.TNFSF13"
## [241] "ATP5MF.PTCD1"    "X.57"            "EPPIN.WFDC6"     "X.58"           
## [245] "X.59"            "X.60"            "X.61"            "X.62"           
## [249] "X.63"            "X.64"            "RNF103.CHMP3"    "X.66"           
## [253] "ARPIN.AP3S2"     "ARPC4.TTLL3"     "X.67"            "X.68"           
## [257] "X.69"            "LY6G6F.LY6G6D"   "X.70"            "CCDC169.SOHLH2" 
## [261] "NT5C1B.RDH14"    "X.71"            "X.72"            "X.73"           
## [265] "TMED7.TICAM2"    "X.74"            "MSANTD3.TMEFF1"  "X.75"           
## [269] "CENPS.CORT"      "X.76"            "X.77"            "TRBV7.4"        
## [273] "CHKB.CPT1B"      "X.78"            "X.79"            "X.80"           
## [277] "X.81"            "X.82"            "X.83"            "CKLF.CMTM1"     
## [281] "ATP6V1G2.DDX39B" "INMT.MINDY4"     "X.85"            "STX16.NPEPL1"   
## [285] "KRTAP5.9"        "X.86"            "SAA2.SAA4"       "ZFP91.CNTF"     
## [289] "X.87"            "MSH5.SAPCD1"     "FXYD6.FXYD2"     "X.88"           
## [293] "X.89"            "X.90"            "X.91"            "X.92"           
## [297] "X.93"            "NEDD8.MDP1"      "TRAV1.1"         "X.94"           
## [301] "X.96"            "X.97"            "KLRC4.KLRK1"     "X.98"           
## [305] "X.99"            "X.100"           "X.101"           "X.102"          
## [309] "X.103"           "X.104"           "TRAV1.2"         "X.107"          
## [313] "X.108"           "X.109"           "X.110"           "X.111"          
## [317] "X.112"           "X.113"           "SLCO1B3.SLCO1B7" "X.114"          
## [321] "X.115"           "X.117"           "X.118"           "X.120"          
## [325] "X.121"           "X.122"           "RPL36A.HNRNPH2"  "X.123"          
## [329] "X.124"           "P2RX5.TAX1BP3"   "X.125"           "X.126"          
## [333] "X.127"           "PPT2.EGFL8"      "X.129"           "X.130"          
## [337] "X.131"           "X.132"           "X.133"           "X.134"          
## [341] "SPECC1L.ADORA2A" "BCL2L2.PABPN1"   "X.137"           "X.138"          
## [345] "PINX1.1"         "X.140"           "X.141"           "X.142"          
## [349] "X.143"           "UBE2F.SCLY"      "X.144"           "FPGT.TNNI3K"    
## [353] "BLOC1S5.TXNDC5"  "X.146"           "POC1B.GALNT4"    "NDUFC2.KCTD14"  
## [357] "X.147"           "ZHX1.C8orf76"    "X.150"           "ST20.MTHFS"     
## [361] "X.151"           "TGIF2.RAB5IF"    "X.152"           "X.153"          
## [365] "X.155"           "X.156"           "X.158"           "X.159"          
## [369] "X.160"           "X.161"           "X.162"           "PMF1.BGLAP"     
## [373] "X.163"           "X.164"           "X.165"           "X.166"          
## [377] "X.169"           "X.170"           "X.171"           "X.172"          
## [381] "X.173"           "X.174"           "X.175"           "X.176"          
## [385] "TEN1.CDK3"       "X.177"           "X.178"           "X.179"          
## [389] "X.180"           "X.181"           "X.182"           "ISY1.RAB43"     
## [393] "X.183"           "X.184"           "X.185"           "X.186"          
## [397] "TMEM256.PLSCR3"  "X.191"           "X.192"           "X.193"          
## [401] "X.195"           "X.196"           "LINC02210.CRHR1" "X.197"          
## [405] "X.198"           "X.200"           "X.201"           "X.202"          
## [409] "X.203"           "X.205"           "CFAP298.TCP10L"  "X.206"          
## [413] "EEF1E1.BLOC1S5"  "X.207"           "X.208"           "X.209"          
## [417] "X.210"           "X.211"           "X.213"           "X.214"          
## [421] "X.215"           "X.216"           "X.217"           "X.218"          
## [425] "X.221"           "X.222"           "X.223"           "X.224"          
## [429] "X.225"           "X.226"           "X.227"           "X.228"          
## [433] "X.229"           "X.230"           "X.231"           "X.232"          
## [437] "X.233"           "X.234"           "X.235"           "X.236"          
## [441] "X.237"           "X.238"           "X.239"           "X.240"          
## [445] "X.241"           "X.242"           "X.244"           "X.245"          
## [449] "X.246"           "X.247"           "X.248"           "X.249"          
## [453] "X.250"           "X.251"           "X.252"           "X.253"          
## [457] "X.254"           "X.255"           "X.256"           "X.257"          
## [461] "X.258"           "X.259"           "X.260"           "X.261"          
## [465] "X.263"           "X.264"           "X.266"           "X.267"          
## [469] "X.268"           "X.269"           "X.270"           "MIA.RAB4B"      
## [473] "X.271"           "X.272"           "X.273"           "X.274"          
## [477] "X.275"           "X.276"           "X.277"           "X.279"          
## [481] "X.280"           "X.281"           "X.283"           "X.285"          
## [485] "X.286"           "X.288"           "ARHGAP19.SLIT1"  "COMMD3.BMI1"    
## [489] "ZNF559.ZNF177"   "X.289"           "TSNAX.DISC1"     "X.290"          
## [493] "X.291"           "X.292"           "BORCS7.ASMT"     "IGHV3.30"       
## [497] "URGCP.MRPS24"    "RPS10.NUDT3"     "TLCD4.RWDD3"     "X.293"

## These genes may be excluded from analysis. Proper gene names
## contain alphanumeric characters only, and start with a letter.

## Warning in CheckSampleIDs(assayData_df): Row names will be ignored. Sample IDs must be in the first column of the data
##   frame.

## Error in .convertPhenoDF(response, type = respType): Regression and categorical data must be a data frame with two columns, sample ID
##   and response, in exactly that order.

super <- AESPCA_pVals(
  object = tt,
  numPCs = 2,
  parallel = FALSE,
  numCores = 8,
  numReps = 2,
  adjustment = "BH")

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

16 Evaluating a log2FC barplot

Figure 2E is now comprised of a plot which shows log2FC values with error bars for selected genes and seeks to show differences between 2.3/uninfected and 2.2/uninfected.

Here is the table Olga used to generate it:

I went looking in the xlsx files produced in 202405 and found that these are the log2FC values and standard errors produced by DESeq2.

It should be noted that in my most recent version of these analyses, these numbers did shift slightly. I am looking into that now.

• Data witout drug

** 2.3 vs Uninfected MØ 2.2 vs Uninfected MØ

16.1 | Gene | Mean | SEM | n | Mean | SEM |n |

|IFI27 | 7.224 | 0.5662 |6 | 2.702 | 0.5669 | 6| |RSAD2 | 6.29 | 0.7312 |6 | 1.623 | 0.7303 | 6| |CCL8 | 6.225 | 0.928 |6 | -0.314| 0.941 | 6| |IFI44L| 5.895 | 0.612 |6 | 2.06 | 0.611 | 6| |OASL | 4.726 | 0.4974 |6 | 1.392 | 0.4973 | 6| |USP18 | 3.644 | 0.483 |6 | 0.999 | 0.4826 | 6| |IDO1 | 7.145 | 1.107 |6 | 1.257 | 1.141 | 6| |IDO2 | 3.935 | 1.3 |6 | 2.557 | 1.341 | 6| |KYNU | 1.07 | 0.2186 |6 | 0.0207| 0.2184 | 6| |AHR | 0.9382 | 0.2236 |6 | 0.5032| 0.2239 | 6| |IL4I1 | 2.593 | 0.4623 |6 | 0.039 | 0.4618 | 6| |SOD2 | 2.76 | 0.349 |6 | 0.4241| 0.3528 | 6| |NOTCH1| 0.7572| 0.275 |6 | 1.495 | 0.2744 | 6| |DLL1 | 0.8268| 0.5285 |6 | 3.455 | 0.5228 | 6| |DLL4 | 1.116 | 0.737 |6 | 4.243 | 0.71 | 6| |HES1 | -0.0183| 0.8599 |6 | 6.536 | 0.7973 | 6| |HEY1 | 0.5533| 0.5789 |6 | 4.181 | 0.6273 | 6|

Ok, I think I found a problem: The NOTCH1 value is actually the adjusted p-value.

• Transporters without drug

** 2.3 vs Uninfected MØ 2.2 vs Uninfected MØ

16.2 | Gene | Mean | SEM | n| Mean | SEM | n|

|ABCB1 | -2.354 | 0.442 | 6| -0.406| 0.431| 6| |ABCG4 | -3.715 | 0.648 | 6| -0.653| 0.630| 6| |ABCB5 | -1.192 | 0.380 | 6| 1.351 | 0.363| 6| |ABCA9 | 1.880 | 0.648 | 6| 3.444 | 0.637| 6| |ABCC2 | 0.454 | 0.321 | 6| 1.818 | 0.314| 6| |AQP2 | -1.191 | 0.529 | 6| 0.745 | 0.514| 6| |AQP3 | -0.940 | 0.402 | 6| 0.431 | 0.395| 6|

• Transporters with drug

** 2.3 vs Uninfected MØ 2.2 vs Uninfected MØ

16.3 |Gene | Mean | SEM | n| Mean | SEM | n |

|ABCB1 | -0.697| 0.349 | 6| -1.255| 0.337 | 6| |ABCG4 | 1.231 | 0.503 | 6| 0.547 | 0.484 | 6| |AQP2 | 0.816 | 0.399 | 6| 0.043 | 0.387 | 6| |AQP3 | -1.286| 0.320 | 6| -1.613| 0.309 | 6| |AQP8 | 0.634 | 0.370 | 6| 0.943 | 0.365 | 6|

Let us now see if I can recapitulate the plot…

nodrug_contrasts <- c("z23nosb_vs_uninf", "z22nosb_vs_uninf")
genes_no_drug <- c("IFI27", "RSAD2", "CCL8", "IFI44L", "OASL", "USP18", "IDO1", "IDO2", "KYNU", "AHR", "IL4I1", "SOD2", "NOTCH1", "DLL1", "DLL4", "HES1", "HEY1")
transporters_no_drug <- c("ABCB1", "ABCG4", "ABCB5", "ABCA9", "ABCC2", "AQP2", "AQP3")
drug_contrasts <- c("z23sb_vs_sb", "z22sb_vs_sb")
transporters_drug <- c("ABCB1", "ABCG4", "AQP2", "AQP3", "AQP8")

These values came out of the data structure called ‘hs_macr_table’

z23nosb_uninf_values <- hs_macr_table[["data"]][["z23nosb_vs_uninf"]]

## Error: object 'hs_macr_table' not found

gene_idx <- z23nosb_uninf_values[["hgnc_symbol"]] %in% genes_no_drug

## Error: object 'z23nosb_uninf_values' not found

nodrug_rows <-  z23nosb_uninf_values[gene_idx, ]

## Error: object 'z23nosb_uninf_values' not found

rownames(nodrug_rows) <- nodrug_rows[["hgnc_symbol"]]

## Error: object 'nodrug_rows' not found

z23_nodrug_values <- nodrug_rows[, c("deseq_logfc", "deseq_lfcse")]

## Error: object 'nodrug_rows' not found

z23_nodrug_values

## Error: object 'z23_nodrug_values' not found

z22nosb_uninf_values <- hs_macr_table[["data"]][["z22nosb_vs_uninf"]]

## Error: object 'hs_macr_table' not found

gene_idx <- z22nosb_uninf_values[["hgnc_symbol"]] %in% genes_no_drug

## Error: object 'z22nosb_uninf_values' not found

nodrug_rows <-  z22nosb_uninf_values[gene_idx, ]

## Error: object 'z22nosb_uninf_values' not found

rownames(nodrug_rows) <- nodrug_rows[["hgnc_symbol"]]

## Error: object 'nodrug_rows' not found

z22_nodrug_values <- nodrug_rows[, c("deseq_logfc", "deseq_lfcse")]

## Error: object 'nodrug_rows' not found

z22_nodrug_values

## Error: object 'z22_nodrug_values' not found

z23_nodrug_values[["state"]] <- "z23_vs_uninfected"

## Error: object 'z23_nodrug_values' not found

z22_nodrug_values[["state"]] <- "z22_vs_uninfected"

## Error: object 'z22_nodrug_values' not found

plot_df <- rbind.data.frame(as.data.frame(z23_nodrug_values), as.data.frame(z22_nodrug_values))

## Error: object 'z23_nodrug_values' not found

plot_df[["gene"]] <- rownames(plot_df)

## Error: object 'plot_df' not found

## I just realized that this is actually just a comparison of z23/z22
## we should just take the adjusted p-values from that contrast for this.
z23_z22_comparison <- hs_macr_table[["data"]][["z23nosb_vs_z22nosb"]]

## Error: object 'hs_macr_table' not found

nodrug_rows <- z23_z22_comparison[gene_idx, ]

## Error: object 'z23_z22_comparison' not found

nodrug_pvalues <- nodrug_rows[, c("deseq_p", "deseq_adjp")]

## Error: object 'nodrug_rows' not found

rownames(nodrug_pvalues) <- nodrug_rows[["hgnc_symbol"]]

## Error: object 'nodrug_rows' not found

nodrug_pvalues

## Error: object 'nodrug_pvalues' not found

ggplot(plot_df, aes(x = gene, y = deseq_logfc, fill = state)) +
  geom_bar(position = position_dodge(), stat = "identity") +
  geom_errorbar(aes(ymin = deseq_logfc - deseq_lfcse,
                    ymax = deseq_logfc + deseq_lfcse),
                width = 0.2, position = position_dodge(0.9)) +
  scale_fill_manual(values = c("#1B9E77", "#7570B3")) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5))

## Error: object 'plot_df' not found

comparison <- c("z23_vs_uninfected", "z22_vs_uninfected")
comparisons <- rep(list(comparison), nrow(plot_df) / 2)

## Error: object 'plot_df' not found

ggplot(plot_df, aes(x = gene, y = deseq_logfc, fill = state, add = deseq_lfcse, facet.by = "state")) +
  geom_bar(position = position_dodge(), stat = "identity") +
  geom_errorbar(aes(ymin = deseq_logfc - deseq_lfcse,
                    ymax = deseq_logfc + deseq_lfcse),
                width = 0.2, position = position_dodge(0.9)) +
  stat_compare_means() +
  stat_compare_means(comparisons = comparisons, label.y = rownames(z23_nodrug_values)) +
  scale_fill_manual(values = c("#1B9E77", "#7570B3")) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5))

## Error: object 'plot_df' not found

Excellent, the values now match up. Now I ust need to figure out why the stupid hgnc IDs got lost… I can see them in the hs_annot data structure, so I must have messed up when I regenered the input to the de. Ok, I got to the same starting point now with identical values. As soon as I did that, I looked at the resulting plot and realized that we are actually just comparing z23 / z22.

Here is why: the plot as it stands is a comparison of the log2FC values of the following two contrasts: z23/uninfected and z22/uninfected; stated differently, this is (z23/uninf)/(z22/uninf) which of course cancels out to just z23/z22.

Therefore it is much more parsimonious to just use the values from z23/z22. I swear I have gone through this exact exercise on so so many occasions in the past it is terrible.

16.4 ggsignificance of the immune modulators

wanted_genes <- c("IFI27", "RSAD2", "CCL8", "IFI44L", "OASL",
                  "USP18", "IDO1", "IDO2", "KYNU", "AHR", "IL4I1",
                  "SOD2", "NOTCH1", "DLL1", "DLL4", "HES1", "HEY1")
modulator_plot <- ggsignif_paired_genes(
  hs_macr, conditions = c("inf_z23", "inf_z22"), genes = wanted_genes,
  excel = "excel/paired_significance_modulators_nodrug.xlsx")

## Removing 9725 low-count genes (11756 remaining).

## transform_counts: Found 2375 values less than 0.

## Warning in transform_counts(count_table, method = transform, ...): NaNs
## produced

## Setting 34274 entries to zero.

## A subset operation de-synced the colors.

## Using Row.names, ensembl_gene_id, ensembl_transcript_id, description, gene_biotype, cds_length, chromosome_name, strand, hgnc_symbol, transcript as id variables
## Using Row.names, ensembl_gene_id, ensembl_transcript_id, description, gene_biotype, cds_length, chromosome_name, strand, hgnc_symbol, transcript as id variables

## Saving to excel/paired_significance_modulators_nodrug.xlsx

pp(file = "images/immune_response_inf_z23_vs_inf_z22_bwplot.pdf", width = 18, height = 9)
modulator_plot$plot

## Warning in wilcox.test.default(c(3.57602516488709, 3.45461601880412,
## 6.92662630290133, : cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(3.14968079679332, 0.467229645986823,
## 5.67587410724857, : cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0, 0, 0.669266872348993, 0, 0, 0, 0, 0, :
## cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(2.43931810389596, 2.30962138477456,
## 2.92880310623182, : cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0.845162885634284, 0.547086425216448,
## 0.619305136636964, : cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
## c(0.983845819547358, : cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
## c(1.74092551199831, : cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
## c(0.459773405747059, : cannot compute exact p-value with ties

dev.off()

## png 
##   2

modulator_plot$plot

## Warning in wilcox.test.default(c(3.57602516488709, 3.45461601880412,
## 6.92662630290133, : cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(3.14968079679332, 0.467229645986823,
## 5.67587410724857, : cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0, 0, 0.669266872348993, 0, 0, 0, 0, 0, :
## cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(2.43931810389596, 2.30962138477456,
## 2.92880310623182, : cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0.845162885634284, 0.547086425216448,
## 0.619305136636964, : cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
## c(0.983845819547358, : cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
## c(1.74092551199831, : cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
## c(0.459773405747059, : cannot compute exact p-value with ties

modulator_plot <- ggsignif_paired_genes(
  hs_macr, conditions = c("inf_sb_z23", "inf_sb_z22"), genes = wanted_genes,
  excel = "excel/paired_significance_modulators_drug.xlsx")

## Removing 9725 low-count genes (11756 remaining).

## transform_counts: Found 2375 values less than 0.

## Warning in transform_counts(count_table, method = transform, ...): NaNs
## produced

## Setting 34274 entries to zero.

## A subset operation de-synced the colors.

## Using Row.names, ensembl_gene_id, ensembl_transcript_id, description, gene_biotype, cds_length, chromosome_name, strand, hgnc_symbol, transcript as id variables
## Using Row.names, ensembl_gene_id, ensembl_transcript_id, description, gene_biotype, cds_length, chromosome_name, strand, hgnc_symbol, transcript as id variables

## Saving to excel/paired_significance_modulators_drug.xlsx

pp(file = "images/immune_response_drug_bwplot.pdf", width = 18, height = 9)
modulator_plot$plot

## Warning in wilcox.test.default(c(2.05034928211145, 1.04459288125667,
## 6.20521657695179, : cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(4.42759260675317, 3.82361053702975,
## 6.84290658530683, : cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0, 0, 2.11732007162309, 2.2425877878219, :
## cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), c(0, :
## cannot compute exact p-value with ties

## Warning: Computation failed in `stat_signif()`.
## Caused by error in `if (p_value < .Machine$double.eps) ...`:
## ! missing value where TRUE/FALSE needed

dev.off()

## png 
##   2

modulator_plot$plot

## Warning in wilcox.test.default(c(2.05034928211145, 1.04459288125667,
## 6.20521657695179, : cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(4.42759260675317, 3.82361053702975,
## 6.84290658530683, : cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0, 0, 2.11732007162309, 2.2425877878219, :
## cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), c(0, :
## cannot compute exact p-value with ties

## Warning: Computation failed in `stat_signif()`.
## Caused by error in `if (p_value < .Machine$double.eps) ...`:
## ! missing value where TRUE/FALSE needed

16.5 ggsignificance of the transporters

## First line is without drug
wanted_genes <- c("ABCB1", "ABCG4", "ABCB5", "AQP2", "AQP3",
                  ## with drug
                  "AQP8")
transporter_plot <- ggsignif_paired_genes(
  hs_macr, conditions = c("inf_z23", "inf_z22"), genes = wanted_genes,
  excel = "excel/transporter_response_nodrug.xlsx")

## Removing 9725 low-count genes (11756 remaining).

## transform_counts: Found 2375 values less than 0.

## Warning in transform_counts(count_table, method = transform, ...): NaNs
## produced

## Setting 34274 entries to zero.

## A subset operation de-synced the colors.

## Using Row.names, ensembl_gene_id, ensembl_transcript_id, description, gene_biotype, cds_length, chromosome_name, strand, hgnc_symbol, transcript as id variables
## Using Row.names, ensembl_gene_id, ensembl_transcript_id, description, gene_biotype, cds_length, chromosome_name, strand, hgnc_symbol, transcript as id variables

## Saving to excel/transporter_response_nodrug.xlsx

transporter_plot$plot

## Warning in wilcox.test.default(c(0.227265042226018, 0, 0, 0, 0, 0,
## 0.342524524890098, : cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0.00392158039284905, 0, 0, 0, 0, 0, 0, :
## cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
## c(0.158465337978994, : cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
## c(0.608106196156396, : cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0.44643072634979, 0, 0, 0, 0.0765472454757495,
## : cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
## c(0.19048433208064, : cannot compute exact p-value with ties

pp(file = "images/transporter_response_nodrug_bwplot.pdf", width = 18, height = 9)
transporter_plot$plot

## Warning in wilcox.test.default(c(0.227265042226018, 0, 0, 0, 0, 0,
## 0.342524524890098, : cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0.00392158039284905, 0, 0, 0, 0, 0, 0, :
## cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
## c(0.158465337978994, : cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
## c(0.608106196156396, : cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0.44643072634979, 0, 0, 0, 0.0765472454757495,
## : cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
## c(0.19048433208064, : cannot compute exact p-value with ties

dev.off()

## png 
##   2

transporter_drug_plot <- ggsignif_paired_genes(
  hs_macr, conditions = c("inf_sb_z23", "inf_sb_z22"), genes = wanted_genes,
  excel = "excel/transporter_response_drug.xlsx")

## Removing 9725 low-count genes (11756 remaining).

## transform_counts: Found 2375 values less than 0.

## Warning in transform_counts(count_table, method = transform, ...): NaNs
## produced

## Setting 34274 entries to zero.

## A subset operation de-synced the colors.

## Using Row.names, ensembl_gene_id, ensembl_transcript_id, description, gene_biotype, cds_length, chromosome_name, strand, hgnc_symbol, transcript as id variables
## Using Row.names, ensembl_gene_id, ensembl_transcript_id, description, gene_biotype, cds_length, chromosome_name, strand, hgnc_symbol, transcript as id variables

## Saving to excel/transporter_response_drug.xlsx

transporter_drug_plot$plot

## Warning in wilcox.test.default(c(0.0177775884891672, 0, 0, 0, 0, 0, 0, 0, :
## cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(1.12902181883524, 0.993926766482496, 0, :
## cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0, 0.00377928281757491, 0, 0, 0, 0,
## 0.181518403937829, : cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), c(0, :
## cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0.524739683612934, 0.430373504939515,
## 0.211690651455837, : cannot compute exact p-value with ties

pp(file = "images/transporter_response_drug_bwplot.pdf", width = 18, height = 9)
transporter_drug_plot$plot

## Warning in wilcox.test.default(c(0.0177775884891672, 0, 0, 0, 0, 0, 0, 0, :
## cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(1.12902181883524, 0.993926766482496, 0, :
## cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0, 0.00377928281757491, 0, 0, 0, 0,
## 0.181518403937829, : cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), c(0, :
## cannot compute exact p-value with ties

## Warning in wilcox.test.default(c(0.524739683612934, 0.430373504939515,
## 0.211690651455837, : cannot compute exact p-value with ties

dev.off()

## png 
##   2

pander::pander(sessionInfo())

## Warning: Your system is mis-configured: '/etc/localtime' is not a symlink

## Warning: It is strongly recommended to set envionment variable TZ to
## 'America/New_York' (or equivalent)

R version 4.5.0 (2025-04-11)

Platform: x86_64-pc-linux-gnu

locale: C

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

other attached packages: rWikiPathways(v.1.28.0), pathwayPCA(v.1.24.0), GSVAdata(v.1.44.0), edgeR(v.4.6.3), ruv(v.0.9.7.1), hpgltools(v.1.2), ggstatsplot(v.0.13.1), enrichplot(v.1.28.4), tidyr(v.1.3.1), tibble(v.3.3.0), UpSetR(v.1.4.0), Heatplus(v.3.16.0), glue(v.1.8.0), ggplot2(v.3.5.2) and ggbreak(v.0.1.6)

loaded via a namespace (and not attached): fs(v.1.6.6), matrixStats(v.1.5.0), bitops(v.1.0-9), blockmodeling(v.1.1.8), doParallel(v.1.0.17), httr(v.1.4.7), RColorBrewer(v.1.1-3), insight(v.1.4.2), numDeriv(v.2016.8-1.1), tools(v.4.5.0), backports(v.1.5.0), R6(v.2.6.1), statsExpressions(v.1.7.1), lazyeval(v.0.2.2), mgcv(v.1.9-3), GetoptLong(v.1.0.5), withr(v.3.0.2), prettyunits(v.1.2.0), gridExtra(v.2.3), preprocessCore(v.1.70.0), cli(v.3.6.5), Biobase(v.2.68.0), logging(v.0.10-108), labeling(v.0.4.3), EBSeq(v.2.6.0), sass(v.0.4.10), mvtnorm(v.1.3-3), genefilter(v.1.90.0), yulab.utils(v.0.2.1), DOSE(v.4.2.0), R.utils(v.2.13.0), dichromat(v.2.0-0.1), limma(v.3.64.3), RSQLite(v.2.4.3), shape(v.1.4.6.1), generics(v.0.1.4), gridGraphics(v.0.5-1), gtools(v.3.9.5), zip(v.2.3.3), dplyr(v.1.1.4), GO.db(v.3.21.0), lars(v.1.3), Matrix(v.1.7-3), S4Vectors(v.0.46.0), abind(v.1.4-8), R.methodsS3(v.1.8.2), lifecycle(v.1.0.4), yaml(v.2.3.10), SummarizedExperiment(v.1.38.1), BiocFileCache(v.2.16.1), gplots(v.3.2.0), qvalue(v.2.40.0), SparseArray(v.1.8.1), paletteer(v.1.6.0), grid(v.4.5.0), blob(v.1.2.4), promises(v.1.3.3), crayon(v.1.5.3), ggtangle(v.0.0.7), lattice(v.0.22-7), cowplot(v.1.2.0), annotate(v.1.86.1), KEGGREST(v.1.48.1), ComplexHeatmap(v.2.24.1), zeallot(v.0.2.0), pillar(v.1.11.0), knitr(v.1.50), varhandle(v.2.0.6), fgsea(v.1.34.2), GenomicRanges(v.1.60.0), rjson(v.0.2.23), boot(v.1.3-31), corpcor(v.1.6.10), codetools(v.0.2-20), fastmatch(v.1.1-6), ggfun(v.0.2.0), Vennerable(v.3.1.0.9000), data.table(v.1.17.8), vctrs(v.0.6.5), png(v.0.1-8), treeio(v.1.32.0), Rdpack(v.2.6.4), testthat(v.3.2.3), gtable(v.0.3.6), rematch2(v.2.1.2), datawizard(v.1.2.0), cachem(v.1.1.0), openxlsx(v.4.2.8), xfun(v.0.53), rbibutils(v.2.3), S4Arrays(v.1.8.1), mime(v.0.13), RcppEigen(v.0.3.4.0.2), correlation(v.0.8.8), ConsensusClusterPlus(v.1.72.0), reformulas(v.0.4.1), survival(v.3.8-3), iterators(v.1.0.14), NOISeq(v.2.52.0), statmod(v.1.5.0), nlme(v.3.1-168), pbkrtest(v.0.5.5), ggtree(v.3.17.1.001), bit64(v.4.6.0-1), filelock(v.1.0.3), progress(v.1.2.3), EnvStats(v.3.1.0), GenomeInfoDb(v.1.44.2), rprojroot(v.2.1.0), bslib(v.0.9.0), KernSmooth(v.2.23-26), colorspace(v.2.1-1), BiocGenerics(v.0.54.0), DBI(v.1.2.3), mnormt(v.2.1.1), DESeq2(v.1.48.1), tidyselect(v.1.2.1), curl(v.7.0.0), bit(v.4.6.0), compiler(v.4.5.0), httr2(v.1.2.1), graph(v.1.86.0), xml2(v.1.4.0), ggdendro(v.0.2.0), desc(v.1.4.3), DelayedArray(v.0.34.1), plotly(v.4.11.0), bayestestR(v.0.17.0), scales(v.1.4.0), caTools(v.1.18.3), psych(v.2.5.6), remaCor(v.0.0.20), RBGL(v.1.84.0), rappdirs(v.0.3.3), stringr(v.1.5.1), digest(v.0.6.37), minqa(v.1.2.8), variancePartition(v.1.38.1), rmarkdown(v.2.29), aod(v.1.3.3), XVector(v.0.48.0), RhpcBLASctl(v.0.23-42), htmltools(v.0.5.8.1), pkgconfig(v.2.0.3), lme4(v.1.1-37), MatrixGenerics(v.1.20.0), dbplyr(v.2.5.0), fastmap(v.1.2.0), GlobalOptions(v.0.1.2), rlang(v.1.1.6), htmlwidgets(v.1.6.4), UCSC.utils(v.1.4.0), shiny(v.1.11.1), farver(v.2.1.2), jquerylib(v.0.1.4), jsonlite(v.2.0.0), BiocParallel(v.1.42.1), GOSemSim(v.2.34.0), R.oo(v.1.27.1), RCurl(v.1.98-1.17), magrittr(v.2.0.3), GenomeInfoDbData(v.1.2.14), ggplotify(v.0.1.2), patchwork(v.1.3.2), parameters(v.0.28.1), Rcpp(v.1.1.0), ape(v.5.8-1), stringi(v.1.8.7), brio(v.1.1.5), MASS(v.7.3-65), DEGreport(v.1.44.0), plyr(v.1.8.9), parallel(v.4.5.0), ggrepel(v.0.9.6), Biostrings(v.2.76.0), splines(v.4.5.0), pander(v.0.6.6), circlize(v.0.4.16), hms(v.1.1.3), locfit(v.1.5-9.12), igraph(v.2.1.4), ggsignif(v.0.6.4), effectsize(v.1.0.1), biomaRt(v.2.64.0), reshape2(v.1.4.4), stats4(v.4.5.0), pkgload(v.1.4.0), rstantools(v.2.5.0), XML(v.3.99-0.19), evaluate(v.1.0.4), RcppParallel(v.5.1.11-1), BiocManager(v.1.30.26), nloptr(v.2.2.1), foreach(v.1.5.2), httpuv(v.1.6.16), purrr(v.1.1.0), reshape(v.0.8.10), clue(v.0.3-66), broom(v.1.0.9), xtable(v.1.8-4), fANCOVA(v.0.6-1), tidytree(v.0.4.6), later(v.1.4.3), viridisLite(v.0.4.2), lmerTest(v.3.1-3), aplot(v.0.2.8), memoise(v.2.0.1), AnnotationDbi(v.1.70.0), IRanges(v.2.42.0), cluster(v.2.1.8.1), sva(v.3.56.0) and GSEABase(v.1.70.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 26f5269f5c71e395abd2c55082004e9884b2170f

## This is hpgltools commit: Wed Sep 10 11:35:36 2025 -0400: 26f5269f5c71e395abd2c55082004e9884b2170f

tmp <- saveme(filename = savefile)

## The savefile is: /lab/singularity/tmrc2_macrophage_deb/202509101136_outputs/savefiles/03differential_expression.rda.xz

## The file does not yet exist.

## The save string is: con <- pipe(paste0('pxz > /lab/singularity/tmrc2_macrophage_deb/202509101136_outputs/savefiles/03differential_expression.rda.xz'), 'wb'); save(list = ls(all.names = TRUE, envir = globalenv()),
##      envir = globalenv(), file = con, compress = FALSE); close(con)

## Error in save(list = ls(all.names = TRUE, envir = globalenv()), envir = globalenv(), : ignoring SIGPIPE signal

tmp <- loadme(filename = savefile)

devtools::load_all(‘~/hpgltools’)