Examining a series of Xenopus laevis early development samples.

1 TODO

1. Ensure I use the new sample_type column.
2. Ensure that the Mag samples are not included.
3. Ensure that the contrasts are performed as per the word document (as shown in the keepers variable below)
4. Make sure I perform all ontology searches and return the xlsx files so that it is obvious which direction and which contrast is being returned.

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

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

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

5 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 
## "#B7B3CA"

sampletype_colors <- get_xlsx_colors_by_column(start_sheet, column = "sample_type")
sampletype_colors

## $D11_Ctrl
##       rgb 
## "#7F9CC1" 
## 
## $D11_Exp
##       rgb 
## "#9B84BF" 
## 
## $V11_Ctrl
##       rgb 
## "#A5B592"

6 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 = "sample_type", colors = sampletype_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 88 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 88 columns.

## The numbers of samples by condition are:

## 
## D11_Ctrl  D11_Exp V11_Ctrl 
##        6        6        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 = "sample_type", colors = sampletype_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 88 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 88 columns.

## The numbers of samples by condition are:

## 
## D11_Ctrl  D11_Exp V11_Ctrl 
##        6        6        6

## The number of samples by batch are:

## 
## FACS  Mag 
##    9    9

facs_se <- subset_se(post_se, subset = "treatment=='FACS'")

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.

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

8 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_Ctrl, D11_Exp, V11_Ctrl
## 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_Ctrl, D11_Exp, V11_Ctrl
## Shapes are defined by t1, t2.

pre_treat <- set_conditions(pre_se, fact = "sample_type", colors = sampletype_colors) |>
  set_batches(fact = "biological_replicate")

## The numbers of samples by condition are:

## 
## D11_Ctrl  D11_Exp V11_Ctrl 
##        6        6        6

## The number of samples by batch are:

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

post_treat <- set_conditions(post_se, fact = "sample_type", colors = sampletype_colors) |>
  set_batches(fact = "biological_replicate")

## The numbers of samples by condition are:

## 
## D11_Ctrl  D11_Exp V11_Ctrl 
##        6        6        6

## 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_Ctrl, D11_Exp, V11_Ctrl
## 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_Ctrl, D11_Exp, V11_Ctrl
## Shapes are defined by R1, R2, R3.

8.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_Ctrl, D11_Exp, V11_Ctrl
## 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_Ctrl, D11_Exp, V11_Ctrl
## Shapes are defined by R1, R2, R3.

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

## calcNormFactors has been renamed to normLibSizes

## 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_Ctrl, D11_Exp, V11_Ctrl
## Shapes are defined by t1, t2.

only_facs_de <- all_pairwise(only_facs)

## D11_Ctrl  D11_Exp V11_Ctrl 
##        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_Ctrl  D11_Exp V11_Ctrl 
##        3        3        3

## calcNormFactors has been renamed to normLibSizes

## conditions
## D11_Ctrl  D11_Exp V11_Ctrl 
##        3        3        3

## conditions
## D11_Ctrl  D11_Exp V11_Ctrl 
##        3        3        3

keepers <- list(
  "D11_Exp_vs_D11_Ctrl" = c("D11_Exp", "D11_Ctrl"),
  "V11_Ctrl_vs_D11_Ctrl" = c("V11_Ctrl", "D11_Ctrl"),
  "D11_Exp_vs_V11_Ctrl" = c("D11_Exp", "V11_Ctrl"))
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-v20260929.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-v20260929.xlsx before writing the tables.

10 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(only_facs)

## The factor D11_Ctrl has 3 rows.

## The factor D11_Exp has 3 rows.

## The factor V11_Ctrl has 3 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_Ctrl", "mean_d11")])[["scatter"]])

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

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

11 Compare the three states

I am going to stop showing pre-deduplication.

conditions(post_treat)

##   D11_exp_FACS D11_exp_FACS_2 D11_exp_FACS_3 D11_exp_MagCar  D11_exp_Mag_2  D11_exp_Mag_3       D11_FACS     D11_FACS_2     D11_FACS_3 
##        D11_Exp        D11_Exp        D11_Exp        D11_Exp        D11_Exp        D11_Exp       D11_Ctrl       D11_Ctrl       D11_Ctrl 
##     D11_MagCar      D11_Mag_2      D11_Mag_3       V11_FACS     V11_FACS_2     V11_FACS_3     V11_MagCar      V11_Mag_2      V11_Mag_3 
##       D11_Ctrl       D11_Ctrl       D11_Ctrl       V11_Ctrl       V11_Ctrl       V11_Ctrl       V11_Ctrl       V11_Ctrl       V11_Ctrl 
## Levels: D11_Ctrl D11_Exp V11_Ctrl

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

## D11_Ctrl  D11_Exp V11_Ctrl 
##        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_Ctrl  D11_Exp V11_Ctrl 
##        6        6        6

