TMRC2 202301: Promastigote (mostly) Differential Expression Analyses.

1 Contrasts

zymodeme_keeper <- list(
    "zymodeme" = c("z23", "z22"))
susceptibility_keepers <- list(
    "resistant_sensitive" = c("resistant", "sensitive"),
    "resistant_ambiguous" = c("resistant", "ambiguous"),
    "sensitive_ambiguous" = c("sensitive", "ambiguous"))

1.1 Zymodeme enzyme gene IDs

Najib read me an email listing off the gene names associated with the zymodeme classification. I took those names and cross referenced them against the Leishmania panamensis gene annotations and found the following:

They are:

1. ALAT: LPAL13_120010900 – alanine aminotransferase
2. ASAT: LPAL13_340013000 – aspartate aminotransferase
3. G6PD: LPAL13_000054100 – glucase-6-phosphate 1-dehydrogenase
4. NH: LPAL13_14006100, LPAL13_180018500 – inosine-guanine nucleoside hydrolase
5. MPI: LPAL13_320022300 (maybe) – mannose phosphate isomerase (I chose phosphomannose isomerase)

Given these 6 gene IDs (NH has two gene IDs associated with it), I can do some looking for specific differences among the various samples.

1.1.1 Expression levels of zymodeme genes

The following creates a colorspace (red to green) heatmap showing the observed expression of these genes in every sample.

my_genes <- c("LPAL13_120010900", "LPAL13_340013000", "LPAL13_000054100",
              "LPAL13_140006100", "LPAL13_180018500", "LPAL13_320022300",
              "other")
my_names <- c("ALAT", "ASAT", "G6PD", "NHv1", "NHv2", "MPI", "other")

zymo_six_genes <- exclude_genes_expt(lp_two_strains, ids = my_genes, method = "keep")

## remove_genes_expt(), before removal, there were 8710 genes, now there are 6.

## There are 84 samples which kept less than 90 percent counts.

## TMRC20001 TMRC20065 TMRC20005 TMRC20066 TMRC20039 TMRC20037 TMRC20038 TMRC20067 
##   0.12101   0.12835   0.13438   0.10804   0.13756   0.11451   0.11446   0.10858 
## TMRC20068 TMRC20041 TMRC20015 TMRC20009 TMRC20010 TMRC20016 TMRC20011 TMRC20012 
##   0.11269   0.12840   0.11072   0.11376   0.10292   0.10451   0.10967   0.11676 
## TMRC20013 TMRC20017 TMRC20014 TMRC20018 TMRC20019 TMRC20070 TMRC20020 TMRC20021 
##   0.11865   0.10594   0.10591   0.11517   0.12027   0.11179   0.11100   0.10764 
## TMRC20022 TMRC20024 TMRC20036 TMRC20069 TMRC20033 TMRC20026 TMRC20031 TMRC20076 
##   0.12952   0.11481   0.12100   0.11631   0.11188   0.13633   0.10013   0.12403 
## TMRC20073 TMRC20055 TMRC20079 TMRC20071 TMRC20078 TMRC20094 TMRC20042 TMRC20058 
##   0.12398   0.13783   0.12639   0.12003   0.13420   0.12021   0.13766   0.11900 
## TMRC20072 TMRC20059 TMRC20048 TMRC20088 TMRC20060 TMRC20077 TMRC20074 TMRC20063 
##   0.14675   0.10984   0.10679   0.12936   0.10882   0.12810   0.12494   0.12239 
## TMRC20053 TMRC20052 TMRC20064 TMRC20075 TMRC20051 TMRC20050 TMRC20049 TMRC20062 
##   0.12145   0.11690   0.12051   0.11261   0.13149   0.11526   0.14480   0.13429 
## TMRC20110 TMRC20080 TMRC20043 TMRC20083 TMRC20054 TMRC20085 TMRC20046 TMRC20089 
##   0.13812   0.12120   0.11306   0.12235   0.12723   0.12316   0.13182   0.11884 
## TMRC20090 TMRC20044 TMRC20105 TMRC20109 TMRC20098 TMRC20096 TMRC20097 TMRC20101 
##   0.11566   0.13691   0.12219   0.11684   0.11771   0.11364   0.11706   0.11784 
## TMRC20092 TMRC20082 TMRC20102 TMRC20099 TMRC20100 TMRC20087 TMRC20104 TMRC20086 
##   0.11465   0.10358   0.11399   0.11888   0.10820   0.12564   0.11723   0.10752 
## TMRC20107 TMRC20081 TMRC20106 TMRC20095 
##   0.09364   0.10335   0.09894   0.06566

