Macrophage response to different strains 202508: Data Set Creation

1 Introduction

This document is intended to create the data structures used to evaluate our TMRC2 samples. In some cases, this includes only those samples starting in 2019; in other instances I am including our previous (2015-2016) samples.

In all cases the processing performed was:

1. Default trimming was performed.
2. Hisat2 was used to map the remaining reads against the Leishmania panamensis genome revision 36.
3. The alignments from hisat2 were used to count reads/gene against the revision 36 annotations with htseq.
4. These alignments were also passed to the pileup functionality of samtools and the vcf/bcf utilities in order to make a matrix of all observed differences between each sample with respect to the reference.
5. The freebayes variant estimation tool was used in addition to #4 to search for variant positions in a more robust fashion.
6. The trimmed reads were passed to kraken2 using a viral database in order to look for samples with potential LRV sequence.
7. An explicit, grep-based search for spliced leader reads was used against all human-derived samples. The results from this were copy/pasted into the sample sheet.

1.1 Multiple datasets

In a couple of important ways the TMRC2 data is much more complex than the TMRC3:

1. It comprises multiple, completely separate queries:
   1. Sequencing the parasite samples
   2. Sequencing a set of human macrophage samples which were infected with specific parasite samples.
2. The parasite transcriptomic samples comprise multiple different types of queries:
   1. Differential expression to look at strain, susceptibility, and clinical outcomes.
   2. Individual variant searches to look for potentially useful SNPs for classification of parasite samples.
3. The human macrophage samples may be used to query both the host and parasite transcriptomes because (at least when not drug treated) there is a tremendous population of parasite reads in them.

1.2 Sample sheet(s)

Our shared online sample sheet is nearly static at the time of this writing (202209), I expect at this point the only likely updates will be to annotate some strains as more or less susceptible to drug treatment.

sample_sheet <- "sample_sheets/macrophage_samples.xlsx"

2 Annotations

Everything which follows depends on the Existing TriTrypDB annotations revision 46, circa 2019. The following block loads a database of these annotations and turns it into a matrix where the rows are genes and columns are all the annotation types provided by TriTrypDB.

The same database was used to create a matrix of orthologous genes between L.panamensis and all of the other species in the TriTrypDB.

The same database of annotations also provides mappings to the set of annotated GO categories for the L.panamensis genome along with gene lengths.

meta <- download_eupath_metadata(webservice = "tritrypdb", overwrite = FALSE)
panamensis_entry <- get_eupath_entry("MHOM", metadata = meta[["valid"]])
panamensis_db <- make_eupath_orgdb(panamensis_entry)
panamensis_pkg <- panamensis_db[["pkgname"]]
package_name <- panamensis_db[["pkgname"]]
if (is.null(panamensis_pkg)) {
  panamensis_pkg <- panamensis_entry[["OrgdbPkg"]]
  package_name <- panamensis_pkg
}

tt <- library(panamensis_pkg, character.only = TRUE)
panamensis_env <- get0(panamensis_pkg)
all_fields <- columns(panamensis_env)
all_lp_annot <- sm(load_orgdb_annotations(
    panamensis_env,
    keytype = "gid",
    fields = c("annot_gene_entrez_id", "annot_gene_name",
               "annot_strand", "annot_chromosome", "annot_cds_length",
               "annot_gene_product")))$genes
## Testing to see just how big the full database is.
## testing <- load_orgdb_annotations(panamensis_pkg, keytype = "gid", fields = "all")

lp_go <- load_orgdb_go(panamensis_pkg)
lp_go <- lp_go[, c("GID", "GO")]
lp_lengths <- all_lp_annot[, c("gid", "annot_cds_length")]
colnames(lp_lengths)  <- c("ID", "length")
all_lp_annot[["annot_gene_product"]] <- tolower(all_lp_annot[["annot_gene_product"]])
orthos <- sm(extract_eupath_orthologs(db = panamensis_pkg))
data_structures <- c(data_structures, "lp_lengths", "lp_go", "all_lp_annot")

all_installed <- rownames(installed.packages())
candidates <- grepl(pattern = "^org.Lpanamensis.MHOM.*v68.eg.db", x = all_installed)
orgdb_pkg_name <- all_installed[candidates]

tt <- library(orgdb_pkg_name, character.only = TRUE)

## Loading required package: AnnotationDbi

## Loading required package: stats4

## Loading required package: BiocGenerics

## Loading required package: generics

## 
## Attaching package: 'generics'

## The following object is masked from 'package:dplyr':
## 
##     explain

## The following objects are masked from 'package:base':
## 
##     as.difftime, as.factor, as.ordered, intersect, is.element, setdiff,
##     setequal, union

## 
## Attaching package: 'BiocGenerics'

## The following object is masked from 'package:dplyr':
## 
##     combine

