1 Introduction
This document is intended to provide an overview of TMRC3 samples which have been sequenced. It includes some plots and analyses showing the relationships among the samples as well as some differential analyses when possible.
2 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()## 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")]3 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).
3.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.
sample_sheet <- "sample_sheets/tmrc3_samples_20210610.xlsx"
hs_expt <- sm(create_expt(sample_sheet,
file_column="hg38100hisatfile", savefile="hs_expt_all.rda",
gene_info=hs_annot))3.2 Minimum coverage sample filtering
I arbitrarily chose 3,000,000 counts as a minimal level of coverage. We may want this to be higher.
hs_valid <- subset_expt(hs_expt, coverage=3000000)## The samples removed (and read coverage) when filtering samples with less than 3e+06 reads are:## TMRC30007 TMRC30010 TMRC30050 TMRC30056 TMRC30124 TMRC30112
## 2588819 52523 808936 2927525 2763416 2410677## subset_expt(): There were 157, now there are 151 samples.4 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.
macr <- set_expt_conditions(hs_expt, fact = "parasitezymodeme") %>%
set_expt_batches(fact = "drug1") %>%
subset_expt(subset = "typeofcells=='Macrophages'")## subset_expt(): There were 157, now there are 14 samples.macr_norm <- sm(normalize_expt(macr, transform="log2", norm="quant", convert="cpm", filter=TRUE))
norm_pca <- plot_pca(macr_norm, plot_labels=FALSE)
pp(file="macrophage_side_experiment/norm_pca.png", image=norm_pca$plot)
plot_3d_pca(norm_pca)## $plot
##
## $file
## [1] "3dpca.html"macr_nb <- sm(normalize_expt(macr, batch="svaseq", filter=TRUE))
macr_nb <- sm(normalize_expt(macr_nb, norm="quant", convert="cpm", transform="log2"))
nb_pca <- plot_pca(macr_nb, plot_labels=FALSE)
pp(file="macrophage_side_experiment/normbatch_pca.png", image=nb_pca$plot)
4.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 rank4.2 Perform DE
##tmp <- normalize_expt(macr, filter=TRUE)
##zy_de <- deseq_pairwise(tmp, model_batch="svaseq")
##zy_ed <- edger_pairwise(tmp, model_batch="svaseq")
zymo_de <- sm(all_pairwise(macr, model_batch="svaseq", filter=TRUE, parallel=FALSE))zymo_table <- sm(combine_de_tables(zymo_de,
excel="macrophage_side_experiment/macrophage_de.xlsx",
parallel=FALSE))pander::pander(sessionInfo())R version 4.0.3 (2020-10-10)
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, parallel, stats, graphics, grDevices, utils, datasets, methods and base
other attached packages: edgeR(v.3.32.1), lme4(v.1.1-26), Matrix(v.1.3-2), BiocParallel(v.1.24.1), variancePartition(v.1.20.0), ruv(v.0.9.7.1), hpgltools(v.1.0), testthat(v.3.0.2), SummarizedExperiment(v.1.20.0), GenomicRanges(v.1.42.0), GenomeInfoDb(v.1.26.2), IRanges(v.2.24.1), S4Vectors(v.0.28.1), MatrixGenerics(v.1.2.1), matrixStats(v.0.58.0), Biobase(v.2.50.0) and BiocGenerics(v.0.36.0)
loaded via a namespace (and not attached): utf8(v.1.1.4), R.utils(v.2.10.1), tidyselect(v.1.1.0), htmlwidgets(v.1.5.3), RSQLite(v.2.2.3), AnnotationDbi(v.1.52.0), grid(v.4.0.3), Rtsne(v.0.15), IHW(v.1.18.0), devtools(v.2.3.2), scatterpie(v.0.1.5), munsell(v.0.5.0), codetools(v.0.2-18), preprocessCore(v.1.52.1), statmod(v.1.4.35), withr(v.2.4.1), colorspace(v.2.0-0), GOSemSim(v.2.16.1), highr(v.0.8), knitr(v.1.31), rstudioapi(v.0.13), Vennerable(v.3.1.0.9000), robustbase(v.0.93-7), DOSE(v.3.16.0), labeling(v.0.4.2), slam(v.0.1-48), lpsymphony(v.1.18.0), GenomeInfoDbData(v.1.2.4), polyclip(v.1.10-0), bit64(v.4.0.5), farver(v.2.0.3), rprojroot(v.2.0.2), downloader(v.0.4), vctrs(v.0.3.6), generics(v.0.1.0), xfun(v.0.21), BiocFileCache(v.1.14.0), fastcluster(v.1.1.25), R6(v.2.5.0), doParallel(v.1.0.16), graphlayouts(v.0.7.1), locfit(v.1.5-9.4), bitops(v.1.0-6), cachem(v.1.0.4), fgsea(v.1.16.0), DelayedArray(v.0.16.1), assertthat(v.0.2.1), scales(v.1.1.1), ggraph(v.2.0.5), enrichplot(v.1.10.2), gtable(v.0.3.0), sva(v.3.38.0), processx(v.3.4.5), tidygraph(v.1.2.0), rlang(v.0.4.10), genefilter(v.1.72.1), lazyeval(v.0.2.2), rtracklayer(v.1.50.0), broom(v.0.7.5), BiocManager(v.1.30.10), yaml(v.2.2.1), reshape2(v.1.4.4), crosstalk(v.1.1.1), GenomicFeatures(v.1.42.1), backports(v.1.2.1), qvalue(v.2.22.0), RBGL(v.1.66.0), clusterProfiler(v.3.18.1), tools(v.4.0.3), usethis(v.2.0.1), ggplot2(v.3.3.3), ellipsis(v.0.3.1), gplots(v.3.1.1), jquerylib(v.0.1.3), RColorBrewer(v.1.1-2), blockmodeling(v.1.0.0), sessioninfo(v.1.1.1), Rcpp(v.1.0.6), plyr(v.1.8.6), progress(v.1.2.2), zlibbioc(v.1.36.0), purrr(v.0.3.4), RCurl(v.1.98-1.2), ps(v.1.5.0), prettyunits(v.1.1.1), openssl(v.1.4.3), viridis(v.0.5.1), cowplot(v.1.1.1), ggrepel(v.0.9.1), colorRamps(v.2.3), fs(v.1.5.0), magrittr(v.2.0.1), data.table(v.1.14.0), DO.db(v.2.9), openxlsx(v.4.2.3), pkgload(v.1.2.0), hms(v.1.0.0), evaluate(v.0.14), xtable(v.1.8-4), pbkrtest(v.0.5-0.1), XML(v.3.99-0.5), gridExtra(v.2.3), compiler(v.4.0.3), biomaRt(v.2.46.3), tibble(v.3.0.6), KernSmooth(v.2.23-18), crayon(v.1.4.1), shadowtext(v.0.0.7), R.oo(v.1.24.0), minqa(v.1.2.4), htmltools(v.0.5.1.1), corpcor(v.1.6.9), mgcv(v.1.8-34), geneplotter(v.1.68.0), tidyr(v.1.1.2), DBI(v.1.1.1), tweenr(v.1.0.1), dbplyr(v.2.1.0), MASS(v.7.3-53.1), rappdirs(v.0.3.3), boot(v.1.3-27), cli(v.2.3.1), R.methodsS3(v.1.8.1), quadprog(v.1.5-8), igraph(v.1.2.6), pkgconfig(v.2.0.3), rvcheck(v.0.1.8), GenomicAlignments(v.1.26.0), plotly(v.4.9.3), xml2(v.1.3.2), foreach(v.1.5.1), annotate(v.1.68.0), bslib(v.0.2.4), XVector(v.0.30.0), EBSeq(v.1.30.0), stringr(v.1.4.0), callr(v.3.5.1), digest(v.0.6.27), graph(v.1.68.0), Biostrings(v.2.58.0), rmarkdown(v.2.7), fastmatch(v.1.1-0), PROPER(v.1.22.0), directlabels(v.2021.1.13), curl(v.4.3), Rsamtools(v.2.6.0), gtools(v.3.8.2), nloptr(v.1.2.2.2), lifecycle(v.1.0.0), nlme(v.3.1-152), jsonlite(v.1.7.2), desc(v.1.2.0), viridisLite(v.0.3.0), askpass(v.1.1), limma(v.3.46.0), fansi(v.0.4.2), pillar(v.1.5.0), lattice(v.0.20-41), DEoptimR(v.1.0-8), fastmap(v.1.1.0), httr(v.1.4.2), pkgbuild(v.1.2.0), survival(v.3.2-7), GO.db(v.3.12.1), glue(v.1.4.2), remotes(v.2.2.0), fdrtool(v.1.2.16), zip(v.2.1.1), iterators(v.1.0.13), pander(v.0.6.3), bit(v.4.0.4), ggforce(v.0.3.2), stringi(v.1.5.3), sass(v.0.3.1), blob(v.1.2.1), DESeq2(v.1.30.1), caTools(v.1.18.1), memoise(v.2.0.0) and dplyr(v.1.0.4)
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 72947fcc6afe09da22d71967059edd84e3063341## This is hpgltools commit: Tue Jun 1 15:57:56 2021 -0400: 72947fcc6afe09da22d71967059edd84e3063341message(paste0("Saving to ", savefile))## Saving to macrophage_side_experiment.rda.xztmp <- sm(saveme(filename=savefile))tmp <- loadme(filename=savefile)