## calcNormFactors has been renamed to normLibSizes

## conditions
## D11_Ctrl  D11_Exp V11_Ctrl 
##        6        6        6

## conditions
## D11_Ctrl  D11_Exp V11_Ctrl 
##        6        6        6

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

## D11_Ctrl  D11_Exp V11_Ctrl 
##        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_Ctrl  D11_Exp V11_Ctrl 
##        6        6        6

## calcNormFactors has been renamed to normLibSizes

## conditions
## D11_Ctrl  D11_Exp V11_Ctrl 
##        6        6        6

## conditions
## D11_Ctrl  D11_Exp V11_Ctrl 
##        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.

11.1 Create tables!

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

## 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_Exp_vs_D11_Ctrl        2131          1231        2122          1243        1368          1496
## 2         V11_Ctrl_vs_D11_Ctrl        3244          2445        3106          2568        2641          2491
## 3 V11_Ctrl_vs_D11_Exp-inverted         981           185        1014           145         603           179

## 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, keepers = keepers,
  excel = glue("excel/treatment_comparisons_sva_table-v{ver}.xlsx"))

## 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_Exp_vs_D11_Ctrl        2347          1586        2340          1653        1858          2040
## 2         V11_Ctrl_vs_D11_Ctrl        2985          2614        2871          2729        2597          2619
## 3 V11_Ctrl_vs_D11_Exp-inverted         942           123        1043           130         519           110

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

11.2 Significant genes

sig_nosva <- extract_significant_genes(
  table_sva, excel = glue("excel/treatment_comparisons_nosva_sig-v{ver}.xlsx"))
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_Exp_vs_D11_Ctrl      2347       1586
## V11_Ctrl_vs_D11_Ctrl     2985       2614
## D11_Exp_vs_V11_Ctrl       942        123

sig_sva <- extract_significant_genes(
  table_sva, excel = glue("excel/treatment_comparisons_sva_sig-v{ver}.xlsx"))
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_Exp_vs_D11_Ctrl      2347       1586
## V11_Ctrl_vs_D11_Ctrl     2985       2614
## D11_Exp_vs_V11_Ctrl       942        123

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

11.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",
                                         excel = "excel/up_cp_v11.xlsx")

## Error in `loadNamespace()` at hpgltools/R/ontology_clusterprofiler.R:962:5:
## ! there is no package called 'org.Xl.eg.db'

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

## Error:
## ! object 'test_up_cp_v11' not found

mf_up_v11_plots[["vol"]]

## Error:
## ! object 'mf_up_v11_plots' not found

mf_up_v11_plots[["tree"]]

## Error:
## ! object 'mf_up_v11_plots' not found

mf_up_v11_plots[["dot"]]

## Error:
## ! object 'mf_up_v11_plots' not found

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

## Error:
## ! object 'test_up_cp_v11' not found

bp_up_v11_plots[["vol"]]

## Error:
## ! object 'bp_up_v11_plots' not found

bp_up_v11_plots[["tree"]]

## Error:
## ! object 'bp_up_v11_plots' not found

bp_up_v11_plots[["dot"]]

## Error:
## ! object 'bp_up_v11_plots' not found

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

## Error:
## ! object 'test_up_cp_v11' not found

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

## Error in `h()`:
## ! error in evaluating the argument 'gse' in selecting a method for function 'plot_topn_gsea': object 'test_up_cp_v11' not found

xl_gsea_v11_plots[[1]]

## Error:
## ! object 'xl_gsea_v11_plots' not found

xl_gsea_v11_plots[[2]]

## Error:
## ! object 'xl_gsea_v11_plots' not found

11.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",
                                           excel = "excel/up_cp_d11.xlsx")

## Error in `loadNamespace()` at hpgltools/R/ontology_clusterprofiler.R:962:5:
## ! there is no package called 'org.Xl.eg.db'

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

## Error:
## ! object 'test_down_cp_v11' not found

mf_down_v11_plots[["vol"]]

## Error:
## ! object 'mf_down_v11_plots' not found

mf_down_v11_plots[["tree"]]

## Error:
## ! object 'mf_down_v11_plots' not found

mf_down_v11_plots[["dot"]]

## Error:
## ! object 'mf_down_v11_plots' not found

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

## Error:
## ! object 'test_down_cp_v11' not found

bp_down_v11_plots[["vol"]]

## Error:
## ! object 'bp_down_v11_plots' not found

bp_down_v11_plots[["tree"]]

## Error:
## ! object 'bp_down_v11_plots' not found

bp_down_v11_plots[["dot"]]

## Error:
## ! object 'bp_down_v11_plots' not found

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

## Error:
## ! object 'test_down_cp_v11' not found

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

## Error in `h()`:
## ! error in evaluating the argument 'gse' in selecting a method for function 'plot_topn_gsea': object 'test_down_cp_v11' not found

