Wellcome to Lina’s RNASeq analyses: 202209

library(hpgltools)

## 
## Attaching package: 'hpgltools'

## The following object is masked from 'package:grDevices':
## 
##     colors

## The following object is masked from 'package:base':
## 
##     tempfile

tt <- sm(devtools::load_all("~/hpgltools"))
knitr::opts_knit$set(progress = TRUE, verbose = TRUE, width = 90,
                     echo = TRUE)
knitr::opts_chunk$set(error = TRUE, fig.width = 8, fig.height = 8,
                      dpi = 96)
old_options <- options(stringsAsFactors = FALSE,
                       knitr.duplicate.label = "allow")
ggplot2::theme_set(ggplot2::theme_bw(base_size = 12))
ver <- "202209"
rundate <- format(Sys.Date(), format = "%Y%m%d")

1 Examining some Wellcome Trust data with Lina

In the moment of writing this text, my brain cannot recall which set of samples these are. I remember primarily that Najib wished to use these to learn more about how I handle the data and through it, to learn some of the methods himself. In terms of the experimental goals and details, that is precious little. I presume though, that the sample sheet will teach me some useful and important details.

These are a group of Cure/Fail Biopsy samples across three visits, some used an antimonial drug while others used miltefosine. All samples except one were done on the same day, but sequenced across a few iterations with a relatively wide range of sequencing depths. There are a couple of samples which have not had their read numbers/mapping rates collected into the sample sheet.

TMRC30126 and TMRC30129

I see that I ran mappings against panamensis, but those do not appear to be in the sample sheet. Do we wish to look at the parasites?

samplesheet <- "sample_sheets/wellcome_trust_samples-v202209.xlsx"

2 Annotation

We take the annotation data from ensembl’s biomart instance. The genome which was used to map the data was hg38 revision 100. My default when using biomart is to load the data from 1 year before the current date.

hs_annot <- load_biomart_annotations(year = "2019")[["annotation"]]

## The biomart annotations file already exists, loading from it.

hs_annot[["transcript"]] <- paste0(rownames(hs_annot), ".", hs_annot[["version"]])
rownames(hs_annot) <- make.names(hs_annot[["ensembl_gene_id"]], unique = TRUE)
tx_gene_map <- hs_annot[, c("transcript", "ensembl_gene_id")]

summary(hs_annot)

##  ensembl_transcript_id  ensembl_gene_id      version      transcript_version    hgnc_symbol        description        gene_biotype   
##  Length   :226754      Length   :226754   Min.   : 1.00   Min.   : 1.000     Length   :226754   Length   :226754   Length   :226754  
##  N.unique :226721      N.unique : 60558   1st Qu.: 6.00   1st Qu.: 1.000     N.unique : 37539   N.unique : 44949   N.unique :    40  
##  N.blank  :     0      N.blank  :     0   Median :12.00   Median : 1.000     N.blank  : 35713   N.blank  :  1871   N.blank  :     0  
##  Min.nchar:    15      Min.nchar:    15   Mean   :10.53   Mean   : 3.074     Min.nchar:     0   Min.nchar:     0   Min.nchar:     3  
##  Max.nchar:    15      Max.nchar:    15   3rd Qu.:15.00   3rd Qu.: 5.000     Max.nchar:    25   Max.nchar:   169   Max.nchar:    34  
##                                           Max.   :26.00   Max.   :17.000                                                             
##                                                                                                                                      
##    cds_length      chromosome_name         strand       start_position       end_position           transcript    
##  Min.   :     3   Length   :226754   Length   :226754   Min.   :      577   Min.   :      647   Length   :226754  
##  1st Qu.:   357   N.unique :    25   N.unique :     2   1st Qu.: 31114489   1st Qu.: 31167025   N.unique :226754  
##  Median :   690   N.blank  :     0   N.blank  :     0   Median : 60378485   Median : 60547984   N.blank  :     0  
##  Mean   :  1134   Min.nchar:     1   Min.nchar:     1   Mean   : 74118835   Mean   : 74206525   Min.nchar:    17  
##  3rd Qu.:  1440   Max.nchar:     2   Max.nchar:     1   3rd Qu.:109375466   3rd Qu.:109441171   Max.nchar:    20  
##  Max.   :107976                                         Max.   :248936581   Max.   :248937043                     
##  NAs    :126871