## The following objects are masked from 'package:hpgltools':
## 
##     conditions, conditions<-, normalize

## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs

## The following objects are masked from 'package:base':
## 
##     Filter, Find, Map, Position, Reduce, anyDuplicated, aperm, append,
##     as.data.frame, basename, cbind, colnames, dirname, do.call,
##     duplicated, eval, evalq, get, grep, grepl, is.unsorted, lapply,
##     mapply, match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
##     rank, rbind, rownames, sapply, saveRDS, table, tapply, unique,
##     unsplit, which.max, which.min

## Loading required package: Biobase

## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.

## Loading required package: IRanges

## Loading required package: S4Vectors

## 
## Attaching package: 'S4Vectors'

## The following objects are masked from 'package:dplyr':
## 
##     first, rename

## The following object is masked from 'package:utils':
## 
##     findMatches

## The following objects are masked from 'package:base':
## 
##     I, expand.grid, unname

## 
## Attaching package: 'IRanges'

## The following objects are masked from 'package:dplyr':
## 
##     collapse, desc, slice

## The following object is masked from 'package:glue':
## 
##     trim

## 
## Attaching package: 'AnnotationDbi'

## The following object is masked from 'package:dplyr':
## 
##     select

## 

panamensis_pkg <- get0(orgdb_pkg_name)
all_fields <- columns(panamensis_pkg)
all_lp_annot <- sm(load_orgdb_annotations(
    panamensis_pkg,
    keytype = "gid",
    fields = c("annot_gene_entrez_id", "annot_gene_name",
               "annot_strand", "annot_chromosome", "annot_cds_length",
               "annot_gene_product")))$genes

lp_go <- load_orgdb_go(panamensis_pkg)

## The chosen keytype was not available.  Using 'GID'.

## This is an orgdb, good.

## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns

lp_go <- lp_go[, c("GID", "GO")]
lp_lengths <- all_lp_annot[, c("gid", "annot_cds_length")]
colnames(lp_lengths)  <- c("ID", "length")
all_lp_annot[["annot_gene_product"]] <- tolower(all_lp_annot[["annot_gene_product"]])
data_structures <- c(data_structures, "lp_lengths", "lp_go", "all_lp_annot")

3 Load a genome

The following block loads the full genome sequence for panamensis. We may use this later to attempt to estimate PCR primers to discern strains.

I am not sure how to increase the number of open files in a container, as a result this does not work.

testing_panamensis <- make_eupath_bsgenome(entry = panamensis_entry)
library(as.character(testing_panamensis), character.only = TRUE)
lp_genome <- get0(as.character(testing_panamensis))
data_structures <- c(data_structures, "lp_genome")

4 Generate Expressionsets and Sample Estimation

The process of sample estimation takes two primary inputs:

1. The sample sheet, which contains all the metadata we currently have on hand, including filenames for the outputs of #3 and #4 above.
2. The gene annotations.

An expressionSet(or summarizedExperiment) is a data structure used in R to examine RNASeq data. It is comprised of annotations, metadata, and expression data. In the case of our processing pipeline, the location of the expression data is provided by the filenames in the metadata.

4.1 Define colors

The following list contains the colors we have chosen to use when plotting the various ways of discerning the data.

color_choices <- list(
  "strain" = list(
    ## "z1.0" = "#333333", ## Changed this to 'braz' to make it easier to find them.
    "z2.0" = "#555555",
    "z3.0" = "#777777",
    "z2.1" = "#874400",
    "z2.2" = "#0000cc",
    "z2.3" = "#cc0000",
    "z2.4" = "#df7000",
    "z3.2" = "#888888",
    "z1.0" = "#cc00cc",
    "z1.5" = "#cc00cc",
    "b2904" = "#cc00cc",
    "unknown" = "#cbcbcb"),
  ## "null" = "#000000"),
  "zymo" = list(
    "none" = "#000000",
    "z22" = "#0000cc",
    "z23" = "#cc0000"),
  "cf" = list(
    "cure" = "#006f00",
    "fail" = "#9dffa0",
    "unknown" = "#cbcbcb",
    "notapplicable" = "#000000"),
  "condition" = list(
    "inf" = "#199c75",
    "inf_sb" = "#d65d00",
    "uninf" = "#6e6ea3",
    "uninf_sb" = "#d83956"),
  "significance" = list(
    "lt0" = "#ffe0e0",
    "lt1" = "#ffa0a0",
    "lt2" = "#f94040",
    "lt4" = "#a00000",
    "gt0" = "#eeccf9",
    "gt1" = "#de8bf9",
    "gt2" = "#ad07e3",
    "gt4" = "#410257"),
  "drug" = list(
    "none" = "#989898",
    "antimony" = "#088b64"),
  "oldnew" = list(
    "previous" = "#2233aa",
    "current" =  "#9c0303"),
  "infectedp" = list(
    "uninfected" = "#676767",
    "infected" = "#ac06e2"),
  "treatment_zymo" = list(
    "inf_sb_z23" = "#E7298A",
    "inf_z23" = "#D95F02",
    "uninf_none" = "#66A61E",
    "uninf_sb_none" = "#E6AB02",
    "inf_z22" = "#1B9E77",
    "inf_sb_z22" = "#7570B3"),
  "donor" = list(
    "d01" = "#8d0000",
    "d02" = "#E6AB02",
    "d09" = "#7570B3",
    "d81" = "#2233aa"),
  "susceptibility" = list(
    "resistant" = "#8563a7",
    "sensitive" = "#8d0000",
    "ambiguous" = "#cbcbcb",
    "unknown" = "#555555"))
