TMRC2 202510: Promastigote (mostly) Differential Expression Analyses.

1 Changelog

202405: Changed excel output directory to match organization scheme in box. Generally this means files go to analyses/transcriptome/{type_of_contrast}/{date}/something_{suffix}.xlsx Where suffix is _table for the full tables and _sig for the significant genes and will include information about whether sva etc was used. 202405: Adding some goseq results.

2 Contrasts

** Note! ** The new definitions of susceptible/resistant are tighter than ever before, as a result there are no longer any ambiguous samples. Thus I removed the ambiguous contrasts in the following block.

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

3 Parasite Ontology data

Just a reminder that in data_structures.Rmd I created lp_go and lp_lengths

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

3.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(lp_two_strains, ids = my_genes, method = "keep")

## Note, I renamed this to subset_genes().

## subset_genes(), before removal, there were 8778 genes, now there are 6.

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

## TMRC20001 TMRC20002 TMRC20065 TMRC20004 TMRC20005 TMRC20066 TMRC20039 TMRC20037 
##   0.11877   0.08774   0.12694   0.11685   0.13185   0.10680   0.13447   0.11226 
## TMRC20038 TMRC20067 TMRC20068 TMRC20041 TMRC20015 TMRC20009 TMRC20010 TMRC20016 
##   0.11215   0.10714   0.11141   0.12526   0.10922   0.11134   0.10190   0.10331 
## TMRC20011 TMRC20012 TMRC20013 TMRC20017 TMRC20014 TMRC20018 TMRC20019 TMRC20070 
##   0.10780   0.11506   0.11721   0.10369   0.10397   0.11309   0.11818   0.11052 
## TMRC20020 TMRC20021 TMRC20022 TMRC20024 TMRC20036 TMRC20069 TMRC20033 TMRC20026 
##   0.10912   0.10639   0.12659   0.11237   0.11854   0.11358   0.10945   0.13413 
## TMRC20031 TMRC20076 TMRC20073 TMRC20055 TMRC20079 TMRC20071 TMRC20078 TMRC20094 
##   0.09875   0.12134   0.12195   0.13396   0.12432   0.11838   0.13113   0.11830 
## TMRC20042 TMRC20058 TMRC20072 TMRC20059 TMRC20048 TMRC20057 TMRC20088 TMRC20056 
##   0.13499   0.11714   0.14327   0.10855   0.10534   0.13253   0.12709   0.13211 
## TMRC20060 TMRC20077 TMRC20074 TMRC20063 TMRC20053 TMRC20052 TMRC20064 TMRC20075 
##   0.10701   0.12627   0.12216   0.12028   0.11890   0.11538   0.11880   0.11144 
## TMRC20051 TMRC20050 TMRC20049 TMRC20062 TMRC20110 TMRC20080 TMRC20043 TMRC20083 
##   0.12924   0.11354   0.14115   0.13212   0.13576   0.11954   0.11182   0.12063 
## TMRC20054 TMRC20085 TMRC20046 TMRC20093 TMRC20089 TMRC20047 TMRC20090 TMRC20044 
##   0.12548   0.12165   0.13005   0.13393   0.11746   0.12230   0.11411   0.13446 
## TMRC20045 TMRC20105 TMRC20108 TMRC20109 TMRC20098 TMRC20096 TMRC20101 TMRC20092 
##   0.12583   0.12061   0.11419   0.11505   0.11619   0.11116   0.11628   0.11304 
## TMRC20082 TMRC20102 TMRC20099 TMRC20100 TMRC20091 TMRC20084 TMRC20087 TMRC20103 
##   0.10253   0.11246   0.11740   0.10696   0.12652   0.11096   0.12368   0.13507 
## TMRC20104 TMRC20086 TMRC20107 TMRC20081 TMRC20095 
##   0.11464   0.10615   0.09249   0.10096   0.06536

strain_norm <- normalize(zymo_six_genes, convert = "rpkm", filter = TRUE, transform = "log2",
                         column = "annot_cds_length")

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

## There appears to be no gene length annotation data, here are the possible columns:

