Examining a series of Xenopus laevis early development samples.

1 Introduction

My understanding of this experiment comes from a conversation from a little over a year ago along with a powerpoint presentation. If I understand/remember them sufficiently, the goal is to understand changes from an approximately 8 (8, I checked) or 16 cell embryo to the cell at 10-14 days (11, I checked). There are two specific segments of this 8 cell embryo which are important for later development. This experiment seeks to follow them during these ~ 2 weeks, I think.

2 My starting material

I have 18 samples; 12 of which are of the D11 segment and 6 of the V11 segment. They are split between groups collected via (I am guessing) a magnetic separation and FACS; the D11 samples are further split in half between those which have a prefix of ‘exp’ and those which do not. My assumption from the presentation is that the exp samples have an additional methionine treatment in the initial 8 cell treatment/application of dye.

My guess therefore is that are looking to make a few observations:

1. Compare/contrast groups of D11/V11 cells.
2. Similarities/differences between FACS and Mag samples.
3. What is the effect of the methionine treatment on the D11 samples?

3 Annotation information

I have been having some ensembl troubles recently, so for now I still just load the gff file for the genome I used. Oh and I just realized I can just download the genbank file and use it for annotations.

I used the gene ID and gene type when counting, so let us just pull those annotations because there is an absurd number of entries in the Xenopus genome.

While I am waiting, I will grab the xenopus genbank file from ensembl. Oh, they do not have laevis, only tropicalis; I guess I will grab that and then get the NCBI genbank file. I already did.

I wonder if these libraries are polyA or riboZero? If so that will likely change the set of annotations I want. I can figure out the answer to this question via IGV, I will do so momentarily.

xl_annot <- load_gff_annotations("reference/xenopus_laevis_v10.1.gff",
                                 id_col = "gene", type = "gene")

## Returning a df with 42 columns and 44457 rows.

rownames(xl_annot) <- make.names(xl_annot[["gene"]], unique = TRUE)

4 Collect preprocessing metadata

The following should read the output logs from fastp/umitoos/hisat/whatever and add the portions of them I think are interesting as new columns to the metadata. I wrote an initial sample sheet for this experiment in the sample_sheets/ directory.

One thing I maybe should change: it does not default to seeking UMIs.

start_sheet <- "sample_sheets/202608_samples.xlsx"

umi_spec <- make_rnaseq_spec(umi = TRUE)
new_meta <- gather_preprocessing_metadata(start_sheet, specification = umi_spec,
                                          species = "xenopus_laevis_v10.1", tag = "gene")

## Did not find the condition column in the sample sheet.

## Filling it in as undefined.

## Did not find the batch column in the sample sheet.

## Filling it in as undefined.

## Checking the state of the condition column.

## Checking the state of the batch column.

## Checking the condition factor.

## Warning in mean.default(input_df[[column]], na.rm = TRUE): argument is not numeric or logical: returning NA
## Warning in mean.default(input_df[[column]], na.rm = TRUE): argument is not numeric or logical: returning NA
## Warning in mean.default(input_df[[column]], na.rm = TRUE): argument is not numeric or logical: returning NA
## Warning in mean.default(input_df[[column]], na.rm = TRUE): argument is not numeric or logical: returning NA
## Warning in mean.default(input_df[[column]], na.rm = TRUE): argument is not numeric or logical: returning NA
## Warning in mean.default(input_df[[column]], na.rm = TRUE): argument is not numeric or logical: returning NA
## Warning in mean.default(input_df[[column]], na.rm = TRUE): argument is not numeric or logical: returning NA
## Warning in mean.default(input_df[[column]], na.rm = TRUE): argument is not numeric or logical: returning NA
## Warning in mean.default(input_df[[column]], na.rm = TRUE): argument is not numeric or logical: returning NA
## Warning in mean.default(input_df[[column]], na.rm = TRUE): argument is not numeric or logical: returning NA
## Warning in mean.default(input_df[[column]], na.rm = TRUE): argument is not numeric or logical: returning NA
## Warning in mean.default(input_df[[column]], na.rm = TRUE): argument is not numeric or logical: returning NA
## Warning in mean.default(input_df[[column]], na.rm = TRUE): argument is not numeric or logical: returning NA
## Warning in mean.default(input_df[[column]], na.rm = TRUE): argument is not numeric or logical: returning NA
## Warning in mean.default(input_df[[column]], na.rm = TRUE): argument is not numeric or logical: returning NA
## Warning in mean.default(input_df[[column]], na.rm = TRUE): argument is not numeric or logical: returning NA
## Warning in mean.default(input_df[[column]], na.rm = TRUE): argument is not numeric or logical: returning NA
## Warning in mean.default(input_df[[column]], na.rm = TRUE): argument is not numeric or logical: returning NA