strain_norm <- normalize_expt(zymo_six_genes, convert="rpkm", filter=TRUE, transform="log2")

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

zymo_heatmap <- plot_sample_heatmap(strain_norm, row_label = my_names)
zymo_heatmap

lp_norm <- normalize_expt(lp_two_strains, filter=TRUE, convert="rpkm",
                          norm="quant", transform="log2")

## Removing 152 low-count genes (8558 remaining).

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

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

zymo_heatmap_all <- plot_sample_heatmap(lp_norm)
zymo_heatmap_all

1.2 Compare to highly expressed, variant genes

I want to compare the above heatmap with one which is comprised of all genes with some ‘significantly high’ expression value and also a not-negligible coefficient of variance.

zymo_high_genes <- normalize_expt(lp_two_strains, filter="cv", cv_min=0.9)

## Removing 5310 low-count genes (3400 remaining).

high_strain_norm <- normalize_expt(zymo_high_genes, convert="rpkm", norm="quant", transform="log2")

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

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

zymo_heatmap <- plot_sample_heatmap(high_strain_norm, row_label = my_names)
zymo_heatmap

I think this plot suggests that the difference between the two primary strains is not really one of a few specific genes, but instead a global pattern.

2 Zymodeme differential expression

2.1 No attempt at batch estimation

two_zymo <- set_expt_conditions(lp_two_strains, fact = "zymodemecategorical") %>%
  subset_expt(subset = "condition!='unknown'")

## 
## z22 z23 
##  43  41

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

zymo_de_nobatch <- all_pairwise(two_zymo, filter = TRUE, model_batch = FALSE)

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

## 
## z22 z23 
##  43  41

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

zymo_table_nobatch <- combine_de_tables(
    zymo_de_nobatch, keepers = zymodeme_keeper,
    rda = glue::glue("rda/zymo_tables_nobatch-v{ver}.rda"),
    excel = glue::glue("excel/zymo_tables_nobatch-v{ver}.xlsx"))

## Deleting the file excel/zymo_tables_nobatch-v202301.xlsx before writing the tables.

## Starting combine_extracted_plots() with do_inverse as: FALSE.

## Plotting volcano plot of the DE results of z23_vs_z22 according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of z23_vs_z22 according to the columns: logFC and FDR using the expressionset colors.

## Plotting volcano plot of the DE results of z23_vs_z22 according to the columns: logFC and adj.P.Val using the expressionset colors.

## Saving de result as zymo_tables_nobatch-v202301 to rda/zymo_tables_nobatch-v202301.rda.

zymo_sig_nobatch <- extract_significant_genes(
    zymo_table_nobatch,
    according_to = "deseq", current_id = "GID", required_id = "GID",
    gmt = glue::glue("gmt/zymodeme_nobatch-v{ver}.gmt"),
    excel = glue::glue("excel/zymo_sig_nobatch_deseq-v{ver}.xlsx"))

## Deleting the file excel/zymo_sig_nobatch_deseq-v202301.xlsx before writing the tables.

## Plotting volcano plot of the DE results of z23_vs_z22 according to the columns: logFC and adj.P.Val using the expressionset colors.

## Number of up IDs in contrast zymodeme: 45.

## Number of down IDs in contrast zymodeme: 85.

2.1.1 Plot DE genes without batch estimation/adjustment

zymo_table_nobatch[["plots"]][["zymodeme"]][["deseq_ma_plots"]][["plot"]]

zymo_table_nobatch[["plots"]][["zymodeme"]][["deseq_vol_plots"]][["plot"]]

Log ratio, mean average plot and volcano plot of the comparison of the two primary zymodeme transcriptomes. When the transcriptomes of the two main strains (43 and 41 samples of z2.3 and z2.1) were compared without any attempt at batch/surrogate estimation with DESeq2, 45 and 85 genes were observed as significantly higher in strain z2.3 and z2.2 respectively using a cutoff of 1.0 logFC and 0.05 FDR adjusted p-value. There remain a large number of genes which are likely significantly different between the two strains, but fall below the 2-fold difference required for ‘significance.’ This follows prior observations that the parasite transcriptomes are constituitively expressed.