## gid, annot_external_db_name, gene_type, annot_aa_sequence_id, annot_annotated_go_component, annot_annotated_go_function, annot_annotated_go_id_component, annot_annotated_go_id_function, annot_annotated_go_id_process, annot_annotated_go_process, annot_anticodon, annot_apollo_link_out, annot_apollo_transcript_description, annot_cds, annot_cds_length, annot_chromosome, annot_coding_end, annot_coding_start, annot_ec_numbers, annot_ec_numbers_derived, annot_end_max, annot_exon_count, annot_external_db_version, annot_five_prime_utr_length, annot_gene_context_end, annot_gene_context_start, annot_gene_end_max, annot_gene_end_max_text, annot_gene_entrez_id, annot_gene_entrez_link, annot_gene_exon_count, annot_gene_hts_noncoding_snps, annot_gene_hts_nonsynonymous_snps, annot_gene_hts_nonsyn_syn_ratio, annot_gene_hts_stop_codon_snps, annot_gene_hts_synonymous_snps, annot_gene_location_text, annot_gene_name, annot_gene_ortholog_number, annot_gene_orthomcl_name, annot_gene_paralog_number, annot_gene_previous_ids, annot_gene_product, annot_gene_start_min, annot_gene_start_min_text, annot_gene_total_hts_snps, annot_gene_transcript_count, annot_gene_type, annot_genomic_sequence_length, annot_genus_species, annot_has_missing_transcripts, annot_interpro_description, annot_interpro_id, annot_isoelectric_point, annot_is_deprecated, annot_is_pseudo, annot_location_text, annot_map_location, annot_mcmc_location, annot_molecular_weight, annot_ncbi_tax_id, annot_orthomcl_link, annot_overview, annot_pfam_description, annot_pfam_id, annot_pirsf_description, annot_pirsf_id, annot_predicted_go_component, annot_predicted_go_function, annot_predicted_go_id_component, annot_predicted_go_id_function, annot_predicted_go_id_process, annot_predicted_go_process, annot_primary_key, annot_prob_map, annot_prob_mcmc, annot_prositeprofiles_description, annot_prositeprofiles_id, annot_protein_length, annot_protein_sequence, annot_protein_source_id, annot_pseudo_string, annot_sequence_database_name, annot_sequence_id, annot_signalp_peptide, annot_smart_description, annot_smart_id, annot_snpoverview, annot_so_id, annot_so_term_definition, annot_so_term_name, annot_so_version, annot_start_min, annot_strand, annot_strand_plus_minus, annot_superfamily_description, annot_superfamily_id, annot_three_prime_utr_length, annot_tigrfam_description, annot_tigrfam_id, annot_tm_count, annot_transcripts_found_per_gene, annot_transcript_index_per_gene, annot_transcript_length, annot_transcript_link, annot_transcript_product, annot_transcript_sequence, annot_trans_found_per_gene_internal, annot_uniprot_ids, annot_uniprot_links

## If one is appropriate, redo the function call with: column='good column'

## Error in hpgl_rpkm(count_table, annotations = annotations, length_column = length_column, : There appears to be no annotation data providing gene length.

zymo_heatmap <- plot_sample_heatmap(strain_norm, row_label = my_names)

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

zymo_heatmap

## Error: object 'zymo_heatmap' not found

lp_norm <- normalize(lp_two_strains, filter = TRUE, convert = "cpm",
                     norm = "quant", transform = "log2", column = "annot_cds_length")

## Removing 142 low-count genes (8636 remaining).

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

zymo_heatmap_all <- plot_sample_heatmap(lp_norm)
zymo_heatmap_all

3.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(lp_two_strains, filter = "cv", cv_min = 0.9)

## Removing 4731 low-count genes (4047 remaining).

high_strain_norm <- normalize(zymo_high_genes, convert = "rpkm",
                              norm = "quant", transform = "log2", column = "annot_cds_length")

## There appears to be no gene length annotation data, here are the possible columns:

## gid, annot_external_db_name, gene_type, annot_aa_sequence_id, annot_annotated_go_component, annot_annotated_go_function, annot_annotated_go_id_component, annot_annotated_go_id_function, annot_annotated_go_id_process, annot_annotated_go_process, annot_anticodon, annot_apollo_link_out, annot_apollo_transcript_description, annot_cds, annot_cds_length, annot_chromosome, annot_coding_end, annot_coding_start, annot_ec_numbers, annot_ec_numbers_derived, annot_end_max, annot_exon_count, annot_external_db_version, annot_five_prime_utr_length, annot_gene_context_end, annot_gene_context_start, annot_gene_end_max, annot_gene_end_max_text, annot_gene_entrez_id, annot_gene_entrez_link, annot_gene_exon_count, annot_gene_hts_noncoding_snps, annot_gene_hts_nonsynonymous_snps, annot_gene_hts_nonsyn_syn_ratio, annot_gene_hts_stop_codon_snps, annot_gene_hts_synonymous_snps, annot_gene_location_text, annot_gene_name, annot_gene_ortholog_number, annot_gene_orthomcl_name, annot_gene_paralog_number, annot_gene_previous_ids, annot_gene_product, annot_gene_start_min, annot_gene_start_min_text, annot_gene_total_hts_snps, annot_gene_transcript_count, annot_gene_type, annot_genomic_sequence_length, annot_genus_species, annot_has_missing_transcripts, annot_interpro_description, annot_interpro_id, annot_isoelectric_point, annot_is_deprecated, annot_is_pseudo, annot_location_text, annot_map_location, annot_mcmc_location, annot_molecular_weight, annot_ncbi_tax_id, annot_orthomcl_link, annot_overview, annot_pfam_description, annot_pfam_id, annot_pirsf_description, annot_pirsf_id, annot_predicted_go_component, annot_predicted_go_function, annot_predicted_go_id_component, annot_predicted_go_id_function, annot_predicted_go_id_process, annot_predicted_go_process, annot_primary_key, annot_prob_map, annot_prob_mcmc, annot_prositeprofiles_description, annot_prositeprofiles_id, annot_protein_length, annot_protein_sequence, annot_protein_source_id, annot_pseudo_string, annot_sequence_database_name, annot_sequence_id, annot_signalp_peptide, annot_smart_description, annot_smart_id, annot_snpoverview, annot_so_id, annot_so_term_definition, annot_so_term_name, annot_so_version, annot_start_min, annot_strand, annot_strand_plus_minus, annot_superfamily_description, annot_superfamily_id, annot_three_prime_utr_length, annot_tigrfam_description, annot_tigrfam_id, annot_tm_count, annot_transcripts_found_per_gene, annot_transcript_index_per_gene, annot_transcript_length, annot_transcript_link, annot_transcript_product, annot_transcript_sequence, annot_trans_found_per_gene_internal, annot_uniprot_ids, annot_uniprot_links

