L. panamensis 202209: TMRC2 Macrophage Differential Expression.

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

1.1 Human data

1.1.1 Primary queries

There is a series of initial questions which make some intuitive sense to me, but these do not necessarily match the set of questions which are intuitive for Olga. 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).

extra <- "z23drugnodrug_vs_z22drugnodrug = (infsbz23 - infz23) - (infsbz22 - infz22)"
new_conditions <- paste0(pData(hs_macrophage)[["macrophagetreatment"]], "_",
                         pData(hs_macrophage)[["macrophagezymodeme"]])
hs_macrophage <- set_expt_conditions(hs_macrophage, fact = new_conditions)

hs_macrophage_de <- all_pairwise(hs_macrophage, model_batch="svaseq", filter=TRUE,
                                 extra_contrasts = extra)

## This DE analysis will perform all pairwise comparisons among:

## 
##   inf_sb_z2.2   inf_sb_z2.3      inf_z2.2      inf_z2.3    uninf_none 
##             6             6             6             6             2 
## uninf_sb_none 
##             2

## This analysis will include surrogate estimates from: svaseq.

## This will pre-filter the input data using normalize_expt's: TRUE argument.

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

## Setting 757 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

## Used reverse contrast for deseq.

## Used reverse contrast for edger.

## Used reverse contrast for deseq.

## Used reverse contrast for basic.

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

tmrc2_human_keepers <- list(
    "z23nosb_vs_uninf" = c("infz23", "uninfnone"),
    "z22nosb_vs_uninf" = c("infz22", "uninfnone"),
    "z23nosb_vs_z22nosb" = c("infz23", "infz22"),
    "z23sb_vs_z22sb" = c("infsbz23", "infsbz22"),
    "z23sb_vs_z23nosb" = c("infsbz23", "infz23"),
    "z22sb_vs_z22nosb" = c("infsbz22", "infz22"),
    "z23sb_vs_sb" = c("infz23", "uninfsbnone"),
    "z22sb_vs_sb" = c("infz22", "uninfsbnone"),
    "z23sb_vs_uninf" = c("infsbz23", "uninfnone"),
    "z22sb_vs_uninf" = c("infsbz22", "uninfnone"),
    "sb_vs_uninf" = c("uninfsbnone", "uninfnone"),
    "extra" = c("z23drugnodrug", "z22drugnodrug"))

1.1.1.2 Write contrast results

Now let us write out the xlsx file containing the above contrasts.

hs_macrophage_table <- combine_de_tables(
    hs_macrophage_de,
    keepers = tmrc2_human_keepers,
    excel=glue::glue("excel/macrophage_human_table-v{ver}.xlsx"))

## Deleting the file excel/macrophage_human_table-v202209.xlsx before writing the tables.

## Warning in combine_single_de_table(li = limma, ed = edger, eb = ebseq, de =
## deseq, : The deseq table seems to be missing.

## Warning in combine_single_de_table(li = limma, ed = edger, eb = ebseq, de =
## deseq, : The basic table seems to be missing.

## Unable to find the table in the set of possible tables.

## The possible tables are: infsbz23_vs_infsbz22, infz22_vs_infsbz22, infz23_vs_infsbz22, uninfnone_vs_infsbz22, uninfsbnone_vs_infsbz22, infz22_vs_infsbz23, infz23_vs_infsbz23, uninfnone_vs_infsbz23, uninfsbnone_vs_infsbz23, infz23_vs_infz22, uninfnone_vs_infz22, uninfsbnone_vs_infz22, uninfnone_vs_infz23, uninfsbnone_vs_infz23, uninfsbnone_vs_uninfnone

hs_macrophage_sig <- extract_significant_genes(
    hs_macrophage_table,
    excel=glue::glue("excel/macrophage_human_sig-v{ver}.xlsx"))

## Deleting the file excel/macrophage_human_sig-v202209.xlsx before writing the tables.

## Unable to find the table in the set of possible tables.

## The possible tables are: infsbz23_vs_infsbz22, infz22_vs_infsbz22, infz23_vs_infsbz22, uninfnone_vs_infsbz22, uninfsbnone_vs_infsbz22, infz22_vs_infsbz23, infz23_vs_infsbz23, uninfnone_vs_infsbz23, uninfsbnone_vs_infsbz23, infz23_vs_infz22, uninfnone_vs_infz22, uninfsbnone_vs_infz22, uninfnone_vs_infz23, uninfsbnone_vs_infz23, uninfsbnone_vs_uninfnone