data_structures <- c(data_structures, "color_choices")

5 Macrophage data

All of the above focused entire on the parasite samples, now let us pull up the macrophage infected samples. This will comprise two datasets, one of the human and one of the parasite.

5.1 Macrophage host data

The metadata for the macrophage samples contains a couple of columns for mapped human and parasite reads. We will therefore use them separately to create two expressionsets, one for each species.

** Note **: I forgot to commit the addition of plot_metadata factors() in the last run of this. In addition, I need to add an explicit month to load_biomart_annotations() or change the function to search a couple more months before it stops trying to find an archive.

hs_annot <- load_biomart_annotations(year = "2020", month = "04", symbol_columns = "hgnc_symbol")

## Using mart: ENSEMBL_MART_ENSEMBL from host: apr2020.archive.ensembl.org.

## Successfully connected to the hsapiens_gene_ensembl database.

## Finished downloading ensembl gene annotations.

## Finished downloading ensembl structure annotations.

## Including symbols, there are 67159 vs the 249740 gene annotations.

## Not dropping haplotype chromosome annotations, set drop_haplotypes = TRUE if this is bad.

## Saving annotations to hsapiens_biomart_annotations.rda.

## Finished save().

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

sanitize_columns <- c("drug", "macrophagetreatment", "macrophagezymodeme")
macr_annot <- hs_annot
rownames(macr_annot) <- gsub(x = rownames(macr_annot),
                             pattern = "^gene:",
                             replacement = "")

hs_macrophage <- create_se(sample_sheet, gene_info = macr_annot,
                           file_column = "hg38100hisatfile") %>%
  set_conditions(fact = "macrophagetreatment") %>%
  set_batches(fact = "macrophagezymodeme") %>%
  sanitize_metadata(columns = sanitize_columns) %>%
  subset_se(nonzero = 12000)

## Reading the sample metadata.

## 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: 69 rows(samples) and 86 columns(metadata fields).

## Matched 21481 annotations and counts.

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

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

## The numbers of samples by condition are:

## 
##      inf   inf_sb    uninf uninf_sb 
##       30       29        5        5

## The number of samples by batch are:

## 
## none z2.2 z2.3 
##   10   30   29

## Recasting the data.frame to DataFrame.

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

## TMRC30162 
##    521145 
## TMRC30162 
##     10208

## Samples removed: 10208

fixed_genenames <- gsub(x = rownames(assay(hs_macrophage)), pattern = "^gene:",
                        replacement = "")
hs_macrophage <- set_se_genenames(hs_macrophage, ids = fixed_genenames)
table(colData(hs_macrophage)$condition)

## 
##      inf   inf_sb    uninf uninf_sb 
##       29       29        5        5

## Added to make a simplified PCA plot.
colData(hs_macrophage)[["experiment"]] <- "macrophage"

## The following 3 lines were copy/pasted to datastructures and should be removed soon.
nostrain <- is.na(colData(hs_macrophage)[["strainid"]])
colData(hs_macrophage)[nostrain, "strainid"] <- "none"

colData(hs_macrophage)[["strain_zymo"]] <- paste0("s", colData(hs_macrophage)[["strainid"]],
                                                "_", colData(hs_macrophage)[["macrophagezymodeme"]])
uninfected <- colData(hs_macrophage)[["strain_zymo"]] == "snone_none"
colData(hs_macrophage)[uninfected, "strain_zymo"] <- "uninfected"

colData(hs_macrophage)[["infectedp"]] <- "infected"
colData(hs_macrophage)[uninfected, "infectedp"] <- "uninfected"

data_structures <- c(data_structures, "hs_macrophage")

5.2 Double-check sample IDs against the sample sheet

1 sample has been excluded from the analysis but is in the sample sheet. I am reasonably certain I know which, but will double-check here.