## If one is appropriate, redo the function call with: column='good column'

## Error in hpgl_rpkm(count_table, annotations = annotations, length_column = length_column, : There appears to be no annotation data providing gene length.

zymo_heatmap <- plot_sample_heatmap(high_strain_norm, row_label = my_names)

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

zymo_heatmap

## Error: object 'zymo_heatmap' not found

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.

4 Zymodeme differential expression

4.1 No attempt at batch estimation

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

## z2.2 z2.3 
##   42   41 
## undefined 
##        83

## Warning: attributes are not identical across measure variables; they will be
## dropped

## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## i Please use the `linewidth` argument instead.
## i The deprecated feature was likely used in the ggsankey package.
##   Please report the issue at <https://github.com/davidsjoberg/ggsankey/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Removing 150 low-count genes (8628 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 4662 entries to zero.

## Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : 
##   contrasts can be applied only to factors with 2 or more levels
## Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : 
##   contrasts can be applied only to factors with 2 or more levels
## Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : 
##   contrasts can be applied only to factors with 2 or more levels
## Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : 
##   contrasts can be applied only to factors with 2 or more levels
## conditions
## z22 z23 
##  42  41 
## Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : 
##   contrasts can be applied only to factors with 2 or more levels
## conditions
## z22 z23 
##  42  41 
## Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : 
##   contrasts can be applied only to factors with 2 or more levels
## conditions
## z22 z23 
##  42  41 
## Error in model.frame.default(object, data, xlev = xlev) : 
##   variable lengths differ (found for 'batch')

## Error in retlst[[meth]]: attempt to select less than one element in get1index

zymo_de_nobatch

## Error: object 'zymo_de_nobatch' not found

## Including the plots causes the rda file to balloon to 3.4Gb in the following invocation.
## Removing them results in... holy crap 2.1Mb
zymo_table_nobatch <- combine_de_tables(
    zymo_de_nobatch, keepers = zymodeme_keeper, label_column = "annot_gene_product",
    rda = glue("rda/zymo_tables_nobatch-v{ver}.rda"),
    excel = glue("{excel_out}/DE_Strain/{ver}/zymo_tables_nobatch-v{ver}.xlsx"))

## Error: object 'zymo_de_nobatch' not found

zymo_table_nobatch

## Error: object 'zymo_table_nobatch' not found

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

## Error: object 'zymo_table_nobatch' not found

zymo_sig_nobatch

## Error: object 'zymo_sig_nobatch' not found

4.1.0.1 Gene ontology comparing the strains

There are too few genes at our current stringencies for a meaningful result.

increased_z22 <- zymo_sig_nobatch[["deseq"]][["downs"]][["zymodeme"]]

## Error: object 'zymo_sig_nobatch' not found

increased_z23 <- zymo_sig_nobatch[["deseq"]][["ups"]][["zymodeme"]]

## Error: object 'zymo_sig_nobatch' not found

z22_goseq <- simple_goseq(increased_z22, go_db = lp_go, length_db = lp_lengths)

## Error: object 'increased_z22' not found

z23_goseq <- simple_goseq(increased_z23, go_db = lp_go, length_db = lp_lengths)

## Error: object 'increased_z23' not found

4.1.1 Plot DE genes without batch estimation/adjustment

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

## Error: object 'zymo_table_nobatch' not found

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

## Error: object 'zymo_table_nobatch' not found

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.

4.2 Attempt SVA estimate

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

## z2.2 z2.3 
##   42   41 
## undefined 
##        83

## Warning: attributes are not identical across measure variables; they will be
## dropped

## Removing 150 low-count genes (8628 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 4662 entries to zero.

## Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : 
##   contrasts can be applied only to factors with 2 or more levels
## Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : 
##   contrasts can be applied only to factors with 2 or more levels
## Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : 
##   contrasts can be applied only to factors with 2 or more levels
## Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : 
##   contrasts can be applied only to factors with 2 or more levels
## conditions
## z22 z23 
##  42  41 
## Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : 
##   contrasts can be applied only to factors with 2 or more levels
## conditions
## z22 z23 
##  42  41 
## Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : 
##   contrasts can be applied only to factors with 2 or more levels
## conditions
## z22 z23 
##  42  41 
## Error in model.frame.default(object, data, xlev = xlev) : 
##   variable lengths differ (found for 'batch')