hs_go <- load_biomart_go()[["go"]]

## The biomart annotations file already exists, loading from it.

hs_length <- hs_annot[, c("ensembl_gene_id", "cds_length")]
colnames(hs_length) <- c("ID", "length")

wellcome_outtime <- create_se(metadata = samplesheet, gene_info = hs_annot,
                             file_column = "hg38_100hisatfile") |>
  set_batches(fact = "drug") |>
  subset_se(nonzero = 15000)

## Reading the sample metadata.

## Did not find the column: sampleid.

## Setting the ID column to the first column.

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

## The sample definitions comprises: 36 rows(samples) and 47 columns(metadata fields).

## Dropping 5 htseq summary rows.

## Matched 21433 annotations and counts out of 21481 count rows and 226754 annotations.

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

## The final summarized experiment has 21481 rows and 47 columns.

## The number of samples by batch are:

## 
## Antimoniate Miltefosine 
##          18          18

## The samples (and read coverage) removed when filtering 15000 non-zero genes are:

## TMRC30129 TMRC30112 
##  33030094   2410677 
## TMRC30129 TMRC30112 
##     14469     14215

## Samples removed: 14469, 14215

outcome_time <- paste0(colData(wellcome_outtime)[["clinicaloutcome"]], "_v",
                       colData(wellcome_outtime)[["visitnumber"]])
wellcome_outtime <- set_conditions(wellcome_outtime, fact = outcome_time)

## The numbers of samples by condition are:

## 
##    Cure_v1    Cure_v2    Cure_v3 Failure_v1 Failure_v2 Failure_v3 
##          5          5          5          6          7          6

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

wellcome_time <- set_conditions(wellcome_outtime, fact = "visitnumber")

## The numbers of samples by condition are:

## 
##  1  2  3 
## 11 12 11

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

colData(wellcome_time)[["condition"]] <- paste0("t", colData(wellcome_time)[["condition"]])

wellcome_outcome <- set_conditions(wellcome_outtime, fact = "clinicaloutcome") |>
  set_batches(fact = "visitnumber")

## The numbers of samples by condition are:

## 
##    Cure Failure 
##      15      19

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

## 
##  1  2  3 
## 11 12 11

colData(wellcome_outcome)[["batch"]] <- paste0("t", colData(wellcome_outcome)[["batch"]])

wellcome_drug <- set_conditions(wellcome_outtime, fact = "drug") |>
  set_batches(fact = "visitnumber")

## The numbers of samples by condition are:

## 
## Antimoniate Miltefosine 
##          18          16

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

## 
##  1  2  3 
## 11 12 11

colData(wellcome_drug)[["batch"]] <- paste0("t", colData(wellcome_drug)[["batch"]])

wellcome_parasite <- set_conditions(wellcome_outtime, fact = "infectingspecie") |>
  set_batches(fact = "visitnumber")

## The numbers of samples by condition are:

## 
## L. V. braziliensis   L. V. panamensis 
##                 11                 23

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

## 
##  1  2  3 
## 11 12 11

colData(wellcome_outcome)[["batch"]] <- paste0("t", colData(wellcome_outcome)[["batch"]])

3 Wellcome plots

3.1 A few global metrics

plot_libsize(wellcome_time)

## Library sizes of 34 samples, 
## ranging from 2,763,416 to 220,494,183.

plot_nonzero(wellcome_time)

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the hpgltools package.
##   Please report the issue to the authors.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

## A non-zero genes plot of 34 samples.
## These samples have an average 47.15248 CPM coverage and 16779 genes observed, ranging from 15293 to
## 17990.

plot_boxplot(wellcome_time)

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

3.2 By visit

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

## Removing 7361 low-count genes (14120 remaining).

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

plot_pca(wellcome_time_norm, plot_labels = FALSE)

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