sample_sheet_ids = c("TMRC30051","TMRC30057","TMRC30059","TMRC30060","TMRC30061","TMRC30062",
                     "TMRC30063","TMRC30064","TMRC30065","TMRC30066","TMRC30067","TMRC30069",
                     "TMRC30117","TMRC30162","TMRC30243","TMRC30244","TMRC30245","TMRC30246",
                     "TMRC30247","TMRC30248","TMRC30249","TMRC30250","TMRC30251","TMRC30252",
                     "TMRC30266","TMRC30267","TMRC30268","TMRC30286","TMRC30326","TMRC30316",
                     "TMRC30317","TMRC30322","TMRC30323","TMRC30328","TMRC30318","TMRC30319",
                     "TMRC30324","TMRC30325","TMRC30320","TMRC30321","TMRC30327","TMRC30312",
                     "TMRC30297","TMRC30298","TMRC30299","TMRC30300","TMRC30295","TMRC30296",
                     "TMRC30303","TMRC30304","TMRC30301","TMRC30302","TMRC30314","TMRC30315",
                     "TMRC30313")
found <- sample_sheet_ids %in% colnames(assay(hs_macrophage))
colnames(assay(hs_macrophage))[!found]

5.3 Subset and create different groupings

all_human <- sanitize_metadata(hs_macrophage, columns = "drug") %>%
  set_conditions(fact = "drug") %>%
  set_batches(fact = "typeofcells")

## Recasting the data.frame to DataFrame.

## The numbers of samples by condition are:

## 
## antimony     none 
##       34       34

## The number of samples by batch are:

## 
## Macrophages        U937 
##          54          14

data_structures <- c(data_structures, "all_human")

## The following 3 lines were copy/pasted to datastructures and should be removed soon.
no_strain_idx <- colData(all_human)[["strainid"]] == "none"
##colData(all_human)[["strainid"]] <- paste0("s", colData(all_human)[["strainid"]],
##                                         "_", colData(all_human)[["macrophagezymodeme"]])
colData(all_human)[no_strain_idx, "strainid"] <- "none"
table(colData(all_human)[["strainid"]])

## 
## 10763 10772 10977 11026 11075 11126 12251 12309 12355 12367  2169  7158  none 
##     2     8     2     2     2     8     7     8     2     7     8     2    10

all_human_types <- set_conditions(all_human, fact = "typeofcells") %>%
  set_batches(fact = "drug")

## The numbers of samples by condition are:

## 
## Macrophages        U937 
##          54          14

## The number of samples by batch are:

## 
## antimony     none 
##       34       34

data_structures <- c(data_structures, "all_human_types")

type_zymo_fact <- paste0(colData(all_human_types)[["condition"]], "_",
                         colData(all_human_types)[["macrophagezymodeme"]])
type_zymo <- set_conditions(all_human_types, fact = type_zymo_fact)

## The numbers of samples by condition are:

## 
## Macrophages_none  Macrophages_z22  Macrophages_z23        U937_none 
##                8               23               23                2 
##         U937_z22         U937_z23 
##                6                6

data_structures <- c(data_structures, "type_zymo")

type_drug_fact <- paste0(colData(all_human_types)[["condition"]], "_",
                         colData(all_human_types)[["drug"]])
type_drug <- set_conditions(all_human_types, fact = type_drug_fact)

## The numbers of samples by condition are:

## 
## Macrophages_antimony     Macrophages_none        U937_antimony 
##                   27                   27                    7 
##            U937_none 
##                    7

data_structures <- c(data_structures, "type_drug")

strain_fact <- colData(all_human_types)[["strainid"]]
table(strain_fact)

## strain_fact
## 10763 10772 10977 11026 11075 11126 12251 12309 12355 12367  2169  7158  none 
##     2     8     2     2     2     8     7     8     2     7     8     2    10

new_conditions <- paste0(colData(hs_macrophage)[["macrophagetreatment"]], "_",
                         colData(hs_macrophage)[["macrophagezymodeme"]])
## Note the sanitize() call is redundant with the addition of sanitize() in the
## datastructures file, but I don't want to wait to rerun that.
hs_macr <- set_conditions(hs_macrophage, fact = new_conditions) %>%
  sanitize_metadata(column = "drug") %>%
  set_se_colors(color_choices[["treatment_zymo"]]) %>%
  subset_se(subset = "typeofcells!='U937'")

## The numbers of samples by condition are:

## 
##    inf_sb_z22    inf_sb_z23       inf_z22       inf_z23    uninf_none 
##            15            14            14            15             5 
## uninf_sb_none 
##             5

## Recasting the data.frame to DataFrame.

data_structures <- c(data_structures, "hs_macr")

ggstats_parasite <- plot_metadata_factors(hs_macr, column = "parasitemappingrate",
                                         type = "ggstats", scale = "log2")
pp(file = "images/ggstats_parasiterate_all_macrophage_drug_treatment.png")
ggstats_parasite
dev.off()