## Error in retlst[[meth]]: attempt to select less than one element in get1index

zymo_de_sva

## Error: object 'zymo_de_sva' not found

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

## Error: object 'zymo_de_sva' not found

zymo_table_sva

## Error: object 'zymo_table_sva' not found

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

## Error: object 'zymo_table_sva' not found

zymo_sig_sva

## Error: object 'zymo_sig_sva' not found

4.2.0.1 Gene ontology comparing the strains

There are too few genes at our current stringencies for a meaningful result.

increased_z22 <- zymo_sig_sva[["deseq"]][["downs"]][["zymodeme"]]

## Error: object 'zymo_sig_sva' not found

increased_z23 <- zymo_sig_sva[["deseq"]][["ups"]][["zymodeme"]]

## Error: object 'zymo_sig_sva' not found

z22_goseq <- simple_goseq(increased_z22, go_db = lp_go, length_db = lp_lengths)

## Error: object 'increased_z22' not found

z23_goseq <- simple_goseq(increased_z23, go_db = lp_go, length_db = lp_lengths)

## Error: object 'increased_z23' not found

4.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"]]

## Error: object 'zymo_table_sva' not found

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

## Error: object 'zymo_table_sva' not found

5 Parasite Susceptibility to Drug (Current)

This susceptibility comparison is using the ‘current’ dataset.

Note again: we no longer have any ambiguous samples, so I commented out a portion of the following block.

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

## resistant sensitive 
##        46        46 
##    cure    fail unknown 
##      40      34      18

## Warning: attributes are not identical across measure variables; they will be
## dropped

## Removing 149 low-count genes (8629 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 5262 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
## resistant sensitive 
##        46        46

## conditions
## resistant sensitive 
##        46        46

## conditions
## resistant sensitive 
##        46        46

sus_de_nobatch

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

## This used a surrogate/batch estimate from: Existing surrogate matrix.

## The primary analysis performed 1 comparisons.

## The logFC agreement among the methods follows:

##                 snstv_vs_r
## basic_vs_deseq      0.9158
## basic_vs_dream      0.9291
## basic_vs_ebseq      0.8898
## basic_vs_edger      0.9147
## basic_vs_limma      0.9427
## basic_vs_noiseq     0.9365
## deseq_vs_dream      0.9853
## deseq_vs_ebseq      0.9493
## deseq_vs_edger      0.9937
## deseq_vs_limma      0.9638
## deseq_vs_noiseq     0.9783
## dream_vs_ebseq      0.9117
## dream_vs_edger      0.9776
## dream_vs_limma      0.9822
## dream_vs_noiseq     0.9569
## ebseq_vs_edger      0.9553
## ebseq_vs_limma      0.8957
## ebseq_vs_noiseq     0.9866
## edger_vs_limma      0.9547
## edger_vs_noiseq     0.9804
## limma_vs_noiseq     0.9405

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

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

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

sus_table_nobatch

## A set of combined differential expression results.

##                             table deseq_sigup deseq_sigdown edger_sigup
## 1 sensitive_vs_resistant-inverted          45           101          43
##   edger_sigdown limma_sigup limma_sigdown
## 1           106          63            91

## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## i Please use tidy evaluation idioms with `aes()`.
## i See also `vignette("ggplot2-in-packages")` for more information.
## i The deprecated feature was likely used in the UpSetR package.
##   Please report the issue to the authors.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

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

## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## i Please use the `linewidth` argument instead.
## i The deprecated feature was likely used in the UpSetR package.
##   Please report the issue to the authors.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Plot describing unique/shared genes in a differential expression table.

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

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

## resistant sensitive 
##        46        46 
##    cure    fail unknown 
##      40      34      18

## Warning: attributes are not identical across measure variables; they will be
## dropped

## Removing 149 low-count genes (8629 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 5262 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
## resistant sensitive 
##        46        46

## conditions
## resistant sensitive 
##        46        46

## conditions
## resistant sensitive 
##        46        46

sus_de_sva

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

## This used a surrogate/batch estimate from: Existing surrogate matrix.

## The primary analysis performed 1 comparisons.

## The logFC agreement among the methods follows:

##                 snstv_vs_r
## basic_vs_deseq      0.9158
## basic_vs_dream      0.9291
## basic_vs_ebseq      0.8898
## basic_vs_edger      0.9147
## basic_vs_limma      0.9427
## basic_vs_noiseq     0.9365
## deseq_vs_dream      0.9853
## deseq_vs_ebseq      0.9493
## deseq_vs_edger      0.9937
## deseq_vs_limma      0.9638
## deseq_vs_noiseq     0.9783
## dream_vs_ebseq      0.9117
## dream_vs_edger      0.9776
## dream_vs_limma      0.9822
## dream_vs_noiseq     0.9569
## ebseq_vs_edger      0.9553
## ebseq_vs_limma      0.8957
## ebseq_vs_noiseq     0.9866
## edger_vs_limma      0.9547
## edger_vs_noiseq     0.9804
## limma_vs_noiseq     0.9405

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

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

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

