1 Introduction
I do not yet know much of the background of these samples. My understanding is that varying titers of vaccine to a balanced set of male and female mice.
2 Annotation data
I chose to use the ~ 2020 mm38_100 genome.
mm_annot <- load_biomart_annotations(species="mmusculus")[["annotation"]]## The biomart annotations file already exists, loading from it.rownames(mm_annot) <- make.names(mm_annot[["ensembl_gene_id"]], unique=TRUE)
drop_tx <- grepl(x=rownames(mm_annot), pattern="\\.")
mm_annot <- mm_annot[!drop_tx, ]
mm_gff <- load_gff_annotations("~/libraries_fs/genome/mm38_100.gff")## Trying attempt: rtracklayer::import.gff3(gff, sequenceRegionsAsSeqinfo = TRUE)## Had a successful gff import with rtracklayer::import.gff3(gff, sequenceRegionsAsSeqinfo = TRUE)## Returning a df with 32 columns and 3899382 rows.mm_gff_idx <- !is.na(mm_gff[["gene_id"]])
mm_gff <- mm_gff[mm_gff_idx, ]
rownames(mm_gff) <- make.names(mm_gff[["gene_id"]], unique=TRUE)
annotations <- merge(mm_annot, mm_gff, by="row.names", all.x=TRUE)
rownames(annotations) <- annotations[["Row.names"]]
annotations[["Row.names"]] <- NULL
rownames(annotations) <- paste0("gene:", rownames(annotations))mm_expt <- create_expt("sample_sheets/initial_metadata_20220221.xlsx", gene_info=annotations)## Reading the sample metadata.## The sample definitions comprises: 12 rows(samples) and 12 columns(metadata fields).## Matched 25661 annotations and counts.## Bringing together the count matrix and gene information.## Some annotations were lost in merging, setting them to 'undefined'.## Saving the expressionset to 'expt.rda'.## The final expressionset has 25760 features and 12 samples.mm_norm <- normalize_expt(mm_expt, transform="log2", convert="cpm", norm="quant", filter=TRUE)## Removing 14542 low-count genes (11218 remaining).plot_pca(mm_norm)$plot
mm_nb <- normalize_expt(mm_expt, transform="log2", norm="quant", filter=TRUE, batch="combat")## Removing 14542 low-count genes (11218 remaining).## Setting 48 low elements to zero.## transform_counts: Found 48 values equal to 0, adding 1 to the matrix.plot_pca(mm_nb)$plot
3 Female only
mm_female <- subset_expt(mm_expt, subset="batch=='female'")## subset_expt(): There were 12, now there are 6 samples.mm_female_norm <- normalize_expt(mm_female, filter=TRUE, transform="log2",
norm="quant", convert="cpm")## Removing 14751 low-count genes (11009 remaining).plot_pca(mm_female_norm)$plot
female_de <- all_pairwise(mm_female, filter=TRUE, model_batch=FALSE)## This DE analysis will perform all pairwise comparisons among:##
## COPS:FliC PBS
## 3 3## This will pre-filter the input data using normalize_expt's: TRUE argument.## Finished running DE analyses, collecting outputs.## Comparing analyses.keeper <- list("vaccination" = c("COPSFliC", "PBS"))
female_table <- combine_de_tables(
female_de,
keeper=keeper,
excel="excel/mm_female_de-v202204.xlsx")## Error in check_xlsx_worksheet(wb, sheet): trying to get slot ".xData" from an object of a basic class ("NULL") with no slotsfemale_sig <- extract_significant_genes(lfc=0.6,
female_table,
excel="excel/mm_female_sig-v202204.xlsx",
according_to="deseq")## Error in extract_significant_genes(lfc = 0.6, female_table, excel = "excel/mm_female_sig-v202204.xlsx", : object 'female_table' not found4 Perform a couple of likely ontology searches
ups <- female_sig[["deseq"]][["ups"]][["vaccination"]]## Error in eval(expr, envir, enclos): object 'female_sig' not foundrownames(ups) <- gsub(x=rownames(ups), pattern="gene:", replacement="")## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'gsub': error in evaluating the argument 'x' in selecting a method for function 'rownames': object 'ups' not founddowns <- female_sig[["deseq"]][["downs"]][["vaccination"]]## Error in eval(expr, envir, enclos): object 'female_sig' not foundrownames(downs) <- gsub(x=rownames(downs), pattern="gene:", replacement="")## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'gsub': error in evaluating the argument 'x' in selecting a method for function 'rownames': object 'downs' not foundup_gp <- simple_gprofiler(ups, species="mmusculus")## Error in simple_gprofiler(ups, species = "mmusculus"): object 'ups' not foundup_gp[["pvalue_plots"]][["bpp_plot_over"]]## Error in eval(expr, envir, enclos): object 'up_gp' not foundup_gp[["pvalue_plots"]][["mfp_plot_over"]]## Error in eval(expr, envir, enclos): object 'up_gp' not foundup_gp[["pvalue_plots"]][["reactome_plot_over"]]## Error in eval(expr, envir, enclos): object 'up_gp' not foundup_gp[["pvalue_plots"]][["kegg_plot_over"]]## Error in eval(expr, envir, enclos): object 'up_gp' not foundup_gp[["pvalue_plots"]][["tf_plot_over"]]## Error in eval(expr, envir, enclos): object 'up_gp' not foundup_gp[["pvalue_plots"]][["hp_plot_over"]]## Error in eval(expr, envir, enclos): object 'up_gp' not founddown_gp <- simple_gprofiler(downs, species="mmusculus")## Error in simple_gprofiler(downs, species = "mmusculus"): object 'downs' not founddown_gp[["pvalue_plots"]][["bpp_plot_over"]]## Error in eval(expr, envir, enclos): object 'down_gp' not founddown_gp[["pvalue_plots"]][["mfp_plot_over"]]## Error in eval(expr, envir, enclos): object 'down_gp' not founddown_gp[["pvalue_plots"]][["reactome_plot_over"]]## Error in eval(expr, envir, enclos): object 'down_gp' not founddown_gp[["pvalue_plots"]][["kegg_plot_over"]]## Error in eval(expr, envir, enclos): object 'down_gp' not founddown_gp[["pvalue_plots"]][["tf_plot_over"]]## Error in eval(expr, envir, enclos): object 'down_gp' not founddown_gp[["pvalue_plots"]][["hp_plot_over"]]## Error in eval(expr, envir, enclos): object 'down_gp' not foundpp(file="images/down_gprofiler_bp.png", image=down_gp[["pvalue_plots"]][["bpp_plot_over"]])## Error in pp(file = "images/down_gprofiler_bp.png", image = down_gp[["pvalue_plots"]][["bpp_plot_over"]]): object 'down_gp' not founddev.off()## png
## 2pp(file="images/down_gprofiler_mf.png", image=down_gp[["pvalue_plots"]][["mfp_plot_over"]])## Error in pp(file = "images/down_gprofiler_mf.png", image = down_gp[["pvalue_plots"]][["mfp_plot_over"]]): object 'down_gp' not founddev.off()## png
## 2pp(file="images/down_gprofiler_reactome.png", image=down_gp[["pvalue_plots"]][["reactome_plot_over"]])## Error in pp(file = "images/down_gprofiler_reactome.png", image = down_gp[["pvalue_plots"]][["reactome_plot_over"]]): object 'down_gp' not founddev.off()## png
## 2pp(file="images/down_gprofiler_kegg.png", image=down_gp[["pvalue_plots"]][["kegg_plot_over"]])## Error in pp(file = "images/down_gprofiler_kegg.png", image = down_gp[["pvalue_plots"]][["kegg_plot_over"]]): object 'down_gp' not founddev.off()## png
## 2pp(file="images/down_gprofiler_tf.png", image=down_gp[["pvalue_plots"]][["tf_plot_over"]], height=16, width=8)## Error in pp(file = "images/down_gprofiler_tf.png", image = down_gp[["pvalue_plots"]][["tf_plot_over"]], : object 'down_gp' not founddev.off()## png
## 2pp(file="images/down_gprofiler_hp.png", image=down_gp[["pvalue_plots"]][["hp_plot_over"]], height=12, width=8)## Error in pp(file = "images/down_gprofiler_hp.png", image = down_gp[["pvalue_plots"]][["hp_plot_over"]], : object 'down_gp' not founddev.off()## png
## 2## down_cp <- simple_clusterprofiler(downs, orgdb="org.Mm.eg.db")pander::pander(sessionInfo())R version 4.1.2 (2021-11-01)
Platform: x86_64-pc-linux-gnu (64-bit)
locale: LC_CTYPE=en_US.UTF-8, LC_NUMERIC=C, LC_TIME=en_US.UTF-8, LC_COLLATE=en_US.UTF-8, LC_MONETARY=en_US.UTF-8, LC_MESSAGES=en_US.UTF-8, LC_PAPER=en_US.UTF-8, LC_NAME=C, LC_ADDRESS=C, LC_TELEPHONE=C, LC_MEASUREMENT=en_US.UTF-8 and LC_IDENTIFICATION=C
attached base packages: stats4, stats, graphics, grDevices, utils, datasets, methods and base
other attached packages: ruv(v.0.9.7.1), hpgltools(v.1.0), testthat(v.3.1.4), reticulate(v.1.25), SummarizedExperiment(v.1.24.0), GenomicRanges(v.1.46.1), GenomeInfoDb(v.1.30.1), IRanges(v.2.28.0), S4Vectors(v.0.32.4), MatrixGenerics(v.1.6.0), matrixStats(v.0.62.0), Biobase(v.2.54.0) and BiocGenerics(v.0.40.0)
loaded via a namespace (and not attached): utf8(v.1.2.2), tidyselect(v.1.1.2), lme4(v.1.1-30), RSQLite(v.2.2.15), AnnotationDbi(v.1.56.2), htmlwidgets(v.1.5.4), grid(v.4.1.2), BiocParallel(v.1.28.3), scatterpie(v.0.1.7), devtools(v.2.4.3), munsell(v.0.5.0), preprocessCore(v.1.56.0), codetools(v.0.2-18), withr(v.2.5.0), colorspace(v.2.0-3), GOSemSim(v.2.20.0), filelock(v.1.0.2), highr(v.0.9), knitr(v.1.39), rstudioapi(v.0.13), DOSE(v.3.20.1), labeling(v.0.4.2), Rdpack(v.2.3.1), GenomeInfoDbData(v.1.2.7), polyclip(v.1.10-0), farver(v.2.1.1), bit64(v.4.0.5), rprojroot(v.2.0.3), downloader(v.0.4), treeio(v.1.18.1), vctrs(v.0.4.1), generics(v.0.1.3), xfun(v.0.31), BiocFileCache(v.2.2.1), R6(v.2.5.1), doParallel(v.1.0.17), graphlayouts(v.0.8.0), locfit(v.1.5-9.6), gridGraphics(v.0.5-1), bitops(v.1.0-7), cachem(v.1.0.6), fgsea(v.1.20.0), DelayedArray(v.0.20.0), assertthat(v.0.2.1), promises(v.1.2.0.1), BiocIO(v.1.4.0), scales(v.1.2.0), ggraph(v.2.0.5), enrichplot(v.1.14.2), gtable(v.0.3.0), sva(v.3.42.0), processx(v.3.7.0), tidygraph(v.1.2.1), rlang(v.1.0.4), genefilter(v.1.76.0), splines(v.4.1.2), rtracklayer(v.1.54.0), lazyeval(v.0.2.2), broom(v.1.0.0), yaml(v.2.3.5), reshape2(v.1.4.4), GenomicFeatures(v.1.46.5), backports(v.1.4.1), httpuv(v.1.6.5), qvalue(v.2.26.0), clusterProfiler(v.4.2.2), tools(v.4.1.2), usethis(v.2.1.6), ggplotify(v.0.1.0), ggplot2(v.3.3.6), ellipsis(v.0.3.2), gplots(v.3.1.3), RColorBrewer(v.1.1-3), jquerylib(v.0.1.4), sessioninfo(v.1.2.2), Rcpp(v.1.0.9), plyr(v.1.8.7), progress(v.1.2.2), zlibbioc(v.1.40.0), purrr(v.0.3.4), RCurl(v.1.98-1.7), ps(v.1.7.1), prettyunits(v.1.1.1), viridis(v.0.6.2), ggrepel(v.0.9.1), fs(v.1.5.2), variancePartition(v.1.24.1), magrittr(v.2.0.3), data.table(v.1.14.2), openxlsx(v.4.2.5), DO.db(v.2.9), pkgload(v.1.3.0), patchwork(v.1.1.1), hms(v.1.1.1), mime(v.0.12), evaluate(v.0.15), xtable(v.1.8-4), pbkrtest(v.0.5.1), RhpcBLASctl(v.0.21-247.1), XML(v.3.99-0.10), gridExtra(v.2.3), compiler(v.4.1.2), biomaRt(v.2.50.3), tibble(v.3.1.7), shadowtext(v.0.1.2), KernSmooth(v.2.23-20), crayon(v.1.5.1), minqa(v.1.2.4), htmltools(v.0.5.3), corpcor(v.1.6.10), ggfun(v.0.0.6), mgcv(v.1.8-40), later(v.1.3.0), geneplotter(v.1.72.0), aplot(v.0.1.6), tidyr(v.1.2.0), DBI(v.1.1.3), tweenr(v.1.0.2), dbplyr(v.2.2.1), MASS(v.7.3-58), rappdirs(v.0.3.3), boot(v.1.3-28), Matrix(v.1.4-1), brio(v.1.1.3), cli(v.3.3.0), rbibutils(v.2.2.8), parallel(v.4.1.2), igraph(v.1.3.3), pkgconfig(v.2.0.3), GenomicAlignments(v.1.30.0), plotly(v.4.10.0), xml2(v.1.3.3), foreach(v.1.5.2), ggtree(v.3.2.1), annotate(v.1.72.0), bslib(v.0.4.0), XVector(v.0.34.0), yulab.utils(v.0.0.5), stringr(v.1.4.0), callr(v.3.7.1), digest(v.0.6.29), Biostrings(v.2.62.0), rmarkdown(v.2.14), fastmatch(v.1.1-3), tidytree(v.0.3.9), edgeR(v.3.36.0), PROPER(v.1.26.0), restfulr(v.0.0.15), curl(v.4.3.2), shiny(v.1.7.1), Rsamtools(v.2.10.0), gtools(v.3.9.3), rjson(v.0.2.21), nloptr(v.2.0.3), lifecycle(v.1.0.1), nlme(v.3.1-158), jsonlite(v.1.8.0), aod(v.1.3.2), desc(v.1.4.1), viridisLite(v.0.4.0), limma(v.3.50.3), fansi(v.1.0.3), pillar(v.1.8.0), lattice(v.0.20-45), KEGGREST(v.1.34.0), fastmap(v.1.1.0), httr(v.1.4.3), pkgbuild(v.1.3.1), survival(v.3.3-1), GO.db(v.3.14.0), glue(v.1.6.2), remotes(v.2.4.2), zip(v.2.2.0), png(v.0.1-7), iterators(v.1.0.14), pander(v.0.6.5), bit(v.4.0.4), ggforce(v.0.3.3), stringi(v.1.7.8), sass(v.0.4.2), blob(v.1.2.3), DESeq2(v.1.34.0), caTools(v.1.18.2), memoise(v.2.0.1), dplyr(v.1.0.9) and ape(v.5.6-2)
message(paste0("This is hpgltools commit: ", get_git_commit()))## If you wish to reproduce this exact build of hpgltools, invoke the following:## > git clone http://github.com/abelew/hpgltools.git## > git reset 86f725d3ee8482b45960e3a4541f1d4ade97a406## This is hpgltools commit: Thu Jul 21 11:08:25 2022 -0400: 86f725d3ee8482b45960e3a4541f1d4ade97a406this_save <- paste0(gsub(pattern="\\.Rmd", replace="", x=rmd_file), "-v", ver, ".rda.xz")
message(paste0("Saving to ", this_save))## Saving to index-v20220725.rda.xztmp <- sm(saveme(filename=this_save))