## png 
##   2

ggstats_parasite

hs_macr_drug_expt <- set_conditions(hs_macr, fact = "drug")

## The numbers of samples by condition are:

## 
## antimony     none 
##       27       27

hs_macr_strain_expt <- set_conditions(hs_macr, fact = "macrophagezymodeme") %>%
  subset_se(subset = "macrophagezymodeme != 'none'")

## The numbers of samples by condition are:

## 
## none  z22  z23 
##    8   23   23

data_structures <- c(data_structures, "hs_macr_strain_expt")

table(colData(hs_macr)[["strainid"]])

## 
## 10763 10772 10977 11026 11075 11126 12251 12309 12355 12367  2169  7158  none 
##     2     6     2     2     2     6     5     6     2     5     6     2     8

Let us see if the sankey plot of these samples looks useful…

ggstats_slreads <- plot_metadata_factors(hs_macrophage, column = "hisatlpsinglemapped",
                                         type = "ggstats", scale = "log2")
pp(file = "images/ggstats_slreads_all_macrophage.png")
ggstats_slreads
dev.off()

## png 
##   2

ggstats_slreads

ggstats_violin <- plot_metadata_factors(hs_macrophage, column = "hisatlpsinglemapped",
                                        scale = "log2")
ggstats_violin

macr_sankey <- plot_meta_sankey(hs_macrophage, color_choices = color_choices,
                                factors = c("oldnew", "drug", "infectedp", "macrophagezymodeme"))
macr_sankey

## A sankey plot describing the metadata of 68 samples,
## including 26 out of 0 nodes and traversing metadata factors:
## oldnew, drug, infectedp, macrophagezymodeme.

Finally, split off the U937 samples.

hs_u937 <- subset_se(hs_macrophage, subset = "typeofcells!='Macrophages'")
data_structures <- c(data_structures, "hs_u937")

5.4 Macrophage parasite data

In the previous block, we used a new invocation of ensembl-derived annotation data, this time we can just use our existing parasite gene annotations.

lp_macrophage <- create_se(
  sample_sheet, file_column = "lpanamensisv36hisatfile", gene_info = all_lp_annot,
  savefile = glue("rda/lp_macrophage-v{ver}.rda"),
  annotation = "org.Lpanamensis.MHOMCOL81L13.v46.eg.db") %>%
  set_conditions(fact = "macrophagezymodeme") %>%
  set_batches(fact = "macrophagetreatment")

## Reading the sample metadata.

## 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: 69 rows(samples) and 86 columns(metadata fields).

## Warning in create_se(sample_sheet, file_column = "lpanamensisv36hisatfile", :
## Some samples were removed when cross referencing the samples against the count
## data.

## Matched 8778 annotations and counts.

## The final summarized experiment has 8778 rows and 86 columns.

## The numbers of samples by condition are:

## 
## none z2.2 z2.3 
##    8   29   29

## The number of samples by batch are:

## 
##      inf   inf_sb    uninf uninf_sb 
##       29       29        4        4

unfilt_written <- write_se(
  lp_macrophage,
  excel = glue("analyses/macrophage_de/{ver}/read_counts/lp_macrophage_reads_unfiltered-v{ver}.xlsx"))

## Writing the first sheet, containing a legend and some summary data.

## Warning in .local(x, row.names, optional, ...): arguments in '...' ignored

## The following samples have less than 5705.7 genes.

##  [1] "TMRC30066" "TMRC30117" "TMRC30244" "TMRC30246" "TMRC30249" "TMRC30266"
##  [7] "TMRC30268" "TMRC30292" "TMRC30300" "TMRC30302" "TMRC30304" "TMRC30309"
## [13] "TMRC30312" "TMRC30313" "TMRC30319" "TMRC30323" "TMRC30325" "TMRC30326"
## [19] "TMRC30327" "TMRC30330" "TMRC30331" "TMRC30332"

## 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.
## 175550 entries are 0.  We are on a log scale, adding 1 to the data.
## 
## Naively calculating coefficient of variation/dispersion with respect to condition.
## 
## Finished calculating dispersion estimates.
## 
## `geom_smooth()` using formula = 'y ~ x'
## This dataset does not support lmer with condition+batch
## 
## Running normalize_se.
## 
## Removing 0 low-count genes (8778 remaining).
## 
## transform_counts: Found 175550 values equal to 0, adding 1 to the matrix.
## 
## `geom_smooth()` using formula = 'y ~ x'
## The factor none has 8 rows.
## 
## The factor z2.2 has 29 rows.
## 
## The factor z2.3 has 29 rows.

lp_macrophage_filt <- subset_se(lp_macrophage, nonzero = 2500) %>%
  semantic_filter(semantic = c("amastin", "gp63", "leishmanolysin"),
                  semantic_column = "annot_gene_product")

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