## Warning in mean.default(sort(x, partial = half + 0L:1L)[half + 0L:1L]): argument is not numeric or logical: returning NA
## Warning in mean.default(sort(x, partial = half + 0L:1L)[half + 0L:1L]): argument is not numeric or logical: returning NA

## Warning in seek_filenames(meta, input_file_spec, new_spec, basedir = basedir, : The filenames are 1entries, but the metadata has 18 entries.

## Writing new metadata to: sample_sheets/202608_samples_modified.xlsx

## Deleting the file sample_sheets/202608_samples_modified.xlsx before writing the tables.

head(new_meta$new_meta)

cell_group_colors <- get_xlsx_colors_by_column(start_sheet, column = "cell_group")
cell_group_colors

## $D11
##       rgb 
## "#8E86AE" 
## 
## $V11
##       rgb 
## "#E16173"

treatment_colors <- get_xlsx_colors_by_column(start_sheet, column = "second_treatment")
treatment_colors

## $methionine
##       rgb 
## "#069A2E" 
## 
## $none
##       rgb 
## "#3465A4"

replicate_colors <- get_xlsx_colors_by_column(start_sheet, column = "biological_replicate")
replicate_colors

## $R1
##       rgb 
## "#DEDCE6" 
## 
## $R2
##       rgb 
## "#B7B3CA" 
## 
## $R3
##       rgb 
## "#8E86AE"

5 Create some data structures

One before and one after deduplication. As of 20260901 this fails because the set of gene annotations has some utterly bizarre entries which confuse featureCounts. Here is an example:

trnar-acg       NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054371;NC_054372;NC_054375;NC_054375;NC_054375;NC_054375;NC_054375;NC_054375;NC_054375;NC_054375;NC_054375;NC_054375;NC_054375;NC_054375;NC_054375;NC_054375;NC_054376;NC_054376;NC_054376;NC_054376;NC_054376;NC_054376;NC_054376;NC_054376;NC_054376;NC_054376;NC_054376;NC_054376;NC_054376;NC_054383;NC_054383;NC_054383;NC_054383       2797755;2801616;2802215;2806099;2942504;2946451;3051317;5498592;3447894;3451745;4742813;4746647;4750495;4773076;4776921;4782481;4786315;4790161;4794009;4797855;4801600;4805438;4809279;4813105;4815460;4819225;4822985;4826749;4830581;4834437;4838270;4842110;4845965;4848857;7149949;7153778;7162533;7166294;7173491;7177342;7181167;7184982;7665012;7680070;7683902;7687731;151383083;130389122;115058601;145423349;145425243;145426617;145428088;145429477;145430087;145431106;145475661;145476095;145485457;145485891;145489412;145489847;45239162;121366482;121367244;121367511;121369517;121372444;121380444;121383079;121385801;121386939;121389021;121389528;45239559;126655978;126659839;126663701;126667563   2797827;2801688;2802287;2806171;2942576;2946523;3051389;5498664;3447966;3451817;4742885;4746719;4750567;4773148;4776993;4782553;4786387;4790233;4794081;4797927;4801672;4805510;4809351;4813177;4815532;4819297;4823057;4826821;4830653;4834509;4838342;4842182;4846037;4848929;7150021;7153850;7162605;7166366;7173562;7177414;7181239;7185054;7665084;7680142;7683974;7687803;151383155;130389194;115058673;145423421;145425315;145426689;145428160;145429549;145430159;145431178;145475733;145476167;145485529;145485963;145489484;145489919;45239234;121366554;121367316;121367583;121369589;121372516;121380516;121383151;121385873;121387011;121389093;121389600;45239631;126656050;126659911;126663773;126667635     -;-;-;-;-;-;-;-;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;-;-;-;-;-;-;-;-;-;-;-;-;+;+;+;+;+;+;+;+;+;+;+;+;+;+;-;+;+;+;+ 5766    0

This is intended to use a single gff entry to list out every arginine tRNA. Hmm, does xenopus use a different codon table? It is listed as -acg which I assume means the anticodon is ACG; but the default codon table has that as threonine.

Either way, this leads to a failure to load the count table until I added an option to make it smrt.

pre_se <- create_se(new_meta[["new_meta"]], file_column = "hisat_count_table", gene_info = xl_annot,
                    savefile = "rda/pre_se.rda") |>
  set_conditions(fact = "cell_group", colors = cell_group_colors) |>
  set_batches(fact = "treatment")

## Reading the sample metadata.

## Checking the state of the condition column.

## Checking the state of the batch column.

## Checking the condition factor.

## The sample definitions comprises: 18 rows(samples) and 87 columns(metadata fields).

## Matched 41731 annotations and counts.

