This document performs a series of diagnostics and differential expression analyses which eventually helped inform the text of Lina Fernanda Giraldo Parra’s and Maria Adelaida Gomez’s paper. In it, they sought to use biopsies over time from people who were observed to cure/fail treatment of cutaneous leishmaniasis. Some of these people were treated with the antimonial ‘glucantime’, others with miltefosine. Biopsies were taken at the beginning, middle, and end of treatment. RNA was taken from the biopsies, made into libraries, and sequenced. The samples were mapped/counted against ensembl hg38 release 100 and the TritrypDB Leishmania panamensis MHOM/COL strain’s genome with hisat2 and salmon. There were very few reads observed from the parasite. As a result (at the time of this writing: 202310), those count tables are not included in the sample sheet (they are now, along with the salmon quantitations against hg38_100), and are not considered in this worksheet.
I collected the available metadata into the sheet ‘all_samples.xlsx’.
The function ‘gather_preprocessing_metadata()’ scans the working directory in which the various preprocessing tasks were performed, trimming, fastqc/fastp, mapping/quantitation, variant searching, etc. It is aware of my default output filenames for stdout/stderr/etc and extracts from them various metrics potentially of interest. I think the sample sheet I copied into this container is the result of the following block, so I am not going to evaluate it here.
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.
## Using mart: ENSEMBL_MART_ENSEMBL from host: Apr2019.archive.ensembl.org.
## Successfully connected to the hsapiens_gene_ensembl database.
## Finished downloading ensembl gene annotations.
## Finished downloading ensembl structure annotations.
## symbol columns is null, pattern matching 'symbol'.
## Including symbols, there are 65071 vs the 229428 gene annotations.
## Dropping haplotype chromosome annotations, set drop_haplotypes = FALSE if this is bad.
## Saving annotations to hsapiens_biomart_annotations.rda.
## Finished save().
## 12 annotation types for 208,460 genes/transcripts downloaded from Apr2019.archive.ensembl.org.
hs_annot <- hs_annot[["annotation"]]
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")]Downloading the Ontology annotation information has a significant chance of failing because it sometimes takes longer than curl’s timeout of 300 seconds. As a result, I may decide to include the saved rda of ontology results in the container, or maybe improve my ontology downloader so that it is no longer susceptible to timeouts…
Actually, I do not think I am using clusterProfiler/gostats/goseq/topgo in this document, perhaps therefore I should just exclude this block? Yeah, it seems like there is a 50/50 chance that this will fail when run within the container, and I am currently not using this set of gene:GO mappings in this document.
The primary datastructure used throughout this document is the expressionSet or summarized Experiment. This document uses the former; but both are created by reading the sample sheet, extracting the filenames of the count data from it, and combining the metadata, annotation data, and counts. My implementation of this process also includes a few extra pieces of information. For example: the colors of various potential conditions in the data.
Before we create the datastructures, I want to have the various color schemes in place. Many of these color choices are taken directly from another project which also looks at host responses.
color_choices <- list(
"cf_visit" = list(
"cure_v1" = "#D95F0E",
"failure_v1" = "#FEC44F",
"cure_v2" = "#DD1C77",
"failure_v2" = "#C994C7",
"cure_v3" = "#3182BD",
"failure_v3" = "#9ECAE1"),
"drug_visit" = list(
"antimoniate_v1" = "#badeef",
"miltefosine_v1" = "#e4a6c1",
"antimoniate_v2" = "#7bc6ea",
"miltefosine_v2" = "#d26895",
"antimoniate_v3" = "#1d91ca",
"miltefosine_v3" = "#d52e75"),
"species_visit" = list(
"lvpanamensis_v1" = "#ffe798",
"lvbraziliensis_v1" = "#b5b5b5",
"lvpanamensis_v2" = "#FFC300",
"lvbraziliensis_v2" = "#888888",
"lvpanamensis_v3" = "#b58b00",
"lvbraziliensis_v3" = "#525252"),
"cf" = list(
"cure" = "#998EC3",
"failure" = "#F1A340"),
"visit" = list(
"v1" = "#33EE33",
"v2" = "#11AA11",
"v3" = "#134413"),
"visitbi" = list(
"v1" = "#BB0000",
"vother" = "#0000BB"),
"species" = list(
"lvpanamensis" = "#FFC300",
"lvbraziliensis" = "#525252"),
"donor" = list(
"d2008" = "#B78415",
"d1029" = "#93752C",
"d1036" = "#7E6EA2",
"d1037" = "#B3499C",
"d1031" = "#BD6332",
"d2002" = "#7D8F31",
"d2009" = "#8E7037",
"d2010" = "#666666",
"d1019" = "#1B9E77",
"d2004" = "#E0A604",
"d2001" = "#CF3F76",
"d2003" = "#A0A811"),
"drug" = list(
"antimoniate" = "#3182AA",
"miltefosine" = "#C994AA"))The set of color choices demonstrates the complexity of the experimental design. We are looking at combinations of time, outcome, species, and drug. Inherent in all of these is an overarching donor effect.
The initial datastructure will use a factor combining the clinical outcome and visit as the experimental condition of interest. We will follow this up by splitting off the various factors of interest.
An aside: when testing the following block within the container, it did not print out the filenames as it was reading them, usually it does. I assume this is not a problem since the data looks fine, but this may warrant further exploration.
The call to sanitize_expt_pData() tries to ensure that Excel does not mess things up with random spaces/capitalization in the metadata.
wellcome_outtime <- create_expt(metadata = samplesheet, gene_info = hs_annot,
file_column = "hg38100hisatfile") %>%
set_expt_batches(fact = "drug") %>%
sanitize_expt_pData(spaces = TRUE)## 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.
## The sample definitions comprises: 36 rows(samples) and 54 columns(metadata fields).
## Warning in file(file, "rt"): cannot open file
## '/mnt/nfs/z1/sw/singularity/CL_biopsies_Colombia/202404201249_outputs/preprocessing/TMRC30097/outputs/02hisat2_hg38_100/hg38_100_sno_gene_gene_id.count.xz':
## No such file or directory
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': error in evaluating the argument 'object' in selecting a method for function 'pData': cannot open the connection
## Error in eval(expr, envir, enclos): object 'wellcome_outtime' not found
The previous block creates the primary datastructure, everything from now on will either:
The following block therefore creates a set of factors which are the concatenation of time+x where x is one of: outcome, drug, species.
outcome_time <- paste0(pData(wellcome_outtime)[["clinicaloutcome"]], "_v",
pData(wellcome_outtime)[["visitnumber"]])## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_outtime' not found
wellcome_outtime <- set_expt_conditions(wellcome_outtime, fact = outcome_time,
colors = color_choices[["cf_visit"]])## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_outtime' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_outtime' not found
## Error in eval(expr, envir, enclos): object 'outcome_time' not found
pData(wellcome_outtime)[["drug_visit"]] <- paste0(pData(wellcome_outtime)[["drug"]], "_",
pData(wellcome_outtime)[["visitnumber"]])## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_outtime' not found
pData(wellcome_outtime)[["species_visit"]] <- paste0(pData(wellcome_outtime)[["infectingspecies"]], "_",
pData(wellcome_outtime)[["visitnumber"]])## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_outtime' not found
First, let us get a quick view of the samples in their native state. There are a few samples which are candidates for removal due to lower coverage.
Random aside: I have some nifty code which tries to autodetect when to color the text in the colored bars white or black. The purple/pink color tricks that code into thinking it is dark enough to be colored white.
## Error in eval(expr, envir, enclos): object 'wellcome_outtime' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'data' in selecting a method for function 'plot_libsize': object 'wellcome_outtime' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'data' in selecting a method for function 'plot_nonzero': object 'wellcome_outtime' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'subset_expt': object 'wellcome_outtime' not found
The following writes out the data in two files, one before and after filtering for samples with less than 15,000 observed genes.
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'exprs': object 'wellcome_outtime' not found
written <- write_expt(wellcome_filtered,
excel = glue("excel/CL_biopsies_filtered_data-v{ver}.xlsx"))## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'exprs': object 'wellcome_filtered' not found
Looking at the nonzero plot, it seems that 15k genes is a reasonable cutoff, which would remove 2 samples. So, for the moment I will split the data up into two sets, pre/post removal.
Lina asked for a CSV copy of the data as CPM.
## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'wellcome_filtered' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'exprs': object 'wellcome_cpm' not found
## Error in eval(expr, p): object 'wellcome_cpm_mtrx' not found
wellcome_scpm <- normalize_expt(wellcome_filtered, filter = TRUE, convert = "cpm", batch = "svaseq")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'wellcome_filtered' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'exprs': object 'wellcome_scpm' not found
wellcome_scpm_write <- write.csv(x = wellcome_cpm_mtrx, file = "excel/CL_biopsies_filtered_cpm_sva.csv")## Error in eval(expr, p): object 'wellcome_cpm_mtrx' not found
We have in place some initial datastructures; let us mix and match them to get a sense of the data. There are a few factors in the experiment that are likely of primary interest: the donor, time, drug used, parasite species/strain, and clinical outcome. There are too many factors for the number of samples to do a full rank model, so I am going to make copies of the expressionset and model each factor separately with the help of sva.
Start with the donor. There are 12 people, most of whom provided three samples. Two of them (d2002 and d2009) had samples which got excluded and so have only 2. I am separating the un/filtered datasets so that we may play with both if we want.
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_outtime' not found
wellcome_donor_filt <- set_expt_conditions(wellcome_filtered, fact = "donor",
colors = color_choices)## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_filtered' not found
In this case, the condition will only be visit number. We have 12 samples for each visit, except v1/v3 which are missing 1 each.
wellcome_time <- set_expt_conditions(wellcome_outtime, fact = "visitnumber",
colors = color_choices[["visit"]])## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_outtime' not found
wellcome_time_filt <- set_expt_conditions(wellcome_filtered, fact = "visitnumber",
colors = color_choices[["visit"]])## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_filtered' not found
Conversely, we can split by clinical outcome. The unfiltered dataset has 15 cures and 21 failures; both filtered samples were failure.
wellcome_outcome <- set_expt_conditions(wellcome_outtime, fact = "clinicaloutcome",
colors = color_choices[["cf"]]) %>%
set_expt_batches(fact = "visitnumber")## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_outtime' not found
wellcome_outcome_filt <- set_expt_conditions(wellcome_filtered, fact = "clinicaloutcome",
colors = color_choices[["cf"]]) %>%
set_expt_batches(fact = "visitnumber")## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_filtered' not found
With respect to drug treatment, we started with 18 of each. The lost samples were both miltefosine.
wellcome_drug <- set_expt_conditions(wellcome_outtime, fact = "drug",
colors = color_choices[["drug"]]) %>%
set_expt_batches(fact = "visitnumber")## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_outtime' not found
wellcome_drug_filt <- set_expt_conditions(wellcome_filtered, fact = "drug",
colors = color_choices[["drug"]]) %>%
set_expt_batches(fact = "visitnumber")## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_filtered' not found
We started with 12 braziliensis and 24 panamensis; one of each was lost.
wellcome_parasite <- set_expt_conditions(wellcome_outtime, fact = "infectingspecies",
colors = color_choices[["species"]]) %>%
set_expt_batches(fact = "visitnumber")## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_outtime' not found
wellcome_parasite_filt <- set_expt_conditions(wellcome_filtered, fact = "infectingspecies",
colors = color_choices[["species"]]) %>%
set_expt_batches(fact = "visitnumber")## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_filtered' not found
Everything I wrote above is visible in the following sankey plot. It is particularly noteworthy that we have no cure braziliensis for miltefosine treatment and only 1 for each of the antominial timepoints.
drug_visit_species_cf <- plot_meta_sankey(
wellcome_filtered, factors = c("drug", "visitnumber", "clinicaloutcome", "infectingspecies"),
color_choices = color_choices)## Error in h(simpleError(msg, call)): error in evaluating the argument 'design' in selecting a method for function 'plot_meta_sankey': object 'wellcome_filtered' not found
## Error in eval(expr, envir, enclos): object 'drug_visit_species_cf' not found
Given the large number of variables in this data, lets start by getting a feeling for the amount of variance in each. Given that we don’t have a full rank with respect to clinical outcome, I assume that trying variance partition with the 4 primary factors will fail, but let us see…
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_filtered' not found
## Error in eval(expr, envir, enclos): object 'varpart_donor' not found
varpart_dvic <- simple_varpart(
wellcome_filtered,
factors = c("drug", "visitnumber", "infectingspecies", "clinicaloutcome"))## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_filtered' not found
## Error in eval(expr, envir, enclos): object 'varpart_dvic' not found
To my eyes, it appears that donor is dominant factor, followed by: visit, species, drug, and outcome. I think this observation may have an effect on the current state of the manuscript - which if I recall properly states that donor is not a significant factor in the data.
Let us ask if there are categories associated with the genes associated with each of these categories and see if they ‘make sense’.
## Error in eval(expr, envir, enclos): object 'varpart_donor' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'tail': object 'varpart_donor' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'top_300_donor_genes' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'rownames': object 'top_300_donor_genes' not found
## Error in eval(expr, envir, enclos): object 'donor_genes_gp' not found
## Error in eval(expr, envir, enclos): object 'donor_genes_gp' not found
So, the top-300 genes with respect to putative donor effect, when passed to gprofiler’s over representation analyses, look like they have lots of variance related to categories that are explicitly important to the general questions we are asking. I think that counts as: Yay!
Now let us look at the other factors in the data and see if anything noteworthy pops out. Given that the dvic datastructure comprises the other factors of interest, we will need to reorder the results with respect to each factor separately.
Visit also is dominated by variance which is explicitly what one would expect: wound healing.
## Error in eval(expr, envir, enclos): object 'varpart_dvic' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'tail': object 'varpart_dvic' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'top_300_visit_genes' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'rownames': object 'top_300_visit_genes' not found
## Error in eval(expr, envir, enclos): object 'visit_genes_gp' not found
I do not have any idea what one should expect vis a vis the two drugs. Looking at the following result, I guess one should not be surprised given that I am unaware of m/any experiments which explicitly compare antimonial treatments with respect to transcriptional profile.
## Error in eval(expr, envir, enclos): object 'varpart_dvic' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'tail': object 'varpart_dvic' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'top_300_drug_genes' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'rownames': object 'top_300_drug_genes' not found
## Error in eval(expr, envir, enclos): object 'drug_genes_gp' not found
## Error in eval(expr, envir, enclos): object 'drug_genes_gp' not found
The violin plot above suggests there is very limited variance with respect to the two species. In addition, there are significantly fewer braziliensis samples which failed treatment (for miltefosine in particular), so I presume variancePartition is going to have trouble extracting genes which are reliable in this context. In addition, I would assume that the human transcriptome has pretty similar responses to the two parasites, exacerbating this confounder.
## Error in eval(expr, envir, enclos): object 'varpart_dvic' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'tail': object 'varpart_dvic' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'top_300_species_genes' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'rownames': object 'top_300_species_genes' not found
## Error in eval(expr, envir, enclos): object 'species_genes_gp' not found
In my mind, clinical outcome is the most generally interesting query. It is also one with limited information in the data, but lots of other studies have been performed in this realm, so it is more likely to have reliable overrepresentation results even with limited information. Thus, when I look at the bar/dot plot the categories look interesting.
The following block provides an example of one of the fun things my gProfiler2 function does: it coerces the gprofiler output to the datastructure used by clusterProfiler and thus may be plotted with the various functions in the ‘enrichplot’ package. It also provides a copy of the cute interactive plots produced by gprofiler.
## Error in eval(expr, envir, enclos): object 'varpart_dvic' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'tail': object 'varpart_dvic' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'top_300_cf_genes' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'rownames': object 'top_300_cf_genes' not found
## Error in eval(expr, envir, enclos): object 'cf_genes_gp' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'dotplot': object 'cf_genes_gp' not found
Here are the PCA results before/after removing the two samples that I called shenanigans on. They are quite similar. The most notable difference is the weirdo failure visit 3 sample on the left goes away.
wellcome_outtime_norm <- normalize_expt(wellcome_outtime, filter = TRUE,
convert = "cpm", transform = "log2", norm = "quant")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'wellcome_outtime' not found
## Error in eval(expr, envir, enclos): object 'wellcome_outtime_norm' not found
## Error in eval(expr, envir, enclos): object 'pre_filter_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'pre_filter_pca' not found
wellcome_outfilt_norm <- normalize_expt(wellcome_filtered, filter = TRUE,
convert = "cpm", transform = "log2", norm = "quant")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'wellcome_filtered' not found
## Error in eval(expr, envir, enclos): object 'wellcome_outfilt_norm' not found
## Error in eval(expr, envir, enclos): object 'post_filter_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'post_filter_pca' not found
Now that we have an initial sense of how the samples relate to each other, let us poke a bit further. In each of the following, I am likely to do one plot before and one after using sva.
I think the pattern of samples by visit is pretty or satisfying in a ‘hey, these things have some internal sense’ way. I guess modifying the counts with surrogates from sva makes it a little bit clearer? It does increase the PC1 variance a bit I suppose.
wellcome_time_norm <- normalize_expt(wellcome_time_filt, norm = "quant", convert = "cpm",
filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'wellcome_time_filt' not found
## Error in eval(expr, envir, enclos): object 'wellcome_time_norm' not found
## Error in eval(expr, envir, enclos): object 'time_norm_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'time_norm_pca' not found
wellcome_time_nb <- normalize_expt(wellcome_time_filt, convert = "cpm",
batch = "svaseq", filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'wellcome_time_filt' not found
## Error in eval(expr, envir, enclos): object 'wellcome_time_nb' not found
## Error in eval(expr, envir, enclos): object 'time_nb_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'time_nb_pca' not found
In contrast, the ability of sva to determine variance which is related to clinical outcome is very limited. The resulting PCA plot is a bit of a mess.
wellcome_outcome_norm <- normalize_expt(wellcome_outcome_filt, norm = "quant", convert = "cpm",
filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'wellcome_outcome_filt' not found
## Error in eval(expr, envir, enclos): object 'wellcome_outcome_norm' not found
## Error in eval(expr, envir, enclos): object 'outcome_norm_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'outcome_norm_pca' not found
wellcome_outcome_nb <- normalize_expt(wellcome_outcome_filt, convert = "cpm",
batch = "svaseq", filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'wellcome_outcome_filt' not found
## Error in eval(expr, envir, enclos): object 'wellcome_outcome_nb' not found
## Error in eval(expr, envir, enclos): object 'outcome_nb_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'outcome_nb_pca' not found
At first glance, it at least appears possible that we can get a handle on differences between the two drug treatment regimes; but I kind of think that is a trick played by our eyes. Glancing at the plot from sva, it appears to disagree with me; at least this recapitulates the variancePartition result I think.
wellcome_drug_norm <- normalize_expt(wellcome_drug_filt, norm = "quant", convert = "cpm",
filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'wellcome_drug_filt' not found
## Error in eval(expr, envir, enclos): object 'wellcome_drug_norm' not found
## Error in eval(expr, envir, enclos): object 'drug_norm_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'drug_norm_pca' not found
wellcome_drug_nb <- normalize_expt(wellcome_drug_filt, convert = "cpm",
batch = "svaseq", filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'wellcome_drug_filt' not found
## Error in eval(expr, envir, enclos): object 'wellcome_drug_nb' not found
## Error in eval(expr, envir, enclos): object 'drug_nb_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'drug_nb_pca' not found
Conversely, it is likely that we want to separate and/or combine the various factors with the visit. In the following block we will use the factor which is a combination of drug and visit.
drugtime <- set_expt_conditions(wellcome_filtered, fact = "drug_visit",
colors = color_choices[["drug_visit"]]) %>%
set_expt_batches(fact = "clinicaloutcome")## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_filtered' not found
drugtime_norm <- normalize_expt(drugtime, norm = "quant", convert = "cpm",
transform = "log2", filter = TRUE)## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'drugtime' not found
## Error in eval(expr, envir, enclos): object 'drugtime_norm' not found
## Error in eval(expr, envir, enclos): object 'drugtime_norm_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'drugtime_norm_pca' not found
drugtime_nb <- normalize_expt(drugtime, batch = "svaseq", convert = "cpm",
transform = "log2", filter = TRUE)## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'drugtime' not found
## Error in eval(expr, envir, enclos): object 'drugtime_nb' not found
## Error in eval(expr, envir, enclos): object 'drugtime_nb_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'drugtime_nb_pca' not found
One interpretation of some queries from reviewers would be to replot the various PCAs with the visits separated from each other. I am not completely sure that is the correct interpretation, but here it is.
drug_v1 <- subset_expt(drugtime, subset = "visitnumber=='v1'") %>%
set_expt_batches(fact = "clinicaloutcome")## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': error in evaluating the argument 'expt' in selecting a method for function 'subset_expt': object 'drugtime' not found
drug_v1_norm <- normalize_expt(drug_v1, norm = "quant", convert = "cpm",
filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'drug_v1' not found
## Error in eval(expr, envir, enclos): object 'drug_v1_norm' not found
## Error in eval(expr, envir, enclos): object 'drug_v1_norm_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'drug_v1_norm_pca' not found
drug_v1_nb <- normalize_expt(drug_v1, batch = "svaseq", convert = "cpm",
filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'drug_v1' not found
## Error in eval(expr, envir, enclos): object 'drug_v1_nb' not found
## Error in eval(expr, envir, enclos): object 'drug_v1_nb_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'drug_v1_nb_pca' not found
I would like to have something useful to say before every block. I don’t.
drug_v2 <- subset_expt(drugtime, subset = "visitnumber=='v2'") %>%
set_expt_batches(fact = "clinicaloutcome")## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': error in evaluating the argument 'expt' in selecting a method for function 'subset_expt': object 'drugtime' not found
drug_v2_norm <- normalize_expt(drug_v2, norm = "quant", convert = "cpm",
filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'drug_v2' not found
## Error in eval(expr, envir, enclos): object 'drug_v2_norm' not found
## Error in eval(expr, envir, enclos): object 'drug_v2_norm_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'drug_v2_norm_pca' not found
drug_v2_nb <- normalize_expt(drug_v2, batch = "svaseq", convert = "cpm",
filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'drug_v2' not found
## Error in eval(expr, envir, enclos): object 'drug_v2_nb' not found
## Error in eval(expr, envir, enclos): object 'drug_v2_nb_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'drug_v2_nb_pca' not found
drug_v3 <- subset_expt(drugtime, subset = "visitnumber=='v3'") %>%
set_expt_batches(fact = "clinicaloutcome")## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': error in evaluating the argument 'expt' in selecting a method for function 'subset_expt': object 'drugtime' not found
drug_v3_norm <- normalize_expt(drug_v3, norm = "quant", convert = "cpm",
filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'drug_v3' not found
## Error in eval(expr, envir, enclos): object 'drug_v3_norm' not found
## Error in eval(expr, envir, enclos): object 'drug_v3_norm_pca' not found
## png
## 2
drug_v3_nb <- normalize_expt(drug_v3, batch = "svaseq", convert = "cpm",
filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'drug_v3' not found
## Error in eval(expr, envir, enclos): object 'drug_v3_nb' not found
## Error in eval(expr, envir, enclos): object 'drug_v3_nb_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'drug_v3_nb_pca' not found
outcome <- set_expt_conditions(wellcome_filtered, fact = "clinicaloutcome",
colors = color_choices[["clinicaloutcome"]]) %>%
set_expt_batches(fact = "drug")## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_filtered' not found
outcome_norm <- normalize_expt(outcome, norm = "quant", convert = "cpm",
transform = "log2", filter = TRUE)## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'outcome' not found
## Error in eval(expr, envir, enclos): object 'outcome_norm' not found
## Error in eval(expr, envir, enclos): object 'outcome_norm_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'outcome_norm_pca' not found
outcome_nb <- normalize_expt(outcome, batch = "svaseq", convert = "cpm",
transform = "log2", filter = TRUE)## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'outcome' not found
## Error in eval(expr, envir, enclos): object 'outcome_nb' not found
## Error in eval(expr, envir, enclos): object 'outcome_nb_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'outcome_nb_pca' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'subset_expt': object 'outcome' not found
outcome_v1_norm <- normalize_expt(outcome_v1, norm = "quant", convert = "cpm",
filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'outcome_v1' not found
## Error in eval(expr, envir, enclos): object 'outcome_v1_norm' not found
## Error in eval(expr, envir, enclos): object 'outcome_v1_norm_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'outcome_v1_norm_pca' not found
outcome_v1_nb <- normalize_expt(outcome_v1, batch = "svaseq", convert = "cpm",
filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'outcome_v1' not found
## Error in eval(expr, envir, enclos): object 'outcome_v1_nb' not found
## Error in eval(expr, envir, enclos): object 'outcome_v1_nb_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'outcome_v1_nb_pca' not found
I would like to have something useful to say before every block. I don’t.
outcome_v2 <- subset_expt(outcome, subset = "visitnumber=='v2'") %>%
set_expt_batches(fact = "clinicaloutcome")## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': error in evaluating the argument 'expt' in selecting a method for function 'subset_expt': object 'outcome' not found
outcome_v2_norm <- normalize_expt(outcome_v2, norm = "quant", convert = "cpm",
filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'outcome_v2' not found
## Error in eval(expr, envir, enclos): object 'outcome_v2_norm' not found
## Error in eval(expr, envir, enclos): object 'outcome_v2_norm_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'outcome_v2_norm_pca' not found
outcome_v2_nb <- normalize_expt(outcome_v2, batch = "svaseq", convert = "cpm",
filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'outcome_v2' not found
## Error in eval(expr, envir, enclos): object 'outcome_v2_nb' not found
## Error in eval(expr, envir, enclos): object 'outcome_v2_nb_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'outcome_v2_nb_pca' not found
outcome_v3 <- subset_expt(outcome, subset = "visitnumber=='v3'") %>%
set_expt_batches(fact = "clinicaloutcome")## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': error in evaluating the argument 'expt' in selecting a method for function 'subset_expt': object 'outcome' not found
outcome_v3_norm <- normalize_expt(outcome_v3, norm = "quant", convert = "cpm",
filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'outcome_v3' not found
## Error in eval(expr, envir, enclos): object 'outcome_v3_norm' not found
## Error in eval(expr, envir, enclos): object 'outcome_v3_norm_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'outcome_v3_norm_pca' not found
outcome_v3_nb <- normalize_expt(outcome_v3, batch = "svaseq", convert = "cpm",
filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'outcome_v3' not found
## Error in eval(expr, envir, enclos): object 'outcome_v3_nb' not found
## Error in eval(expr, envir, enclos): object 'outcome_v3_nb_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'outcome_v3_nb_pca' not found
Repeat the above process, but this time using clinical outcome and time combined.
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_filtered' not found
wellcome_outtime_norm <- normalize_expt(wellcome_outtime, norm = "quant", convert = "cpm",
filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'wellcome_outtime' not found
## Error in eval(expr, envir, enclos): object 'wellcome_outtime_norm' not found
## Error in eval(expr, envir, enclos): object 'outtime_norm_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'outtime_norm_pca' not found
wellcome_outtime_nb <- normalize_expt(wellcome_outtime, convert = "cpm",
batch = "svaseq", filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'wellcome_outtime' not found
## Error in eval(expr, envir, enclos): object 'wellcome_outtime_nb' not found
## Error in eval(expr, envir, enclos): object 'outtime_nb_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'outtime_nb_pca' not found
outtime_v1 <- subset_expt(wellcome_outtime, subset = "visitnumber=='v1'") %>%
set_expt_batches(fact = "drug")## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': error in evaluating the argument 'expt' in selecting a method for function 'subset_expt': object 'wellcome_outtime' not found
outtime_v1_norm <- normalize_expt(outtime_v1, norm = "quant", convert = "cpm",
filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'outtime_v1' not found
## Error in eval(expr, envir, enclos): object 'outtime_v1_norm' not found
## Error in eval(expr, envir, enclos): object 'outtime_v1_norm_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'outtime_v1_norm_pca' not found
outtime_v1_nb <- normalize_expt(outtime_v1, batch = "svaseq", convert = "cpm",
filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'outtime_v1' not found
## Error in eval(expr, envir, enclos): object 'outtime_v1_nb' not found
## Error in eval(expr, envir, enclos): object 'outtime_v1_nb_pca' not found
outtime_v2 <- subset_expt(wellcome_outtime, subset = "visitnumber=='v2'") %>%
set_expt_batches(fact = "drug")## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': error in evaluating the argument 'expt' in selecting a method for function 'subset_expt': object 'wellcome_outtime' not found
outtime_v2_norm <- normalize_expt(outtime_v2, norm = "quant", convert = "cpm",
filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'outtime_v2' not found
## Error in eval(expr, envir, enclos): object 'outtime_v2_norm' not found
## Error in eval(expr, envir, enclos): object 'outtime_v2_norm_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'outtime_v2_norm_pca' not found
outtime_v2_nb <- normalize_expt(outtime_v2, batch = "svaseq", convert = "cpm",
filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'outtime_v2' not found
## Error in eval(expr, envir, enclos): object 'outtime_v2_nb' not found
## Error in eval(expr, envir, enclos): object 'outtime_v2_nb_pca' not found
outtime_v3 <- subset_expt(wellcome_outtime, subset = "visitnumber=='v3'") %>%
set_expt_batches(fact = "drug")## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': error in evaluating the argument 'expt' in selecting a method for function 'subset_expt': object 'wellcome_outtime' not found
outtime_v3_norm <- normalize_expt(outtime_v3, norm = "quant", convert = "cpm",
filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'outtime_v3' not found
## Error in eval(expr, envir, enclos): object 'outtime_v3_norm' not found
## Error in eval(expr, envir, enclos): object 'outtime_v3_norm_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'outtime_v3_norm_pca' not found
outtime_v3_nb <- normalize_expt(outtime_v3, batch = "svaseq", convert = "cpm",
filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'outtime_v3' not found
## Error in eval(expr, envir, enclos): object 'outtime_v3_nb' not found
## Error in eval(expr, envir, enclos): object 'outtime_v3_nb_pca' not found
speciestime <- set_expt_conditions(wellcome_outtime, fact = "species_visit",
colors = color_choices[["species_visit"]]) %>%
set_expt_batches(fact = "clinicaloutcome")## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_outtime' not found
speciestime_norm <- normalize_expt(speciestime, norm = "quant", convert = "cpm",
transform = "log2", filter = TRUE)## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'speciestime' not found
## Error in eval(expr, envir, enclos): object 'speciestime_norm' not found
## Error in eval(expr, envir, enclos): object 'speciestime_norm_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'speciestime_norm_pca' not found
speciestime_nb <- normalize_expt(speciestime, batch = "svaseq", convert = "cpm",
transform = "log2", filter = TRUE)## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'speciestime' not found
## Error in eval(expr, envir, enclos): object 'speciestime_nb' not found
## Error in eval(expr, envir, enclos): object 'speciestime_nb_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'speciestime_nb_pca' not found
species_v1 <- subset_expt(speciestime, subset = "visitnumber=='v1'") %>%
set_expt_batches(fact = "clinicaloutcome")## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': error in evaluating the argument 'expt' in selecting a method for function 'subset_expt': object 'speciestime' not found
species_v1_norm <- normalize_expt(species_v1, norm = "quant", convert = "cpm",
filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'species_v1' not found
## Error in eval(expr, envir, enclos): object 'species_v1_norm' not found
## Error in eval(expr, envir, enclos): object 'species_v1_norm_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'species_v1_norm_pca' not found
species_v1_nb <- normalize_expt(species_v1, batch = "svaseq", convert = "cpm",
filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'species_v1' not found
## Error in eval(expr, envir, enclos): object 'species_v1_nb' not found
## Error in eval(expr, envir, enclos): object 'species_v1_nb_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'species_v1_nb_pca' not found
species_v2 <- subset_expt(speciestime, subset = "visitnumber=='v2'") %>%
set_expt_batches(fact = "clinicaloutcome")## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': error in evaluating the argument 'expt' in selecting a method for function 'subset_expt': object 'speciestime' not found
species_v2_norm <- normalize_expt(species_v2, norm = "quant", convert = "cpm",
filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'species_v2' not found
## Error in eval(expr, envir, enclos): object 'species_v2_norm' not found
## Error in eval(expr, envir, enclos): object 'species_v2_norm_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'species_v2_norm_pca' not found
species_v2_nb <- normalize_expt(species_v2, batch = "svaseq", convert = "cpm",
filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'species_v2' not found
## Error in eval(expr, envir, enclos): object 'species_v2_nb' not found
## Error in eval(expr, envir, enclos): object 'species_v2_nb_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'species_v2_nb_pca' not found
species_v3 <- subset_expt(speciestime, subset = "visitnumber=='v3'") %>%
set_expt_batches(fact = "clinicaloutcome")## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': error in evaluating the argument 'expt' in selecting a method for function 'subset_expt': object 'speciestime' not found
species_v3_norm <- normalize_expt(species_v3, norm = "quant", convert = "cpm",
filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'species_v3' not found
## Error in eval(expr, envir, enclos): object 'species_v3_norm' not found
## Error in eval(expr, envir, enclos): object 'species_v3_norm_pca' not found
## png
## 2
species_v3_nb <- normalize_expt(species_v3, batch = "svaseq", convert = "cpm",
filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'species_v3' not found
## Error in eval(expr, envir, enclos): object 'species_v3_nb' not found
## Error in eval(expr, envir, enclos): object 'species_v3_nb_pca' not found
## png
## 2
## Error in eval(expr, envir, enclos): object 'species_v3_nb_pca' not found
wellcome_parasite_norm <- normalize_expt(wellcome_parasite_filt, norm = "quant", convert = "cpm",
filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'wellcome_parasite_filt' not found
## Error in eval(expr, envir, enclos): object 'wellcome_parasite_norm' not found
wellcome_parasite_nb <- normalize_expt(wellcome_parasite_filt, norm = "quant", convert = "cpm",
batch = "svaseq", filter = TRUE, transform = "log2")## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'normalize_expt': object 'wellcome_parasite_filt' not found
## Error in eval(expr, envir, enclos): object 'wellcome_parasite_nb' not found
Ok, so that is a big pile of pictures, now let us perform some DE analyses and see what pops out.
Given the above variance partition and PCA plots, we can surmise that some metadata factors are much more likely to give interesting results than others (Drug far more than C/F, for example). With that in mind, let us perform the various possible DE analyses with each factor and see what comes out. Each of these runs of all pairwise will be performed once with SVA estimates in the model, and once with the drug treatment as the batch factor.
One thing to recall, my sanitizer removes punctuation from concatenated factors; so ‘cure_v2’ turns into ‘curev2’.
The following list names each contrast and the two element vector following provides each numerator and denominator.
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"))Start with the combined factor of {outcome}_{visit}. Given the PCA, I expect to see a little bit of information with batch in the model, oddly less information with SVA.
One of the review comments seemed to me to suggest that using the variancePartition dream method might be useful. I am not certain if that is currently enabled. I definitely got it working, though. Also note that filter in this context is on a gene basis, not sample.
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_filtered' not found
## Error in eval(expr, envir, enclos): object 'outtime_sva_de' not found
outtime_sva_table <- combine_de_tables(
outtime_sva_de,
keepers = outtime_keepers,
excel = glue("excel/CL_biopsies_outtime_table_sva-v{ver}.xlsx"))## Error in eval(expr, envir, enclos): object 'outtime_sva_de' not found
## Error in eval(expr, envir, enclos): object 'outtime_sva_table' not found
outtime_sva_sig <- extract_significant_genes(
outtime_sva_table,
excel = glue("excel/CL_biopsies_outtime_sig_sva-v{ver}.xlsx"))## Error in eval(expr, envir, enclos): object 'outtime_sva_table' not found
## Error in eval(expr, envir, enclos): object 'outtime_sva_sig' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_filtered' not found
## Error in eval(expr, envir, enclos): object 'outtime_batch_de' not found
outtime_batch_table <- combine_de_tables(
outtime_batch_de,
keepers = outtime_keepers,
excel = glue("excel/CL_biopsies_outtime_table_batch-v{ver}.xlsx"))## Error in eval(expr, envir, enclos): object 'outtime_batch_de' not found
## Error in eval(expr, envir, enclos): object 'outtime_batch_table' not found
outtime_batch_sig <- extract_significant_genes(
outtime_batch_table,
excel = glue("excel/CL_biopsies_outtime_sig_batch-v{ver}.xlsx"))## Error in eval(expr, envir, enclos): object 'outtime_batch_table' not found
## Error in eval(expr, envir, enclos): object 'outtime_batch_sig' not found
There was a request some time ago to provide a power analysis. I arbitrarily chose this dataset in order to invoke PROPER and provide the resulting plots…
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'exprs': object 'outtime_batch_table' not found
## Error in eval(expr, envir, enclos): object 'outtime_batch_proper' not found
## Error in eval(expr, envir, enclos): object 'outtime_batch_proper' not found
## Error in eval(expr, envir, enclos): object 'outtime_batch_proper' not found
## Error in eval(expr, envir, enclos): object 'outtime_batch_proper' not found
## Error in eval(expr, envir, enclos): object 'outtime_batch_proper' not found
## Error in eval(expr, envir, enclos): object 'outtime_batch_proper' not found
## Error in eval(expr, envir, enclos): object 'outtime_batch_proper' not found
Now let us compare the time points, with and without SVA.
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_time_filt' not found
## Error in eval(expr, envir, enclos): object 'time_sva_de' not found
time_sva_table <- combine_de_tables(
time_sva_de,
keepers = time_keepers,
excel = glue("excel/CL_biopsies_time_table_sva-v{ver}.xlsx"))## Error in eval(expr, envir, enclos): object 'time_sva_de' not found
## Error in eval(expr, envir, enclos): object 'time_sva_table' not found
time_sva_sig <- extract_significant_genes(
time_sva_table,
excel = glue("excel/CL_biopsies_time_sig_sva-v{ver}.xlsx"))## Error in eval(expr, envir, enclos): object 'time_sva_table' not found
## Error in eval(expr, envir, enclos): object 'time_sva_sig' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_time_filt' not found
## Error in eval(expr, envir, enclos): object 'time_batch_de' not found
time_batch_table <- combine_de_tables(
time_batch_de,
keepers = time_keepers,
excel = glue("excel/CL_biopsies_time_table_batch-v{ver}.xlsx"))## Error in eval(expr, envir, enclos): object 'time_batch_de' not found
## Error in eval(expr, envir, enclos): object 'time_batch_table' not found
time_batch_sig <- extract_significant_genes(
time_batch_table,
excel = glue("excel/CL_biopsies_time_batch_sva-v{ver}.xlsx"))## Error in eval(expr, envir, enclos): object 'time_batch_table' not found
## Error in eval(expr, envir, enclos): object 'time_batch_sig' not found
I neglected to set a contrast name/value for this.
The infecting strain I think should prove one of the more interesting comparisons.
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_parasite_filt' not found
## Error in eval(expr, envir, enclos): object 'parasite_sva_de' not found
parasite_sva_table <- combine_de_tables(
parasite_sva_de,
excel = glue("excel/CL_biopsies_parasite_table_sva-v{ver}.xlsx"))## Error in eval(expr, envir, enclos): object 'parasite_sva_de' not found
## Error in eval(expr, envir, enclos): object 'parasite_sva_table' not found
parasite_sva_sig <- extract_significant_genes(
parasite_sva_table,
excel = glue("excel/CL_biopsies_parasite_sig_sva-v{ver}.xlsx"))## Error in eval(expr, envir, enclos): object 'parasite_sva_table' not found
## Error in eval(expr, envir, enclos): object 'parasite_sva_sig' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_parasite_filt' not found
## Error in eval(expr, envir, enclos): object 'parasite_batch_de' not found
parasite_batch_table <- combine_de_tables(
parasite_batch_de,
excel = glue("excel/CL_biopsies_parasite_table_batch-v{ver}.xlsx"))## Error in eval(expr, envir, enclos): object 'parasite_batch_de' not found
## Error in eval(expr, envir, enclos): object 'parasite_batch_table' not found
parasite_batch_sig <- extract_significant_genes(
parasite_batch_table,
excel = glue("excel/CL_biopsies_parasite_sig_sva-v{ver}.xlsx"))## Error in eval(expr, envir, enclos): object 'parasite_batch_table' not found
## Error in eval(expr, envir, enclos): object 'parasite_batch_sig' not found
As stated above, clinical outcome by itself does not have much information content in this dataset.
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_outcome_filt' not found
## Error in eval(expr, envir, enclos): object 'outcome_sva_de' not found
outcome_sva_table <- combine_de_tables(
outcome_sva_de,
excel = glue("excel/CL_biopsies_outcome_table_sva-v{ver}.xlsx"))## Error in eval(expr, envir, enclos): object 'outcome_sva_de' not found
## Error in eval(expr, envir, enclos): object 'outcome_sva_table' not found
outcome_sva_sig <- extract_significant_genes(
outcome_sva_table,
excel = glue("excel/CL_biopsies_outcome_sig_sva-v{ver}.xlsx"))## Error in eval(expr, envir, enclos): object 'outcome_sva_table' not found
## Error in eval(expr, envir, enclos): object 'outcome_sva_sig' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_outcome_filt' not found
## Error in eval(expr, envir, enclos): object 'outcome_batch_de' not found
outcome_batch_table <- combine_de_tables(
outcome_batch_de,
excel = glue("excel/CL_biopsies_outcome_table_batch-v{ver}.xlsx"))## Error in eval(expr, envir, enclos): object 'outcome_batch_de' not found
## Error in eval(expr, envir, enclos): object 'outcome_batch_table' not found
outcome_batch_sig <- extract_significant_genes(
outcome_batch_table,
excel = glue("excel/CL_biopsies_outcome_sig_sva-v{ver}.xlsx"))## Error in eval(expr, envir, enclos): object 'outcome_batch_table' not found
## Error in eval(expr, envir, enclos): object 'outcome_batch_sig' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_drug_filt' not found
## Error in eval(expr, envir, enclos): object 'drug_sva_de' not found
drug_sva_table <- combine_de_tables(
drug_sva_de,
excel = glue("excel/CL_biopsies_drug_table_sva-v{ver}.xlsx"))## Error in eval(expr, envir, enclos): object 'drug_sva_de' not found
## Error in eval(expr, envir, enclos): object 'drug_sva_table' not found
drug_sva_sig <- extract_significant_genes(
drug_sva_table,
excel = glue("excel/CL_biopsies_drug_sig_sva-v{ver}.xlsx"))## Error in eval(expr, envir, enclos): object 'drug_sva_table' not found
## Error in eval(expr, envir, enclos): object 'drug_sva_sig' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': object 'wellcome_drug_filt' not found
## Error in eval(expr, envir, enclos): object 'drug_batch_de' not found
drug_batch_table <- combine_de_tables(
drug_batch_de,
excel = glue("excel/CL_biopsies_drug_table_batch-v{ver}.xlsx"))## Error in eval(expr, envir, enclos): object 'drug_batch_de' not found
## Error in eval(expr, envir, enclos): object 'drug_batch_table' not found
drug_batch_sig <- extract_significant_genes(
drug_batch_table,
excel = glue("excel/CL_biopsies_drug_sig_sva-v{ver}.xlsx"))## Error in eval(expr, envir, enclos): object 'drug_batch_table' not found
## Error in eval(expr, envir, enclos): object 'drug_batch_sig' not found
This requires setting up the contrasts of interest. I am going to assume cross species/cross time are not interesting, and that v3/v1 and v3/v2 are not interesting.
Note to self, do not forget that I sanitize the data by removing punctuation!
speciestime_keepers <- list(
"pan_v2v1" = c("lvpanamensisv2", "lvpanamensisv1"),
"pan_v3v1" = c("lvpanamensisv3", "lvpanamensisv1"),
"braz_v2v1" = c("lvbraziliensisv2", "lvbraziliensisv1"),
"braz_v3v1" = c("lvbraziliensisv3", "lvbraziliensisv1"),
"v1_panbraz" = c("lvbraziliensisv1", "lvpanamensisv1"),
"v2_panbraz" = c("lvbraziliensisv2", "lvpanamensisv2"),
"v3_panbraz" = c("lvbraziliensisv3", "lvpanamensisv3"))## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': object 'speciestime' not found
## Error in eval(expr, envir, enclos): object 'speciestime_sva_de' not found
speciestime_sva_table <- combine_de_tables(
speciestime_sva_de, keepers = speciestime_keepers,
excel = glue("excel/CL_biopsies_speciestime_table_sva-v{ver}.xlsx"))## Error in eval(expr, envir, enclos): object 'speciestime_sva_de' not found
## Error in eval(expr, envir, enclos): object 'speciestime_sva_table' not found
speciestime_sva_sig <- extract_significant_genes(
speciestime_sva_table,
excel = glue("excel/CL_biopsies_speciestime_sig_sva-v{ver}.xlsx"))## Error in eval(expr, envir, enclos): object 'speciestime_sva_table' not found
## Error in eval(expr, envir, enclos): object 'speciestime_sva_sig' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': object 'speciestime' not found
## Error in eval(expr, envir, enclos): object 'speciestime_sva_de' not found
speciestime_sva_table <- combine_de_tables(
speciestime_sva_de, keepers = speciestime_keepers,
excel = glue("excel/CL_biopsies_speciestime_table_sva-v{ver}.xlsx"))## Error in eval(expr, envir, enclos): object 'speciestime_sva_de' not found
## Error in eval(expr, envir, enclos): object 'speciestime_sva_table' not found
speciestime_sva_sig <- extract_significant_genes(
speciestime_sva_table,
excel = glue("excel/CL_biopsies_speciestime_sig_sva-v{ver}.xlsx"))## Error in eval(expr, envir, enclos): object 'speciestime_sva_table' not found
## Error in eval(expr, envir, enclos): object 'speciestime_sva_sig' not found
Same rules here as strain+time
drugtime_keepers <- list(
"milt_v2v1" = c("miltefosinev2", "miltefosinev1"),
"milt_v3v1" = c("miltefosinev3", "miltefosinev1"),
"gluc_v2v1" = c("antimoniatev2", "antimoniatev1"),
"gluc_v3v1" = c("antimoniatev3", "antimoniatev1"),
"v1_miltgluc" = c("miltefosinev1", "antimoniatev1"),
"v2_miltgluc" = c("miltefosinev2", "antimoniatev2"),
"v3_miltgluc" = c("miltefosinev3", "antimoniatev3"))## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': object 'drugtime' not found
## Error in eval(expr, envir, enclos): object 'drugtime_sva_de' not found
drugtime_sva_table <- combine_de_tables(
drugtime_sva_de, keepers = drugtime_keepers,
excel = glue("excel/CL_biopsies_drugtime_table_sva-v{ver}.xlsx"))## Error in eval(expr, envir, enclos): object 'drugtime_sva_de' not found
## Error in eval(expr, envir, enclos): object 'drugtime_sva_table' not found
drugtime_sva_sig <- extract_significant_genes(
drugtime_sva_table,
excel = glue("excel/CL_biopsies_drugtime_sig_sva-v{ver}.xlsx"))## Error in eval(expr, envir, enclos): object 'drugtime_sva_table' not found
## Error in eval(expr, envir, enclos): object 'drugtime_sva_sig' not found
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': object 'drugtime' not found
## Error in eval(expr, envir, enclos): object 'drugtime_sva_de' not found
drugtime_sva_table <- combine_de_tables(
drugtime_sva_de, keepers = drugtime_keepers,
excel = glue("excel/CL_biopsies_drugtime_table_sva-v{ver}.xlsx"))## Error in eval(expr, envir, enclos): object 'drugtime_sva_de' not found
## Error in eval(expr, envir, enclos): object 'drugtime_sva_table' not found
drugtime_sva_sig <- extract_significant_genes(
drugtime_sva_table,
excel = glue("excel/CL_biopsies_drugtime_sig_sva-v{ver}.xlsx"))## Error in eval(expr, envir, enclos): object 'drugtime_sva_table' not found
## Error in eval(expr, envir, enclos): object 'drugtime_sva_sig' not found
Given the results in the previous blocks, let us use gProfiler2 in look for over represented categories/ontologies in the suspect databases: GO/reactome/kegg/TF/wikipath/etc.
I think for the moment I will do this entirely with the sva results.
I have a series of wrapper functions for gProfiler2/goseq/clusterProfiler/topGO/GOstats. Some of them are smart enough to take as input the result from extract_significant_genes().
There is also a function which writes pretty-ified ontology results, but it works on a single-table basis.
## Error in eval(expr, envir, enclos): object 'outtime_sva_sig' not found
## Error in eval(expr, envir, enclos): object 'outtime_gprofiler' not found
## Error in eval(expr, envir, enclos): object 'outtime_gprofiler' not found
## Error in eval(expr, envir, enclos): object 'outtime_gprofiler' not found
## Note, I am writing a set of methods for write_gprofiler_data so that it is
## aware of the result from all_gprofiler and can handle either all contrasts
## or a set of desired contrasts.
write_all <- function(result) {
for (n in names(result)) {
name <- glue("excel/gprofiler_{n}.xlsx")
res <- result[[n]]
written <- try(write_gprofiler_data(res, excel = name))
}
}
written <- write_all(outtime_gprofiler)## Error in eval(expr, envir, enclos): object 'outtime_gprofiler' not found
## Error in eval(expr, envir, enclos): object 'time_sva_sig' not found
## Error in eval(expr, envir, enclos): object 'visit_gprofiler' not found
## Error in eval(expr, envir, enclos): object 'parasite_sva_sig' not found
## Error in eval(expr, envir, enclos): object 'parasite_gprofiler' not found
## Error in eval(expr, envir, enclos): object 'outcome_sva_sig' not found
## Error in eval(expr, envir, enclos): object 'outcome_gprofiler' not found
## Error in eval(expr, envir, enclos): object 'drug_sva_sig' not found
## Error in eval(expr, envir, enclos): object 'drug_gprofiler' not found
## Error in eval(expr, envir, enclos): object 'speciestime_sva_sig' not found
## Error in eval(expr, envir, enclos): object 'speciestime_gprofiler' not found
From my perspective, it seems that the most interesting GSVA comparison is to look for scores changing between the cure and fail samples. The following block therefore passes the raw expression data and performs the gene set variance analysis of it against the mSigDB C2 set of categories by default. We can pretty easily change the mSigDB release/category set; except I would need download the dataset to the container to use a newer revision and I am not sure about potential licensing restrictions; so this will just use the publicly available signatures in the ‘GSVAdata’ package.
## Error in eval(expr, envir, enclos): object 'wellcome_outcome_filt' not found
## Error in eval(expr, envir, enclos): object 'wt_gsva' not found
## Error in eval(expr, envir, enclos): object 'wt_gsva' not found
## Error in eval(expr, envir, enclos): object 'wt_gsva_sig' not found
Hmm, did I get confused, is C7 what I thought it was? Yeah, it is the immunologic signature sets. C2 is the larger set of experiments.
## Error in eval(expr, envir, enclos): object 'wellcome_filtered' not found
wellcome_gsva_c2_sig <- get_sig_gsva_categories(
wellcome_gsva_c2,
excel = "excel/CL_biopsies_gsva_outcome_c2.xlsx")## Error in eval(expr, envir, enclos): object 'wellcome_gsva_c2' not found
The following blocks may be used to save and reload the state of the data.
R version 4.3.1 (2023-06-16)
Platform: x86_64-conda-linux-gnu (64-bit)
locale: C
attached base packages: stats4, stats, graphics, grDevices, utils, datasets, methods and base
other attached packages: BiocParallel(v.1.36.0), variancePartition(v.1.32.5), hpgltools(v.1.0), Matrix(v.1.6-5), glue(v.1.7.0), SummarizedExperiment(v.1.32.0), GenomicRanges(v.1.54.1), GenomeInfoDb(v.1.38.8), IRanges(v.2.36.0), S4Vectors(v.0.40.2), MatrixGenerics(v.1.14.0), matrixStats(v.1.3.0), Biobase(v.2.62.0) and BiocGenerics(v.0.48.1)
loaded via a namespace (and not attached): splines(v.4.3.1), later(v.1.3.2), BiocIO(v.1.12.0), ggplotify(v.0.1.2), bitops(v.1.0-7), filelock(v.1.0.3), tibble(v.3.2.1), polyclip(v.1.10-6), graph(v.1.80.0), XML(v.3.99-0.16.1), lifecycle(v.1.0.4), Rdpack(v.2.6), doParallel(v.1.0.17), lattice(v.0.22-6), MASS(v.7.3-60.0.1), backports(v.1.4.1), magrittr(v.2.0.3), openxlsx(v.4.2.5.2), limma(v.3.58.1), plotly(v.4.10.4), sass(v.0.4.9), rmarkdown(v.2.26), jquerylib(v.0.1.4), yaml(v.2.3.8), httpuv(v.1.6.15), zip(v.2.3.1), RColorBrewer(v.1.1-3), cowplot(v.1.1.3), DBI(v.1.2.2), minqa(v.1.2.6), lubridate(v.1.9.3), abind(v.1.4-5), zlibbioc(v.1.48.2), EnvStats(v.2.8.1), purrr(v.1.0.2), ggraph(v.2.2.1), RCurl(v.1.98-1.14), yulab.utils(v.0.1.4), tweenr(v.2.0.3), rappdirs(v.0.3.3), GenomeInfoDbData(v.1.2.11), enrichplot(v.1.22.0), ggrepel(v.0.9.5), pbkrtest(v.0.5.2), tidytree(v.0.4.6), annotate(v.1.80.0), codetools(v.0.2-20), DelayedArray(v.0.28.0), ggforce(v.0.4.2), DOSE(v.3.28.2), xml2(v.1.3.6), tidyselect(v.1.2.1), aplot(v.0.2.2), farver(v.2.1.1), viridis(v.0.6.5), lme4(v.1.1-35.3), BiocFileCache(v.2.10.2), GenomicAlignments(v.1.38.2), jsonlite(v.1.8.8), tidygraph(v.1.3.1), iterators(v.1.0.14), foreach(v.1.5.2), tools(v.4.3.1), progress(v.1.2.3), treeio(v.1.26.0), Rcpp(v.1.0.12), gridExtra(v.2.3), SparseArray(v.1.2.4), xfun(v.0.43), qvalue(v.2.34.0), dplyr(v.1.1.4), withr(v.3.0.0), numDeriv(v.2016.8-1.1), BiocManager(v.1.30.22), fastmap(v.1.1.1), boot(v.1.3-30), fansi(v.1.0.6), caTools(v.1.18.2), digest(v.0.6.35), gridGraphics(v.0.5-1), timechange(v.0.3.0), R6(v.2.5.1), mime(v.0.12), colorspace(v.2.1-0), GO.db(v.3.18.0), gtools(v.3.9.5), biomaRt(v.2.58.2), RSQLite(v.2.3.6), RhpcBLASctl(v.0.23-42), utf8(v.1.2.4), tidyr(v.1.3.1), generics(v.0.1.3), data.table(v.1.15.4), corpcor(v.1.6.10), rtracklayer(v.1.62.0), graphlayouts(v.1.1.1), prettyunits(v.1.2.0), httr(v.1.4.7), htmlwidgets(v.1.6.4), S4Arrays(v.1.2.1), scatterpie(v.0.2.2), pkgconfig(v.2.0.3), gtable(v.0.3.4), blob(v.1.2.4), XVector(v.0.42.0), remaCor(v.0.0.18), shadowtext(v.0.1.3), htmltools(v.0.5.8.1), fgsea(v.1.28.0), GSEABase(v.1.64.0), scales(v.1.3.0), png(v.0.1-8), fANCOVA(v.0.6-1), ggfun(v.0.1.4), knitr(v.1.46), reshape2(v.1.4.4), rjson(v.0.2.21), nlme(v.3.1-164), curl(v.5.2.1), nloptr(v.2.0.3), cachem(v.1.0.8), stringr(v.1.5.1), KernSmooth(v.2.23-22), parallel(v.4.3.1), HDO.db(v.0.99.1), AnnotationDbi(v.1.64.1), restfulr(v.0.0.15), pillar(v.1.9.0), grid(v.4.3.1), vctrs(v.0.6.5), gplots(v.3.1.3.1), promises(v.1.3.0), dbplyr(v.2.3.4), xtable(v.1.8-4), evaluate(v.0.23), mvtnorm(v.1.2-4), cli(v.3.6.2), compiler(v.4.3.1), Rsamtools(v.2.18.0), rlang(v.1.1.3), crayon(v.1.5.2), plyr(v.1.8.9), fs(v.1.6.3), pander(v.0.6.5), stringi(v.1.8.3), viridisLite(v.0.4.2), lmerTest(v.3.1-3), munsell(v.0.5.1), Biostrings(v.2.70.3), lazyeval(v.0.2.2), aod(v.1.3.3), GOSemSim(v.2.28.1), patchwork(v.1.2.0), hms(v.1.1.3), bit64(v.4.0.5), ggplot2(v.3.5.0), KEGGREST(v.1.42.0), statmod(v.1.5.0), shiny(v.1.8.1.1), rbibutils(v.2.2.16), igraph(v.2.0.3), broom(v.1.0.5), memoise(v.2.0.1), bslib(v.0.7.0), ggtree(v.3.10.1), fastmatch(v.1.1-4), bit(v.4.0.5) and ape(v.5.8)
## If you wish to reproduce this exact build of hpgltools, invoke the following:
## > git clone http://github.com/abelew/hpgltools.git
## > git reset 24f917df4a7665fa4d6050160cc8117a26eab4d5
## This is hpgltools commit: Thu Apr 18 16:11:22 2024 -0400: 24f917df4a7665fa4d6050160cc8117a26eab4d5
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))## Error in save(list = ls(all.names = TRUE, envir = globalenv()), envir = globalenv(), : ignoring SIGPIPE signal