wellcome_time_varpart <- simple_varpart(wellcome_time_norm)

## The model of ~ condition + batch has 4 levels and rank 4

wellcome_time_varpart$partition_plot

wellcome_time_nb <- normalize(wellcome_time, norm = "quant", convert = "cpm",
                                   batch = "svaseq", filter = TRUE, transform = "log2")

## Warning in normalize_se(object, ...): Quantile normalization and sva do not always play well together.

## Removing 7361 low-count genes (14120 remaining).

## transform_counts: Found 569 values less than 0.

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

plot_pca(wellcome_time_nb, plot_labels = FALSE)

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

wellcome_time_nb_varpart <- simple_varpart(wellcome_time_nb)

## The model of ~ condition + batch has 4 levels and rank 4

wellcome_time_nb_varpart$partition_plot

3.3 Cure/Fail

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

## Removing 7361 low-count genes (14120 remaining).

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

plot_pca(wellcome_outcome_norm, plot_labels = FALSE)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by Cure, Failure
## Shapes are defined by tt1, tt2, tt3.

wellcome_outcome_varpart <- simple_varpart(wellcome_outcome_norm)

## The model of ~ condition + batch has 4 levels and rank 4

wellcome_outcome_varpart$partition_plot

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

## Removing 7361 low-count genes (14120 remaining).

## transform_counts: Found 819 values less than 0.

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

plot_pca(wellcome_outcome_nb, plot_labels = FALSE)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by Cure, Failure
## Shapes are defined by tt1, tt2, tt3.

wellcome_outcome_nb_varpart <- simple_varpart(wellcome_outcome_nb)

## The model of ~ condition + batch has 4 levels and rank 4

wellcome_outcome_nb_varpart$partition_plot

3.4 Drug

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

## Removing 7361 low-count genes (14120 remaining).

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

plot_pca(wellcome_drug_norm, plot_labels = FALSE)

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

wellcome_drug_varpart <- simple_varpart(wellcome_drug_norm)

## The model of ~ condition + batch has 4 levels and rank 4

wellcome_drug_varpart$partition_plot

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

## Removing 7361 low-count genes (14120 remaining).

## transform_counts: Found 868 values less than 0.

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

plot_pca(wellcome_drug_nb, plot_labels = FALSE)

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

wellcome_drug_nb_varpart <- simple_varpart(wellcome_drug_nb)

## The model of ~ condition + batch has 4 levels and rank 4

wellcome_drug_nb_varpart$partition_plot

3.5 Outcome and time

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

## Removing 7361 low-count genes (14120 remaining).

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

plot_pca(wellcome_outtime_norm, plot_labels = FALSE)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by Cure_v1, Cure_v2, Cure_v3, Failure_v1, Failure_v2, Failure_v3
## Shapes are defined by Antimoniate, Miltefosine.

wellcome_outtime_varpart <- simple_varpart(wellcome_outtime_norm)

## The model of ~ condition + batch has 7 levels and rank 7

wellcome_outtime_varpart$partition_plot

wellcome_outtime_nb <- normalize(wellcome_outtime, norm = "quant", convert = "cpm",
                                   batch = "svaseq", filter = TRUE, transform = "log2")

## Warning in normalize_se(object, ...): Quantile normalization and sva do not always play well together.

## Removing 7361 low-count genes (14120 remaining).

## transform_counts: Found 673 values less than 0.

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

