TMRC2 Infecting macrophages: 202203

1 TODO

1. Compare these samples to previous in-vitro L.major infected macrophages.
2. Look at shared/unique genes with respect to antimonial treatment.
3. Send along rpkm counts/tables.
4. Perform enrichment analyses of high/low genes vs. DE genes.
5. Compare to MSigDB and/or GSVA

2 Introduction

3 Annotation

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

hs_annot <- load_biomart_annotations(year="2020", month="jan")

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

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

4 Sample Estimation

I used two mapping methods for this data, hisat2 and salmon. Most analyses use hisat2, which is a more traditional map-and-count method. In contrast, salmon uses what may be thought of as a indexed voting method (so that multi-matches are discounted and the votes split among all matches). Salmon also required a pre-existing database of known transcripts (though later versions may actually use mapping from things like hisat), while hisat uses the genome and a database of known transcripts (and optionally can search for splicing junctions to find new transcripts).

4.1 Generate expressionsets

The first thing to note is the large range in coverage. There are multiple samples with coverage which is too low to use. These will be removed shortly.

sanitize_columns <- c("donor", "drug", "condition", "batch", "macrophagetreatment",
                      "macrophageinfectionzymodeme")
macr <- create_expt("sample_sheets/tmrc2_macrophage_samples_202203.xlsx",
                    file_column = "hg38100hisatfile",
                    savefile = "tmrc2_macrophage.rda",
                    gene_info = hs_annot) %>%
  exclude_genes_expt(column = "gene_biotype", method = "keep",
                     pattern = "protein_coding") %>%
  subset_expt(nonzero = 11000) %>%
  set_expt_factors(columns = sanitize_columns) %>%
  set_expt_factors(columns=sanitize_columns, class="factor") %>%
  set_expt_batches(fact="macrophagetreatment") %>%
  set_expt_conditions(fact="macrophageinfectionzymodeme")

## Reading the sample metadata.

## The sample definitions comprises: 28 rows(samples) and 67 columns(metadata fields).

## Matched 21447 annotations and counts.

## Bringing together the count matrix and gene information.

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

## The final expressionset has 21481 rows and 28 columns.

## Before removal, there were 21481 genes, now there are 19923.

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

## named numeric(0)

## subset_expt(): There were 28, now there are 28 samples.

5 Macrophages

These samples are rather different from all of the others. The following section is therefore written primarily in response to a separate set of emails from Olga and Maria Adelaida; here is a snippet:

Dear all, about the samples corresponding to infected macrophages with three sensitive (2.2) and three resistant (2.3) clinical strains of L. (V.) panamensis, I send you the results of parasite burden by detection of 7SLRNA. I think these results are interesting, but the sample size is very small. Doctor Najib or Trey could you please send me the quality data and PCA analysis of these samples?

and

Hi Doctor, thank you. These samples corresponding to primary macrophages infected with clinical strains 2.2 (n=3) and 2.3 (n = 3). These information is in the file: TMRC project 3: excel Host TMRC3 v1.1, rows 137 to 150.

Thus I added 3 columns to the end of the sample sheet which seek to include this information. The first is ‘drug’ and encodes both the infection state (no for the two controls and yes for everything else), the second is zymodeme which I took from the tmrc2 sample sheet, the last is drug, which is either no or sb.

table(pData(macr)[["condition"]])

## 
##    none unknown    z2.2    z2.3 
##       4       2      10      12

table(pData(macr)[["batch"]])

## 
##      inf   inf_sb    uninf uninf_sb 
##       12       12        2        2

table(pData(macr)[["macrophagetreatment"]])

## 
##      inf   inf_sb    uninf uninf_sb 
##       12       12        2        2

table(pData(macr)[["macrophageinfectionzymodeme"]])

## 
##    none unknown    z2.2    z2.3 
##       4       2      10      12

macr_libsize <- plot_libsize(macr, yscale=TRUE)
pp(file="macrophage_side_experiment/macr_libsize.png", image=macr_libsize$plot)

## Warning in pp(file = "macrophage_side_experiment/macr_libsize.png", image =
## macr_libsize$plot): There is no device to shut down.

macr_nonzero <- plot_nonzero(macr)
pp(file="macrophage_side_experiment/macr_nonzero.png", image=macr_nonzero$plot)

## Warning in pp(file = "macrophage_side_experiment/macr_nonzero.png", image =
## macr_nonzero$plot): There is no device to shut down.