When the same data was plotted via a volcano plot, the relatively small range of fold changes compared to the large range of adjusted p-values is visible.

2.2 Attempt SVA estimate

zymo_de_sva <- all_pairwise(two_zymo, filter = TRUE, model_batch = "svaseq")

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

## 
## z22 z23 
##  43  41

## 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 (8558 remaining).

## Setting 431 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

zymo_table_sva <- combine_de_tables(
    zymo_de_sva, keepers = zymodeme_keeper,
    rda = glue::glue("rda/zymo_tables_sva-v{ver}.rda"),
    excel = glue::glue("excel/zymo_tables_sva-v{ver}.xlsx"))

## Deleting the file excel/zymo_tables_sva-v202301.xlsx before writing the tables.

## Starting combine_extracted_plots() with do_inverse as: FALSE.

## Plotting volcano plot of the DE results of z23_vs_z22 according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of z23_vs_z22 according to the columns: logFC and FDR using the expressionset colors.

## Plotting volcano plot of the DE results of z23_vs_z22 according to the columns: logFC and adj.P.Val using the expressionset colors.

## Saving de result as zymo_tables_sva-v202301 to rda/zymo_tables_sva-v202301.rda.

zymo_sig_sva <- extract_significant_genes(
    zymo_table_sva,
    according_to = "deseq",
    current_id = "GID", required_id = "GID",
    gmt = glue::glue("gmt/zymodeme_sva-v{ver}.gmt"),
    excel = glue::glue("excel/zymo_sig_sva-v{ver}.xlsx"))

## Deleting the file excel/zymo_sig_sva-v202301.xlsx before writing the tables.

## Plotting volcano plot of the DE results of z23_vs_z22 according to the columns: logFC and adj.P.Val using the expressionset colors.

## Number of up IDs in contrast zymodeme: 45.

## Number of down IDs in contrast zymodeme: 82.

2.2.1 Plot zymodeme DE genes with sva batch estimation/adjustment

When estimates from SVA were included in the statistical model used by EdgeR, DESeq2, and limma; a nearly identical view of the data emerged. I think this shows with a high degree of confidence, that sva is not having a significant effect on this dataset.

zymo_table_sva[["plots"]][["zymodeme"]][["deseq_ma_plots"]][["plot"]]

zymo_table_sva[["plots"]][["zymodeme"]][["deseq_vol_plots"]][["plot"]]

3 Parasite Susceptibility to Drug (Current)

This susceptibility comparison is using the ‘current’ dataset.

sus_de_nobatch <- all_pairwise(lp_susceptibility, filter = TRUE, model_batch = FALSE)

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

## 
## ambiguous resistant sensitive   unknown 
##        13        12        52        24

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

sus_table_nobatch <- combine_de_tables(
    sus_de_nobatch, keepers = susceptibility_keepers,
    excel = glue::glue("excel/sus_tables_nobatch-v{ver}.xlsx"))

## Deleting the file excel/sus_tables_nobatch-v202301.xlsx before writing the tables.

## Starting combine_extracted_plots() with do_inverse as: TRUE.

## Plotting volcano plot of the DE results of sensitive_vs_resistant according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_resistant according to the columns: logFC and FDR using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_resistant according to the columns: logFC and adj.P.Val using the expressionset colors.

## Starting combine_extracted_plots() with do_inverse as: FALSE.

## Plotting volcano plot of the DE results of resistant_vs_ambiguous according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of resistant_vs_ambiguous according to the columns: logFC and FDR using the expressionset colors.

## Plotting volcano plot of the DE results of resistant_vs_ambiguous according to the columns: logFC and adj.P.Val using the expressionset colors.

## Starting combine_extracted_plots() with do_inverse as: FALSE.

## Plotting volcano plot of the DE results of sensitive_vs_ambiguous according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_ambiguous according to the columns: logFC and FDR using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_ambiguous according to the columns: logFC and adj.P.Val using the expressionset colors.

sus_sig_nobatch <- extract_significant_genes(
    sus_table_nobatch,
    excel = glue::glue("excel/sus_sig_nobatch-v{ver}.xlsx"))

## Deleting the file excel/sus_sig_nobatch-v202301.xlsx before writing the tables.