sus_table_sva

## A set of combined differential expression results.

##                             table deseq_sigup deseq_sigdown edger_sigup
## 1 sensitive_vs_resistant-inverted          45           101          43
##   edger_sigdown limma_sigup limma_sigdown
## 1           106          63            91

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

## Plot describing unique/shared genes in a differential expression table.

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

## A set of genes deemed significant according to deseq.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##                     deseq_up deseq_down
## resistant_sensitive       45        101

## 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_se(lp_susceptibility, subset = "condition!='ambiguous'") %>%
##   subset_se(subset = "condition!='unknown'")
## no_ambiguous_de_sva <- all_pairwise(no_ambiguous, filter = TRUE, model_batch = "svaseq")
## no_ambiguous_de_sva
## 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("excel/no_ambiguous_tables_sva-v{ver}.xlsx"))
## no_ambiguous_table_sva
## no_ambiguous_sig_sva <- extract_significant_genes(
##     no_ambiguous_table_sva, according_to = "deseq",
##     excel = glue("excel/no_ambiguous_sig_sva-v{ver}.xlsx"))
## no_ambiguous_sig_sva

5.0.0.1 Gene ontology comparing the susceptibility

sus_sig_sva

## A set of genes deemed significant according to deseq.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##                     deseq_up deseq_down
## resistant_sensitive       45        101

increased_resistant <- sus_sig_sva[["deseq"]][["ups"]][["resistant_sensitive"]]
increased_sensitive <- sus_sig_sva[["deseq"]][["downs"]][["resistant_sensitive"]]
resistant_goseq <- simple_goseq(increased_resistant, go_db = lp_go, length_db = lp_lengths)

## Found 13 go_db genes and 45 length_db genes out of 45.

sensitive_goseq <- simple_goseq(increased_sensitive, go_db = lp_go, length_db = lp_lengths)

## Found 24 go_db genes and 101 length_db genes out of 101.

5.0.1 Plot Susceptibility DE genes with sva batch estimation/adjustment

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

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

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

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

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.

6 Comparing DE results from strain/sensitivity

## zymo_table_sva[["plots"]][["zymodeme"]][["deseq_ma_plots"]][["plot"]]
zy_df <- zymo_table_sva[["data"]][["zymodeme"]]

## Error: object 'zymo_table_sva' not found

sus_df <- sus_table_sva[["data"]][["resistant_sensitive"]]

both_df <- merge(zy_df, sus_df, by = "row.names")

## Error: object 'zy_df' not found

plot_df <- both_df[, c("deseq_logfc.x", "deseq_logfc.y")]

## Error: object 'both_df' not found

rownames(plot_df) <- both_df[["Row.names"]]

## Error: object 'both_df' not found

colnames(plot_df) <- c("z23_vs_z22", "sensitive_vs_resistant")

## Error: object 'plot_df' not found

compare <- plot_linear_scatter(plot_df)

## Error: object 'plot_df' not found

pp(file = "images/compare_sus_zy.png")
compare$scatter

## Error in compare$scatter: object of type 'closure' is not subsettable

dev.off()

## png 
##   2

compare$scatter

## Error in compare$scatter: object of type 'closure' is not subsettable

compare$cor

## Error in compare$cor: object of type 'closure' is not subsettable

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

## ambiguous resistant sensitive   unknown 
##         5        12        30        45 
##    cure    fail unknown 
##      40      34      18

## Warning: attributes are not identical across measure variables; they will be
## dropped

## Removing 149 low-count genes (8629 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 5262 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
## ambiguous resistant sensitive   unknown 
##         5        12        30        45

## conditions
## ambiguous resistant sensitive   unknown 
##         5        12        30        45

## conditions
## ambiguous resistant sensitive   unknown 
##         5        12        30        45

sushist_de_nobatch

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

## This used a surrogate/batch estimate from: Existing surrogate matrix.

## The primary analysis performed 6 comparisons.

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

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

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

sushist_table_nobatch

## A set of combined differential expression results.

##                             table deseq_sigup deseq_sigdown edger_sigup
## 1 sensitive_vs_resistant-inverted          46           125          45
##   edger_sigdown limma_sigup limma_sigdown
## 1           125          62            84

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

## Plot describing unique/shared genes in a differential expression table.

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

## A set of genes deemed significant according to limma, edger, deseq, ebseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##                     limma_up limma_down edger_up edger_down deseq_up deseq_down
## resistant_sensitive       62         84       45        125       46        125
##                     ebseq_up ebseq_down basic_up basic_down
## resistant_sensitive       36         61       49         88

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

## ambiguous resistant sensitive   unknown 
##         5        12        30        45 
##    cure    fail unknown 
##      40      34      18

## Warning: attributes are not identical across measure variables; they will be
## dropped

## Removing 149 low-count genes (8629 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 5262 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
## ambiguous resistant sensitive   unknown 
##         5        12        30        45