## Some annotations were lost in merging, setting them to 'undefined'.

## The final summarized experiment has 41896 rows and 87 columns.

## The numbers of samples by condition are:

## 
## D11 V11 
##  12   6

## The number of samples by batch are:

## 
## FACS  Mag 
##    9    9

post_se <- create_se(new_meta[["new_meta"]], file_column = "umi_dedup_output_count",
                     gene_info = xl_annot, savefile = "rda/post_se.rda") |>
  set_conditions(fact = "cell_group", colors = cell_group_colors) |>
  set_batches(fact = "treatment")

## Reading the sample metadata.

## Checking the state of the condition column.

## Checking the state of the batch column.

## Checking the condition factor.

## The sample definitions comprises: 18 rows(samples) and 87 columns(metadata fields).

## Matched 41731 annotations and counts.

## Some annotations were lost in merging, setting them to 'undefined'.

## The final summarized experiment has 41896 rows and 87 columns.

## The numbers of samples by condition are:

## 
## D11 V11 
##  12   6

## The number of samples by batch are:

## 
## FACS  Mag 
##    9    9

combined <- paste0(conditions(pre_se), "_", colData(pre_se)[["second_treatment"]])
colData(pre_se)[["combined_treatment"]] <- combined
colData(post_se)[["combined_treatment"]] <- combined

NEB kit uses much lower Mg+2 concentration for the fragmentation, allowing a greater degree of flexibility when choosing fragmentation conditions. Less PCR and 2 fewer cleanup steps which should help lower the degree of duplication – this agrees well with the similarities of the pre/post quantreads observed.

Note to self, if I load the rda, I will need to reset the conditions/batches.

6 Reads pre/post

plot_quantreads(pre_se)

## Library sizes of 18 samples, 
## ranging from 14,231,005 to 19,233,956.

plot_quantreads(post_se)

## Library sizes of 18 samples, 
## ranging from 9,448,270 to 16,339,297.

plot_boxplot(pre_se)

## 254838 entries are 0.  We are on a log scale, adding 1 to the data.

## Plot describing the gene distribution from a dataset.

plot_boxplot(post_se)

## 254721 entries are 0.  We are on a log scale, adding 1 to the data.
## Plot describing the gene distribution from a dataset.

7 Distribution pre/post

pre_norm <- normalize(pre_se, transform = "log2", convert = "cpm", norm = "quant", filter = TRUE) |>
  set_batches(fact = "tech_replicate")

## Removing 22343 low-count genes (19553 remaining).

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

## The number of samples by batch are:

## 
## t1 t2 
##  6 12

plot_pca(pre_norm)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by D11, V11
## Shapes are defined by t1, t2.

post_norm <- normalize(post_se, transform = "log2", convert = "cpm", norm = "quant", filter = TRUE) |>
  set_batches(fact = "tech_replicate")

## Removing 22342 low-count genes (19554 remaining).

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

## The number of samples by batch are:

## 
## t1 t2 
##  6 12

plot_pca(post_norm)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by D11, V11
## Shapes are defined by t1, t2.

pre_treat <- set_conditions(pre_se, fact = "combined_treatment") |>
  set_batches(fact = "biological_replicate")

## The numbers of samples by condition are:

## 
## D11_methionine       D11_none       V11_none 
##              6              6              6

## Setting colors with no pre-defined colors, using the Dark2 palette.
## The number of samples by batch are:

## 
## R1 R2 R3 
##  6  6  6

post_treat <- set_conditions(post_se, fact = "combined_treatment") |>
  set_batches(fact = "biological_replicate")

## The numbers of samples by condition are:

## 
## D11_methionine       D11_none       V11_none 
##              6              6              6

## Setting colors with no pre-defined colors, using the Dark2 palette.
## The number of samples by batch are:

## 
## R1 R2 R3 
##  6  6  6

pre_treat_norm <- normalize(pre_treat, transform = "log2", convert = "cpm",
                            norm = "quant", filter = TRUE)

## Removing 22343 low-count genes (19553 remaining).

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

post_norm <- normalize(post_treat, transform = "log2", convert = "cpm",
                       norm = "quant", filter = TRUE)

## Removing 22342 low-count genes (19554 remaining).

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

plot_pca(pre_treat_norm)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by D11_methionine, D11_none, V11_none
## Shapes are defined by R1, R2, R3.

plot_pca(post_norm)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by D11_methionine, D11_none, V11_none
## Shapes are defined by R1, R2, R3.

7.1 See if sva has an effect

pre_treat_nb <- normalize(pre_treat, transform = "log2", convert = "cpm",
                          batch = "svaseq", filter = TRUE)

## Removing 22343 low-count genes (19553 remaining).