## TMRC30066 TMRC30117 TMRC30244 TMRC30246 TMRC30266 TMRC30268 TMRC30304 TMRC30309 
##      3080      1147      1662      2834       822      3444       289       188 
## TMRC30312 TMRC30313 TMRC30319 TMRC30323 TMRC30325 TMRC30326 TMRC30327 TMRC30330 
##        76        96       374        84       356       375       129       181 
## TMRC30066 TMRC30117 TMRC30244 TMRC30246 TMRC30266 TMRC30268 TMRC30304 TMRC30309 
##      1890       888      1135      1796       649      1915       207       166 
## TMRC30312 TMRC30313 TMRC30319 TMRC30323 TMRC30325 TMRC30326 TMRC30327 TMRC30330 
##        76        84       270        74       279       303       123       135

## Samples removed: 1890, 888, 1135, 1796, 649, 1915, 207, 166, 76, 84, 270, 74, 279, 303, 123, 135
## semantic_expt_filter(): Removed 68 genes.

data_structures <- c(data_structures, "lp_macrophage", "lp_macrophage_filt")

filt_written <- write_se(lp_macrophage_filt,
  excel = glue("analyses/macrophage_de/{ver}/read_counts/lp_macrophage_reads_filtered-v{ver}.xlsx"))

## Writing the first sheet, containing a legend and some summary data.

## Warning in .local(x, row.names, optional, ...): arguments in '...' ignored

## The following samples have less than 5661.5 genes.

## [1] "TMRC30249" "TMRC30292" "TMRC30300" "TMRC30302" "TMRC30331" "TMRC30332"

## 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.
## 44583 entries are 0.  We are on a log scale, adding 1 to the data.
## 
## Naively calculating coefficient of variation/dispersion with respect to condition.
## 
## Finished calculating dispersion estimates.
## 
## `geom_smooth()` using formula = 'y ~ x'
## Running normalize_se.
## 
## Removing 0 low-count genes (8710 remaining).
## 
## transform_counts: Found 44583 values equal to 0, adding 1 to the matrix.
## 
## `geom_smooth()` using formula = 'y ~ x'
## The factor z2.2 has 21 rows.
## 
## The factor z2.3 has 29 rows.

lp_macrophage <- lp_macrophage_filt

lp_macrophage_nosb <- subset_se(lp_macrophage, subset="batch!='inf_sb'")
lp_nosb_write <- write_se(
  lp_macrophage_nosb,
  excel = glue("analyses/macrophage_de/{ver}/read_counts/lp_macrophage_nosb_reads-v{ver}.xlsx"))

## Writing the first sheet, containing a legend and some summary data.

## Warning in .local(x, row.names, optional, ...): arguments in '...' ignored

## 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.
## 6396 entries are 0.  We are on a log scale, adding 1 to the data.
## 
## Naively calculating coefficient of variation/dispersion with respect to condition.
## 
## Finished calculating dispersion estimates.
## 
## `geom_smooth()` using formula = 'y ~ x'
## The dataset has a minimal or missing set of conditions/batches.
## 
## Running normalize_se.
## 
## Removing 119 low-count genes (8591 remaining).
## 
## transform_counts: Found 3163 values equal to 0, adding 1 to the matrix.
## 
## `geom_smooth()` using formula = 'y ~ x'
## The factor z2.2 has 14 rows.
## 
## The factor z2.3 has 15 rows.

data_structures <- c(data_structures, "lp_macrophage_nosb")

spec <- make_rnaseq_spec()
test <- gather_preprocessing_metadata(sample_sheet, specification = spec)

## Did not find the condition column in the sample sheet.
## Filling it in as undefined.
## Did not find the batch column in the sample sheet.
## Filling it in as undefined.
## Checking the state of the condition column.
## Checking the state of the batch column.
## Checking the condition factor.

## Warning in gather_preprocessing_metadata(sample_sheet, specification = spec):
## Column: trimomatic_input already exists, replacing it.

## Warning in gather_preprocessing_metadata(sample_sheet, specification = spec):
## Column: trimomatic_output already exists, replacing it.

## Warning in gather_preprocessing_metadata(sample_sheet, specification = spec):
## Column: trimomatic_percent already exists, replacing it.

## Warning in gather_preprocessing_metadata(sample_sheet, specification = spec):
## Column: fastqc_pct_gc already exists, replacing it.

## Warning in gather_preprocessing_metadata(sample_sheet, specification = spec):
## Column: hisat_genome_single_concordant already exists, replacing it.

## Warning in gather_preprocessing_metadata(sample_sheet, specification = spec):
## Column: hisat_genome_multi_concordant already exists, replacing it.

## Warning in gather_preprocessing_metadata(sample_sheet, specification = spec):
## Column: hisat_genome_single_all already exists, replacing it.