macr_norm <- sm(normalize_expt(macr, transform="log2", norm="quant", convert="cpm", filter=TRUE))
norm_pca <- plot_pca(macr_norm, plot_labels=FALSE)
norm_pca$plot

pp(file="macrophage_side_experiment/pca_colors_zymodeme.png", image=norm_pca$plot)

## Warning in pp(file = "macrophage_side_experiment/pca_colors_zymodeme.png", :
## There is no device to shut down.

treatment_macr <- set_expt_conditions(macr, fact="macrophagetreatment") %>%
  set_expt_batches(fact="macrophageinfectionzymodeme")
treatment_norm <- sm(normalize_expt(treatment_macr, transform="log2", norm="quant", convert="cpm", filter=TRUE))
treatment_pca <- plot_pca(treatment_norm, plot_labels=FALSE)
treatment_pca$plot
pp(file="macrophage_side_experiment/pca_colors_treatment.png", image=treatment_pca$plot)

## Warning in pp(file = "macrophage_side_experiment/pca_colors_treatment.png", :
## There is no device to shut down.

plot_3d_pca(norm_pca)$plot

macr_nb <- sm(normalize_expt(macr, norm = "quant", convert = "cpm",
                             transform = "log2", batch = "svaseq", filter = TRUE))
nb_pca <- plot_pca(macr_nb, plot_labels=FALSE)
pp(file="macrophage_side_experiment/normbatch_pca.png", image=nb_pca$plot)

## Warning in pp(file = "macrophage_side_experiment/normbatch_pca.png", image =
## nb_pca$plot): There is no device to shut down.

5.1 Write the data

macr_written <- sm(write_expt(macr, excel="macrophage_side_experiment/macrophage_expt.xlsx"))

## Error in checkModelStatus(fit, showWarnings = showWarnings, colinearityCutoff = colinearityCutoff,  : 
##   The variables specified in this model are redundant,
## so the design matrix is not full rank

## Error in checkModelStatus(fit, showWarnings = showWarnings, colinearityCutoff = colinearityCutoff,  : 
##   The variables specified in this model are redundant,
## so the design matrix is not full rank

5.2 Perform DE

There were a couple of contrasts explicitly requested for this data:

• Macrophages infected with 2.3 without SbV vs macrophagues infected with 2.2 without SbV
• Macrophages infected with 2.3 without SbV vs macrophages uninfected (M0)

Unfortunately, I think to really answer these questions we will require more replicates of a few of these conditions, most notably the uninfected samples.

treatment_macr <- set_expt_conditions(macr, fact="macrophagetreatment") %>%
  set_expt_batches(fact="macrophageinfectionzymodeme")
##tmp <- normalize_expt(macr, filter=TRUE)
##zy_de <- deseq_pairwise(tmp, model_batch="svaseq")
##zy_ed <- edger_pairwise(tmp, model_batch="svaseq")
combined_condition <- paste0(pData(macr)[["macrophagetreatment"]], "_",
                             pData(macr)[["macrophageinfectionzymodeme"]])
pData(macr[["expressionset"]])[["combined"]] <- combined_condition
macr <- set_expt_conditions(macr, fact = "combined")
zymo_de <- all_pairwise(macr, model_batch = "svaseq", filter = TRUE,
                        do_ebseq = FALSE)

## batch_counts: Before batch/surrogate estimation, 697 entries are x==0: 0%.

## Plotting a PCA before surrogate/batch inclusion.

## Using svaseq to visualize before/after batch inclusion.

## Performing a test normalization with: raw

## Removing 0 low-count genes (11425 remaining).

## batch_counts: Before batch/surrogate estimation, 697 entries are x==0: 0%.

## batch_counts: Before batch/surrogate estimation, 15145 entries are 0<x<1: 5%.

## Setting 785 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

keepers <- list(
    "z22_uninf" = c("infz22", "uninfnone"),
    "z23_uninf" = c("infz23", "uninfnone"),
    "z22_sbuninf" = c("infsbz22", "uninfsbnone"),
    "z23_sbuninf" = c("infsbz23", "uninfsbnone"),
    "zymo_nodrug" = c("infz23", "infz22"),
    "z23_uninf" = c("infz23", "uninfnozymo"))
zymo_table <- combine_de_tables(
    zymo_de, keepers=keepers,
    excel="macrophage_side_experiment/macrophage_de.xlsx")