## conditions
## ambiguous resistant sensitive   unknown 
##         5        12        30        45

## conditions
## ambiguous resistant sensitive   unknown 
##         5        12        30        45

sushist_de_sva

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

## This used a surrogate/batch estimate from: Existing surrogate matrix.

## The primary analysis performed 6 comparisons.

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

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

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

sushist_table_sva

## A set of combined differential expression results.

##                             table deseq_sigup deseq_sigdown edger_sigup
## 1 sensitive_vs_resistant-inverted          46           125          45
##   edger_sigdown limma_sigup limma_sigdown
## 1           125          62            84

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

## Plot describing unique/shared genes in a differential expression table.

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

## A set of genes deemed significant according to deseq.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##                     deseq_up deseq_down
## resistant_sensitive       46        125

8 Cure/Fail association

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

## cure fail 
##   40   34 
## resistant sensitive 
##        36        38

## Warning: attributes are not identical across measure variables; they will be
## dropped

## Removing 154 low-count genes (8624 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 3817 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
## cure fail 
##   40   34

## conditions
## cure fail 
##   40   34

## conditions
## cure fail 
##   40   34

cf_de_nobatch

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

## This used a surrogate/batch estimate from: Existing surrogate matrix.

## The primary analysis performed 1 comparisons.

## The logFC agreement among the methods follows:

##                 fail_vs_cr
## basic_vs_deseq      0.8091
## basic_vs_dream      0.8638
## basic_vs_ebseq      0.8554
## basic_vs_edger      0.8363
## basic_vs_limma      0.8765
## basic_vs_noiseq     0.7507
## deseq_vs_dream      0.6499
## deseq_vs_ebseq      0.7856
## deseq_vs_edger      0.9930
## deseq_vs_limma      0.6458
## deseq_vs_noiseq     0.9068
## dream_vs_ebseq      0.8367
## dream_vs_edger      0.6746
## dream_vs_limma      0.9847
## dream_vs_noiseq     0.5225
## ebseq_vs_edger      0.8143
## ebseq_vs_limma      0.8234
## ebseq_vs_noiseq     0.6591
## edger_vs_limma      0.6698
## edger_vs_noiseq     0.9300
## limma_vs_noiseq     0.5183

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

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

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

cf_table_nobatch

## A set of combined differential expression results.

##          table deseq_sigup deseq_sigdown edger_sigup edger_sigdown limma_sigup
## 1 fail_vs_cure           0             9           1            17           0
##   limma_sigdown
## 1             0

## Only fail_vs_cure_down has information, cannot create an UpSet.

## Plot describing unique/shared genes in a differential expression table.

## NULL

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

## A set of genes deemed significant according to limma, edger, deseq, ebseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##              limma_up limma_down edger_up edger_down deseq_up deseq_down
## fail_vs_cure        0          0        1         17        0          9
##              ebseq_up ebseq_down basic_up basic_down
## fail_vs_cure        0          0        0          0

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

## cure fail 
##   40   34 
## resistant sensitive 
##        36        38

## Warning: attributes are not identical across measure variables; they will be
## dropped

## Removing 154 low-count genes (8624 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 3817 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
## cure fail 
##   40   34

## conditions
## cure fail 
##   40   34

## conditions
## cure fail 
##   40   34

cf_de

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

## This used a surrogate/batch estimate from: Existing surrogate matrix.

## The primary analysis performed 1 comparisons.

## The logFC agreement among the methods follows:

##                 fail_vs_cr
## basic_vs_deseq      0.8091
## basic_vs_dream      0.8638
## basic_vs_ebseq      0.8554
## basic_vs_edger      0.8363
## basic_vs_limma      0.8765
## basic_vs_noiseq     0.7507
## deseq_vs_dream      0.6499
## deseq_vs_ebseq      0.7856
## deseq_vs_edger      0.9930
## deseq_vs_limma      0.6458
## deseq_vs_noiseq     0.9068
## dream_vs_ebseq      0.8367
## dream_vs_edger      0.6746
## dream_vs_limma      0.9847
## dream_vs_noiseq     0.5225
## ebseq_vs_edger      0.8143
## ebseq_vs_limma      0.8234
## ebseq_vs_noiseq     0.6591
## edger_vs_limma      0.6698
## edger_vs_noiseq     0.9300
## limma_vs_noiseq     0.5183

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

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

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

cf_table

## A set of combined differential expression results.

##          table deseq_sigup deseq_sigdown edger_sigup edger_sigdown limma_sigup
## 1 fail_vs_cure           0             9           1            17           0
##   limma_sigdown
## 1             0

## Only fail_vs_cure_down has information, cannot create an UpSet.

## Plot describing unique/shared genes in a differential expression table.

## NULL

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

## A set of genes deemed significant according to limma, edger, deseq, ebseq, basic.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##              limma_up limma_down edger_up edger_down deseq_up deseq_down
## fail_vs_cure        0          0        1         17        0          9
##              ebseq_up ebseq_down basic_up basic_down
## fail_vs_cure        0          0        0          0

I am not going to mess with GO searches for this.