## Plotting volcano plot of the DE results of sensitive_vs_resistant according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of resistant_vs_ambiguous according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_ambiguous according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_resistant according to the columns: logFC and FDR using the expressionset colors.

## Plotting volcano plot of the DE results of resistant_vs_ambiguous according to the columns: logFC and FDR using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_ambiguous according to the columns: logFC and FDR using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_resistant according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of resistant_vs_ambiguous according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_ambiguous according to the columns: logFC and adj.P.Val using the expressionset colors.

sus_de_sva <- all_pairwise(lp_susceptibility, filter = TRUE, model_batch = "svaseq")

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

## 
## ambiguous resistant sensitive   unknown 
##        13        12        52        24

## 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 (8576 remaining).

## Setting 386 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

sus_table_sva <- combine_de_tables(
    sus_de_sva, keepers = susceptibility_keepers,
    excel = glue::glue("excel/sus_tables_sva-v{ver}.xlsx"))

## Deleting the file excel/sus_tables_sva-v202301.xlsx before writing the tables.

## Starting combine_extracted_plots() with do_inverse as: TRUE.

## Plotting volcano plot of the DE results of sensitive_vs_resistant according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_resistant according to the columns: logFC and FDR using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_resistant according to the columns: logFC and adj.P.Val using the expressionset colors.

## Starting combine_extracted_plots() with do_inverse as: FALSE.

## Plotting volcano plot of the DE results of resistant_vs_ambiguous according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of resistant_vs_ambiguous according to the columns: logFC and FDR using the expressionset colors.

## Plotting volcano plot of the DE results of resistant_vs_ambiguous according to the columns: logFC and adj.P.Val using the expressionset colors.

## Starting combine_extracted_plots() with do_inverse as: FALSE.

## Plotting volcano plot of the DE results of sensitive_vs_ambiguous according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_ambiguous according to the columns: logFC and FDR using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_ambiguous according to the columns: logFC and adj.P.Val using the expressionset colors.

sus_sig_sva <- extract_significant_genes(
    sus_table_sva, according_to = "deseq",
    excel = glue::glue("excel/sus_sig_sva-v{ver}.xlsx"))

## Deleting the file excel/sus_sig_sva-v202301.xlsx before writing the tables.

## Plotting volcano plot of the DE results of sensitive_vs_resistant according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of resistant_vs_ambiguous according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_ambiguous according to the columns: logFC and adj.P.Val using the expressionset colors.

## To get a more true sense of sensitive vs resistant with sva, we kind of need to get rid of the
## unknown samples and perhaps the ambiguous.
no_ambiguous <- subset_expt(lp_susceptibility, subset="condition!='ambiguous'") %>%
  subset_expt(subset="condition!='unknown'")

## subset_expt(): There were 101, now there are 88 samples.

## subset_expt(): There were 88, now there are 64 samples.

no_ambiguous_de_sva <- all_pairwise(no_ambiguous, filter = TRUE, model_batch = "svaseq")

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

## 
## resistant sensitive 
##        12        52

## 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 (8559 remaining).

## Setting 213 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

## Let us see if my keeper code will fail hard or soft with extra contrasts...
no_ambiguous_table_sva <- combine_de_tables(
    no_ambiguous_de_sva, keepers = susceptibility_keepers,
    excel = glue::glue("excel/no_ambiguous_tables_sva-v{ver}.xlsx"))

## Deleting the file excel/no_ambiguous_tables_sva-v202301.xlsx before writing the tables.

## Starting combine_extracted_plots() with do_inverse as: TRUE.

## Plotting volcano plot of the DE results of sensitive_vs_resistant according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_resistant according to the columns: logFC and FDR using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_resistant according to the columns: logFC and adj.P.Val using the expressionset colors.

## Warning in extract_keepers_lst(extracted, keepers, legend[["table_names"]], :
## FOUND NEITHER resistant_vs_ambiguous NOR ambiguous_vs_resistant!

no_ambiguous_sig_sva <- extract_significant_genes(
    no_ambiguous_table_sva, according_to = "deseq",
    excel = glue::glue("excel/no_ambiguous_sig_sva-v{ver}.xlsx"))

## Deleting the file excel/no_ambiguous_sig_sva-v202301.xlsx before writing the tables.
## Plotting volcano plot of the DE results of sensitive_vs_resistant according to the columns: logFC and adj.P.Val using the expressionset colors.