## transform_counts: Found 1832 values less than 0.

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

post_treat_nb <- normalize(post_treat, transform = "log2", convert = "cpm",
                           batch = "svaseq", filter = TRUE)

## Removing 22342 low-count genes (19554 remaining).

## transform_counts: Found 2059 values less than 0.

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

plot_pca(pre_treat_nb)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by D11_methionine, D11_none, V11_none
## Shapes are defined by R1, R2, R3.

plot_pca(post_treat_nb)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by D11_methionine, D11_none, V11_none
## Shapes are defined by R1, R2, R3.

8 Other comparisons of interest

Compare the three states with only the FACS

only_facs <- subset_se(post_treat, subset = 'treatment=="FACS"') |>
  set_batches(fact = "tech_replicate")

## The number of samples by batch are:

## 
## t1 t2 
##  3  6

only_norm <- normalize(only_facs, transform = "log2", convert = "cpm", filter = TRUE, norm = "tmm")

## Removing 23367 low-count genes (18529 remaining).

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

plot_pca(only_norm)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by D11_methionine, D11_none, V11_none
## Shapes are defined by t1, t2.

only_facs_de <- all_pairwise(only_facs)

## D11_methionine       D11_none       V11_none 
##              3              3              3 
## t1 t2 
##  3  6

## Basic step 0/3: Filtering data.

## 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 9860 entries to zero.

## converting counts to integer mode

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## This contrast put the denominator first.
## This contrast put the denominator first.
## This contrast put the denominator first.

## conditions
## D11_methionine       D11_none       V11_none 
##              3              3              3

## conditions
## D11_methionine       D11_none       V11_none 
##              3              3              3

## conditions
## D11_methionine       D11_none       V11_none 
##              3              3              3

keepers <- list(
  "d11m_vs_d11" = c("D11_methionine", "D11_none"),
  "v11_vs_d11" = c("V11_none", "D11_none"),
  "v11_vs_d11m" = c("V11_none", "D11_methionine"))
only_facs_tables <- combine_de_tables(only_facs_de, keepers = keepers, excel = glue("excel/treatment_comparisons_only_facs_tables-v{ver}.xlsx"))

## Deleting the file excel/treatment_comparisons_only_facs_tables-v20260924.xlsx before writing the tables.

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

only_facs_sig <- extract_significant_genes(only_facs_tables, excel = glue("excel/treatment_comparisons_only_facs_sig-v{ver}.xlsx"))

## Deleting the file excel/treatment_comparisons_only_facs_sig-v20260924.xlsx before writing the tables.

9 Compare to proteomics result

Likely plot from calculate_aucc of scatter of conditions vs each other.

peptide_values <- read_metadata("exogenous_data/FACS_cell_fate_Px_processed.xlsx")
d11_columns <- grepl(x = colnames(peptide_values), pattern = "d11_ctrl")
met_columns <- grepl(x = colnames(peptide_values), pattern = "d11_met")
v11_columns <- grepl(x = colnames(peptide_values), pattern = "v11_ctrl")
rownames(peptide_values) = make.names(peptide_values[["genesymbol"]], unique = TRUE)

peptide_values[["mean_d11"]] <- rowMeans(peptide_values[, d11_columns])
peptide_values[["mean_met"]] <- rowMeans(peptide_values[, met_columns])
peptide_values[["mean_v11"]] <- rowMeans(peptide_values[, v11_columns])

rna_means <- mean_by_factor(post_treat_nb)

## The factor D11_methionine has 6 rows.

## The factor D11_none has 6 rows.

## The factor V11_none has 6 rows.

combined <- merge(rna_means[["medians"]], peptide_values, by = "row.names")
rownames(combined) <- combined[["Row.names"]]
combined[["Row.names"]] <- NULL

pp(file = "images/rna_vs_proteomics_d11.png",
   image = plot_linear_scatter(combined[, c("D11_none", "mean_d11")])[["scatter"]])

pp(file = "images/rna_vs_proteomics_v11.png",
   image = plot_linear_scatter(combined[, c("V11_none", "mean_v11")])[["scatter"]])

pp(file = "images/rna_vs_proteomics_d11met.png",
   image = plot_linear_scatter(combined[, c("D11_methionine", "mean_met")])[["scatter"]])

10 Compare the three states

I am going to stop showing pre-deduplication.

de_nosva <- all_pairwise(post_treat, filter = TRUE)  ## batch in model using technical replicate

## D11_methionine       D11_none       V11_none 
##              6              6              6 
## R1 R2 R3 
##  6  6  6

## Removing 22342 low-count genes (19554 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 31757 entries to zero.

## converting counts to integer mode

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## This contrast put the denominator first.
## This contrast put the denominator first.
## This contrast put the denominator first.