## Warning in gather_preprocessing_metadata(sample_sheet, specification = spec):
## Column: hisat_genome_multi_all already exists, replacing it.

## Warning in gather_preprocessing_metadata(sample_sheet, specification = spec):
## Column: hisat_count_table already exists, replacing it.

## preprocessing/TMRC30051/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30057/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30059/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30060/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30061/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30062/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30063/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30064/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30065/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30066/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30067/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30069/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30117/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30162/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30243/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30244/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30245/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30246/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30247/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30248/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30249/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30250/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30251/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30252/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30266/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30267/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30268/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30286/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30291/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30292/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30293/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30294/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30295/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30296/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30297/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30298/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30299/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30300/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30301/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30302/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30303/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30304/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30305/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30306/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30307/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30308/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30309/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30310/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30311/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30312/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30313/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30314/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30315/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30316/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30317/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30318/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30319/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30320/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30321/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30322/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30323/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30324/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30325/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30326/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30327/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30328/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30330/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30331/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz
## preprocessing/TMRC30332/outputs/*hisat*_*/*_*genome*_gene_ID_fcounts.csv.xz

## Writing new metadata to: sample_sheets/macrophage_samples_modified.xlsx

6 Save all data structures into one rda

found_idx <- data_structures %in% ls()
if (sum(!found_idx) > 0) {
  not_found <- data_structures[!found_idx]
  warning("Some datastructures were not generated: ", toString(not_found), ".")
  data_structures <- data_structures[found_idx]
}
save(list = data_structures, file = glue("rda/tmrc2_data_structures-v{ver}.rda"))

pander::pander(sessionInfo())

## Warning: Your system is mis-configured: '/etc/localtime' is not a symlink

## Warning: It is strongly recommended to set envionment variable TZ to
## 'America/New_York' (or equivalent)

R version 4.5.0 (2025-04-11)

Platform: x86_64-pc-linux-gnu

locale: C

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

other attached packages: ruv(v.0.9.7.1), BiocParallel(v.1.42.1), variancePartition(v.1.38.1), org.Lpanamensis.MHOMCOL81L13.v68.eg.db(v.2024.09), AnnotationDbi(v.1.70.0), IRanges(v.2.42.0), S4Vectors(v.0.46.0), Biobase(v.2.68.0), BiocGenerics(v.0.54.0), generics(v.0.1.4), dplyr(v.1.1.4), Heatplus(v.3.16.0), ggplot2(v.3.5.2), glue(v.1.8.0) and hpgltools(v.1.2)