8.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"]]

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

9 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?

In my container image, this dataset is not currently loaded, so turning this off.

## I just fixed this in the datasets Rmd, but until that propagates just set it manually
annotation(lp_se) <- annotation(lp_macrophage)
tmrc2_macrophage_norm <- normalize(lp_macrophage, transform="log2", convert="cpm",
                                   norm="quant", filter=TRUE)

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

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

all_tmrc2 <- hpgltools:::combine_se(lp_se, lp_macrophage)

all_nosb <- all_tmrc2
colData(all_nosb)[["stage"]] <- "promastigote"
na_idx <- is.na(colData(all_nosb)[["macrophagetreatment"]])
colData(all_nosb)[na_idx, "macrophagetreatment"] <- "undefined"
all_nosb <- subset_se(all_nosb, subset = "macrophagetreatment!='inf_sb'")
ama_idx <- colData(all_nosb)[["macrophagetreatment"]] == "inf"
colData(all_nosb)[ama_idx, "stage" ] <- "amastigote"
colData(all_nosb)[["batch"]] <- colData(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_se_batches(all_nosb, fact = "condition") %>%
  set_se_conditions(fact = "stage")

## Error in set_se_conditions(., fact = "stage"): could not find function "set_se_conditions"

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

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

## z2.1 z2.2 z2.3 z2.4 
##    7   56   56    2 
##   amastigote promastigote 
##           29           92

## Warning: attributes are not identical across measure variables; they will be
## dropped

## Removing 94 low-count genes (8684 remaining).

## Warning in RColorBrewer::brewer.pal(12, "Dark2"): n too large, allowed maximum for palette Dark2 is 8
## Returning the palette you asked for with that many colors

## Potentially check over the experimental design, there appear to be missing values.

## Warning in plot_pca(data = mtrx, design = design, state = state, plot_colors =
## plot_colors, : There are NA values in the component data.  This can lead to
## weird plotting errors.

## Warning in RColorBrewer::brewer.pal(12, "Dark2"): n too large, allowed maximum for palette Dark2 is 8
## Returning the palette you asked for with that many colors

## Potentially check over the experimental design, there appear to be missing values.

## Warning in plot_pca(data = mtrx, design = design, state = state, plot_colors =
## plot_colors, : There are NA values in the component data.  This can lead to
## weird plotting errors.

## 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 13046 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
## z21 z22 z23 z24 
##   7  56  56   2

## conditions
## z21 z22 z23 z24 
##   7  56  56   2

## conditions
## z21 z22 z23 z24 
##   7  56  56   2

pro_ama_table <- combine_de_tables(
  pro_ama,
  excel = glue("{excel_out}/DE_promastigote_amastigote/{ver}/pro_vs_ama_table-v{ver}.xlsx"))

## Error in names(all_colors) <- rownames(colData(exp)): attempt to set an attribute on NULL

pro_ama_sig <- extract_significant_genes(
    pro_ama_table,
    excel = glue("{excel_out}/DE_promastigote_amastigote/{ver}/pro_vs_ama_sig-v{ver}.xlsx"))

## Error: object 'pro_ama_table' not found

9.0.0.1 Gene ontology comparing the life stages

increased_promastigote <- pro_ama_sig[["deseq"]][["ups"]][["promastigote_vs_amastigote"]]

## Error: object 'pro_ama_sig' not found

increased_amastigote <- pro_ama_sig[["deseq"]][["downs"]][["promastigote_vs_amastigote"]]

## Error: object 'pro_ama_sig' not found

promastigote_goseq <- simple_goseq(increased_promastigote, go_db = lp_go, length_db = lp_lengths)

## Error: object 'increased_promastigote' not found

promastigote_goseq

## Error: object 'promastigote_goseq' not found

amastigote_goseq <- simple_goseq(increased_amastigote, go_db = lp_go, length_db = lp_lengths)

## Error: object 'increased_amastigote' not found

amastigote_goseq

## Error: object 'amastigote_goseq' not found

## silly, topgo wants the gene id column to be 'ID', I should fix this.
colnames(lp_go) <- c("ID", "GO")
promastigote_topgo <- simple_topgo(increased_promastigote, go_db = lp_go)

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

enrichplot::dotplot(promastigote_topgo$enrich_results$bp)

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

amastigote_topgo <- simple_topgo(increased_amastigote, go_db = lp_go)

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

enrichplot::dotplot(amastigote_topgo$enrich_results$bp)

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

9.0.1 Plot promastigote/amastigote DE genes

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

## Error: object 'pro_ama_table' not found

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

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: edgeR(v.4.6.3), ruv(v.0.9.7.1), hpgltools(v.1.2), testthat(v.3.2.3), glue(v.1.8.0) and Heatplus(v.3.16.0)

loaded via a namespace (and not attached): fs(v.1.6.6), matrixStats(v.1.5.0), bitops(v.1.0-9), enrichplot(v.1.28.4), blockmodeling(v.1.1.8), devtools(v.2.4.6), httr(v.1.4.7), RColorBrewer(v.1.1-3), numDeriv(v.2016.8-1.1), tools(v.4.5.0), backports(v.1.5.0), R6(v.2.6.1), lazyeval(v.0.2.2), mgcv(v.1.9-3), withr(v.3.0.2), gridExtra(v.2.3), preprocessCore(v.1.70.0), cli(v.3.6.5), Biobase(v.2.68.0), labeling(v.0.4.3), EBSeq(v.2.6.0), sass(v.0.4.10), robustbase(v.0.99-6), mvtnorm(v.1.3-3), S7(v.0.2.0), genefilter(v.1.90.0), Rsamtools(v.2.24.0), systemfonts(v.1.3.1), yulab.utils(v.0.2.1), gson(v.0.1.0), DOSE(v.4.2.0), R.utils(v.2.13.0), dichromat(v.2.0-0.1), sessioninfo(v.1.2.3), limma(v.3.64.3), rstudioapi(v.0.17.1), RSQLite(v.2.4.3), generics(v.0.1.4), gridGraphics(v.0.5-1), BiocIO(v.1.18.0), gtools(v.3.9.5), zip(v.2.3.3), dplyr(v.1.1.4), GO.db(v.3.21.0), 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), gplots(v.3.2.0), qvalue(v.2.40.0), SparseArray(v.1.8.1), 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), GenomicFeatures(v.1.60.0), annotate(v.1.86.1), KEGGREST(v.1.48.1), 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), ggiraph(v.0.9.2), ggfun(v.0.2.0), fontLiberation(v.0.1.0), Vennerable(v.3.1.0.9000), data.table(v.1.17.8), remotes(v.2.5.0), vctrs(v.0.6.5), png(v.0.1-8), treeio(v.1.32.0), Rdpack(v.2.6.4), gtable(v.0.3.6), 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), reformulas(v.0.4.1), survival(v.3.8-3), NOISeq(v.2.52.0), iterators(v.1.0.14), statmod(v.1.5.0), ellipsis(v.0.3.2), nlme(v.3.1-168), pbkrtest(v.0.5.5), ggtree(v.3.99.0), usethis(v.3.2.1), bit64(v.4.6.0-1), fontquiver(v.0.2.1), EnvStats(v.3.1.0), UpSetR(v.1.4.0), GenomeInfoDb(v.1.44.2), rprojroot(v.2.1.1), bslib(v.0.9.0), KernSmooth(v.2.23-26), BiocGenerics(v.0.54.0), DBI(v.1.2.3), DESeq2(v.1.48.1), tidyselect(v.1.2.1), bit(v.4.6.0), compiler(v.4.5.0), curl(v.7.0.0), graph(v.1.86.0), desc(v.1.4.3), fontBitstreamVera(v.0.1.1), DelayedArray(v.0.34.1), plotly(v.4.11.0), rtracklayer(v.1.68.0), scales(v.1.4.0), caTools(v.1.18.3), DEoptimR(v.1.1-4), 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), ggsankey(v.0.0.99999), 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), fastmap(v.1.2.0), 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.4), GenomeInfoDbData(v.1.2.14), ggplotify(v.0.1.2), patchwork(v.1.3.2), Rcpp(v.1.1.0), ape(v.5.8-1), gdtools(v.0.4.4), stringi(v.1.8.7), brio(v.1.1.5), MASS(v.7.3-65), plyr(v.1.8.9), pkgbuild(v.1.4.8), parallel(v.4.5.0), ggrepel(v.0.9.6), forcats(v.1.0.0), Biostrings(v.2.76.0), splines(v.4.5.0), pander(v.0.6.6), locfit(v.1.5-9.12), igraph(v.2.1.4), fastcluster(v.1.3.0), reshape2(v.1.4.4), stats4(v.4.5.0), pkgload(v.1.4.1), XML(v.3.99-0.19), evaluate(v.1.0.4), BiocManager(v.1.30.26), nloptr(v.2.2.1), PROPER(v.1.40.0), foreach(v.1.5.2), httpuv(v.1.6.16), tidyr(v.1.3.1), purrr(v.1.1.0), ggplot2(v.4.0.0), broom(v.1.0.10), xtable(v.1.8-4), restfulr(v.0.0.16), fANCOVA(v.0.6-1), tidytree(v.0.4.6), later(v.1.4.3), viridisLite(v.0.4.2), tibble(v.3.3.0), lmerTest(v.3.1-3), clusterProfiler(v.4.16.0), aplot(v.0.2.8), memoise(v.2.0.1), AnnotationDbi(v.1.70.0), GenomicAlignments(v.1.44.0), IRanges(v.2.42.0), sva(v.3.56.0) and GSEABase(v.1.70.0)

message(paste0("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 159b347956b8e87168310c6d2f8096e8f512c924

## This is hpgltools commit: Mon Oct 13 17:35:31 2025 -0400: 159b347956b8e87168310c6d2f8096e8f512c924

## message(paste0("Saving to ", savefile))
## tmp <- sm(saveme(filename = savefile))

tmp <- loadme(filename = savefile)