plot_pca(wellcome_outtime_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 Cure_v1, Cure_v2, Cure_v3, Failure_v1, Failure_v2, Failure_v3
## Shapes are defined by Antimoniate, Miltefosine.

wellcome_outtime_nb_varpart <- simple_varpart(wellcome_outtime_nb)

## The model of ~ condition + batch has 7 levels and rank 7

wellcome_outtime_nb_varpart$partition_plot

3.6 Parasite Species

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

## Removing 7361 low-count genes (14120 remaining).

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

plot_pca(wellcome_parasite_norm, plot_labels=FALSE)$plot

wellcome_parasite_varpart <- simple_varpart(wellcome_parasite_norm)

## The model of ~ condition + batch has 4 levels and rank 4

wellcome_parasite_varpart$partition_plot

wellcome_parasite_nb <- normalize(wellcome_parasite, norm = "quant", convert = "cpm",
                                   batch = "svaseq", filter = TRUE, transform = "log2")

## Warning in normalize_se(object, ...): Quantile normalization and sva do not always play well together.

## Removing 7361 low-count genes (14120 remaining).

## transform_counts: Found 649 values less than 0.

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

plot_pca(wellcome_parasite_nb)$plot

wellcome_parasite_nb_varpart <- simple_varpart(wellcome_parasite_nb)

## The model of ~ condition + batch has 4 levels and rank 4

wellcome_parasite_nb_varpart$partition_plot

4 DE analyses

4.1 Outcome and time

TODO: Also do this without sva.

outtime_de <- all_pairwise(wellcome_outtime, model_fstring = "~ 0 + condition", filter = TRUE)

##    Cure_v1    Cure_v2    Cure_v3 Failure_v1 Failure_v2 Failure_v3 
##          5          5          5          6          7          6

## Removing 7361 low-count genes (14120 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 23321 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.
## This contrast put the denominator first.
## This contrast put the denominator first.
## This contrast put the denominator first.
## This contrast put the denominator first.
## This contrast put the denominator first.
## This contrast put the denominator first.
## This contrast put the denominator first.
## This contrast put the denominator first.
## This contrast put the denominator first.
## This contrast put the denominator first.
## This contrast put the denominator first.
## This contrast put the denominator first.

## conditions
##    Cure_v1    Cure_v2    Cure_v3 Failure_v1 Failure_v2 Failure_v3 
##          5          5          5          6          7          6

## calcNormFactors has been renamed to normLibSizes

## conditions
##    Cure_v1    Cure_v2    Cure_v3 Failure_v1 Failure_v2 Failure_v3 
##          5          5          5          6          7          6

## conditions
##    Cure_v1    Cure_v2    Cure_v3 Failure_v1 Failure_v2 Failure_v3 
##          5          5          5          6          7          6

outtime_keepers <- list(
    "curev1v2" = c("Curev2", "Curev1"),
    "curev1v3" = c("Curev3", "Curev1"),
    "curev2v3" = c("Curev3", "Curev2"),
    "failv1v2" = c("Failurev2", "Failurev1"),
    "failv1v3" = c("Failurev3", "Failurev1"),
    "failv2v3" = c("Failurev3", "Failurev2"),
    "cfv1" = c("Failurev1", "Curev1"),
    "cfv2" = c("Failurev2", "Curev2"),
    "cfv3" = c("Failurev3", "Curev3"))
outtime_table <- combine_de_tables(
    outtime_de,
    keepers = outtime_keepers,
    excel = glue::glue("excel/wellcome_outtime_table_sva-v{ver}.xlsx"))

## Deleting the file excel/wellcome_outtime_table_sva-v202209.xlsx before writing the tables.

## The keepers has no elements in the coefficients.

## Here are the keepers: Curev2, Curev1, Curev3, Curev1, Curev3, Curev2, Failurev2, Failurev1, Failurev3, Failurev1, Failurev3, Failurev2, Failurev1, Curev1, Failurev2, Curev2, Failurev3, Curev3

## Here are the coefficients: Failure_v3, Failure_v2, Failure_v3, Failure_v1, Failure_v2, Failure_v1, Failure_v3, Cure_v3, Failure_v2, Cure_v3, Failure_v1, Cure_v3, Failure_v3, Cure_v2, Failure_v2, Cure_v2, Failure_v1, Cure_v2, Cure_v3, Cure_v2, Failure_v3, Cure_v1, Failure_v2, Cure_v1, Failure_v1, Cure_v1, Cure_v3, Cure_v1, Cure_v2, Cure_v1

## Error in `extract_keepers()`:
## ! Unable to find the set of contrasts to keep, fix this and try again.

outtime_sig <- extract_significant_genes(
    outtime_table,
    excel = glue::glue("excel/wellcome_outtime_sig_sva-v{ver}.xlsx"))

## Deleting the file excel/wellcome_outtime_sig_sva-v202209.xlsx before writing the tables.

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

outtime_batch_de <- all_pairwise(wellcome_outtime, filter = TRUE)

##    Cure_v1    Cure_v2    Cure_v3 Failure_v1 Failure_v2 Failure_v3 
##          5          5          5          6          7          6 
## Antimoniate Miltefosine 
##          18          16

## Removing 7361 low-count genes (14120 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 23321 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.
## This contrast put the denominator first.
## This contrast put the denominator first.
## This contrast put the denominator first.
## This contrast put the denominator first.
## This contrast put the denominator first.
## This contrast put the denominator first.
## This contrast put the denominator first.
## This contrast put the denominator first.
## This contrast put the denominator first.
## This contrast put the denominator first.
## This contrast put the denominator first.
## This contrast put the denominator first.

## conditions
##    Cure_v1    Cure_v2    Cure_v3 Failure_v1 Failure_v2 Failure_v3 
##          5          5          5          6          7          6

## calcNormFactors has been renamed to normLibSizes

## conditions
##    Cure_v1    Cure_v2    Cure_v3 Failure_v1 Failure_v2 Failure_v3 
##          5          5          5          6          7          6

## conditions
##    Cure_v1    Cure_v2    Cure_v3 Failure_v1 Failure_v2 Failure_v3 
##          5          5          5          6          7          6

outtime_batch_table <- combine_de_tables(
    outtime_batch_de,
    keepers = outtime_keepers,
    excel = glue::glue("excel/wellcome_outtime_table_batch-v{ver}.xlsx"))

## Deleting the file excel/wellcome_outtime_table_batch-v202209.xlsx before writing the tables.

## The keepers has no elements in the coefficients.

## Here are the keepers: Curev2, Curev1, Curev3, Curev1, Curev3, Curev2, Failurev2, Failurev1, Failurev3, Failurev1, Failurev3, Failurev2, Failurev1, Curev1, Failurev2, Curev2, Failurev3, Curev3

## Here are the coefficients: Failure_v3, Failure_v2, Failure_v3, Failure_v1, Failure_v2, Failure_v1, Failure_v3, Cure_v3, Failure_v2, Cure_v3, Failure_v1, Cure_v3, Failure_v3, Cure_v2, Failure_v2, Cure_v2, Failure_v1, Cure_v2, Cure_v3, Cure_v2, Failure_v3, Cure_v1, Failure_v2, Cure_v1, Failure_v1, Cure_v1, Cure_v3, Cure_v1, Cure_v2, Cure_v1

## Error in `extract_keepers()`:
## ! Unable to find the set of contrasts to keep, fix this and try again.

outtime_sig <- extract_significant_genes(
    outtime_batch_table,
    excel = glue::glue("excel/wellcome_outtime_sig_batch-v{ver}.xlsx"))

## Deleting the file excel/wellcome_outtime_sig_batch-v202209.xlsx before writing the tables.

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

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

## t1 t2 t3 
## 11 12 11

## Removing 7361 low-count genes (14120 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 23321 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
## t1 t2 t3 
## 11 12 11

## calcNormFactors has been renamed to normLibSizes

## conditions
## t1 t2 t3 
## 11 12 11

## conditions
## t1 t2 t3 
## 11 12 11

time_keepers <- list(
    "v1v2" = c("t2", "t1"),
    "v1v3" = c("t3", "t1"),
    "v2v3" = c("t3", "t2"))
time_sva_table <- combine_de_tables(
    time_de,
    keepers = time_keepers,
    excel = glue::glue("excel/wellcome_time_table_sva-v{ver}.xlsx"))

## Deleting the file excel/wellcome_time_table_sva-v202209.xlsx before writing the tables.

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

time_sva_sig <- extract_significant_genes(
    time_sva_table,
    excel = glue::glue("excel/wellcome_time_sig_sva-v{ver}.xlsx"))

## Deleting the file excel/wellcome_time_sig_sva-v202209.xlsx before writing the tables.

4.2 Outcome

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

## L. V. braziliensis   L. V. panamensis 
##                 11                 23

## Removing 7361 low-count genes (14120 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 23321 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
## LVbraziliensis   LVpanamensis 
##             11             23

## calcNormFactors has been renamed to normLibSizes

## conditions
## LVbraziliensis   LVpanamensis 
##             11             23

## conditions
## LVbraziliensis   LVpanamensis 
##             11             23

parasite_sva_table <- combine_de_tables(
    parasite_sva_de,
    excel = glue::glue("excel/wellcome_parasite_table_sva-v{ver}.xlsx"))

## Deleting the file excel/wellcome_parasite_table_sva-v202209.xlsx before writing the tables.

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

parasite_sva_sig <- extract_significant_genes(
    parasite_sva_table,
    excel = glue::glue("excel/wellcome_parasite_sig_sva-v{ver}.xlsx"))

## Deleting the file excel/wellcome_parasite_sig_sva-v202209.xlsx before writing the tables.

parasite_batch_de <- all_pairwise(wellcome_parasite, filter = TRUE)

## L. V. braziliensis   L. V. panamensis 
##                 11                 23 
##  1  2  3 
## 11 12 11

## Removing 7361 low-count genes (14120 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 23321 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
## LVbraziliensis   LVpanamensis 
##             11             23

## calcNormFactors has been renamed to normLibSizes

## conditions
## LVbraziliensis   LVpanamensis 
##             11             23

## conditions
## LVbraziliensis   LVpanamensis 
##             11             23

parasite_batch_table <- combine_de_tables(
    parasite_batch_de,
    excel = glue::glue("excel/wellcome_parasite_table_batch-v{ver}.xlsx"))

## Deleting the file excel/wellcome_parasite_table_batch-v202209.xlsx before writing the tables.

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

parasite_batch_sig <- extract_significant_genes(
    parasite_batch_table,
    excel = glue::glue("excel/wellcome_parasite_sig_sva-v{ver}.xlsx"))

## Deleting the file excel/wellcome_parasite_sig_sva-v202209.xlsx before writing the tables.

4.3 Parasite

## wellcome_outcome_norm <- normalize_expt(wellcome_outcome, norm = "quant", convert = "cpm",

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

##    Cure Failure 
##      15      19

## Removing 7361 low-count genes (14120 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 23321 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
##    Cure Failure 
##      15      19

## calcNormFactors has been renamed to normLibSizes

## conditions
##    Cure Failure 
##      15      19

## conditions
##    Cure Failure 
##      15      19

outcome_sva_table <- combine_de_tables(
    outcome_sva_de,
    excel = glue::glue("excel/wellcome_outcome_table_sva-v{ver}.xlsx"))

## Deleting the file excel/wellcome_outcome_table_sva-v202209.xlsx before writing the tables.

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

outcome_sva_sig <- extract_significant_genes(
    outcome_sva_table,
    excel = glue::glue("excel/wellcome_outcome_sig_sva-v{ver}.xlsx"))

## Deleting the file excel/wellcome_outcome_sig_sva-v202209.xlsx before writing the tables.

outcome_batch_de <- all_pairwise(wellcome_outcome, filter = TRUE)

##    Cure Failure 
##      15      19 
## tt1 tt2 tt3 
##  11  12  11

## Removing 7361 low-count genes (14120 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 23321 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
##    Cure Failure 
##      15      19

## calcNormFactors has been renamed to normLibSizes

## conditions
##    Cure Failure 
##      15      19

## conditions
##    Cure Failure 
##      15      19

outcome_batch_table <- combine_de_tables(
    outcome_batch_de,
    excel = glue::glue("excel/wellcome_outcome_table_batch-v{ver}.xlsx"))

## Deleting the file excel/wellcome_outcome_table_batch-v202209.xlsx before writing the tables.

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

outcome_batch_sig <- extract_significant_genes(
    outcome_batch_table,
    excel = glue::glue("excel/wellcome_outcome_sig_sva-v{ver}.xlsx"))

## Deleting the file excel/wellcome_outcome_sig_sva-v202209.xlsx before writing the tables.

wellcome_gsva_c2 <- simple_gsva(wellcome_outcome, signature_category = "c2")

## gsva requires the annotation field to be filled in. Setting it to orgdb given.

## Converting the rownames() of the SummarizedExperiment to ENTREZID.

## 1666 ENSEMBL ID's didn't have a matching ENTEREZ ID. Dropping them now.

## Before conversion, the dataset has 21481 entries.

## After conversion, the SummarizedExperiment has 19750 entries.

## ℹ GSVA version 2.6.14

## ℹ Searching for rows with constant values

## ! 1791 rows with constant values throughout the columns

## ℹ Calculating ssGSEA scores for 3000 gene sets

## Calculating ssGSEA scores ■■■■■                             15% | ETA:  7sCalculating ssGSEA scores ■■■■■■                            18% | ETA:  9sCalculating ssGSEA scores ■■■■■■■                           21% | ETA:  8sCalculating ssGSEA scores ■■■■■■■■                          24% | ETA:  8sCalculating ssGSEA scores ■■■■■■■■■                         26% | ETA:  7sCalculating ssGSEA scores ■■■■■■■■■■                        29% | ETA:  7sCalculating ssGSEA scores ■■■■■■■■■■■                       32% | ETA:  6sCalculating ssGSEA scores ■■■■■■■■■■■■                      35% | ETA:  6sCalculating ssGSEA scores ■■■■■■■■■■■■                      38% | ETA:  6sCalculating ssGSEA scores ■■■■■■■■■■■■■                     41% | ETA:  5sCalculating ssGSEA scores ■■■■■■■■■■■■■■                    44% | ETA:  5sCalculating ssGSEA scores ■■■■■■■■■■■■■■■                   47% | ETA:  5sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■                  50% | ETA:  4sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■                 53% | ETA:  4sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■                56% | ETA:  4sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■■               59% | ETA:  4sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■■■              62% | ETA:  3sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■■■              65% | ETA:  3sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■■■■             68% | ETA:  3sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■■■■■            71% | ETA:  3sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■■■■■■           74% | ETA:  2sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■■■■■■■          76% | ETA:  2sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■■■■■■■■         79% | ETA:  2sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■■■■■■■■■        82% | ETA:  1sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■■■■■■■■■■       85% | ETA:  1sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■■■■■■■■■■       88% | ETA:  1sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■■■■■■■■■■■      91% | ETA:  1sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■     94% | ETA:  1sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■    97% | ETA:  0s                                                                          ℹ Normalizing ssGSEA scores
## ✔ Calculations finished
## Reading the sample metadata.
## Did not find the column: sampleid.
## The rownames do not appear numeric, using them.
## Checking the state of the condition column.
## Checking the state of the batch column.
## Checking the condition factor.
## The sample definitions comprises: 34 rows(samples) and 50 columns(metadata fields).

## Warning in merge_counts_annotations(gene_info, all_count_tables, tx_gene_map): Even after changing the rownames in gene info, they do not
## match the count table.

## Even after changing the rownames in gene info, they do not match the count table.
## Here are the first few rownames from the count tables:
## 1, 2, 3, 4, 5, 6
## Here are the first few rownames from the gene information table:
## NAKAMURA_CANCER_MICROENVIRONMENT_UP, NAKAMURA_CANCER_MICROENVIRONMENT_DN, WEST_ADRENOCORTICAL_TUMOR_MARKERS_UP, WEST_ADRENOCORTICAL_TUMOR_MARKERS_DN, WINTER_HYPOXIA_UP, WINTER_HYPOXIA_DN

## Error in `merge.data.table()`:
## ! The following columns listed in `by` are missing from x: [row.names]

wellcome_gsva_c2_sig <- get_sig_gsva_categories(
    wellcome_gsva_c2,
    excel = "excel/wellcome_gsva_c2.xlsx")

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

wellcome_gsva_c7 <- simple_gsva(wellcome_outcome, signature_category = "c7")

## gsva requires the annotation field to be filled in. Setting it to orgdb given.
## Converting the rownames() of the SummarizedExperiment to ENTREZID.
## 1666 ENSEMBL ID's didn't have a matching ENTEREZ ID. Dropping them now.
## Before conversion, the dataset has 21481 entries.
## After conversion, the SummarizedExperiment has 19750 entries.
## ℹ GSVA version 2.6.14
## ℹ Searching for rows with constant values
## ! 1791 rows with constant values throughout the columns
## ℹ Calculating ssGSEA scores for 3000 gene sets
## Calculating ssGSEA scores ■■■■■                             15% | ETA:  6sCalculating ssGSEA scores ■■■■■■                            18% | ETA:  6sCalculating ssGSEA scores ■■■■■■■                           21% | ETA:  6sCalculating ssGSEA scores ■■■■■■■■                          24% | ETA:  6sCalculating ssGSEA scores ■■■■■■■■■                         26% | ETA:  6sCalculating ssGSEA scores ■■■■■■■■■■                        29% | ETA:  5sCalculating ssGSEA scores ■■■■■■■■■■■                       32% | ETA:  5sCalculating ssGSEA scores ■■■■■■■■■■■■                      35% | ETA:  5sCalculating ssGSEA scores ■■■■■■■■■■■■                      38% | ETA:  5sCalculating ssGSEA scores ■■■■■■■■■■■■■                     41% | ETA:  4sCalculating ssGSEA scores ■■■■■■■■■■■■■■                    44% | ETA:  4sCalculating ssGSEA scores ■■■■■■■■■■■■■■■                   47% | ETA:  4sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■                  50% | ETA:  4sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■                 53% | ETA:  4sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■                56% | ETA:  3sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■■               59% | ETA:  3sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■■■              62% | ETA:  3sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■■■              65% | ETA:  3sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■■■■             68% | ETA:  2sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■■■■■            71% | ETA:  2sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■■■■■■           74% | ETA:  2sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■■■■■■■          76% | ETA:  2sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■■■■■■■■         79% | ETA:  2sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■■■■■■■■■        82% | ETA:  1sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■■■■■■■■■■       85% | ETA:  1sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■■■■■■■■■■       88% | ETA:  1sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■■■■■■■■■■■      91% | ETA:  1sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■     94% | ETA:  0sCalculating ssGSEA scores ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■    97% | ETA:  0s                                                                          ℹ Normalizing ssGSEA scores
## ✔ Calculations finished
## Reading the sample metadata.
## Did not find the column: sampleid.
## The rownames do not appear numeric, using them.
## Checking the state of the condition column.
## Checking the state of the batch column.
## Checking the condition factor.
## The sample definitions comprises: 34 rows(samples) and 50 columns(metadata fields).

## Warning in merge_counts_annotations(gene_info, all_count_tables, tx_gene_map): Even after changing the rownames in gene info, they do not
## match the count table.

## Even after changing the rownames in gene info, they do not match the count table.
## Here are the first few rownames from the count tables:
## 1, 2, 3, 4, 5, 6
## Here are the first few rownames from the gene information table:
## NAKAMURA_CANCER_MICROENVIRONMENT_UP, NAKAMURA_CANCER_MICROENVIRONMENT_DN, WEST_ADRENOCORTICAL_TUMOR_MARKERS_UP, WEST_ADRENOCORTICAL_TUMOR_MARKERS_DN, WINTER_HYPOXIA_UP, WINTER_HYPOXIA_DN

## Error in `merge.data.table()`:
## ! The following columns listed in `by` are missing from x: [row.names]

wellcome_gsva_c7_sig <- get_sig_gsva_categories(
    wellcome_gsva_c7,
    excel = "excel/wellcome_gsva_c7.xlsx")

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

outtime_gprofiler <- all_gprofiler(outtime_sig)

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

outtime_gprofiler$curev1v2_up$pvalue_plots$bpp_plot_over

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

outtime_gprofiler$failv1v2_up$pvalue_plots$kegg_plot_over

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

outtime_gprofiler$cfv1_up$pvalue_plots$reactome_plot_over

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

outtime_gprofiler$cfv1_down$pvalue_plots$reactome_plot_over

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