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.

umi_spec <- make_rnaseq_spec(umi = TRUE)
new_meta <- gather_preprocessing_metadata("sample_sheets/202608_samples.xlsx", 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)

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

## Setting colors with no pre-defined colors, using the Dark2 palette.

## The number of samples by batch are:

## 
##   FACS    Mag MagCar 
##      9      6      3

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

## Setting colors with no pre-defined colors, using the Dark2 palette.

## The number of samples by batch are:

## 
##   FACS    Mag MagCar 
##      9      6      3

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

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

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.

7.1 See if sva has an effect

pre_nb <- normalize(pre_se, transform = "log2", convert = "cpm", batch = "svaseq", filter = TRUE) |>
  set_batches(fact = "tech_replicate")

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

## transform_counts: Found 1086 values less than 0.

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

## The number of samples by batch are:

## 
## t1 t2 
##  6 12

post_nb <- normalize(post_se, transform = "log2", convert = "cpm", batch = "svaseq", filter = TRUE) |>
  set_batches(fact = "tech_replicate")

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

## transform_counts: Found 1090 values less than 0.

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

## The number of samples by batch are:

## 
## t1 t2 
##  6 12

plot_pca(pre_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, V11
## Shapes are defined by t1, t2.

plot_pca(post_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, V11
## Shapes are defined by t1, t2.

hmm not as much as I might have guessed given the original distributions.

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

## D11 V11 
##  12   6 
##   FACS    Mag MagCar 
##      9      6      3

## Removing 22343 low-count genes (19553 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 32022 entries to zero.

## converting counts to integer mode

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## This contrast put the denominator first.

## conditions
## D11 V11 
##  12   6

## conditions
## D11 V11 
##  12   6

## conditions
## D11 V11 
##  12   6

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

## D11 V11 
##  12   6

## Removing 22343 low-count genes (19553 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 32022 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.

## conditions
## D11 V11 
##  12   6

## conditions
## D11 V11 
##  12   6

## conditions
## D11 V11 
##  12   6

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

## D11 V11 
##  12   6 
##   FACS    Mag MagCar 
##      9      6      3

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

## conditions
## D11 V11 
##  12   6

## conditions
## D11 V11 
##  12   6

## conditions
## D11 V11 
##  12   6

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

## D11 V11 
##  12   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.

## conditions
## D11 V11 
##  12   6

## conditions
## D11 V11 
##  12   6

## conditions
## D11 V11 
##  12   6

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

## The logFC agreement among the methods follows:

##                 V11_vs_D11
## basic_vs_deseq   0.8671891
## basic_vs_dream   0.9211830
## basic_vs_ebseq   0.8264844
## basic_vs_edger   0.8663105
## basic_vs_limma   0.9326426
## basic_vs_noiseq  0.9181156
## deseq_vs_dream   0.9562242
## deseq_vs_ebseq   0.9811104
## deseq_vs_edger   0.9991253
## deseq_vs_limma   0.9508950
## deseq_vs_noiseq  0.9861404
## dream_vs_ebseq   0.9282299
## dream_vs_edger   0.9550018
## dream_vs_limma   0.9951704
## dream_vs_noiseq  0.9625600
## ebseq_vs_edger   0.9862286
## ebseq_vs_limma   0.9226252
## ebseq_vs_noiseq  0.9596426
## edger_vs_limma   0.9496427
## edger_vs_noiseq  0.9849323
## limma_vs_noiseq  0.9623787

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

## The logFC agreement among the methods follows:

##                 V11_vs_D11
## basic_vs_deseq   0.7631688
## basic_vs_dream   0.7792545
## basic_vs_ebseq   0.8264844
## basic_vs_edger   0.7520345
## basic_vs_limma   0.7931846
## basic_vs_noiseq  0.9181156
## deseq_vs_dream   0.9520734
## deseq_vs_ebseq   0.9105133
## deseq_vs_edger   0.9669265
## deseq_vs_limma   0.9241889
## deseq_vs_noiseq  0.8966109
## dream_vs_ebseq   0.9040798
## dream_vs_edger   0.9619298
## dream_vs_limma   0.9715965
## dream_vs_noiseq  0.9028118
## ebseq_vs_edger   0.9231723
## ebseq_vs_limma   0.8916599
## ebseq_vs_noiseq  0.9596426
## edger_vs_limma   0.9325085
## edger_vs_noiseq  0.8911288
## limma_vs_noiseq  0.8992606

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

## The logFC agreement among the methods follows:

##                 V11_vs_D11
## basic_vs_deseq   0.8675544
## basic_vs_dream   0.9282444
## basic_vs_ebseq   0.8305458
## basic_vs_edger   0.8723469
## basic_vs_limma   0.9393038
## basic_vs_noiseq  0.9257728
## deseq_vs_dream   0.9505560
## deseq_vs_ebseq   0.9795649
## deseq_vs_edger   0.9956103
## deseq_vs_limma   0.9462991
## deseq_vs_noiseq  0.9803102
## dream_vs_ebseq   0.9265426
## dream_vs_edger   0.9549568
## dream_vs_limma   0.9952848
## dream_vs_noiseq  0.9630444
## ebseq_vs_edger   0.9863539
## ebseq_vs_limma   0.9222900
## ebseq_vs_noiseq  0.9575381
## edger_vs_limma   0.9504886
## edger_vs_noiseq  0.9839447
## limma_vs_noiseq  0.9637417

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

## The logFC agreement among the methods follows:

##                 V11_vs_D11
## basic_vs_deseq   0.8113203
## basic_vs_dream   0.8341480
## basic_vs_ebseq   0.8305458
## basic_vs_edger   0.8244049
## basic_vs_limma   0.8457657
## basic_vs_noiseq  0.9257728
## deseq_vs_dream   0.9546893
## deseq_vs_ebseq   0.9193795
## deseq_vs_edger   0.9767537
## deseq_vs_limma   0.9329747
## deseq_vs_noiseq  0.9161984
## dream_vs_ebseq   0.9208076
## dream_vs_edger   0.9752972
## dream_vs_limma   0.9780335
## dream_vs_noiseq  0.9312580
## ebseq_vs_edger   0.9383801
## ebseq_vs_limma   0.9109500
## ebseq_vs_noiseq  0.9575381
## edger_vs_limma   0.9524538
## edger_vs_noiseq  0.9311229
## limma_vs_noiseq  0.9290768

7.2 Create tables!

pre_table_nosva <- combine_de_tables(pre_de_nosva, excel = glue("excel/pre_nosva_table-v{ver}.xlsx"))

## Deleting the file excel/pre_nosva_table-v20260903.xlsx before writing the tables.

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

pre_table_nosva

## A set of combined differential expression results.

##        table deseq_sigup deseq_sigdown edger_sigup edger_sigdown limma_sigup
## 1 V11_vs_D11        1290          1767        1359          1724        1810
##   limma_sigdown
## 1          1069

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

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

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

pre_table_sva <- combine_de_tables(pre_de_sva, excel = glue("excel/pre_sva_table-v{ver}.xlsx"))

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

pre_table_sva

## A set of combined differential expression results.

##        table deseq_sigup deseq_sigdown edger_sigup edger_sigdown limma_sigup
## 1 V11_vs_D11         381          1287         377          1323         225
##   limma_sigdown
## 1           689

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

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

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

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

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

post_table_nosva

## A set of combined differential expression results.

##        table deseq_sigup deseq_sigdown edger_sigup edger_sigdown limma_sigup
## 1 V11_vs_D11        1322          1750        1325          1761        1780
##   limma_sigdown
## 1          1074

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

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

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

## Deleting the file excel/post_sva_table-v20260903.xlsx before writing the tables.

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

post_table_sva

## A set of combined differential expression results.

##        table deseq_sigup deseq_sigdown edger_sigup edger_sigdown limma_sigup
## 1 V11_vs_D11         860          1392         843          1442         460
##   limma_sigdown
## 1           810

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

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

7.3 Significant genes

pre_sig_nosva <- extract_significant_genes(pre_table_sva, excel = glue("excel/pre_nosva_sig-v{ver}.xlsx"))

## Deleting the file excel/pre_nosva_sig-v20260903.xlsx before writing the tables.

pre_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
## V11_vs_D11      381       1287

pre_sig_sva <- extract_significant_genes(pre_table_sva, excel = glue("excel/pre_sva_sig-v{ver}.xlsx"))

## Deleting the file excel/pre_sva_sig-v20260903.xlsx before writing the tables.

pre_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
## V11_vs_D11      381       1287

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

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

post_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
## V11_vs_D11      860       1392

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

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

post_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
## V11_vs_D11      860       1392

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

sig_df <- post_sig_sva[["deseq"]][["ups"]][[1]]
table_df <- post_table_sva[["data"]][[1]]
test_cp <- simple_clusterprofiler(sig_df, table_df, orgdb = "org.Xl.eg.db", orgdb_from = "SYMBOL")

## There are 749 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.31% 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, : 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...

## 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.09% 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, : 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: 26505,10096,9235,9491,7485,79085

## --> return NULL...

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

## preparing geneSet collections...

## --> Expected input gene ID: 64216,10555,441531,2195,136853,79602

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

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

mf_plots <- plot_enrichresult(test_cp[["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.
## ! # 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_plots <- plot_enrichresult(test_cp[["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.

cc_plots <- plot_enrichresult(test_cp[["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.

xl_gsea_plots <- plot_topn_gsea(test_cp[["go_data"]][["GO_gse"]])

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.0), 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), 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))