## Unable to find the table in the set of possible tables.

## The possible tables are: infsbz23_vs_infsbz22, infz22_vs_infsbz22, infz23_vs_infsbz22, uninfnone_vs_infsbz22, uninfsbnone_vs_infsbz22, infz22_vs_infsbz23, infz23_vs_infsbz23, uninfnone_vs_infsbz23, uninfsbnone_vs_infsbz23, infz23_vs_infz22, uninfnone_vs_infz22, uninfsbnone_vs_infz22, uninfnone_vs_infz23, uninfsbnone_vs_infz23, uninfsbnone_vs_uninfnone

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

z23nosb_vs_uninf_volcano <- plot_volcano_de(
    table = hs_macrophage_table[["data"]][["z23nosb_vs_uninf"]],
    fc_col = "deseq_logfc", p_col = "deseq_adjp",
    shapes_by_state = FALSE, color_by = "fc",  label = 10, label_column = "hgncsymbol")

## Error in (function(x, i, exact) if (is.matrix(i)) as.matrix(x)[[i]] else .subset2(x, : object 'table_column' not found

plotly::ggplotly(z23nosb_vs_uninf_volcano$plot)

## Error in plotly::ggplotly(z23nosb_vs_uninf_volcano$plot): object 'z23nosb_vs_uninf_volcano' not found

z22nosb_vs_uninf_volcano <- plot_volcano_de(
    table = hs_macrophage_table[["data"]][["z22nosb_vs_uninf"]],
    fc_col = "deseq_logfc", p_col = "deseq_adjp",
    shapes_by_state = FALSE, color_by = "fc",  label = 10, label_column = "hgncsymbol")

## Error in (function(x, i, exact) if (is.matrix(i)) as.matrix(x)[[i]] else .subset2(x, : object 'table_column' not found

plotly::ggplotly(z22nosb_vs_uninf_volcano$plot)

## Error in plotly::ggplotly(z22nosb_vs_uninf_volcano$plot): object 'z22nosb_vs_uninf_volcano' not found

hs_macrophage_table[["plots"]][["z23nosb_vs_uninf"]][["deseq_vol_plots"]][["plot"]] +
  xlim(-10, 25) +
  ylim(0, 40)

## Warning: Removed 1 rows containing missing values (geom_point).

hs_macrophage_table[["plots"]][["z22nosb_vs_uninf"]][["deseq_vol_plots"]][["plot"]] +
  xlim(-10, 25) +
  ylim(0, 40)

hs_macrophage_table[["plots"]][["z23nosb_vs_z22nosb"]][["deseq_vol_plots"]][["plot"]] +
  xlim(-10, 10) +
  ylim(0, 60)

hs_macrophage_table[["plots"]][["z23sb_vs_z22sb"]][["deseq_vol_plots"]][["plot"]] +
  xlim(-10, 10) +
  ylim(0, 60)

shared <- Vennerable::Venn(list("drug" = rownames(hs_macrophage_sig[["deseq"]][["ups"]][["z23sb_vs_z22sb"]]),
                                "nodrug" = rownames(hs_macrophage_sig[["deseq"]][["ups"]][["z23nosb_vs_z22nosb"]])))
pp(file="images/drug_nodrug_venn_up.png")
Vennerable::plot(shared)
dev.off()

## png 
##   2

Vennerable::plot(shared)

shared <- Vennerable::Venn(list("drug" = rownames(hs_macrophage_sig[["deseq"]][["downs"]][["z23sb_vs_z22sb"]]),
                                "nodrug" = rownames(hs_macrophage_sig[["deseq"]][["downs"]][["z23nosb_vs_z22nosb"]])))
pp(file="images/drug_nodrug_venn_down.png")
Vennerable::plot(shared)
dev.off()

## png 
##   2

hs_macrophage_table[["plots"]][["z23sb_vs_z23nosb"]][["deseq_vol_plots"]][["plot"]] +
  xlim(-8, 8) +
  ylim(0, 210)