## conditions
## D11_methionine       D11_none       V11_none 
##              6              6              6

## conditions
## D11_methionine       D11_none       V11_none 
##              6              6              6

## conditions
## D11_methionine       D11_none       V11_none 
##              6              6              6

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

## D11_methionine       D11_none       V11_none 
##              6              6              6

## Removing 22342 low-count genes (19554 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 31757 entries to zero.

## This received a matrix of SVs.

## converting counts to integer mode

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## This contrast put the denominator first.
## This contrast put the denominator first.
## This contrast put the denominator first.

## conditions
## D11_methionine       D11_none       V11_none 
##              6              6              6

## conditions
## D11_methionine       D11_none       V11_none 
##              6              6              6

## conditions
## D11_methionine       D11_none       V11_none 
##              6              6              6

de_nosva

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

de_sva

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

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

## The primary analysis performed 3 comparisons.

10.1 Create tables!

table_nosva <- combine_de_tables(de_nosva, excel = glue("excel/post_nosva_table-v{ver}.xlsx"))

## Deleting the file excel/post_nosva_table-v20260924.xlsx before writing the tables.

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

table_nosva

## A set of combined differential expression results.

##                        table deseq_sigup deseq_sigdown edger_sigup edger_sigdown limma_sigup limma_sigdown
## 1 D11_none_vs_D11_methionine        1231          2131        1243          2122        1496          1368
## 2 V11_none_vs_D11_methionine         185           981         145          1014         179           603
## 3       V11_none_vs_D11_none        3244          2445        3106          2568        2641          2491

## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## ℹ The deprecated feature was likely used in the UpSetR package.
##   Please report the issue at <https://github.com/hms-dbmi/UpSetR/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## ℹ The deprecated feature was likely used in the UpSetR package.
##   Please report the issue at <https://github.com/hms-dbmi/UpSetR/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

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

table_sva <- combine_de_tables(de_sva, excel = glue("excel/post_sva_table-v{ver}.xlsx"))

## Deleting the file excel/post_sva_table-v20260924.xlsx before writing the tables.
## Looking for subscript invalid names, end of extract_keepers.

table_sva

## A set of combined differential expression results.

##                        table deseq_sigup deseq_sigdown edger_sigup edger_sigdown limma_sigup limma_sigdown
## 1 D11_none_vs_D11_methionine        1586          2347        1653          2339        2040          1858
## 2 V11_none_vs_D11_methionine         123           942         130          1043         110           519
## 3       V11_none_vs_D11_none        2985          2614        2871          2729        2597          2619

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

10.2 Significant genes

sig_nosva <- extract_significant_genes(table_sva, excel = glue("excel/post_nosva_sig-v{ver}.xlsx"))

## Deleting the file excel/post_nosva_sig-v20260924.xlsx before writing the tables.

sig_nosva

## 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
## D11_none_vs_D11_methionine     1586       2347
## V11_none_vs_D11_methionine      123        942
## V11_none_vs_D11_none           2985       2614

sig_sva <- extract_significant_genes(table_sva, excel = glue("excel/post_sva_sig-v{ver}.xlsx"))

## Deleting the file excel/post_sva_sig-v20260924.xlsx before writing the tables.

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
## D11_none_vs_D11_methionine     1586       2347
## V11_none_vs_D11_methionine      123        942
## V11_none_vs_D11_none           2985       2614

10.3 Ontologies

I have two favorite tools for seeking out significant over representation: gProfiler2 and clusterProfiler. The latter depends on the xenopus annotation package ‘org.Xl.eg.db’ I will also need to make sure that the IDs I chose match it. I do not think gProfiler2 has Xenopus laevis, but does have tropicalis. I may be able to map genes across for that? I may give it a shot and see what happens.

10.3.1 Increased in the V11 samples vs. D11 untreated

comparison_table <- "V11_none_vs_D11_none"
sig_df_up_v11 <- sig_sva[["deseq"]][["ups"]][[comparison_table]]
table_df_v11 <- table_sva[["data"]][[1]]
test_up_cp_v11 <- simple_clusterprofiler(sig_df_up_v11, table_df_v11,
                                         orgdb = "org.Xl.eg.db", orgdb_from = "SYMBOL")

## There are 2534 genes deemed significant out of 17253.

## There are 1 entries which match the search string: Xenopus laevis.

## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).

## preparing geneSet collections...