3.0.1 Plot zymodeme DE genes with sva batch estimation/adjustment

sus_table_nobatch[["plots"]][["resistant_sensitive"]][["deseq_ma_plots"]][["plot"]]

sus_table_nobatch[["plots"]][["resistant_sensitive"]][["deseq_vol_plots"]][["plot"]]

sus_table_sva[["plots"]][["resistant_sensitive"]][["deseq_ma_plots"]][["plot"]]

sus_table_sva[["plots"]][["resistant_sensitive"]][["deseq_vol_plots"]][["plot"]]

no_ambiguous_table_sva[["plots"]][["resistant_sensitive"]][["deseq_ma_plots"]][["plot"]]

no_ambiguous_table_sva[["plots"]][["resistant_sensitive"]][["deseq_vol_plots"]][["plot"]]

Given that resistance/sensitivity tends to be correlated with strain, one might expect similar results. One caveat in this context though: there are fewer strains with resistance/sensitivity definitions. This when the analysis was repeated without the ambiguous/unknown samples, a few more genes were observed as significant.

4 Parasite Susceptibility to Drug (Historical)

This susceptibility comparison is using the historical dataset.

sushist_de_nobatch <- all_pairwise(lp_susceptibility_historical, filter = TRUE, model_batch = FALSE)

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

## 
## ambiguous resistant sensitive   unknown 
##         5        12        29        55

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

sushist_table_nobatch <- combine_de_tables(
    sushist_de_nobatch, keepers = susceptibility_keepers,
    excel = glue::glue("excel/sushist_tables_nobatch-v{ver}.xlsx"))

## Deleting the file excel/sushist_tables_nobatch-v202301.xlsx before writing the tables.

## Starting combine_extracted_plots() with do_inverse as: TRUE.

## Plotting volcano plot of the DE results of sensitive_vs_resistant according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_resistant according to the columns: logFC and FDR using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_resistant according to the columns: logFC and adj.P.Val using the expressionset colors.

## Starting combine_extracted_plots() with do_inverse as: FALSE.

## Plotting volcano plot of the DE results of resistant_vs_ambiguous according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of resistant_vs_ambiguous according to the columns: logFC and FDR using the expressionset colors.

## Plotting volcano plot of the DE results of resistant_vs_ambiguous according to the columns: logFC and adj.P.Val using the expressionset colors.

## Starting combine_extracted_plots() with do_inverse as: FALSE.

## Plotting volcano plot of the DE results of sensitive_vs_ambiguous according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_ambiguous according to the columns: logFC and FDR using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_ambiguous according to the columns: logFC and adj.P.Val using the expressionset colors.

sushist_sig_nobatch <- extract_significant_genes(
    sushist_table_nobatch,
    excel = glue::glue("excel/sushist_sig_nobatch-v{ver}.xlsx"))

## Deleting the file excel/sushist_sig_nobatch-v202301.xlsx before writing the tables.

## Plotting volcano plot of the DE results of sensitive_vs_resistant according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of resistant_vs_ambiguous according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_ambiguous according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_resistant according to the columns: logFC and FDR using the expressionset colors.

## Plotting volcano plot of the DE results of resistant_vs_ambiguous according to the columns: logFC and FDR using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_ambiguous according to the columns: logFC and FDR using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_resistant according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of resistant_vs_ambiguous according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_ambiguous according to the columns: logFC and adj.P.Val using the expressionset colors.

sushist_de_sva <- all_pairwise(lp_susceptibility_historical, filter = TRUE, model_batch = "svaseq")

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

## 
## ambiguous resistant sensitive   unknown 
##         5        12        29        55

## 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 (8576 remaining).

## Setting 374 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

sushist_table_sva <- combine_de_tables(
    sushist_de_sva, keepers = susceptibility_keepers,
    excel = glue::glue("excel/sushist_tables_sva-v{ver}.xlsx"))

## Deleting the file excel/sushist_tables_sva-v202301.xlsx before writing the tables.

## Starting combine_extracted_plots() with do_inverse as: TRUE.

## Plotting volcano plot of the DE results of sensitive_vs_resistant according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_resistant according to the columns: logFC and FDR using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_resistant according to the columns: logFC and adj.P.Val using the expressionset colors.

## Starting combine_extracted_plots() with do_inverse as: FALSE.