hs_macrophage_table[["plots"]][["z22sb_vs_z22nosb"]][["deseq_vol_plots"]][["plot"]] +
  xlim(-8, 8) +
  ylim(0, 210)

hs_macrophage_table[["plots"]][["z23sb_vs_sb"]][["deseq_vol_plots"]][["plot"]] +
  xlim(-10, 28) +
  ylim(0, 140)

hs_macrophage_table[["plots"]][["z22sb_vs_sb"]][["deseq_vol_plots"]][["plot"]] +
  xlim(-10, 28) +
  ylim(0, 140)

shared <- Vennerable::Venn(list("z23" = rownames(hs_macrophage_sig[["deseq"]][["ups"]][["z23sb_vs_sb"]]),
                                "z22" = rownames(hs_macrophage_sig[["deseq"]][["ups"]][["z22sb_vs_sb"]])))
pp(file="images/z23_z22_drug_venn_up.png")
Vennerable::plot(shared)
dev.off()

## png 
##   2

Vennerable::plot(shared)

shared <- Vennerable::Venn(list("z23" = rownames(hs_macrophage_sig[["deseq"]][["downs"]][["z23sb_vs_sb"]]),
                                "z22" = rownames(hs_macrophage_sig[["deseq"]][["downs"]][["z22sb_vs_sb"]])))
pp(file="images/z23_z22_drug_venn_down.png")
Vennerable::plot(shared)
dev.off()

## png 
##   2

Vennerable::plot(shared)

1.2 Parasite

lp_macrophage_de <- all_pairwise(lp_macrophage,
                                 model_batch="svaseq", filter=TRUE)

## This DE analysis will perform all pairwise comparisons among:

## 
## z2.2 z2.3 
##    5    6

## This analysis will include surrogate estimates from: svaseq.

## This will pre-filter the input data using normalize_expt's: TRUE argument.

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

## Setting 45 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

tmrc2_parasite_keepers <- list(
    "z23nosb_vs_z22nosb" = c("z23", "z22"))
lp_macrophage_table <- combine_de_tables(
  lp_macrophage_de, keepers = tmrc2_parasite_keepers,
  excel=glue::glue("excel/macrophage_parasite_infection_de-v{ver}.xlsx"))

## Deleting the file excel/macrophage_parasite_infection_de-v202209.xlsx before writing the tables.

lp_macrophage_sig <- extract_significant_genes(
    lp_macrophage_table,
    excel=glue::glue("excel/macrophage_parasite_sig-v{ver}.xlsx"))

## Deleting the file excel/macrophage_parasite_sig-v202209.xlsx before writing the tables.

pp(file="images/lp_macrophage_z23_z22.png", image=lp_macrophage_table[["plots"]][["z23nosb_vs_z22nosb"]][["deseq_vol_plots"]][["plot"]])

## Warning in pp(file = "images/lp_macrophage_z23_z22.png", image =
## lp_macrophage_table[["plots"]][["z23nosb_vs_z22nosb"]][["deseq_vol_plots"]]
## [["plot"]]): There is no device to shut down.

up_genes <- lp_macrophage_sig[["deseq"]][["ups"]][[1]]
dim(up_genes)

## [1] 40 58

down_genes <- lp_macrophage_sig[["deseq"]][["downs"]][[1]]
dim(down_genes)

## [1] 75 58

2 Over representation searches

all_gp <- all_gprofiler(hs_macrophage_sig)

## Performing gProfiler GO search of 670 genes against hsapiens.

## Performing gProfiler KEGG search of 670 genes against hsapiens.

## Performing gProfiler REAC search of 670 genes against hsapiens.

## Performing gProfiler MI search of 670 genes against hsapiens.

## Performing gProfiler TF search of 670 genes against hsapiens.

## Performing gProfiler CORUM search of 670 genes against hsapiens.

## Performing gProfiler HP search of 670 genes against hsapiens.

## Performing gProfiler GO search of 386 genes against hsapiens.

## Performing gProfiler KEGG search of 386 genes against hsapiens.

## Performing gProfiler REAC search of 386 genes against hsapiens.

## Performing gProfiler MI search of 386 genes against hsapiens.

## Performing gProfiler TF search of 386 genes against hsapiens.

## Performing gProfiler CORUM search of 386 genes against hsapiens.

## Performing gProfiler HP search of 386 genes against hsapiens.