loaded via a namespace (and not attached): fs(v.1.6.6), matrixStats(v.1.5.0), bitops(v.1.0-9), doParallel(v.1.0.17), httr(v.1.4.7), RColorBrewer(v.1.1-3), insight(v.1.4.1), numDeriv(v.2016.8-1.1), tools(v.4.5.0), backports(v.1.5.0), R6(v.2.6.1), statsExpressions(v.1.7.1), lazyeval(v.0.2.2), mgcv(v.1.9-3), withr(v.3.0.2), gridExtra(v.2.3), prettyunits(v.1.2.0), cli(v.3.6.5), performance(v.0.15.1), labeling(v.0.4.3), sass(v.0.4.10), prismatic(v.1.1.2), BWStest(v.0.2.3), mvtnorm(v.1.3-3), readr(v.2.1.5), genefilter(v.1.90.0), pbapply(v.1.7-4), Rsamtools(v.2.24.0), yulab.utils(v.0.2.1), DOSE(v.4.2.0), R.utils(v.2.13.0), dichromat(v.2.0-0.1), limma(v.3.64.3), RSQLite(v.2.4.3), BiocIO(v.1.18.0), vroom(v.1.6.5), gtools(v.3.9.5), zip(v.2.3.3), GO.db(v.3.21.0), Matrix(v.1.7-3), abind(v.1.4-8), R.methodsS3(v.1.8.2), lifecycle(v.1.0.4), yaml(v.2.3.10), edgeR(v.4.6.3), SummarizedExperiment(v.1.38.1), gplots(v.3.2.0), qvalue(v.2.40.0), SparseArray(v.1.8.1), BiocFileCache(v.2.16.1), Rtsne(v.0.17), paletteer(v.1.6.0), grid(v.4.5.0), blob(v.1.2.4), promises(v.1.3.3), crayon(v.1.5.3), lattice(v.0.22-7), cowplot(v.1.2.0), GenomicFeatures(v.1.60.0), annotate(v.1.86.1), KEGGREST(v.1.48.1), zeallot(v.0.2.0), pillar(v.1.11.0), knitr(v.1.50), varhandle(v.2.0.6), fgsea(v.1.34.2), GenomicRanges(v.1.60.0), rjson(v.0.2.23), boot(v.1.3-31), corpcor(v.1.6.10), kSamples(v.1.2-12), codetools(v.0.2-20), fastmatch(v.1.1-6), data.table(v.1.17.8), vctrs(v.0.6.5), png(v.0.1-8), Rdpack(v.2.6.4), gtable(v.0.3.6), rematch2(v.2.1.2), datawizard(v.1.2.0), cachem(v.1.1.0), xfun(v.0.53), openxlsx(v.4.2.8), rbibutils(v.2.3), S4Arrays(v.1.8.1), mime(v.0.13), correlation(v.0.8.8), reformulas(v.0.4.1), coda(v.0.19-4.1), survival(v.3.8-3), iterators(v.1.0.14), statmod(v.1.5.0), gmp(v.0.7-5), directlabels(v.2025.6.24), nlme(v.3.1-168), pbkrtest(v.0.5.5), bit64(v.4.6.0-1), EnvStats(v.3.1.0), progress(v.1.2.3), filelock(v.1.0.3), GenomeInfoDb(v.1.44.2), bslib(v.0.9.0), KernSmooth(v.2.23-26), DBI(v.1.2.3), DESeq2(v.1.48.1), tidyselect(v.1.2.1), bit(v.4.6.0), compiler(v.4.5.0), curl(v.7.0.0), httr2(v.1.2.1), graph(v.1.86.0), xml2(v.1.4.0), DelayedArray(v.0.34.1), plotly(v.4.11.0), bayestestR(v.0.17.0), rtracklayer(v.1.68.0), scales(v.1.4.0), caTools(v.1.18.3), remaCor(v.0.0.20), quadprog(v.1.5-8), multcompView(v.0.1-10), rappdirs(v.0.3.3), stringr(v.1.5.1), digest(v.0.6.37), ggsankey(v.0.0.99999), minqa(v.1.2.8), rmarkdown(v.2.29), aod(v.1.3.3), XVector(v.0.48.0), RhpcBLASctl(v.0.23-42), htmltools(v.0.5.8.1), pkgconfig(v.2.0.3), lme4(v.1.1-37), MatrixGenerics(v.1.20.0), dbplyr(v.2.5.0), fastmap(v.1.2.0), rlang(v.1.1.6), htmlwidgets(v.1.6.4), UCSC.utils(v.1.4.0), shiny(v.1.11.1), SuppDists(v.1.1-9.9), farver(v.2.1.2), jquerylib(v.0.1.4), jsonlite(v.2.0.0), GOSemSim(v.2.34.0), R.oo(v.1.27.1), RCurl(v.1.98-1.17), magrittr(v.2.0.3), GenomeInfoDbData(v.1.2.14), patchwork(v.1.3.2), parameters(v.0.28.1), Rcpp(v.1.1.0), stringi(v.1.8.7), MASS(v.7.3-65), plyr(v.1.8.9), parallel(v.4.5.0), ggrepel(v.0.9.6), Biostrings(v.2.76.0), PMCMRplus(v.1.9.12), splines(v.4.5.0), pander(v.0.6.6), hms(v.1.1.3), locfit(v.1.5-9.12), fastcluster(v.1.3.0), ggsignif(v.0.6.4), effectsize(v.1.0.1), reshape2(v.1.4.4), biomaRt(v.2.64.0), rstantools(v.2.4.0), XML(v.3.99-0.19), evaluate(v.1.0.4), RcppParallel(v.5.1.11-1), tzdb(v.0.5.0), nloptr(v.2.2.1), foreach(v.1.5.2), httpuv(v.1.6.16), MatrixModels(v.0.5-4), BayesFactor(v.0.9.12-4.7), tidyr(v.1.3.1), purrr(v.1.1.0), broom(v.1.0.9), xtable(v.1.8-4), restfulr(v.0.0.16), Rmpfr(v.1.1-1), fANCOVA(v.0.6-1), later(v.1.4.3), viridisLite(v.0.4.2), tibble(v.3.3.0), lmerTest(v.3.1-3), ggstatsplot(v.0.13.1), memoise(v.2.0.1), GenomicAlignments(v.1.44.0), sva(v.3.56.0) and GSEABase(v.1.70.0)

message("This is hpgltools commit: ", get_git_commit())

## If you wish to reproduce this exact build of hpgltools, invoke the following:

## > git clone http://github.com/abelew/hpgltools.git

## > git reset 69244d8a6c50468fcd448838c4ad80c01e44ba37

## This is hpgltools commit: Sun Aug 31 21:11:04 2025 -0400: 69244d8a6c50468fcd448838c4ad80c01e44ba37

message("Saving to ", savefile)

## Saving to 01datasets.rda.xz

# tmp <- sm(saveme(filename = savefile))

tmp <- loadme(filename = savefile)