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)## 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")## Deleting the file excel/mm_female_de-v202204.xlsx before writing the tables.female_sig <- extract_significant_genes(lfc=0.6,
female_table,
excel="excel/mm_female_sig-v202204.xlsx",
according_to="deseq")## Deleting the file excel/mm_female_sig-v202204.xlsx before writing the tables.4 Perform a couple of likely ontology searches
ups <- female_sig[["deseq"]][["ups"]][["vaccination"]]
rownames(ups) <- gsub(x=rownames(ups), pattern="gene:", replacement="")
downs <- female_sig[["deseq"]][["downs"]][["vaccination"]]
rownames(downs) <- gsub(x=rownames(downs), pattern="gene:", replacement="")
up_gp <- simple_gprofiler(ups, species="mmusculus")## Performing gProfiler GO search of 130 genes against mmusculus.## GO search found 3 hits.## Performing gProfiler KEGG search of 130 genes against mmusculus.## KEGG search found 2 hits.## Performing gProfiler REAC search of 130 genes against mmusculus.## REAC search found 1 hits.## Performing gProfiler MI search of 130 genes against mmusculus.## MI search found 0 hits.## Performing gProfiler TF search of 130 genes against mmusculus.## TF search found 0 hits.## Performing gProfiler CORUM search of 130 genes against mmusculus.## CORUM search found 0 hits.## Performing gProfiler HP search of 130 genes against mmusculus.## HP search found 1 hits.up_gp[["pvalue_plots"]][["kegg_plot_over"]]## NULLdown_gp <- simple_gprofiler(downs, species="mmusculus")## Performing gProfiler GO search of 369 genes against mmusculus.## GO search found 198 hits.## Performing gProfiler KEGG search of 369 genes against mmusculus.## KEGG search found 6 hits.## Performing gProfiler REAC search of 369 genes against mmusculus.## REAC search found 23 hits.## Performing gProfiler MI search of 369 genes against mmusculus.## MI search found 0 hits.## Performing gProfiler TF search of 369 genes against mmusculus.## TF search found 77 hits.## Performing gProfiler CORUM search of 369 genes against mmusculus.## CORUM search found 0 hits.## Performing gProfiler HP search of 369 genes against mmusculus.## HP search found 28 hits.pp(file="images/down_gprofiler_bp.png", image=down_gp[["pvalue_plots"]][["bpp_plot_over"]])## Warning in pp(file = "images/down_gprofiler_bp.png", image =
## down_gp[["pvalue_plots"]][["bpp_plot_over"]]): There is no device to shut down.dev.off()## png
## 2pp(file="images/down_gprofiler_mf.png", image=down_gp[["pvalue_plots"]][["mfp_plot_over"]])## Warning in pp(file = "images/down_gprofiler_mf.png", image =
## down_gp[["pvalue_plots"]][["mfp_plot_over"]]): There is no device to shut down.dev.off()## png
## 2pp(file="images/down_gprofiler_reactome.png", image=down_gp[["pvalue_plots"]][["reactome_plot_over"]])## Warning in pp(file = "images/down_gprofiler_reactome.png", image =
## down_gp[["pvalue_plots"]][["reactome_plot_over"]]): There is no device to shut
## down.dev.off()## png
## 2pp(file="images/down_gprofiler_kegg.png", image=down_gp[["pvalue_plots"]][["kegg_plot_over"]])## Warning in pp(file = "images/down_gprofiler_kegg.png", image =
## down_gp[["pvalue_plots"]][["kegg_plot_over"]]): There is no device to shut down.dev.off()## png
## 2pp(file="images/down_gprofiler_tf.png", image=down_gp[["pvalue_plots"]][["tf_plot_over"]], height=16, width=8)## Warning in pp(file = "images/down_gprofiler_tf.png", image =
## down_gp[["pvalue_plots"]][["tf_plot_over"]], : There is no device to shut down.dev.off()## png
## 2pp(file="images/down_gprofiler_hp.png", image=down_gp[["pvalue_plots"]][["hp_plot_over"]], height=12, width=8)## Warning in pp(file = "images/down_gprofiler_hp.png", image =
## down_gp[["pvalue_plots"]][["hp_plot_over"]], : There is no device to shut down.dev.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): rappdirs(v.0.3.3), rtracklayer(v.1.54.0), R.methodsS3(v.1.8.1), tidyr(v.1.2.0), ggplot2(v.3.3.6), bit64(v.4.0.5), knitr(v.1.39), R.utils(v.2.11.0), DelayedArray(v.0.20.0), data.table(v.1.14.2), KEGGREST(v.1.34.0), RCurl(v.1.98-1.6), doParallel(v.1.0.17), generics(v.0.1.2), GenomicFeatures(v.1.46.5), preprocessCore(v.1.56.0), callr(v.3.7.0), RhpcBLASctl(v.0.21-247.1), usethis(v.2.1.6), RSQLite(v.2.2.14), shadowtext(v.0.1.2), bit(v.4.0.4), enrichplot(v.1.14.2), xml2(v.1.3.3), httpuv(v.1.6.5), assertthat(v.0.2.1), viridis(v.0.6.2), xfun(v.0.31), hms(v.1.1.1), jquerylib(v.0.1.4), evaluate(v.0.15), promises(v.1.2.0.1), IHW(v.1.22.0), DEoptimR(v.1.0-11), fansi(v.1.0.3), restfulr(v.0.0.13), progress(v.1.2.2), caTools(v.1.18.2), dbplyr(v.2.1.1), igraph(v.1.3.1), DBI(v.1.1.2), geneplotter(v.1.72.0), htmlwidgets(v.1.5.4), purrr(v.0.3.4), ellipsis(v.0.3.2), dplyr(v.1.0.9), backports(v.1.4.1), annotate(v.1.72.0), aod(v.1.3.2), biomaRt(v.2.50.3), vctrs(v.0.4.1), remotes(v.2.4.2), cachem(v.1.0.6), withr(v.2.5.0), ggforce(v.0.3.3), robustbase(v.0.95-0), GenomicAlignments(v.1.30.0), treeio(v.1.18.1), fdrtool(v.1.2.17), prettyunits(v.1.1.1), DOSE(v.3.20.1), ape(v.5.6-2), lazyeval(v.0.2.2), crayon(v.1.5.1), genefilter(v.1.76.0), edgeR(v.3.36.0), pkgconfig(v.2.0.3), slam(v.0.1-50), labeling(v.0.4.2), tweenr(v.1.0.2), nlme(v.3.1-157), pkgload(v.1.2.4), devtools(v.2.4.3), rlang(v.1.0.2), lifecycle(v.1.0.1), downloader(v.0.4), filelock(v.1.0.2), BiocFileCache(v.2.2.1), rprojroot(v.2.0.3), polyclip(v.1.10-0), graph(v.1.72.0), Matrix(v.1.4-1), aplot(v.0.1.5), lpsymphony(v.1.22.0), boot(v.1.3-28), processx(v.3.5.3), png(v.0.1-7), viridisLite(v.0.4.0), rjson(v.0.2.21), bitops(v.1.0-7), R.oo(v.1.24.0), KernSmooth(v.2.23-20), pander(v.0.6.5), Biostrings(v.2.62.0), blob(v.1.2.3), stringr(v.1.4.0), qvalue(v.2.26.0), gProfileR(v.0.7.0), gridGraphics(v.0.5-1), scales(v.1.2.0), memoise(v.2.0.1), magrittr(v.2.0.3), plyr(v.1.8.7), gplots(v.3.1.3), zlibbioc(v.1.40.0), compiler(v.4.1.2), scatterpie(v.0.1.7), BiocIO(v.1.4.0), RColorBrewer(v.1.1-3), lme4(v.1.1-29), DESeq2(v.1.34.0), Rsamtools(v.2.10.0), cli(v.3.3.0), XVector(v.0.34.0), patchwork(v.1.1.1), ps(v.1.7.0), MASS(v.7.3-57), mgcv(v.1.8-40), tidyselect(v.1.1.2), stringi(v.1.7.6), highr(v.0.9), yaml(v.2.3.5), GOSemSim(v.2.20.0), locfit(v.1.5-9.5), ggrepel(v.0.9.1), grid(v.4.1.2), sass(v.0.4.1), fastmatch(v.1.1-3), tools(v.4.1.2), parallel(v.4.1.2), rstudioapi(v.0.13), foreach(v.1.5.2), gridExtra(v.2.3), farver(v.2.1.0), ggraph(v.2.0.5), digest(v.0.6.29), shiny(v.1.7.1), Rcpp(v.1.0.8.3), broom(v.0.8.0), later(v.1.3.0), httr(v.1.4.3), AnnotationDbi(v.1.56.2), Rdpack(v.2.3), colorspace(v.2.0-3), brio(v.1.1.3), XML(v.3.99-0.9), fs(v.1.5.2), splines(v.4.1.2), RBGL(v.1.70.0), yulab.utils(v.0.0.4), PROPER(v.1.26.0), tidytree(v.0.3.9), graphlayouts(v.0.8.0), ggplotify(v.0.1.0), plotly(v.4.10.0), sessioninfo(v.1.2.2), xtable(v.1.8-4), jsonlite(v.1.8.0), nloptr(v.2.0.3), ggtree(v.3.2.1), tidygraph(v.1.2.1), corpcor(v.1.6.10), ggfun(v.0.0.6), Vennerable(v.3.1.0.9000), R6(v.2.5.1), pillar(v.1.7.0), htmltools(v.0.5.2), mime(v.0.12), glue(v.1.6.2), fastmap(v.1.1.0), minqa(v.1.2.4), clusterProfiler(v.4.2.2), BiocParallel(v.1.28.3), codetools(v.0.2-18), fgsea(v.1.20.0), pkgbuild(v.1.3.1), utf8(v.1.2.2), lattice(v.0.20-45), bslib(v.0.3.1), tibble(v.3.1.7), sva(v.3.42.0), pbkrtest(v.0.5.1), curl(v.4.3.2), gtools(v.3.9.2.1), zip(v.2.2.0), GO.db(v.3.14.0), openxlsx(v.4.2.5), survival(v.3.3-1), limma(v.3.50.3), rmarkdown(v.2.14), desc(v.1.4.1), munsell(v.0.5.0), DO.db(v.2.9), GenomeInfoDbData(v.1.2.7), iterators(v.1.0.14), variancePartition(v.1.24.1), reshape2(v.1.4.4), gtable(v.0.3.0) and rbibutils(v.2.2.8)
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 12cc39929e5afa444a4ddf2f7b7a0472713666e3## This is hpgltools commit: Wed Jun 22 15:35:44 2022 -0400: 12cc39929e5afa444a4ddf2f7b7a0472713666e3this_save <- paste0(gsub(pattern="\\.Rmd", replace="", x=rmd_file), "-v", ver, ".rda.xz")
message(paste0("Saving to ", this_save))## Saving to index-v20220705.rda.xztmp <- sm(saveme(filename=this_save))