## Performing gProfiler GO search of 435 genes against hsapiens.

## Performing gProfiler KEGG search of 435 genes against hsapiens.

## Performing gProfiler REAC search of 435 genes against hsapiens.

## Performing gProfiler MI search of 435 genes against hsapiens.

## Performing gProfiler TF search of 435 genes against hsapiens.

## Performing gProfiler CORUM search of 435 genes against hsapiens.

## Performing gProfiler HP search of 435 genes against hsapiens.

## Performing gProfiler GO search of 7 genes against hsapiens.

## Performing gProfiler KEGG search of 7 genes against hsapiens.

## Performing gProfiler REAC search of 7 genes against hsapiens.

## Performing gProfiler MI search of 7 genes against hsapiens.

## Performing gProfiler TF search of 7 genes against hsapiens.

## Performing gProfiler CORUM search of 7 genes against hsapiens.

## Performing gProfiler HP search of 7 genes against hsapiens.

## Performing gProfiler GO search of 529 genes against hsapiens.

## Performing gProfiler KEGG search of 529 genes against hsapiens.

## Performing gProfiler REAC search of 529 genes against hsapiens.

## Performing gProfiler MI search of 529 genes against hsapiens.

## Performing gProfiler TF search of 529 genes against hsapiens.

## Performing gProfiler CORUM search of 529 genes against hsapiens.

## Performing gProfiler HP search of 529 genes against hsapiens.

## Performing gProfiler GO search of 730 genes against hsapiens.

## Performing gProfiler KEGG search of 730 genes against hsapiens.

## Performing gProfiler REAC search of 730 genes against hsapiens.

## Performing gProfiler MI search of 730 genes against hsapiens.

## Performing gProfiler TF search of 730 genes against hsapiens.

## Performing gProfiler CORUM search of 730 genes against hsapiens.

## Performing gProfiler HP search of 730 genes against hsapiens.

## Performing gProfiler GO search of 484 genes against hsapiens.

## Performing gProfiler KEGG search of 484 genes against hsapiens.

## Performing gProfiler REAC search of 484 genes against hsapiens.

## Performing gProfiler MI search of 484 genes against hsapiens.

## Performing gProfiler TF search of 484 genes against hsapiens.

## Performing gProfiler CORUM search of 484 genes against hsapiens.

## Performing gProfiler HP search of 484 genes against hsapiens.

## Performing gProfiler GO search of 422 genes against hsapiens.

## Performing gProfiler KEGG search of 422 genes against hsapiens.

## Performing gProfiler REAC search of 422 genes against hsapiens.

## Performing gProfiler MI search of 422 genes against hsapiens.

## Performing gProfiler TF search of 422 genes against hsapiens.

## Performing gProfiler CORUM search of 422 genes against hsapiens.

## Performing gProfiler HP search of 422 genes against hsapiens.

## Performing gProfiler GO search of 555 genes against hsapiens.

## Performing gProfiler KEGG search of 555 genes against hsapiens.

## Performing gProfiler REAC search of 555 genes against hsapiens.

## Performing gProfiler MI search of 555 genes against hsapiens.

## Performing gProfiler TF search of 555 genes against hsapiens.

## Performing gProfiler CORUM search of 555 genes against hsapiens.

## Performing gProfiler HP search of 555 genes against hsapiens.

## Performing gProfiler GO search of 780 genes against hsapiens.

## Performing gProfiler KEGG search of 780 genes against hsapiens.

## Performing gProfiler REAC search of 780 genes against hsapiens.

## Performing gProfiler MI search of 780 genes against hsapiens.

## Performing gProfiler TF search of 780 genes against hsapiens.

## Performing gProfiler CORUM search of 780 genes against hsapiens.

## Performing gProfiler HP search of 780 genes against hsapiens.

## Performing gProfiler GO search of 454 genes against hsapiens.

## Performing gProfiler KEGG search of 454 genes against hsapiens.

## Performing gProfiler REAC search of 454 genes against hsapiens.

## Performing gProfiler MI search of 454 genes against hsapiens.

## Performing gProfiler TF search of 454 genes against hsapiens.

## Performing gProfiler CORUM search of 454 genes against hsapiens.

## Performing gProfiler HP search of 454 genes against hsapiens.