pp(file="images/z22_vs_uninfected_ma.png", image=zymo_table$plots$z22_uninf$deseq_ma_plots$plot)

## Warning in pp(file = "images/z22_vs_uninfected_ma.png", image =
## zymo_table$plots$z22_uninf$deseq_ma_plots$plot): There is no device to shut
## down.

pp(file="images/z22_vs_sb_uninfected_ma.png", image=zymo_table$plots$z22_sbuninf$deseq_ma_plots$plot)

## Warning in pp(file = "images/z22_vs_sb_uninfected_ma.png", image =
## zymo_table$plots$z22_sbuninf$deseq_ma_plots$plot): There is no device to shut
## down.

pp(file="images/z23_vs_uninfected_ma.png", image=zymo_table$plots$z23_uninf$deseq_ma_plots$plot)

## Warning in pp(file = "images/z23_vs_uninfected_ma.png", image =
## zymo_table$plots$z23_uninf$deseq_ma_plots$plot): There is no device to shut
## down.

pp(file="images/z23_vs_sb_uninfected_ma.png", image=zymo_table$plots$z23_sbuninf$deseq_ma_plots$plot)

## Warning in pp(file = "images/z23_vs_sb_uninfected_ma.png", image =
## zymo_table$plots$z23_sbuninf$deseq_ma_plots$plot): There is no device to shut
## down.

zymo_sig <- extract_significant_genes(zymo_table, according_to="deseq")

z22_ma <- plot_ma_de(zymo_table[["data"]][["z22_uninf"]], expr_col="deseq_basemean",
                 fc_col="deseq_logfc", p_col="deseq_adjp")
z22_ma$plot
z22_uninf_up_sig <- zymo_sig[["deseq"]][["ups"]][["z22_uninf"]]
dim(z22_uninf_up_sig)

## [1] 398  47

z22_uninf_down_sig <- zymo_sig[["deseq"]][["downs"]][["z22_uninf"]]
dim(z22_uninf_down_sig)

## [1]  5 47

z22_sbuninf_up_sig <- zymo_sig[["deseq"]][["ups"]][["z22_sbuninf"]]
dim(z22_sbuninf_up_sig)

## [1] 11 47

z22_sbuninf_down_sig <- zymo_sig[["deseq"]][["downs"]][["z22_sbuninf"]]
dim(z22_sbuninf_down_sig)

## [1]  0 47

z23_ma <- plot_ma_de(zymo_table[["data"]][["z23_uninf"]], expr_col="deseq_basemean",
                     fc_col="deseq_logfc", p_col="deseq_adjp")
z23_ma$plot
z23_uninf_up_sig <- zymo_sig[["deseq"]][["ups"]][["z23_uninf"]]
dim(z23_uninf_up_sig)

## [1] 615  47

z23_uninf_down_sig <- zymo_sig[["deseq"]][["downs"]][["z23_uninf"]]
dim(z23_uninf_down_sig)

## [1] 358  47

z23_sbuninf_up_sig <- zymo_sig[["deseq"]][["ups"]][["z23_sbuninf"]]
dim(z23_sbuninf_up_sig)

## [1] 529  47

z23_sbuninf_down_sig <- zymo_sig[["deseq"]][["downs"]][["z23_sbuninf"]]
dim(z23_sbuninf_down_sig)

## [1] 333  47

5.3 Compare groups of genes

z22_z23_nosb_groups <- list(
    "z22_up" = rownames(z22_uninf_up_sig),
    "z23_up" = rownames(z23_uninf_up_sig))
z22_z23_nosb_venn <- Vennerable::Venn(Sets = z22_z23_nosb_groups)
pp(file="images/z22_z23_nosb_venn.png")

## NULL

z22_z23_nosb_venn_plot <- Vennerable::plot(z22_z23_nosb_venn, doWeights = FALSE)
dev.off()

## png 
##   2

annot <- fData(macr)
z23_nosb_unique <- z22_z23_nosb_venn@IntersectionSets[["01"]]
wanted <- rownames(annot) %in% z23_nosb_unique
z23_wanted <- annot[wanted, ]
written <- write_xlsx(z23_wanted, excel="excel/z23_nosb_unique_genes.xlsx")

