1 Introduction
This is a very preliminary look into how we might combine our library of human infected samples with a set of (mostly macrophage) samples which I downloaded from the SRA. The current data is comprised of ~ 370 samples salmon mapped against hg38 revision 91.
2 Annotation data for download human data.
This short document should gather the annotation data required for playing with the data downloaded from SRA and friends.
hs_gff <- "reference/hg38_91.gff"
hs_gff_annot <- sm(load_gff_annotations(hs_gff, id_col="gene_id"))
hs_annot <- sm(load_biomart_annotations())$annotation
hs_annot$ID <- hs_annot$geneID
rownames(hs_annot) <- make.names(hs_annot[["ensembl_gene_id"]], unique=TRUE)
dim(hs_annot)## [1] 197995 12hs_size_dist <- plot_histogram(hs_annot[["cds_length"]])
hs_size_dist +
ggplot2::scale_x_continuous(limits=c(0, 20000))## Warning: Removed 103681 rows containing non-finite values (stat_bin).## Warning: Removed 103681 rows containing non-finite values (stat_density).## Warning: Removed 2 rows containing missing values (geom_bar).
hs_go_biomart <- sm(load_biomart_go())
hs_lengths <- hs_annot[, c("ensembl_gene_id", "cds_length")]
colnames(hs_lengths) <- c("ID", "width")
rownames(hs_annot) <- make.names(
paste0(hs_annot[["ensembl_transcript_id"]], ".",
hs_annot[["transcript_version"]]),
unique=TRUE)
hs_tx_gene <- hs_annot[, c("ensembl_gene_id", "ensembl_transcript_id")]
hs_tx_gene[["id"]] <- rownames(hs_tx_gene)
hs_tx_gene <- hs_tx_gene[, c("id", "ensembl_gene_id")]
new_hs_annot <- hs_annot
rownames(new_hs_annot) <- make.names(hs_annot[["ensembl_gene_id"]], unique=TRUE)hs_expt <- create_expt("sample_sheets/all_samples_202104.xlsx",
gene_info=new_hs_annot,
tx_gene_map=hs_tx_gene,
file_column="salmonfile")## Reading the sample metadata.## The sample definitions comprises: 369 rows(samples) and 62 columns(metadata fields).## Reading count tables.## Using the transcript to gene mapping.## Reading salmon data with tximport.## Finished reading count data.## Matched 19629 annotations and counts.## Bringing together the count matrix and gene information.## The mapped IDs are not the rownames of your gene information, changing them now.## Some annotations were lost in merging, setting them to 'undefined'.## Saving the expressionset to 'expt.rda'.## The final expressionset has 19629 rows and 369 columns.cds_entries <- fData(hs_expt)
cds_entries <- cds_entries[["gene_biotype"]] == "protein_coding"
hs_cds_expt <- hs_expt
hs_cds_expt$expressionset <- hs_cds_expt$expressionset[cds_entries, ]
new_cds_entries <- fData(hs_cds_expt)3 Initial plots
Najib loves his PCA!
plot_nonzero(hs_cds_expt)$plot## Warning: ggrepel: 360 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
hs_cds_norm <- normalize_expt(hs_cds_expt, transform = "log2", convert = "cpm",
norm = "quant", filter = TRUE)## Removing 0 low-count genes (19201 remaining).## transform_counts: Found 6662 values equal to 0, adding 1 to the matrix.hs_cds_pca <- plot_pca(hs_cds_norm)## plot labels was not set and there are more than 100 samples, disabling it.hs_cds_pca$plot## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
hs_cds_nb <- normalize_expt(hs_cds_expt, transform = "log2", convert = "cpm",
filter = "simple", batch = "svaseq")## Removing 54 low-count genes (19147 remaining).## batch_counts: Before batch/surrogate estimation, 1921423 entries are x==0: 27%.## batch_counts: Before batch/surrogate estimation, 965762 entries are 0<x<1: 14%.## Setting 526138 low elements to zero.## transform_counts: Found 526138 values equal to 0, adding 1 to the matrix.hs_cds_nb_pca <- plot_pca(hs_cds_nb)## plot labels was not set and there are more than 100 samples, disabling it.hs_cds_nb_pca$plot
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: parallel, stats, graphics, grDevices, utils, datasets, methods and base
other attached packages: ruv(v.0.9.7.1), hpgltools(v.1.0), testthat(v.3.0.2), R6(v.2.5.0), Biobase(v.2.50.0) and BiocGenerics(v.0.36.0)
loaded via a namespace (and not attached): utf8(v.1.1.4), tidyselect(v.1.1.0), lme4(v.1.1-26), RSQLite(v.2.2.3), AnnotationDbi(v.1.52.0), grid(v.4.0.3), BiocParallel(v.1.24.1), devtools(v.2.3.2), scatterpie(v.0.1.5), munsell(v.0.5.0), preprocessCore(v.1.52.1), codetools(v.0.2-18), 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), stats4(v.4.0.3), DOSE(v.3.16.0), MatrixGenerics(v.1.2.1), labeling(v.0.4.2), tximport(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), doParallel(v.1.0.16), GenomeInfoDb(v.1.26.2), 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), splines(v.4.0.3), 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), GenomicFeatures(v.1.42.1), backports(v.1.2.1), qvalue(v.2.22.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), 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), S4Vectors(v.0.28.1), SummarizedExperiment(v.1.20.0), ggrepel(v.0.9.1), colorRamps(v.2.3), fs(v.1.5.0), variancePartition(v.1.20.0), magrittr(v.2.0.1), data.table(v.1.14.0), DO.db(v.2.9), openxlsx(v.4.2.3), matrixStats(v.0.58.0), 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), IRanges(v.2.24.1), 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), minqa(v.1.2.4), htmltools(v.0.5.1.1), corpcor(v.1.6.9), mgcv(v.1.8-34), 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), readr(v.1.4.0), Matrix(v.1.3-2), cli(v.2.3.1), igraph(v.1.2.6), GenomicRanges(v.1.42.0), pkgconfig(v.2.0.3), rvcheck(v.0.1.8), GenomicAlignments(v.1.26.0), 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), stringr(v.1.4.0), callr(v.3.5.1), digest(v.0.6.27), Biostrings(v.2.58.0), rmarkdown(v.2.7), fastmatch(v.1.1-0), edgeR(v.3.32.1), PROPER(v.1.22.0), 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), 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), 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), 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 ad4239cab7311278f24a38d1b488513292f5f835## This is hpgltools commit: Mon Apr 26 14:24:06 2021 -0400: ad4239cab7311278f24a38d1b488513292f5f835this_save <- paste0(gsub(pattern="\\.Rmd", replace="", x=rmd_file), "-v", ver, ".rda.xz")
message(paste0("Saving to ", this_save))## Saving to 02_estimation_202104-v202104.rda.xztmp <- sm(saveme(filename=this_save))4 Estimates!
hs_expt <- set_expt_conditions(hs_expt, fact="infectstate")
hs_expt <- set_expt_batches(hs_expt, fact="studypmid")
plot_libsize(hs_expt)$plot
hs_norm <- normalize_expt(hs_expt, transform="log2", convert="cpm",
norm="quant", filter="simple")## Removing 81 low-count genes (19548 remaining).## transform_counts: Found 7380 values equal to 0, adding 1 to the matrix.hs_cor <- plot_corheat(hs_norm)
hs_pca <- plot_pca(hs_norm)## Potentially check over the experimental design, there appear to be missing values.## Warning in plot_pca(hs_norm): There are NA values in the component data. This
## can lead to weird plotting errors.## plot labels was not set and there are more than 100 samples, disabling it.hs_pca$plot## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure## Warning: Removed 267 rows containing missing values (geom_point).
hs_nb <- normalize_expt(hs_expt, transform="log2", convert="cpm",
norm="quant", filter="simple", batch="svaseq")## Warning in normalize_expt(hs_expt, transform = "log2", convert = "cpm", :
## Quantile normalization and sva do not always play well together.## Removing 81 low-count genes (19548 remaining).## batch_counts: Before batch/surrogate estimation, 7380 entries are x==0: 0%.## batch_counts: Before batch/surrogate estimation, 2972809 entries are 0<x<1: 41%.## Setting 320780 low elements to zero.## transform_counts: Found 320780 values equal to 0, adding 1 to the matrix.hs_nb_pca <- plot_pca(hs_nb)## Potentially check over the experimental design, there appear to be missing values.## Warning in plot_pca(hs_nb): There are NA values in the component data. This can
## lead to weird plotting errors.## plot labels was not set and there are more than 100 samples, disabling it.hs_nb_pca$plot## Warning: Removed 267 rows containing missing values (geom_point).
hs_expt <- set_expt_conditions(hs_expt, fact="expttime")
hs_norm <- normalize_expt(hs_expt, transform="log2",
norm="quant", filter="simple", batch="svaseq")## Warning in normalize_expt(hs_expt, transform = "log2", norm = "quant", filter =
## "simple", : Quantile normalization and sva do not always play well together.## Removing 81 low-count genes (19548 remaining).## batch_counts: Before batch/surrogate estimation, 7380 entries are x==0: 0%.## batch_counts: Before batch/surrogate estimation, 1634088 entries are 0<x<1: 23%.## Error in solve.default(t(mod) %*% mod) :
## Lapack routine dgesv: system is exactly singular: U[44,44] = 0## Warning in all_adjusters(count_table, design = design, estimate_type = method, :
## It is highly likely that the underlying reason for this error is too many 0's in
## the dataset, please try doing a filtering of the data and retry.## Error in solve.default(t(mod) %*% mod) :
## Lapack routine dgesv: system is exactly singular: U[44,44] = 0## Warning in do_batch(count_table, method = batch, expt_design = expt_design, :
## The batch_counts call failed. Returning non-batch reduced data.## transform_counts: Found 7380 values equal to 0, adding 1 to the matrix.hs_pca <- plot_pca(hs_norm)## Potentially check over the experimental design, there appear to be missing values.## Warning in plot_pca(hs_norm): There are NA values in the component data. This
## can lead to weird plotting errors.## plot labels was not set and there are more than 100 samples, disabling it.hs_pca$plot## Warning: Removed 267 rows containing missing values (geom_point).
pander::pander(sessionInfo())
message(paste0("This is hpgltools commit: ", get_git_commit()))
this_save <- paste0(gsub(pattern="\\.Rmd", replace="", x=rmd_file), "-v", ver, ".rda.xz")
message(paste0("Saving to ", this_save))
tmp <- sm(saveme(filename=this_save))