## Performing gProfiler GO search of 986 genes against hsapiens.

## Performing gProfiler KEGG search of 986 genes against hsapiens.

## Performing gProfiler REAC search of 986 genes against hsapiens.

## Performing gProfiler MI search of 986 genes against hsapiens.

## Performing gProfiler TF search of 986 genes against hsapiens.

## Performing gProfiler CORUM search of 986 genes against hsapiens.

## Performing gProfiler HP search of 986 genes against hsapiens.

## Performing gProfiler GO search of 1081 genes against hsapiens.

## Performing gProfiler KEGG search of 1081 genes against hsapiens.

## Performing gProfiler REAC search of 1081 genes against hsapiens.

## Performing gProfiler MI search of 1081 genes against hsapiens.

## Performing gProfiler TF search of 1081 genes against hsapiens.

## Performing gProfiler CORUM search of 1081 genes against hsapiens.

## Performing gProfiler HP search of 1081 genes against hsapiens.

## Performing gProfiler GO search of 701 genes against hsapiens.

## Performing gProfiler KEGG search of 701 genes against hsapiens.

## Performing gProfiler REAC search of 701 genes against hsapiens.

## Performing gProfiler MI search of 701 genes against hsapiens.

## Performing gProfiler TF search of 701 genes against hsapiens.

## Performing gProfiler CORUM search of 701 genes against hsapiens.

## Performing gProfiler HP search of 701 genes against hsapiens.

## Performing gProfiler GO search of 1025 genes against hsapiens.

## Performing gProfiler KEGG search of 1025 genes against hsapiens.

## Performing gProfiler REAC search of 1025 genes against hsapiens.

## Performing gProfiler MI search of 1025 genes against hsapiens.

## Performing gProfiler TF search of 1025 genes against hsapiens.

## Performing gProfiler CORUM search of 1025 genes against hsapiens.

## Performing gProfiler HP search of 1025 genes against hsapiens.

## Performing gProfiler GO search of 427 genes against hsapiens.

## Performing gProfiler KEGG search of 427 genes against hsapiens.

## Performing gProfiler REAC search of 427 genes against hsapiens.

## Performing gProfiler MI search of 427 genes against hsapiens.

## Performing gProfiler TF search of 427 genes against hsapiens.

## Performing gProfiler CORUM search of 427 genes against hsapiens.

## Performing gProfiler HP search of 427 genes against hsapiens.

## Performing gProfiler GO search of 910 genes against hsapiens.

## Performing gProfiler KEGG search of 910 genes against hsapiens.

## Performing gProfiler REAC search of 910 genes against hsapiens.

## Performing gProfiler MI search of 910 genes against hsapiens.

## Performing gProfiler TF search of 910 genes against hsapiens.

## Performing gProfiler CORUM search of 910 genes against hsapiens.

## Performing gProfiler HP search of 910 genes against hsapiens.

## Performing gProfiler GO search of 983 genes against hsapiens.

## Performing gProfiler KEGG search of 983 genes against hsapiens.

## Performing gProfiler REAC search of 983 genes against hsapiens.

## Performing gProfiler MI search of 983 genes against hsapiens.

## Performing gProfiler TF search of 983 genes against hsapiens.

## Performing gProfiler CORUM search of 983 genes against hsapiens.

## Performing gProfiler HP search of 983 genes against hsapiens.

## Performing gProfiler GO search of 572 genes against hsapiens.

## Performing gProfiler KEGG search of 572 genes against hsapiens.

## Performing gProfiler REAC search of 572 genes against hsapiens.

## Performing gProfiler MI search of 572 genes against hsapiens.

## Performing gProfiler TF search of 572 genes against hsapiens.

## Performing gProfiler CORUM search of 572 genes against hsapiens.

## Performing gProfiler HP search of 572 genes against hsapiens.

## Performing gProfiler GO search of 644 genes against hsapiens.

## Performing gProfiler KEGG search of 644 genes against hsapiens.

## Performing gProfiler REAC search of 644 genes against hsapiens.

## Performing gProfiler MI search of 644 genes against hsapiens.

## Performing gProfiler TF search of 644 genes against hsapiens.

## Performing gProfiler CORUM search of 644 genes against hsapiens.

## Performing gProfiler HP search of 644 genes against hsapiens.