z22_nosb_unique <- z22_z23_nosb_venn@IntersectionSets[["10"]]
wanted <- rownames(annot) %in% z22_nosb_unique
z22_wanted <- annot[wanted, ]
written <- write_xlsx(z22_wanted, excel="excel/z22_nosb_unique_genes.xlsx")

z22_nosb_shared <- z22_z23_nosb_venn@IntersectionSets[["11"]]
wanted <- rownames(annot) %in% z22_nosb_shared
shared_wanted <- annot[wanted, ]
written <- write_xlsx(shared_wanted, excel="excel/z22_z23_nosb_shared_genes.xlsx")

z22_unique_gp <- simple_gprofiler(rownames(z22_wanted))

## Performing gProfiler GO search of 311 genes against hsapiens.

## GO search found 70 hits.

## Performing gProfiler KEGG search of 311 genes against hsapiens.

## KEGG search found 3 hits.

## Performing gProfiler REAC search of 311 genes against hsapiens.

## REAC search found 13 hits.

## Performing gProfiler MI search of 311 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 311 genes against hsapiens.

## TF search found 32 hits.

## Performing gProfiler CORUM search of 311 genes against hsapiens.

## CORUM search found 1 hits.

## Performing gProfiler HP search of 311 genes against hsapiens.

## HP search found 0 hits.

pp(file="images/z22_unique_bp.png", image=z22_unique_gp[["pvalue_plots"]][["bpp_plot_over"]], height=9, width=6)

## Warning in pp(file = "images/z22_unique_bp.png", image =
## z22_unique_gp[["pvalue_plots"]][["bpp_plot_over"]], : There is no device to shut
## down.

pp(file="images/z22_unique_kegg.png", image=z22_unique_gp[["pvalue_plots"]][["kegg_plot_over"]])

## Warning in pp(file = "images/z22_unique_kegg.png", image =
## z22_unique_gp[["pvalue_plots"]][["kegg_plot_over"]]): There is no device to shut
## down.

pp(file="images/z22_unique_reac.png", image=z22_unique_gp[["pvalue_plots"]][["reactome_plot_over"]])

## Warning in pp(file = "images/z22_unique_reac.png", image =
## z22_unique_gp[["pvalue_plots"]][["reactome_plot_over"]]): There is no device to
## shut down.

z23_unique_gp <- simple_gprofiler(rownames(z23_wanted))

## Performing gProfiler GO search of 528 genes against hsapiens.

## GO search found 468 hits.

## Performing gProfiler KEGG search of 528 genes against hsapiens.

## KEGG search found 21 hits.

## Performing gProfiler REAC search of 528 genes against hsapiens.

## REAC search found 31 hits.

## Performing gProfiler MI search of 528 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 528 genes against hsapiens.

## TF search found 92 hits.

## Performing gProfiler CORUM search of 528 genes against hsapiens.

## CORUM search found 2 hits.

## Performing gProfiler HP search of 528 genes against hsapiens.

## HP search found 0 hits.

pp(file="images/z23_unique_bp.png", image=z23_unique_gp[["pvalue_plots"]][["bpp_plot_over"]], height=9, width=6)

## Warning in pp(file = "images/z23_unique_bp.png", image =
## z23_unique_gp[["pvalue_plots"]][["bpp_plot_over"]], : There is no device to shut
## down.

pp(file="images/z23_unique_kegg.png", image=z23_unique_gp[["pvalue_plots"]][["kegg_plot_over"]])

## Warning in pp(file = "images/z23_unique_kegg.png", image =
## z23_unique_gp[["pvalue_plots"]][["kegg_plot_over"]]): There is no device to shut
## down.

pp(file="images/z23_unique_reac.png", image=z23_unique_gp[["pvalue_plots"]][["reactome_plot_over"]])

## Warning in pp(file = "images/z23_unique_reac.png", image =
## z23_unique_gp[["pvalue_plots"]][["reactome_plot_over"]]): There is no device to
## shut down.

pp(file="images/z23_unique_tf.png", image=z23_unique_gp[["pvalue_plots"]][["tf_plot_over"]])

## Warning in pp(file = "images/z23_unique_tf.png", image =
## z23_unique_gp[["pvalue_plots"]][["tf_plot_over"]]): There is no device to shut
## down.

z22_z23_shared_gp <- simple_gprofiler(rownames(shared_wanted))

## Performing gProfiler GO search of 87 genes against hsapiens.

## GO search found 54 hits.

## Performing gProfiler KEGG search of 87 genes against hsapiens.