## Plotting volcano plot of the DE results of resistant_vs_ambiguous according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of resistant_vs_ambiguous according to the columns: logFC and FDR using the expressionset colors.

## Plotting volcano plot of the DE results of resistant_vs_ambiguous according to the columns: logFC and adj.P.Val using the expressionset colors.

## Starting combine_extracted_plots() with do_inverse as: FALSE.

## Plotting volcano plot of the DE results of sensitive_vs_ambiguous according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_ambiguous according to the columns: logFC and FDR using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_ambiguous according to the columns: logFC and adj.P.Val using the expressionset colors.

sushist_sig_sva <- extract_significant_genes(
    sushist_table_sva, according_to = "deseq",
    excel = glue::glue("excel/sushist_sig_sva-v{ver}.xlsx"))

## Deleting the file excel/sushist_sig_sva-v202301.xlsx before writing the tables.

## Plotting volcano plot of the DE results of sensitive_vs_resistant according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of resistant_vs_ambiguous according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of sensitive_vs_ambiguous according to the columns: logFC and adj.P.Val using the expressionset colors.

5 Cure/Fail association

##cf_nb_input <- subset_expt(cf_expt, subset="condition!='unknown'")
cf_de_nobatch <- all_pairwise(lp_cf_known, filter = TRUE, model_batch = FALSE)

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

## 
## cure fail 
##   38   38

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

cf_table_nobatch <- combine_de_tables(cf_de_nobatch, excel = glue::glue("excel/cf_tables_nobatch-v{ver}.xlsx"))

## Deleting the file excel/cf_tables_nobatch-v202301.xlsx before writing the tables.

## Starting combine_extracted_plots() with do_inverse as: FALSE.

## Plotting volcano plot of the DE results of fail_vs_cure according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of fail_vs_cure according to the columns: logFC and FDR using the expressionset colors.

## Plotting volcano plot of the DE results of fail_vs_cure according to the columns: logFC and adj.P.Val using the expressionset colors.

cf_sig_nobatch <- extract_significant_genes(cf_table_nobatch, excel = glue::glue("excel/cf_sig_nobatch-v{ver}.xlsx"))

## Deleting the file excel/cf_sig_nobatch-v202301.xlsx before writing the tables.
## Plotting volcano plot of the DE results of fail_vs_cure according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of fail_vs_cure according to the columns: logFC and FDR using the expressionset colors.

## Plotting volcano plot of the DE results of fail_vs_cure according to the columns: logFC and adj.P.Val using the expressionset colors.

cf_de <- all_pairwise(lp_cf_known, filter = TRUE, model_batch = "svaseq")

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

## 
## cure fail 
##   38   38

## 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 (8548 remaining).

## Setting 117 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

cf_table <- combine_de_tables(cf_de, excel = glue::glue("excel/cf_tables-v{ver}.xlsx"))

## Deleting the file excel/cf_tables-v202301.xlsx before writing the tables.

## Starting combine_extracted_plots() with do_inverse as: FALSE.

## Plotting volcano plot of the DE results of fail_vs_cure according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of fail_vs_cure according to the columns: logFC and FDR using the expressionset colors.

## Plotting volcano plot of the DE results of fail_vs_cure according to the columns: logFC and adj.P.Val using the expressionset colors.

cf_sig <- extract_significant_genes(cf_table, excel = glue::glue("excel/cf_sig-v{ver}.xlsx"))

## Deleting the file excel/cf_sig-v202301.xlsx before writing the tables.
## Plotting volcano plot of the DE results of fail_vs_cure according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of fail_vs_cure according to the columns: logFC and FDR using the expressionset colors.

## Plotting volcano plot of the DE results of fail_vs_cure according to the columns: logFC and adj.P.Val using the expressionset colors.

5.1 Cure/Fail DE plots

It is not surprising that few or no genes are deemed significantly differentially expressed across samples which were taken from cure or fail patients.

cf_table_nobatch[["plots"]][["fail_vs_cure"]][["deseq_ma_plots"]][["plot"]]

dev <- pp(file = "images/cf_ma.png")
cf_table[["plots"]][["fail_vs_cure"]][["deseq_ma_plots"]][["plot"]]
closed <- dev.off()
cf_table[["plots"]][["fail_vs_cure"]][["deseq_ma_plots"]][["plot"]]

6 Combining the macrophage infected amastigotes with in-vitro promastigotes