## Performing gProfiler GO search of 362 genes against hsapiens.

## Performing gProfiler KEGG search of 362 genes against hsapiens.

## Performing gProfiler REAC search of 362 genes against hsapiens.

## Performing gProfiler MI search of 362 genes against hsapiens.

## Performing gProfiler TF search of 362 genes against hsapiens.

## Performing gProfiler CORUM search of 362 genes against hsapiens.

## Performing gProfiler HP search of 362 genes against hsapiens.

## Performing gProfiler GO search of 503 genes against hsapiens.

## Performing gProfiler KEGG search of 503 genes against hsapiens.

## Performing gProfiler REAC search of 503 genes against hsapiens.

## Performing gProfiler MI search of 503 genes against hsapiens.

## Performing gProfiler TF search of 503 genes against hsapiens.

## Performing gProfiler CORUM search of 503 genes against hsapiens.

## Performing gProfiler HP search of 503 genes against hsapiens.

pp(file="images/z23_uninf_reactome_up.png", image=all_gp[["z23nosb_vs_uninf_up"]][["pvalue_plots"]][["reactome_plot_over"]], height=12, width=9)

## Error in all_gp[["z23nosb_vs_uninf_up"]][["pvalue_plots"]][["reactome_plot_over"]]: subscript out of bounds

all_gp[["z23nosb_vs_uninf_up"]][["pvalue_plots"]][["kegg_plot_over"]]

## Error in all_gp[["z23nosb_vs_uninf_up"]][["pvalue_plots"]][["kegg_plot_over"]]: subscript out of bounds

all_gp[["z23nosb_vs_uninf_up"]][["pvalue_plots"]][["mfp_plot_over"]]

## Error in all_gp[["z23nosb_vs_uninf_up"]][["pvalue_plots"]][["mfp_plot_over"]]: subscript out of bounds

all_gp[["z23nosb_vs_uninf_up"]][["pvalue_plots"]][["tf_plot_over"]]

## Error in all_gp[["z23nosb_vs_uninf_up"]][["pvalue_plots"]][["tf_plot_over"]]: subscript out of bounds

pp(file="images/z22_uninf_reactome_up.png", image=all_gp[["z22nosb_vs_uninf_up"]][["pvalue_plots"]][["reactome_plot_over"]], height=12, width=9)

## Error in all_gp[["z22nosb_vs_uninf_up"]][["pvalue_plots"]][["reactome_plot_over"]]: subscript out of bounds

all_gp[["z22nosb_vs_uninf_up"]][["pvalue_plots"]][["kegg_plot_over"]]

## Error in all_gp[["z22nosb_vs_uninf_up"]][["pvalue_plots"]][["kegg_plot_over"]]: subscript out of bounds

all_gp[["z22nosb_vs_uninf_up"]][["pvalue_plots"]][["mfp_plot_over"]]

## Error in all_gp[["z22nosb_vs_uninf_up"]][["pvalue_plots"]][["mfp_plot_over"]]: subscript out of bounds

all_gp[["z22nosb_vs_uninf_up"]][["pvalue_plots"]][["tf_plot_over"]]

## Error in all_gp[["z22nosb_vs_uninf_up"]][["pvalue_plots"]][["tf_plot_over"]]: subscript out of bounds

all_gp[["z23sb_vs_z22sb_up"]][["pvalue_plots"]][["reactome_plot_over"]]

## Error in all_gp[["z23sb_vs_z22sb_up"]][["pvalue_plots"]][["reactome_plot_over"]]: subscript out of bounds

all_gp[["z23sb_vs_z22sb_up"]][["pvalue_plots"]][["kegg_plot_over"]]

## Error in all_gp[["z23sb_vs_z22sb_up"]][["pvalue_plots"]][["kegg_plot_over"]]: subscript out of bounds

all_gp[["z23sb_vs_z22sb_up"]][["pvalue_plots"]][["mfp_plot_over"]]

## Error in all_gp[["z23sb_vs_z22sb_up"]][["pvalue_plots"]][["mfp_plot_over"]]: subscript out of bounds

all_gp[["z23sb_vs_z22sb_up"]][["pvalue_plots"]][["tf_plot_over"]]

## Error in all_gp[["z23sb_vs_z22sb_up"]][["pvalue_plots"]][["tf_plot_over"]]: subscript out of bounds