## KEGG search found 8 hits.

## Performing gProfiler REAC search of 87 genes against hsapiens.

## REAC search found 10 hits.

## Performing gProfiler MI search of 87 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 87 genes against hsapiens.

## TF search found 9 hits.

## Performing gProfiler CORUM search of 87 genes against hsapiens.

## CORUM search found 2 hits.

## Performing gProfiler HP search of 87 genes against hsapiens.

## HP search found 0 hits.

pp(file="images/z23_shared_bp.png", image=z22_z23_shared_gp[["pvalue_plots"]][["bpp_plot_over"]], height=9, width=6)

## Warning in pp(file = "images/z23_shared_bp.png", image =
## z22_z23_shared_gp[["pvalue_plots"]][["bpp_plot_over"]], : There is no device to
## shut down.

pp(file="images/z23_shared_kegg.png", image=z22_z23_shared_gp[["pvalue_plots"]][["kegg_plot_over"]])

## Warning in pp(file = "images/z23_shared_kegg.png", image =
## z22_z23_shared_gp[["pvalue_plots"]][["kegg_plot_over"]]): There is no device to
## shut down.

pp(file="images/z23_shared_reac.png", image=z22_z23_shared_gp[["pvalue_plots"]][["reactome_plot_over"]])

## Warning in pp(file = "images/z23_shared_reac.png", image =
## z22_z23_shared_gp[["pvalue_plots"]][["reactome_plot_over"]]): There is no device
## to shut down.

pp(file="images/z23_shared_tf.png", image=z22_z23_shared_gp[["pvalue_plots"]][["tf_plot_over"]])

## Warning in pp(file = "images/z23_shared_tf.png", image =
## z22_z23_shared_gp[["pvalue_plots"]][["tf_plot_over"]]): There is no device to
## shut down.

z22_uninf_up_gp <- simple_gprofiler(z22_uninf_up_sig)

## Performing gProfiler GO search of 398 genes against hsapiens.

## GO search found 117 hits.

## Performing gProfiler KEGG search of 398 genes against hsapiens.

## KEGG search found 8 hits.

## Performing gProfiler REAC search of 398 genes against hsapiens.

## REAC search found 18 hits.

## Performing gProfiler MI search of 398 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 398 genes against hsapiens.

## TF search found 51 hits.

## Performing gProfiler CORUM search of 398 genes against hsapiens.

## CORUM search found 1 hits.

## Performing gProfiler HP search of 398 genes against hsapiens.

## HP search found 1 hits.

pp(file="images/z22_uninf_up_gp_bp.png", image=z22_uninf_up_gp$pvalue_plots$bpp_plot_over, height=20, width=10)

## Warning in pp(file = "images/z22_uninf_up_gp_bp.png", image =
## z22_uninf_up_gp$pvalue_plots$bpp_plot_over, : There is no device to shut down.

z22_uninf_down_gp <- simple_gprofiler(z22_uninf_down_sig)

## Performing gProfiler GO search of 5 genes against hsapiens.

## GO search found 0 hits.

## Performing gProfiler KEGG search of 5 genes against hsapiens.

## KEGG search found 0 hits.

## Performing gProfiler REAC search of 5 genes against hsapiens.

## REAC search found 0 hits.

## Performing gProfiler MI search of 5 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 5 genes against hsapiens.

## TF search found 0 hits.

## Performing gProfiler CORUM search of 5 genes against hsapiens.

## CORUM search found 0 hits.

## Performing gProfiler HP search of 5 genes against hsapiens.

## HP search found 0 hits.

z22_sbuninf_up_gp <- simple_gprofiler(z22_sbuninf_up_sig)

## Performing gProfiler GO search of 11 genes against hsapiens.

## GO search found 0 hits.

## Performing gProfiler KEGG search of 11 genes against hsapiens.

## KEGG search found 0 hits.

## Performing gProfiler REAC search of 11 genes against hsapiens.

## REAC search found 0 hits.

## Performing gProfiler MI search of 11 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 11 genes against hsapiens.

## TF search found 2 hits.

## Performing gProfiler CORUM search of 11 genes against hsapiens.

## CORUM search found 0 hits.

## Performing gProfiler HP search of 11 genes against hsapiens.

## HP search found 0 hits.

z22_sbuninf_down_gp <- simple_gprofiler(z22_sbuninf_down_sig)

## Performing gProfiler GO search of 0 genes against hsapiens.