One query we have not yet addressed: what are the similarities and differences among the strains used to infect the macrophage samples and the promastigote samples used in the TMRC2 parasite data?

tmrc2_macrophage_norm <- normalize_expt(lp_macrophage, transform="log2", convert="cpm",
                                        norm="quant", filter=TRUE)

## Removing 169 low-count genes (8541 remaining).

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

all_tmrc2 <- combine_expts(lp_expt, lp_macrophage)

## 
##   b2904 unknown    z1.0    z1.5    z2.0    z2.1    z2.2    z2.3    z2.4    z3.0 
##       1       2       1       1       1       7      54      50       2       1 
##    z3.2 
##       1

all_nosb <- all_tmrc2
pData(all_nosb)[["stage"]] <- "promastigote"
na_idx <- is.na(pData(all_nosb)[["macrophagetreatment"]])
pData(all_nosb)[na_idx, "macrophagetreatment"] <- "undefined"
all_nosb <- subset_expt(all_nosb, subset="macrophagetreatment!='inf_sb'")

## subset_expt(): There were 121, now there are 119 samples.

ama_idx <- pData(all_nosb)[["macrophagetreatment"]] == "inf"
pData(all_nosb)[ama_idx, "stage" ] <- "amastigote"
pData(all_nosb)[["batch"]] <- pData(all_nosb)[["stage"]]

I think the above picture is sort of the opposite of what we want to compare in a DE analysis for this set of data, e.g. we want to compare promastigotes from amastigotes?

all_nosb <- set_expt_batches(all_nosb, fact="condition") %>%
  set_expt_conditions(fact="stage")

## 
##   b2904 unknown    z1.0    z1.5    z2.0    z2.1    z2.2    z2.3    z2.4    z3.0 
##       1       2       1       1       1       7      54      48       2       1 
##    z3.2 
##       1 
## 
##   amastigote promastigote 
##           18          101

two_zymo <- subset_expt(all_nosb, subset="zymodemecategorical=='z22'|zymodemecategorical=='z23'|zymodemecategorical=='unknown'")

## subset_expt(): There were 119, now there are 86 samples.

pro_ama <- all_pairwise(all_nosb, filter=TRUE, model_batch="svaseq")

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

## 
##   amastigote promastigote 
##           18          101

## 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 (8585 remaining).

## Setting 696 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

pro_ama_table <- combine_de_tables(
    pro_ama,
    excel = glue::glue("excel/tmrc2_pro_vs_ama_table-v{ver}.xlsx"))

## Deleting the file excel/tmrc2_pro_vs_ama_table-v202301.xlsx before writing the tables.

## Starting combine_extracted_plots() with do_inverse as: FALSE.

## Plotting volcano plot of the DE results of promastigote_vs_amastigote according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of promastigote_vs_amastigote according to the columns: logFC and FDR using the expressionset colors.

## Plotting volcano plot of the DE results of promastigote_vs_amastigote according to the columns: logFC and adj.P.Val using the expressionset colors.

pro_ama_sig <- extract_significant_genes(
    pro_ama_table,
    excel = glue::glue("excel/tmrc2_pro_vs_ama_sig-v{ver}.xlsx"))

## Deleting the file excel/tmrc2_pro_vs_ama_sig-v202301.xlsx before writing the tables.
## Plotting volcano plot of the DE results of promastigote_vs_amastigote according to the columns: logFC and adj.P.Val using the expressionset colors.

## Plotting volcano plot of the DE results of promastigote_vs_amastigote according to the columns: logFC and FDR using the expressionset colors.

## Plotting volcano plot of the DE results of promastigote_vs_amastigote according to the columns: logFC and adj.P.Val using the expressionset colors.

6.0.1 Plot promastigote/amastigote DE genes

pro_ama_table[["plots"]][["promastigote_vs_amastigote"]][["deseq_ma_plots"]][["plot"]]

I am a little surprised by this plot, I somewhat expected there to be few genes which passed the 2-fold difference demarcation line.

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 911e7d4beebdc73267ec4be631a305574289efd3

## This is hpgltools commit: Tue Jan 17 10:36:44 2023 -0500: 911e7d4beebdc73267ec4be631a305574289efd3

## Saving to tmrc2_differential_expression_202301.rda.xz

tmp <- loadme(filename = savefile)