all_gp[["z23sb_vs_z22sb_down"]][["pvalue_plots"]][["reactome_plot_over"]]

## Error in all_gp[["z23sb_vs_z22sb_down"]][["pvalue_plots"]][["reactome_plot_over"]]: subscript out of bounds

all_gp[["z23sb_vs_z23nosb_up"]][["pvalue_plots"]][["reactome_plot_over"]]

## Error in all_gp[["z23sb_vs_z23nosb_up"]][["pvalue_plots"]][["reactome_plot_over"]]: subscript out of bounds

all_gp[["z23sb_vs_z23nosb_up"]][["pvalue_plots"]][["kegg_plot_over"]]

## Error in all_gp[["z23sb_vs_z23nosb_up"]][["pvalue_plots"]][["kegg_plot_over"]]: subscript out of bounds

all_gp[["z23sb_vs_z23nosb_up"]][["pvalue_plots"]][["mfp_plot_over"]]

## Error in all_gp[["z23sb_vs_z23nosb_up"]][["pvalue_plots"]][["mfp_plot_over"]]: subscript out of bounds

all_gp[["z23sb_vs_z23nosb_up"]][["pvalue_plots"]][["tf_plot_over"]]

## Error in all_gp[["z23sb_vs_z23nosb_up"]][["pvalue_plots"]][["tf_plot_over"]]: subscript out of bounds

all_gp[["z22sb_vs_z22nosb_up"]][["pvalue_plots"]][["reactome_plot_over"]]

## Error in all_gp[["z22sb_vs_z22nosb_up"]][["pvalue_plots"]][["reactome_plot_over"]]: subscript out of bounds

all_gp[["z22sb_vs_z22nosb_up"]][["pvalue_plots"]][["kegg_plot_over"]]

## Error in all_gp[["z22sb_vs_z22nosb_up"]][["pvalue_plots"]][["kegg_plot_over"]]: subscript out of bounds

all_gp[["z22sb_vs_z22nosb_up"]][["pvalue_plots"]][["mfp_plot_over"]]

## Error in all_gp[["z22sb_vs_z22nosb_up"]][["pvalue_plots"]][["mfp_plot_over"]]: subscript out of bounds

all_gp[["z22sb_vs_z22nosb_up"]][["pvalue_plots"]][["tf_plot_over"]]

## Error in all_gp[["z22sb_vs_z22nosb_up"]][["pvalue_plots"]][["tf_plot_over"]]: subscript out of bounds

up_goseq <- simple_goseq(up_genes, go_db=lp_go, length_db=lp_lengths)

## Found 11 go_db genes and 40 length_db genes out of 40.

## Testing that go categories are defined.

## Removing undefined categories.

## Gathering synonyms.

## Gathering category definitions.

## The score column is null, defaulting to score.

## Possible columns are:

## [1] "category"                 "over_represented_pvalue" 
## [3] "under_represented_pvalue" "numDEInCat"              
## [5] "numInCat"                 "term"                    
## [7] "ontology"                 "qvalue"

## The score column is null, defaulting to score.
## Possible columns are:

## [1] "category"                 "over_represented_pvalue" 
## [3] "under_represented_pvalue" "numDEInCat"              
## [5] "numInCat"                 "term"                    
## [7] "ontology"                 "qvalue"

## The score column is null, defaulting to score.
## Possible columns are:

## [1] "category"                 "over_represented_pvalue" 
## [3] "under_represented_pvalue" "numDEInCat"              
## [5] "numInCat"                 "term"                    
## [7] "ontology"                 "qvalue"

## View categories over represented in the 2.3 samples
up_goseq$pvalue_plots$bpp_plot_over

down_goseq <- simple_goseq(down_genes, go_db=lp_go, length_db=lp_lengths)

## Found 21 go_db genes and 75 length_db genes out of 75.

## Testing that go categories are defined.

## Removing undefined categories.

## Gathering synonyms.

## Gathering category definitions.

## The score column is null, defaulting to score.

## Possible columns are:

## [1] "category"                 "over_represented_pvalue" 
## [3] "under_represented_pvalue" "numDEInCat"              
## [5] "numInCat"                 "term"                    
## [7] "ontology"                 "qvalue"