## GO search found 0 hits.

## Performing gProfiler KEGG search of 0 genes against hsapiens.

## KEGG search found 0 hits.

## Performing gProfiler REAC search of 0 genes against hsapiens.

## REAC search found 0 hits.

## Performing gProfiler MI search of 0 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 0 genes against hsapiens.

## TF search found 0 hits.

## Performing gProfiler CORUM search of 0 genes against hsapiens.

## CORUM search found 0 hits.

## Performing gProfiler HP search of 0 genes against hsapiens.

## HP search found 0 hits.

z23_uninf_up_gp <- simple_gprofiler(z23_uninf_up_sig)

## Performing gProfiler GO search of 615 genes against hsapiens.

## GO search found 484 hits.

## Performing gProfiler KEGG search of 615 genes against hsapiens.

## KEGG search found 29 hits.

## Performing gProfiler REAC search of 615 genes against hsapiens.

## REAC search found 47 hits.

## Performing gProfiler MI search of 615 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 615 genes against hsapiens.

## TF search found 102 hits.

## Performing gProfiler CORUM search of 615 genes against hsapiens.

## CORUM search found 3 hits.

## Performing gProfiler HP search of 615 genes against hsapiens.

## HP search found 7 hits.

pp(file="images/z23_uninf_up_gp_bp.png", image=z23_uninf_up_gp$pvalue_plots$bpp_plot_over, height=20, width=10)

## Warning in pp(file = "images/z23_uninf_up_gp_bp.png", image =
## z23_uninf_up_gp$pvalue_plots$bpp_plot_over, : There is no device to shut down.

6 Parasite reads

sanitize_columns <- c("donor", "drug", "condition", "batch", "macrophagetreatment",
                      "macrophageinfectionzymodeme")
library(org.Lpanamensis.MHOMCOL81L13.v46.eg.db)

## Loading required package: AnnotationDbi

## 

pan_db <- org.Lpanamensis.MHOMCOL81L13.v46.eg.db
all_lp_annot <- sm(load_orgdb_annotations(
    pan_db,
    keytype = "gid",
    fields = c("annot_gene_entrez_id", "annot_gene_name",
               "annot_strand", "annot_chromosome", "annot_cds_length",
               "annot_gene_product")))$genes
macr_lp <- create_expt(metadata="sample_sheets/tmrc2_macrophage_samples_202203.xlsx",
                       gene_info=all_lp_annot, file_column = "lpanamensisv36hisatfile")

## Reading the sample metadata.

## The sample definitions comprises: 28 rows(samples) and 67 columns(metadata fields).

## Warning in file(file, "rt"): cannot open file '/mnt/cbcb/fs01_abelew/
## cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/
## preprocessing/TMRC30051/outputs/hisat2_lpanamensis_v36/
## r1_trimmed.count_lpanamensis_v36_sno_gene_ID.count.xz': No such file or
## directory

## Error in file(file, "rt"): cannot open the connection

  exclude_genes_expt(column = "gene_biotype", method = "keep",
                     pattern = "protein_coding") %>%
  subset_expt(subset="typeofcells=='Macrophages'") %>%
  subset_expt(nonzero = 11000) %>%
  set_expt_factors(columns = sanitize_columns) %>%
  set_expt_factors(columns=sanitize_columns, class="factor") %>%
  set_expt_batches(fact="macrophagetreatment") %>%
  set_expt_conditions(fact="macrophageinfectionzymodeme")

## 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': error in evaluating the argument 'object' in selecting a method for function 'fData': error in evaluating the argument 'object' in selecting a method for function 'fData': error in evaluating the argument 'object' in selecting a method for function 'sampleNames': error in evaluating the argument 'object' in selecting a method for function 'sampleNames': argument "expt" is missing, with no default

if (!isTRUE(get0("skip_load"))) {
  pander::pander(sessionInfo())
  message(paste0("This is hpgltools commit: ", get_git_commit()))
  message(paste0("Saving to ", savefile))
  tmp <- sm(saveme(filename=savefile))
}

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

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

## > git reset c07ff314c871617d0f3c45257d3c46e61270ed3e

## This is hpgltools commit: Thu Mar 17 11:25:56 2022 -0400: c07ff314c871617d0f3c45257d3c46e61270ed3e

## Saving to tmrc2_macrophage_202203.rda.xz

tmp <- loadme(filename=savefile)