1 Introduction
2 Preprocessing
I passed these samples to my default preprocessing pipeline on the umiacs cluster with the following command:
cd preprocessing
start=$(pwd)
for i in $(/bin/ls -d *); do
cd $i
cyoa --task pipe --method prnaseq --species hd38_100 --gff_type gene \
--stranded reverse --gff_tag ID \
--input $(/bin/ls *.fastq.gz | tr '\n' ':' | sed 's/:$//g')
cd $start
done3 Genome Annotation
I will extract annotations from ensembl, since it is my favorite. My load_biomart_annotation() function attempts to find an archive server from one year ago, this may be made explicit if one wishes.
One important caveat, when creating the expressionset data structure, the gene IDs must match and the gff usually uses IDs which are a specific combination of columns from the ensembl table. It usually takes a moment to figure out how the two data sources line up; as a result sometimes it is easier to just extract annotations directly from the gff file.
hs_annot <- load_biomart_annotations(species="hsapiens", year="2020")## The biomart annotations file already exists, loading from it.4 Metadata collection
I often just copy/paste the relevant fields from the output files created by trimomatic/htseq/hisat/etc. However, I did make a fun function: gather_preprocessing_metadata() which in theory at least can fill out a metadata xlsx file.
written <- gather_preprocessing_metadata("sample_sheets/all_samples.xlsx")I noticed a couple problems with my function:
- It does not like dates in the MM/DD/YYYY format. There are a couple reasons for which this may have happened.
- It did not save the condition/batch columns I created, so I copy/pasted them from the previous sheet.
5 Initial Data Structure
Let us create an expression set of the metadata and annotations. Note that the
annot <- hs_annot[["annotation"]]
rownames(annot) <- make.names(paste0("gene:", annot[["ensembl_gene_id"]]), unique=TRUE)
rownames(annot) <- gsub(pattern="\\.", replacement=":", x=rownames(annot))
hs_expt <- create_expt(metadata="sample_sheets/all_samples_modified.xlsx",
gene_info=annot, file_column="hisatcounttable") %>%
set_expt_conditions(fact="treatment") %>%
set_expt_batches(fact="replicatenumber") %>%
exclude_genes_expt(column="gene_biotype", method="keep", pattern="protein_coding")## 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.## The sample definitions comprises: 18 rows(samples) and 24 columns(metadata fields).## Matched 21444 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 21481 features and 18 samples.## remove_genes_expt(), before removal, there were 21481 genes, now there are 19909.written <- write_expt(hs_expt, excel="excel/hs_expt.xlsx")## Deleting the file excel/hs_expt.xlsx before writing the tables.## Writing the first sheet, containing a legend and some summary data.## `geom_smooth()` using formula 'y ~ x'## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following object is masked from 'package:S4Vectors':
##
## expand##
## Total:107 s## `geom_smooth()` using formula 'y ~ x'##
## Total:77 s6 Initial plots
Let us look briefly at a couple of likely plots of the data.
plot_libsize(hs_expt)$plot
plot_nonzero(hs_expt)$plot## 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.
hs_expt_norm <- normalize_expt(hs_expt, filter=TRUE, convert="cpm",
norm="quant", transform="log2")## Removing 7890 low-count genes (12019 remaining).plot_pca(hs_expt_norm)$plot
hs_nbsva <- normalize_expt(hs_expt, filter=TRUE, convert="cpm",
batch="svaseq", transform="log2")## Removing 7890 low-count genes (12019 remaining).## Setting 13 low elements to zero.## transform_counts: Found 13 values equal to 0, adding 1 to the matrix.plot_pca(hs_nbsva)$plot
pander::pander(sessionInfo())R version 4.2.1 (2022-06-23)
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: splines, stats4, stats, graphics, grDevices, utils, datasets, methods and base
other attached packages: ruv(v.0.9.7.1), lme4(v.1.1-30), Matrix(v.1.4-1), BiocParallel(v.1.30.3), variancePartition(v.1.26.0), hpgltools(v.1.0), testthat(v.3.1.4), reticulate(v.1.25), SummarizedExperiment(v.1.26.1), GenomicRanges(v.1.48.0), GenomeInfoDb(v.1.32.2), IRanges(v.2.30.0), S4Vectors(v.0.34.0), MatrixGenerics(v.1.8.1), matrixStats(v.0.62.0), Biobase(v.2.56.0) and BiocGenerics(v.0.42.0)
loaded via a namespace (and not attached): rappdirs(v.0.3.3), rtracklayer(v.1.56.1), tidyr(v.1.2.0), ggplot2(v.3.3.6), bit64(v.4.0.5), knitr(v.1.39), DelayedArray(v.0.22.0), data.table(v.1.14.2), KEGGREST(v.1.36.3), RCurl(v.1.98-1.7), doParallel(v.1.0.17), generics(v.0.1.3), GenomicFeatures(v.1.48.3), preprocessCore(v.1.58.0), callr(v.3.7.1), RhpcBLASctl(v.0.21-247.1), usethis(v.2.1.6), RSQLite(v.2.2.15), shadowtext(v.0.1.2), bit(v.4.0.4), enrichplot(v.1.16.1), 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), fansi(v.1.0.3), restfulr(v.0.0.15), progress(v.1.2.2), caTools(v.1.18.2), dbplyr(v.2.2.1), igraph(v.1.3.4), DBI(v.1.1.3), 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.74.0), aod(v.1.3.2), biomaRt(v.2.52.0), 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), GenomicAlignments(v.1.32.1), treeio(v.1.20.1), prettyunits(v.1.1.1), DOSE(v.3.22.0), ape(v.5.6-2), lazyeval(v.0.2.2), crayon(v.1.5.1), genefilter(v.1.78.0), edgeR(v.3.38.1), pkgconfig(v.2.0.3), labeling(v.0.4.2), tweenr(v.1.0.2), nlme(v.3.1-158), pkgload(v.1.3.0), devtools(v.2.4.4), rlang(v.1.0.4), lifecycle(v.1.0.1), miniUI(v.0.1.1.1), downloader(v.0.4), filelock(v.1.0.2), BiocFileCache(v.2.4.0), directlabels(v.2021.1.13), rprojroot(v.2.0.3), polyclip(v.1.10-0), aplot(v.0.1.6), boot(v.1.3-28), processx(v.3.7.0), png(v.0.1-7), viridisLite(v.0.4.0), rjson(v.0.2.21), bitops(v.1.0-7), KernSmooth(v.2.23-20), pander(v.0.6.5), Biostrings(v.2.64.0), blob(v.1.2.3), stringr(v.1.4.0), qvalue(v.2.28.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.42.0), compiler(v.4.2.1), scatterpie(v.0.1.7), BiocIO(v.1.6.0), RColorBrewer(v.1.1-3), Rsamtools(v.2.12.0), cli(v.3.3.0), XVector(v.0.36.0), urlchecker(v.1.0.1), patchwork(v.1.1.1), ps(v.1.7.1), MASS(v.7.3-58), mgcv(v.1.8-40), tidyselect(v.1.1.2), stringi(v.1.7.8), highr(v.0.9), yaml(v.2.3.5), GOSemSim(v.2.22.0), locfit(v.1.5-9.6), ggrepel(v.0.9.1), grid(v.4.2.1), sass(v.0.4.2), fastmatch(v.1.1-3), tools(v.4.2.1), parallel(v.4.2.1), rstudioapi(v.0.13), foreach(v.1.5.2), gridExtra(v.2.3), farver(v.2.1.1), Rtsne(v.0.16), ggraph(v.2.0.5), digest(v.0.6.29), shiny(v.1.7.2), quadprog(v.1.5-8), Rcpp(v.1.0.9), broom(v.1.0.0), later(v.1.3.0), httr(v.1.4.3), AnnotationDbi(v.1.58.0), Rdpack(v.2.4), colorspace(v.2.0-3), brio(v.1.1.3), XML(v.3.99-0.10), fs(v.1.5.2), yulab.utils(v.0.0.5), PROPER(v.1.28.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.4.1), tidygraph(v.1.2.1), corpcor(v.1.6.10), ggfun(v.0.0.6), R6(v.2.5.1), profvis(v.0.3.7), pillar(v.1.8.0), htmltools(v.0.5.3), mime(v.0.12), glue(v.1.6.2), fastmap(v.1.1.0), minqa(v.1.2.4), clusterProfiler(v.4.4.4), codetools(v.0.2-18), fgsea(v.1.22.0), pkgbuild(v.1.3.1), utf8(v.1.2.2), lattice(v.0.20-45), bslib(v.0.4.0), tibble(v.3.1.8), sva(v.3.44.0), pbkrtest(v.0.5.1), curl(v.4.3.2), gtools(v.3.9.3), zip(v.2.2.0), GO.db(v.3.15.0), openxlsx(v.4.2.5), survival(v.3.3-1), limma(v.3.52.2), rmarkdown(v.2.14), desc(v.1.4.1), munsell(v.0.5.0), DO.db(v.2.9), fastcluster(v.1.2.3), GenomeInfoDbData(v.1.2.8), iterators(v.1.0.14), 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 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-v202208.rda.xztmp <- sm(saveme(filename=this_save))