## GSEA analysis...

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (26.32% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize = minSize, : There were 66 pathways for which P-values were not
## calculated properly due to unbalanced (positive and negative) gene-level statistic values. For such pathways pval, padj, NES, log2err are set
## to NA. You can try to increase the value of the argument nPermSimple (for example set it nPermSimple = 10000)

## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize = minSize, : For some of the pathways the P-values were likely
## overestimated. For such pathways log2err is set to NA.

## leading edge analysis...

## done...

## Reading KEGG annotation online: "https://rest.kegg.jp/link/xla/pathway"...

## Reading KEGG annotation online: "https://rest.kegg.jp/list/pathway/xla"...

## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).

## preparing geneSet collections...

## GSEA analysis...

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (27.16% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : For some of the pathways the P-values were likely overestimated. For such pathways log2err is set to NA.

## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize = minSize, : For some pathways, in reality P-values are less than
## 1e-10. You can set the `eps` argument to zero for better estimation.

## leading edge analysis...

## done...

## --> No gene can be mapped....

## --> Expected input gene ID: 57409,4283,6139,79364,3608,10040

## --> return NULL...

## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).

## preparing geneSet collections...

## --> Expected input gene ID: 79083,23314,84310,152195,84154,23034

## --> No gene can be mapped....

## Error in (function (cl, name, valueClass)  : 
##   'organism' is not a slot in class "NULL"

mf_up_v11_plots <- plot_enrichresult(test_up_cp_v11[["go_data"]][["MF_enrich"]])

## Warning in (function (model, data, ...) : Arguments in `...` must be used.
## ✖ Problematic argument:
## • by = "Count"
## ℹ Did you misspell an argument name?

## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.

## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.

mf_up_v11_plots[["vol"]]

mf_up_v11_plots[["tree"]]

mf_up_v11_plots[["dot"]]

bp_up_v11_plots <- plot_enrichresult(test_up_cp_v11[["go_data"]][["BP_enrich"]])

## Warning in (function (model, data, ...) : Arguments in `...` must be used.
## ✖ Problematic argument:
## • by = "Count"
## ℹ Did you misspell an argument name?

## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.

bp_up_v11_plots[["vol"]]

bp_up_v11_plots[["tree"]]

bp_up_v11_plots[["dot"]]

cc_up_v11_plots <- plot_enrichresult(test_up_cp_v11[["go_data"]][["CC_enrich"]])

## Warning in (function (model, data, ...) : Arguments in `...` must be used.
## ✖ Problematic argument:
## • by = "Count"
## ℹ Did you misspell an argument name?

## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.

xl_gsea_v11_plots <- plot_topn_gsea(test_up_cp_v11[["go_data"]][["GO_gse"]])
xl_gsea_v11_plots[[1]]

xl_gsea_v11_plots[[2]]

10.3.2 Increased in the D11 samples vs V11

sig_df_down_v11 <- sig_sva[["deseq"]][["downs"]][[comparison_table]]
test_down_cp_v11 <- simple_clusterprofiler(sig_df_down_v11, table_df_v11,
                                           orgdb = "org.Xl.eg.db", orgdb_from = "SYMBOL")

## There are 2235 genes deemed significant out of 17253.

## There are 1 entries which match the search string: Xenopus laevis.

## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).

## preparing geneSet collections...

## GSEA analysis...

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (26.32% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize = minSize, : There were 53 pathways for which P-values were not
## calculated properly due to unbalanced (positive and negative) gene-level statistic values. For such pathways pval, padj, NES, log2err are set
## to NA. You can try to increase the value of the argument nPermSimple (for example set it nPermSimple = 10000)

## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize = minSize, : For some of the pathways the P-values were likely
## overestimated. For such pathways log2err is set to NA.

## leading edge analysis...

## done...

## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).

## preparing geneSet collections...

## GSEA analysis...

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (27.16% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : For some of the pathways the P-values were likely overestimated. For such pathways log2err is set to NA.

## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize = minSize, : For some pathways, in reality P-values are less than
## 1e-10. You can set the `eps` argument to zero for better estimation.

## leading edge analysis...

## done...

## --> No gene can be mapped....

## --> Expected input gene ID: 10525,51053,90480,22882,1579,4489

## --> return NULL...

## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).

## preparing geneSet collections...

## --> Expected input gene ID: 3950,6907,554235,51534,7203,153562

## --> No gene can be mapped....

## Error in (function (cl, name, valueClass)  : 
##   'organism' is not a slot in class "NULL"

mf_down_v11_plots <- plot_enrichresult(test_down_cp_v11[["go_data"]][["MF_enrich"]])

## Warning in (function (model, data, ...) : Arguments in `...` must be used.
## ✖ Problematic argument:
## • by = "Count"
## ℹ Did you misspell an argument name?

## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.

## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.

mf_down_v11_plots[["vol"]]

mf_down_v11_plots[["tree"]]

mf_down_v11_plots[["dot"]]

bp_down_v11_plots <- plot_enrichresult(test_down_cp_v11[["go_data"]][["BP_enrich"]])

## Warning in (function (model, data, ...) : Arguments in `...` must be used.
## ✖ Problematic argument:
## • by = "Count"
## ℹ Did you misspell an argument name?

## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.

bp_down_v11_plots[["vol"]]

bp_down_v11_plots[["tree"]]

bp_down_v11_plots[["dot"]]

cc_down_v11_plots <- plot_enrichresult(test_down_cp_v11[["go_data"]][["CC_enrich"]])

## Warning in (function (model, data, ...) : Arguments in `...` must be used.
## ✖ Problematic argument:
## • by = "Count"
## ℹ Did you misspell an argument name?

## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
## ! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.

xl_gsea_down_v11_plots <- plot_topn_gsea(test_down_cp_v11[["go_data"]][["GO_gse"]])
xl_gsea_down_v11_plots[[1]]

xl_gsea_down_v11_plots[[2]]

pander::pander(sessionInfo())

R version 4.5.3 (2026-03-11)

Platform: x86_64-pc-linux-gnu

locale: LC_CTYPE=en_US.UTF-8, LC_NUMERIC=C, LC_TIME=en_US.UTF-8, LC_COLLATE=en_US.UTF-8, LC_MONETARY=en_US.UTF-8, LC_MESSAGES=en_US.UTF-8, LC_PAPER=en_US.UTF-8, LC_NAME=C, LC_ADDRESS=C, LC_TELEPHONE=C, LC_MEASUREMENT=en_US.UTF-8 and LC_IDENTIFICATION=C

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

other attached packages: GSVAdata(v.1.46.0), edgeR(v.4.8.2), ruv(v.0.9.7.2), hpgltools(v.2026.03), testthat(v.3.3.2) and reticulate(v.1.47.0)

loaded via a namespace (and not attached): fs(v.2.1.0), matrixStats(v.1.5.0), bitops(v.1.1-0), enrichplot(v.1.30.5), blockmodeling(v.1.1.8), devtools(v.2.5.2), httr(v.1.4.9), RColorBrewer(v.1.1-3), numDeriv(v.2016.8-1.1), tools(v.4.5.3), backports(v.1.5.1), R6(v.2.6.1), lazyeval(v.0.2.3), mgcv(v.1.9-4), withr(v.3.0.3), gridExtra(v.2.3.1), preprocessCore(v.1.72.0), cli(v.3.6.6), Biobase(v.2.70.0), scatterpie(v.0.2.6), EBSeq(v.2.8.0), labeling(v.0.4.3), sass(v.0.4.10), robustbase(v.0.99-7), mvtnorm(v.1.4-2), S7(v.0.2.2), readr(v.2.2.0), genefilter(v.1.92.0), ggarchery(v.0.4.4), Rsamtools(v.2.26.0), systemfonts(v.1.3.2), yulab.utils(v.0.2.5), ggupset(v.0.4.1), gson(v.0.2.1), DOSE(v.4.4.0), R.utils(v.2.13.0), dichromat(v.2.0-1), sessioninfo(v.1.2.4), limma(v.3.66.0), rstudioapi(v.0.19.0), RSQLite(v.3.53.3), BiocIO(v.1.20.0), generics(v.0.1.4), gridGraphics(v.0.5-1), vroom(v.1.7.1), gtools(v.3.9.5), zip(v.3.0.2), dplyr(v.1.2.1), GO.db(v.3.22.0), Matrix(v.1.7-6), S4Vectors(v.0.48.1), abind(v.1.4-8), R.methodsS3(v.1.8.2), lifecycle(v.1.0.5), yaml(v.2.3.12), SummarizedExperiment(v.1.40.0), gplots(v.3.3.0), qvalue(v.2.42.0), SparseArray(v.1.10.10), grid(v.4.5.3), blob(v.1.3.0), promises(v.1.5.0), crayon(v.1.5.3), ggtangle(v.0.1.2), lattice(v.0.23-1), cowplot(v.1.2.0), cigarillo(v.1.0.0), GenomicFeatures(v.1.62.0), annotate(v.1.88.0), KEGGREST(v.1.50.0), pillar(v.1.11.1), knitr(v.1.51), varhandle(v.2.0.6), fgsea(v.1.36.2), GenomicRanges(v.1.62.1), rjson(v.0.2.23), boot(v.1.3-32), corpcor(v.1.6.10), codetools(v.0.2-20), fastmatch(v.1.1-8), glue(v.1.8.1), ggiraph(v.0.9.6), ggfun(v.0.2.1), fontLiberation(v.0.1.0), Vennerable(v.3.1.0.9000), data.table(v.1.18.6.1), vctrs(v.0.7.3), png(v.0.1-9), treeio(v.1.34.0), Rdpack(v.2.6.6), gtable(v.0.3.6), cachem(v.1.1.0), openxlsx(v.4.2.9), xfun(v.0.60), rbibutils(v.2.4.1), S4Arrays(v.1.10.1), mime(v.0.13), RcppEigen(v.0.3.4.0.2), Seqinfo(v.1.0.0), reformulas(v.0.4.4), survival(v.3.8-11), NOISeq(v.2.54.0), iterators(v.1.0.14), org.Xl.eg.db(v.3.22.0), statmod(v.1.5.2), ellipsis(v.0.3.3), nlme(v.3.1-168), pbkrtest(v.0.5.5), ggtree(v.4.0.5), usethis(v.3.2.1), bit64(v.4.8.6), fontquiver(v.0.2.1), EnvStats(v.3.1.0), UpSetR(v.1.4.1), rprojroot(v.2.1.1), bslib(v.0.12.0), KernSmooth(v.2.23-26), otel(v.0.2.0), BiocGenerics(v.0.56.0), DBI(v.1.3.0), DESeq2(v.1.50.2), tidyselect(v.1.2.1), bit(v.4.6.0), compiler(v.4.5.3), curl(v.8.0.0), graph(v.1.88.1), desc(v.1.4.3), fontBitstreamVera(v.0.1.1), DelayedArray(v.0.36.1), plotly(v.4.12.1), rtracklayer(v.1.70.1), scales(v.1.4.0), caTools(v.1.18.4), DEoptimR(v.1.2-1), remaCor(v.0.0.20), RBGL(v.1.86.0), rappdirs(v.0.3.4), stringr(v.1.6.0), digest(v.0.6.39), minqa(v.1.2.8), variancePartition(v.1.40.2), rmarkdown(v.2.32), aod(v.1.3.3), XVector(v.0.50.0), RhpcBLASctl(v.0.23-42), htmltools(v.0.5.9), pkgconfig(v.2.0.3), lme4(v.2.0-6), MatrixGenerics(v.1.22.0), fastmap(v.1.2.0), rlang(v.1.3.0), htmlwidgets(v.1.6.4), shiny(v.1.14.0), farver(v.2.1.2), jquerylib(v.0.1.4), jsonlite(v.2.0.0), BiocParallel(v.1.44.0), GOSemSim(v.2.36.0), R.oo(v.1.27.1), RCurl(v.1.98-1.20), magrittr(v.2.0.5), ggplotify(v.0.1.3), patchwork(v.1.3.2), Rcpp(v.1.1.2), ape(v.5.8-1), ggnewscale(v.0.5.2), gdtools(v.0.5.1), stringi(v.1.8.9), brio(v.1.1.5), MASS(v.7.3-66), plyr(v.1.8.9), pkgbuild(v.1.4.8), parallel(v.4.5.3), ggrepel(v.0.9.8), Biostrings(v.2.78.0), splines(v.4.5.3), pander(v.0.6.6), hms(v.1.1.4), locfit(v.1.5-9.12), igraph(v.2.3.3), enrichit(v.0.2.1), reshape2(v.1.4.5), restez(v.2.1.5), stats4(v.4.5.3), pkgload(v.1.5.3), XML(v.3.99-0.24), evaluate(v.1.0.5), BiocManager(v.1.30.27), tzdb(v.0.5.0), nloptr(v.2.2.1), PROPER(v.1.42.0), foreach(v.1.5.2), tweenr(v.2.0.3), httpuv(v.1.6.17), tidyr(v.1.3.2), purrr(v.1.2.2), polyclip(v.1.10-7), ggplot2(v.4.0.3), ggforce(v.0.5.0), broom(v.1.0.13), xtable(v.1.8-8), restfulr(v.0.0.17), fANCOVA(v.0.6-1), tidytree(v.0.4.8), tidydr(v.0.0.6), later(v.1.4.8), viridisLite(v.0.4.3), tibble(v.3.3.1), lmerTest(v.3.2-1), clusterProfiler(v.4.18.4), aplot(v.0.3.1), GenomicAlignments(v.1.46.0), memoise(v.2.0.1), AnnotationDbi(v.1.72.0), IRanges(v.2.44.0), cluster(v.2.1.8.2), sva(v.3.58.0) and GSEABase(v.1.72.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 6dab42248ffc95c07d0cd3ea4f3e1e2c21c912de

## This is hpgltools commit: Mon Aug 31 11:10:05 2026 -0400: 6dab42248ffc95c07d0cd3ea4f3e1e2c21c912de

#this_save <- paste0(gsub(pattern = "\\.Rmd", replace = "", x = rmd_file), "-v", ver, ".rda.xz")
#message("Saving to ", this_save)
#tmp <- sm(saveme(filename = this_save))