## The score column is null, defaulting to score.
## Possible columns are:

## [1] "category"                 "over_represented_pvalue" 
## [3] "under_represented_pvalue" "numDEInCat"              
## [5] "numInCat"                 "term"                    
## [7] "ontology"                 "qvalue"

## The score column is null, defaulting to score.
## Possible columns are:

## [1] "category"                 "over_represented_pvalue" 
## [3] "under_represented_pvalue" "numDEInCat"              
## [5] "numInCat"                 "term"                    
## [7] "ontology"                 "qvalue"

## View categories over represented in the 2.2 samples
down_goseq$pvalue_plots$bpp_plot_over

3 GSVA

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

## subset_expt(): There were 28, now there are 26 samples.

## subset_expt(): There were 26, now there are 24 samples.

hs_gsva_c2 <- simple_gsva(hs_infected)

## Converting the rownames() of the expressionset to ENTREZID.

## 1622 ENSEMBL ID's didn't have a matching ENTEREZ ID. Dropping them now.

## Before conversion, the expressionset has 21481 entries.

## After conversion, the expressionset has 20013 entries.

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

## Converting the rownames() of the expressionset to ENTREZID.

## 1622 ENSEMBL ID's didn't have a matching ENTEREZ ID. Dropping them now.

## Before conversion, the expressionset has 21481 entries.

## After conversion, the expressionset has 20013 entries.

hs_gsva_c2_sig <- get_sig_gsva_categories(hs_gsva_c2, excel = "excel/hs_macrophage_gsva_c2_sig.xlsx")

## Starting limma pairwise comparison.

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

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Choosing the non-intercept containing model.

## Assuming this data is similar to a micro array and not performign voom.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

## The factor inf_sb_z2.2 has 6 rows.

## The factor inf_sb_z2.3 has 6 rows.

## The factor inf_z2.2 has 6 rows.

## The factor inf_z2.3 has 6 rows.

## Testing each factor against the others.

## Scoring inf_sb_z2.2 against everything else.

## Scoring inf_sb_z2.3 against everything else.

## Scoring inf_z2.2 against everything else.

## Scoring inf_z2.3 against everything else.

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

hs_gsva_c2_sig$raw_plot

hs_gsva_c7_sig <- get_sig_gsva_categories(hs_gsva_c7, excel = "excel/hs_macrophage_gsva_c7_sig.xlsx")

## Starting limma pairwise comparison.

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

## Using the libsize from expt$libsize.

## Limma step 1/6: choosing model.

## Choosing the non-intercept containing model.

## Assuming this data is similar to a micro array and not performign voom.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

## The factor inf_sb_z2.2 has 6 rows.

## The factor inf_sb_z2.3 has 6 rows.

## The factor inf_z2.2 has 6 rows.

## The factor inf_z2.3 has 6 rows.

## Testing each factor against the others.

## Scoring inf_sb_z2.2 against everything else.

## Scoring inf_sb_z2.3 against everything else.

## Scoring inf_z2.2 against everything else.

## Scoring inf_z2.3 against everything else.

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

hs_gsva_c7_sig$raw_plot

4 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)

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)

expt <- subset_expt(hs_macrophage, subset="macrophagetreatment!='uninf'") %>%
  subset_expt(subset="macrophagetreatment!='uninf_sb'")
expt <- set_expt_conditions(expt, fact="macrophagezymodeme")

symbol_vector <- fData(expt)[[symbol_column]]
names(symbol_vector) <- rownames(fData(expt))
symbol_df <- as.data.frame(symbol_vector)

assay_df <- merge(symbol_df, as.data.frame(exprs(expt)), 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(pData(expt))
factor_df[["SampleID"]] <- rownames(factor_df)
factor_df <- dplyr::select(factor_df, "SampleID", everything())
factor_df <- factor_df[, c("SampleID", factors)]

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

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

if (!isTRUE(get0("skip_load"))) {
  pander::pander(sessionInfo())
  message(paste0("This is hpgltools commit: ", get_git_commit()))
  message(paste0("Saving to ", savefile))
  tmp <- sm(saveme(filename = savefile))
}

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

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

## > git reset

## This is hpgltools commit:

## Saving to tmrc2_macrophage_visualization_202209.rda.xz

tmp <- loadme(filename = savefile)