xl_gsea_down_v11_plots[[1]]

## Error:
## ! object 'xl_gsea_down_v11_plots' not found

xl_gsea_down_v11_plots[[2]]

## Error:
## ! object 'xl_gsea_down_v11_plots' not found

pander::pander(sessionInfo())

R version 4.6.1 (2026-06-24)

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: edgeR(v.4.10.5), 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.32.0), 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.6.1), 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.74.0), cli(v.3.6.6), Biobase(v.2.72.0), scatterpie(v.0.2.6), EBSeq(v.2.10.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.94.0), Rsamtools(v.2.28.0), systemfonts(v.1.3.2), yulab.utils(v.0.2.5), gson(v.0.2.1), R.utils(v.2.13.0), DOSE(v.4.6.0), dichromat(v.2.0-1), sessioninfo(v.1.2.4), limma(v.3.68.5), rstudioapi(v.0.19.0), RSQLite(v.3.53.3), BiocIO(v.1.22.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.23.1), Matrix(v.1.7-6), S4Vectors(v.0.50.3), abind(v.1.4-8), R.methodsS3(v.1.8.2), lifecycle(v.1.0.5), yaml(v.2.3.12), SummarizedExperiment(v.1.42.0), gplots(v.3.3.0), qvalue(v.2.44.0), SparseArray(v.1.12.3), grid(v.4.6.1), blob(v.1.3.0), promises(v.1.5.0), crayon(v.1.5.3), ggtangle(v.0.1.3), lattice(v.0.23-1), cigarillo(v.1.2.1), GenomicFeatures(v.1.64.0), annotate(v.1.90.0), KEGGREST(v.1.52.2), pillar(v.1.11.1), knitr(v.1.52), varhandle(v.2.0.6), GenomicRanges(v.1.64.0), rjson(v.0.2.23), boot(v.1.3-32), corpcor(v.1.6.10), codetools(v.0.2-20), 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.36.1), Rdpack(v.2.6.6), gtable(v.0.3.6), cachem(v.1.1.0), openxlsx(v.4.2.9), xfun(v.0.61), rbibutils(v.2.4.1), S4Arrays(v.1.12.1), mime(v.0.13), RcppEigen(v.0.3.4.0.2), Seqinfo(v.1.2.0), reformulas(v.0.4.4), survival(v.3.8-12), aisdk(v.1.4.12), NOISeq(v.2.56.0), iterators(v.1.0.14), statmod(v.1.5.2), ellipsis(v.0.3.3), nlme(v.3.1-171), pbkrtest(v.0.5.5), ggtree(v.4.2.0), usethis(v.3.2.2), 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-27), otel(v.0.2.0), BiocGenerics(v.0.58.1), DBI(v.1.3.0), DESeq2(v.1.52.0), tidyselect(v.1.2.1), processx(v.3.9.0), bit(v.4.6.0), compiler(v.4.6.1), curl(v.8.0.0), graph(v.1.90.0), httr2(v.1.3.0), desc(v.1.4.3), fontBitstreamVera(v.0.1.1), DelayedArray(v.0.38.2), plotly(v.4.12.1), rtracklayer(v.1.72.0), scales(v.1.4.0), caTools(v.1.18.4), DEoptimR(v.1.2-2), remaCor(v.0.0.20), RBGL(v.1.88.0), callr(v.3.8.0), rappdirs(v.0.3.4), stringr(v.1.6.0), digest(v.0.6.39), minqa(v.1.2.8), variancePartition(v.1.42.0), rmarkdown(v.2.32), aod(v.1.3.3), XVector(v.0.52.0), RhpcBLASctl(v.0.23-42), htmltools(v.0.5.9), pkgconfig(v.2.0.3), lme4(v.2.0-6), MatrixGenerics(v.1.24.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.46.0), R.oo(v.1.27.1), GOSemSim(v.2.38.3), 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.6.1), ggrepel(v.0.9.8), Biostrings(v.2.80.2), splines(v.4.6.1), pander(v.0.6.6), hms(v.1.1.4), locfit(v.1.5-9.12), ps(v.1.9.3), igraph(v.2.3.3), enrichit(v.0.2.5), reshape2(v.1.4.5), restez(v.2.1.5), stats4(v.4.6.1), pkgload(v.1.5.3), XML(v.3.99-0.25), evaluate(v.1.0.5), BiocManager(v.1.30.27), tzdb(v.0.5.0), nloptr(v.2.2.1), PROPER(v.1.44.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.20.0), aplot(v.0.3.2), GenomicAlignments(v.1.48.0), memoise(v.2.0.1), AnnotationDbi(v.1.74.0), IRanges(v.2.46.0), cluster(v.2.1.8.3), sva(v.3.60.0) and GSEABase(v.1.74.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))