This document is intended to create the data structures used to evaluate our TMRC2 samples. In some cases, this includes only those samples starting in 2019; in other instances I am including our previous (2015-2016) samples.
In all cases the processing performed was:
I am thinking that this meeting will bring Maria Adelaida fully back into the analyses of the parasite data, and therefore may focus primarily on the goals rather than the analyses?
In a couple of important ways the TMRC2 data is much more complex than the TMRC3:
Everything which follows depends on the Existing TriTrypDB annotations revision 46, circa 2019. The following block loads a database of these annotations and turns it into a matrix where the rows are genes and columns are all the annotation types provided by TriTrypDB.
The same database was used to create a matrix of orthologous genes between L.panamensis and all of the other species in the TriTrypDB.
The same database of annotations also provides mappings to the set of annotated GO categories for the L.panamensis genome along with gene lengths.
## Loading taxonomy and species database to cross reference against the download.
## Working on 1/88: org.Tvivax.Y486.v68.eg.db.
## Working on 2/88: org.Tcongolense.IL3000.v68.eg.db.
## Working on 3/88: org.Laethiopica.L147.v68.eg.db.
## Working on 4/88: org.Ltropica.L590.v68.eg.db.
## Working on 5/88: org.Tcruzi.Tula.cl2.v68.eg.db.
## Working on 6/88: org.Lpanamensis.MHOMCOL81L13.v68.eg.db.
## Working on 7/88: org.Lbraziliensis.MHOMBR75M2903.v68.eg.db.
## Working on 8/88: org.Tcruzi.Dm28c.2014.v68.eg.db.
## Working on 9/88: org.Tbrucei.brucei.TREU927.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species/strain.
## Found more than one taxonomy ID match, returning the first match.
## Working on 10/88: org.Lmajor.Friedlin.v68.eg.db.
## Found 279 candidate genera matching Leishmania
## Found an exact match for the combination genus/species not strain for Leishmania major.
## Found more than one taxonomy ID match, returning the first match.
## Working on 11/88: org.Tcruzi.CL.Brener.v68.eg.db.
## Working on 12/88: org.Tcruzi.Esmeraldo.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species/strain.
## Found more than one taxonomy ID match, returning the first match.
## Working on 13/88: org.Lbraziliensis.MHOMBR75M2904.v68.eg.db.
## Working on 14/88: org.Trangeli.SC58.v68.eg.db.
## Working on 15/88: org.Linfantum.JPCM5.v68.eg.db.
## Working on 16/88: org.Tbrucei.gambiense.DAL972.v68.eg.db.
## Working on 17/88: org.Lmajor.LV39c5.v68.eg.db.
## Working on 18/88: org.Lmajor.SD.75.1.v68.eg.db.
## Working on 19/88: org.Tcruzi.JR.cl.4.v68.eg.db.
## Working on 20/88: org.Lmexicana.MHOMGT2001U1103.v68.eg.db.
## Working on 21/88: org.Ldonovani.BPK282A1.v68.eg.db.
## Working on 22/88: org.Adeanei.Cavalho.ATCC.PRA.265.v68.eg.db.
## Found 4 candidate genera matching Angomonas
## Found an exact match for the combination genus/species not strain for Angomonas deanei.
## Found more than one taxonomy ID match, returning the first match.
## Working on 23/88: org.Bayalai.B08.376.v68.eg.db.
## Found 20 candidate genera matching Blechomonas
## Found an exact match for the combination genus/species not strain for Blechomonas ayalai.
## Found more than one taxonomy ID match, returning the first match.
## Working on 24/88: org.Bnonstop.P57.v68.eg.db.
## Found 28 candidate genera matching Blastocrithidia
## Found a genus, but not species for Blastocrithidia nonstop P57, not adding taxon ID number.
## Working on 25/88: org.Bsaltans.Lake.Konstanz.v68.eg.db.
## Found 28 candidate genera matching Bodo
## Found an exact match for the combination genus/species not strain for Bodo saltans.
## Working on 26/88: org.Cfasciculata.Cf.Cl.v68.eg.db.
## Found 53 candidate genera matching Crithidia
## Found an exact match for the combination genus/species not strain for Crithidia fasciculata.
## Working on 27/88: org.Emonterogeii.LV88.v68.eg.db.
## Found 13 candidate genera matching Endotrypanum
## Found an exact match for the combination genus/species not strain for Endotrypanum monterogeii.
## Found more than one taxonomy ID match, returning the first match.
## Working on 28/88: org.Lamazonensis.MHOMBR71973M2269.v68.eg.db.
## Found 279 candidate genera matching Leishmania
## Found an exact match for the combination genus/species not strain for Leishmania amazonensis.
## Found more than one taxonomy ID match, returning the first match.
## Working on 29/88: org.Lamazonensis.PH8.v68.eg.db.
## Found 279 candidate genera matching Leishmania
## Found an exact match for the combination genus/species not strain for Leishmania amazonensis.
## Found more than one taxonomy ID match, returning the first match.
## Working on 30/88: org.Larabica.LEM1108.v68.eg.db.
## Found 279 candidate genera matching Leishmania
## Found an exact match for the combination genus/species not strain for Leishmania arabica.
## Working on 31/88: org.Lbraziliensis.MHOMBR75M2904.2019.v68.eg.db.
## Found 279 candidate genera matching Leishmania
## Found an exact match for the combination genus/species not strain for Leishmania braziliensis.
## Found more than one taxonomy ID match, returning the first match.
## Working on 32/88: org.Ldonovani.BHU.1220.v68.eg.db.
## Found 279 candidate genera matching Leishmania
## Found an exact match for the combination genus/species not strain for Leishmania donovani.
## Found more than one taxonomy ID match, returning the first match.
## Working on 33/88: org.Ldonovani.CL.SL.v68.eg.db.
## Found 279 candidate genera matching Leishmania
## Found an exact match for the combination genus/species not strain for Leishmania donovani.
## Found more than one taxonomy ID match, returning the first match.
## Working on 34/88: org.Ldonovani.HU3.v68.eg.db.
## Found 279 candidate genera matching Leishmania
## Found an exact match for the combination genus/species not strain for Leishmania donovani.
## Found more than one taxonomy ID match, returning the first match.
## Working on 35/88: org.Ldonovani.LV9.v68.eg.db.
## Found 279 candidate genera matching Leishmania
## Found an exact match for the combination genus/species not strain for Leishmania donovani.
## Found more than one taxonomy ID match, returning the first match.
## Working on 36/88: org.Lenriettii.LEM3045.v68.eg.db.
## Found 279 candidate genera matching Leishmania
## Found an exact match for the combination genus/species not strain for Leishmania enriettii.
## Found more than one taxonomy ID match, returning the first match.
## Working on 37/88: org.Lenriettii.MCAVBR2001CUR178.v68.eg.db.
## Found 279 candidate genera matching Leishmania
## Found an exact match for the combination genus/species not strain for Leishmania enriettii.
## Found more than one taxonomy ID match, returning the first match.
## Working on 38/88: org.Lgerbilli.LEM452.v68.eg.db.
## Found 279 candidate genera matching Leishmania
## Found an exact match for the combination genus/species not strain for Leishmania gerbilli.
## Found more than one taxonomy ID match, returning the first match.
## Working on 39/88: org.Lmajor.Friedlin.2021.v68.eg.db.
## Found 279 candidate genera matching Leishmania
## Found an exact match for the combination genus/species not strain for Leishmania major.
## Found more than one taxonomy ID match, returning the first match.
## Working on 40/88: org.Lmartiniquensis.LEM2494.v68.eg.db.
## Found 279 candidate genera matching Leishmania
## Found an exact match for the combination genus/species not strain for Leishmania martiniquensis.
## Found more than one taxonomy ID match, returning the first match.
## Working on 41/88: org.Lmartiniquensis.MHOMTH2012LSCM1.v68.eg.db.
## Found 279 candidate genera matching Leishmania
## Found an exact match for the combination genus/species not strain for Leishmania martiniquensis.
## Found more than one taxonomy ID match, returning the first match.
## Working on 42/88: org.Lorientalis.MHOMTH2014LSCM4.v68.eg.db.
## Found 279 candidate genera matching Leishmania
## Found an exact match for the combination genus/species not strain for Leishmania orientalis.
## Found more than one taxonomy ID match, returning the first match.
## Working on 43/88: org.Lpanamensis.MHOMPA94PSC.1.v68.eg.db.
## Found 279 candidate genera matching Leishmania
## Found an exact match for the combination genus/species not strain for Leishmania panamensis.
## Found more than one taxonomy ID match, returning the first match.
## Working on 44/88: org.Lpyrrhocoris.H10.v68.eg.db.
## Found 73 candidate genera matching Leptomonas
## Found an exact match for the combination genus/species not strain for Leptomonas pyrrhocoris.
## Found more than one taxonomy ID match, returning the first match.
## Working on 45/88: org.Lseymouri.ATCC.30220.v68.eg.db.
## Found 73 candidate genera matching Leptomonas
## Found an exact match for the combination genus/species not strain for Leptomonas seymouri.
## Working on 46/88: org.Lsp.Ghana.MHOMGH2012GH5.v68.eg.db.
## Found 279 candidate genera matching Leishmania
## Found an exact match for the combination genus/species not strain for Leishmania sp..
## Working on 47/88: org.Lsp.Namibia.MPRONA1975252LV425.v68.eg.db.
## Found 279 candidate genera matching Leishmania
## Found an exact match for the combination genus/species not strain for Leishmania sp..
## Working on 48/88: org.Ltarentolae.Parrot.TarII.v68.eg.db.
## Found 279 candidate genera matching Leishmania
## Found an exact match for the combination genus/species not strain for Leishmania tarentolae.
## Found more than one taxonomy ID match, returning the first match.
## Working on 49/88: org.Ltarentolae.Parrot.Tar.II.2019.v68.eg.db.
## Found 279 candidate genera matching Leishmania
## Found an exact match for the combination genus/species not strain for Leishmania tarentolae.
## Found more than one taxonomy ID match, returning the first match.
## Working on 50/88: org.Lturanica.LEM423.v68.eg.db.
## Found 279 candidate genera matching Leishmania
## Found an exact match for the combination genus/species not strain for Leishmania turanica.
## Working on 51/88: org.Pconfusum.CUL13.v68.eg.db.
## Found 4 candidate genera matching Paratrypanosoma
## Found an exact match for the combination genus/species not strain for Paratrypanosoma confusum.
## Found more than one taxonomy ID match, returning the first match.
## Working on 52/88: org.Phertigi.MCOEPA1965C119.v68.eg.db.
## Found 3 candidate genera matching Porcisia
## Found an exact match for the combination genus/species not strain for Porcisia hertigi.
## Working on 53/88: org.Tbrucei.EATRO1125.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma brucei.
## Found more than one taxonomy ID match, returning the first match.
## Working on 54/88: org.Tbrucei.Lister.427.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma brucei.
## Found more than one taxonomy ID match, returning the first match.
## Working on 55/88: org.Tbrucei.Lister.427.2018.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma brucei.
## Found more than one taxonomy ID match, returning the first match.
## Working on 56/88: org.Tcongolense.IL3000.2019.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma congolense.
## Working on 57/88: org.Tcongolense.Tc1148.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma congolense.
## Working on 58/88: org.Tcruzi.231.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma cruzi.
## Working on 59/88: org.Tcruzi.Berenice.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma cruzi.
## Working on 60/88: org.Tcruzi.Brazil.A4.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma cruzi.
## Working on 61/88: org.Tcruzi.Bug2148.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma cruzi.
## Working on 62/88: org.Tcruzi.CL.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma cruzi.
## Working on 63/88: org.Tcruzi.CL.Brener.Esmeraldo.like.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma cruzi.
## Working on 64/88: org.Tcruzi.CL.Brener.Non.Esmeraldo.like.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma cruzi.
## Working on 65/88: org.Tcruzi.Dm28c.2017.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma cruzi.
## Working on 66/88: org.Tcruzi.Dm28c.2018.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma cruzi.
## Working on 67/88: org.Tcruzi.G.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma cruzi.
## Working on 68/88: org.Tcruzi.S11.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma cruzi.
## Working on 69/88: org.Tcruzi.S15.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma cruzi.
## Working on 70/88: org.Tcruzi.S154a.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma cruzi.
## Working on 71/88: org.Tcruzi.S162a.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma cruzi.
## Working on 72/88: org.Tcruzi.S23b.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma cruzi.
## Working on 73/88: org.Tcruzi.S44a.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma cruzi.
## Working on 74/88: org.Tcruzi.S92a.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma cruzi.
## Working on 75/88: org.Tcruzi.Sylvio.X101.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma cruzi.
## Working on 76/88: org.Tcruzi.Sylvio.X101.2012.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma cruzi.
## Working on 77/88: org.Tcruzi.TCC.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma cruzi.
## Working on 78/88: org.Tcruzi.Y.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma cruzi.
## Working on 79/88: org.Tcruzi.Y.C6.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma cruzi.
## Working on 80/88: org.Tcruzi.Ycl2.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma cruzi.
## Working on 81/88: org.Tcruzi.Ycl4.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma cruzi.
## Working on 82/88: org.Tcruzi.Ycl6.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma cruzi.
## Working on 83/88: org.Tcruzi.marinkellei.B7.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma cruzi.
## Working on 84/88: org.Tequiperdum.OVI.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma equiperdum.
## Found more than one taxonomy ID match, returning the first match.
## Working on 85/88: org.Tevansi.STIB.805.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma evansi.
## Found more than one taxonomy ID match, returning the first match.
## Working on 86/88: org.Tgrayi.ANR4.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma grayi.
## Working on 87/88: org.Tmelophagium.St.Kilda.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma melophagium.
## Working on 88/88: org.Ttheileri.isolate.Edinburgh.v68.eg.db.
## Found 549 candidate genera matching Trypanosoma
## Found an exact match for the combination genus/species not strain for Trypanosoma theileri.
## Found the following hits: Leishmania panamensis MHOM/COL/81/L13, Leishmania braziliensis MHOM/BR/75/M2903, Leishmania braziliensis MHOM/BR/75/M2904, Leishmania mexicana MHOM/GT/2001/U1103, Leishmania sp. Ghana MHOM/GH/2012/GH5, choosing the first.
## Using: Leishmania panamensis MHOM/COL/81/L13.
## org.Lpanamensis.MHOMCOL81L13.v68.eg.db is already installed and a copy should be found at: /data/renv/library/R-4.3/x86_64-conda-linux-gnu/org.Lpanamensis.MHOMCOL81L13.v68.eg.db/extdata/org.Lpanamensis.MHOMCOL81L13.v68.eg.sqlite.
panamensis_pkg <- panamensis_db[["pkgname"]]
package_name <- panamensis_db[["pkgname"]]
if (is.null(panamensis_pkg)) {
panamensis_pkg <- panamensis_entry[["OrgdbPkg"]]
package_name <- panamensis_pkg
}
tt <- library(panamensis_pkg, character.only = TRUE)
## Loading required package: AnnotationDbi
##
## Attaching package: 'AnnotationDbi'
## The following object is masked from 'package:dplyr':
##
## select
##
panamensis_env <- get0(panamensis_pkg)
all_fields <- columns(panamensis_env)
all_lp_annot <- sm(load_orgdb_annotations(
panamensis_env,
keytype = "gid",
fields = c("annot_gene_entrez_id", "annot_gene_name",
"annot_strand", "annot_chromosome", "annot_cds_length",
"annot_gene_product")))$genes
## Testing to see just how big the full database is.
## testing <- load_orgdb_annotations(panamensis_pkg, keytype = "gid", fields = "all")
lp_go <- load_orgdb_go(panamensis_pkg)
lp_go <- lp_go[, c("GID", "GO")]
lp_lengths <- all_lp_annot[, c("gid", "annot_cds_length")]
colnames(lp_lengths) <- c("ID", "length")
all_lp_annot[["annot_gene_product"]] <- tolower(all_lp_annot[["annot_gene_product"]])
orthos <- sm(extract_eupath_orthologs(db = panamensis_pkg))
data_structures <- c(data_structures, "lp_lengths", "lp_go", "all_lp_annot")
The following block loads the full genome sequence for panamensis. We may use this later to attempt to estimate PCR primers to discern strains.
I am not sure how to increase the number of open files in a container, as a result this does not work.
The process of sample estimation takes two primary inputs:
An expressionSet(or summarizedExperiment) is a data structure used in R to examine RNASeq data. It is comprised of annotations, metadata, and expression data. In the case of our processing pipeline, the location of the expression data is provided by the filenames in the metadata.
The following list contains the colors we have chosen to use when plotting the various ways of discerning the data.
color_choices <- list(
"strain" = list(
## "z1.0" = "#333333", ## Changed this to 'braz' to make it easier to find them.
"z2.0" = "#555555",
"z3.0" = "#777777",
"z2.1" = "#874400",
"z2.2" = "#0000cc",
"z2.3" = "#cc0000",
"z2.4" = "#df7000",
"z3.2" = "#888888",
"z1.0" = "#cc00cc",
"z1.5" = "#cc00cc",
"b2904" = "#cc00cc",
"unknown" = "#cbcbcb"),
## "null" = "#000000"),
"zymo" = list(
"z22" = "#0000cc",
"z23" = "#cc0000"),
"cf" = list(
"cure" = "#006f00",
"fail" = "#9dffa0",
"unknown" = "#cbcbcb",
"notapplicable" = "#000000"),
"condition" = list(
"inf" = "#199c75",
"inf_sb" = "#d65d00",
"uninf" = "#6e6ea3",
"uninf_sb" = "#d83956"),
"significance" = list(
"lt0" = "#ffe0e0",
"lt1" = "#ffa0a0",
"lt2" = "#f94040",
"lt4" = "#a00000",
"gt0" = "#eeccf9",
"gt1" = "#de8bf9",
"gt2" = "#ad07e3",
"gt4" = "#410257"),
"drug" = list(
"none" = "#989898",
"antimony" = "#088b64"),
"oldnew" = list(
"previous" = "#2233aa",
"current" = "#9c0303"),
"infectedp" = list(
"uninfected" = "#676767",
"infected" = "#ac06e2"),
"treatment_zymo" = list(
"infsb_z23" = "#E7298A",
"inf_z23" = "#D95F02",
"uninf_none" = "#66A61E",
"uninfsb_none" = "#E6AB02",
"inf_z22" = "#1B9E77",
"infsb_z22" = "#7570B3"),
"susceptibility" = list(
"resistant" = "#8563a7",
"sensitive" = "#8d0000",
"ambiguous" = "#cbcbcb",
"unknown" = "#555555"))
data_structures <- c(data_structures, "color_choices")
All of the above focused entire on the parasite samples, now let us pull up the macrophage infected samples. This will comprise two datasets, one of the human and one of the parasite.
The metadata for the macrophage samples contains a couple of columns for mapped human and parasite reads. We will therefore use them separately to create two expressionsets, one for each species.
** Note **: I forgot to commit the addition of plot_metadata factors() in the last run of this. In addition, I need to add an explicit month to load_biomart_annotations() or change the function to search a couple more months before it stops trying to find an archive.
## Using mart: ENSEMBL_MART_ENSEMBL from host: apr2020.archive.ensembl.org.
## Successfully connected to the hsapiens_gene_ensembl database.
## Finished downloading ensembl gene annotations.
## Finished downloading ensembl structure annotations.
## symbol columns is null, pattern matching 'symbol'.
## Including symbols, there are 68503 vs the 249740 gene annotations.
## Not dropping haplotype chromosome annotations, set drop_haplotypes = TRUE if this is bad.
## Saving annotations to hsapiens_biomart_annotations.rda.
## Finished save().
hs_annot <- hs_annot[["annotation"]]
hs_annot[["transcript"]] <- paste0(rownames(hs_annot), ".", hs_annot[["transcript_version"]])
rownames(hs_annot) <- make.names(hs_annot[["ensembl_gene_id"]], unique = TRUE)
rownames(hs_annot) <- paste0("gene:", rownames(hs_annot))
tx_gene_map <- hs_annot[, c("transcript", "ensembl_gene_id")]
sanitize_columns <- c("drug", "macrophagetreatment", "macrophagezymodeme")
macr_annot <- hs_annot
rownames(macr_annot) <- gsub(x = rownames(macr_annot),
pattern = "^gene:",
replacement = "")
hs_macrophage <- create_expt(sample_sheet, gene_info = macr_annot,
file_column = "hg38100hisatfile") %>%
set_expt_conditions(fact = "macrophagetreatment") %>%
set_expt_batches(fact = "macrophagezymodeme") %>%
sanitize_expt_pData(columns = sanitize_columns) %>%
subset_expt(nonzero = 12000)
## 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: 69 rows(samples) and 83 columns(metadata fields).
## Matched 21481 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 69 samples.
## The numbers of samples by condition are:
##
## inf inf_sb uninf uninf_sb
## 30 29 5 5
## The number of samples by batch are:
##
## none z2.2 z2.3
## 10 30 29
## The samples (and read coverage) removed when filtering 12000 non-zero genes are:
## TMRC30162
## 10208
## by number of genes.
## subset_expt(): There were 69, now there are 68 samples.
fixed_genenames <- gsub(x = rownames(exprs(hs_macrophage)), pattern = "^gene:",
replacement = "")
hs_macrophage <- set_expt_genenames(hs_macrophage, ids = fixed_genenames)
table(pData(hs_macrophage)$condition)
##
## inf inf_sb uninf uninf_sb
## 29 29 5 5
## The following 3 lines were copy/pasted to datastructures and should be removed soon.
nostrain <- is.na(pData(hs_macrophage)[["strainid"]])
pData(hs_macrophage)[nostrain, "strainid"] <- "none"
pData(hs_macrophage)[["strain_zymo"]] <- paste0("s", pData(hs_macrophage)[["strainid"]],
"_", pData(hs_macrophage)[["macrophagezymodeme"]])
uninfected <- pData(hs_macrophage)[["strain_zymo"]] == "snone_none"
pData(hs_macrophage)[uninfected, "strain_zymo"] <- "uninfected"
pData(hs_macrophage)[["infectedp"]] <- "infected"
pData(hs_macrophage)[uninfected, "infectedp"] <- "uninfected"
data_structures <- c(data_structures, "hs_macrophage")
1 sample has been excluded from the analysis but is in the sample sheet. I am reasonably certain I know which, but will double-check here.
sample_sheet_ids = c("TMRC30051","TMRC30057","TMRC30059","TMRC30060","TMRC30061","TMRC30062",
"TMRC30063","TMRC30064","TMRC30065","TMRC30066","TMRC30067","TMRC30069",
"TMRC30117","TMRC30162","TMRC30243","TMRC30244","TMRC30245","TMRC30246",
"TMRC30247","TMRC30248","TMRC30249","TMRC30250","TMRC30251","TMRC30252",
"TMRC30266","TMRC30267","TMRC30268","TMRC30286","TMRC30326","TMRC30316",
"TMRC30317","TMRC30322","TMRC30323","TMRC30328","TMRC30318","TMRC30319",
"TMRC30324","TMRC30325","TMRC30320","TMRC30321","TMRC30327","TMRC30312",
"TMRC30297","TMRC30298","TMRC30299","TMRC30300","TMRC30295","TMRC30296",
"TMRC30303","TMRC30304","TMRC30301","TMRC30302","TMRC30314","TMRC30315",
"TMRC30313")
found <- sample_sheet_ids %in% colnames(exprs(hs_macr))
colnames(exprs(hs_macr))[!found]
all_human <- sanitize_expt_pData(hs_macrophage, columns = "drug") %>%
set_expt_conditions(fact = "drug") %>%
set_expt_batches(fact = "typeofcells")
## The numbers of samples by condition are:
##
## antimony none
## 34 34
## The number of samples by batch are:
##
## Macrophages U937
## 54 14
data_structures <- c(data_structures, "all_human")
## The following 3 lines were copy/pasted to datastructures and should be removed soon.
no_strain_idx <- pData(all_human)[["strainid"]] == "none"
##pData(all_human)[["strainid"]] <- paste0("s", pData(all_human)[["strainid"]],
## "_", pData(all_human)[["macrophagezymodeme"]])
pData(all_human)[no_strain_idx, "strainid"] <- "none"
table(pData(all_human)[["strainid"]])
##
## 10763 10772 10977 11026 11075 11126 12251 12309 12355 12367 2169 7158 none
## 2 8 2 2 2 8 7 8 2 7 8 2 10
all_human_types <- set_expt_conditions(all_human, fact = "typeofcells") %>%
set_expt_batches(fact = "drug")
## The numbers of samples by condition are:
##
## Macrophages U937
## 54 14
## The number of samples by batch are:
##
## antimony none
## 34 34
data_structures <- c(data_structures, "all_human_types")
type_zymo_fact <- paste0(pData(all_human_types)[["condition"]], "_",
pData(all_human_types)[["macrophagezymodeme"]])
type_zymo <- set_expt_conditions(all_human_types, fact = type_zymo_fact)
## The numbers of samples by condition are:
##
## Macrophages_none Macrophages_z22 Macrophages_z23 U937_none
## 8 23 23 2
## U937_z22 U937_z23
## 6 6
data_structures <- c(data_structures, "type_zymo")
type_drug_fact <- paste0(pData(all_human_types)[["condition"]], "_",
pData(all_human_types)[["drug"]])
type_drug <- set_expt_conditions(all_human_types, fact = type_drug_fact)
## The numbers of samples by condition are:
##
## Macrophages_antimony Macrophages_none U937_antimony
## 27 27 7
## U937_none
## 7
data_structures <- c(data_structures, "type_drug")
strain_fact <- pData(all_human_types)[["strainid"]]
table(strain_fact)
## strain_fact
## 10763 10772 10977 11026 11075 11126 12251 12309 12355 12367 2169 7158 none
## 2 8 2 2 2 8 7 8 2 7 8 2 10
new_conditions <- paste0(pData(hs_macrophage)[["macrophagetreatment"]], "_",
pData(hs_macrophage)[["macrophagezymodeme"]])
## Note the sanitize() call is redundant with the addition of sanitize() in the
## datastructures file, but I don't want to wait to rerun that.
hs_macr <- set_expt_conditions(hs_macrophage, fact = new_conditions) %>%
sanitize_expt_pData(column = "drug") %>%
set_expt_colors(color_choices[["treatment_zymo"]]) %>%
subset_expt(subset = "typeofcells!='U937'")
## The numbers of samples by condition are:
##
## inf_z22 inf_z23 infsb_z22 infsb_z23 uninf_none uninfsb_none
## 14 15 15 14 5 5
## The samples excluded are: TMRC30309, TMRC30293, TMRC30294, TMRC30291, TMRC30292, TMRC30307, TMRC30308, TMRC30310, TMRC30331, TMRC30311, TMRC30332, TMRC30305, TMRC30306, TMRC30330.
## subset_expt(): There were 68, now there are 54 samples.
data_structures <- c(data_structures, "hs_macr")
ggstats_parasite <- plot_metadata_factors(hs_macr, column = "parasitemappingrate",
type = "ggstats", scale = "log2")
pp(file = "images/ggstats_parasiterate_all_macrophage_drug_treatment.png")
ggstats_parasite
dev.off()
## png
## 2
## The numbers of samples by condition are:
##
## antimony none
## 27 27
hs_macr_strain_expt <- set_expt_conditions(hs_macr, fact = "macrophagezymodeme") %>%
subset_expt(subset = "macrophagezymodeme != 'none'")
## The numbers of samples by condition are:
##
## none z22 z23
## 8 23 23
## The samples excluded are: TMRC30059, TMRC30060, TMRC30266, TMRC30268, TMRC30326, TMRC30327, TMRC30312, TMRC30313.
## subset_expt(): There were 54, now there are 46 samples.
##
## 10763 10772 10977 11026 11075 11126 12251 12309 12355 12367 2169 7158 none
## 2 6 2 2 2 6 5 6 2 5 6 2 8
Let us see if the sankey plot of these samples looks useful…
ggstats_slreads <- plot_metadata_factors(hs_macrophage, column = "hisatlpsinglemapped",
type = "ggstats", scale = "log2")
pp(file = "images/ggstats_slreads_all_macrophage.png")
ggstats_slreads
dev.off()
## png
## 2
ggstats_violin <- plot_metadata_factors(hs_macrophage, column = "hisatlpsinglemapped",
scale = "log2")
ggstats_violin
macr_sankey <- plot_meta_sankey(hs_macrophage, color_choices = color_choices,
factors = c("oldnew", "drug", "infectedp", "macrophagezymodeme"))
macr_sankey
## A sankey plot describing the metadata of 68 samples,
## including 26 out of 0 nodes and traversing metadata factors:
## .
Finally, split off the U937 samples.
## The samples excluded are: TMRC30051, TMRC30057, TMRC30059, TMRC30060, TMRC30061, TMRC30062, TMRC30063, TMRC30064, TMRC30065, TMRC30066, TMRC30067, TMRC30069, TMRC30117, TMRC30243, TMRC30244, TMRC30245, TMRC30246, TMRC30247, TMRC30248, TMRC30249, TMRC30250, TMRC30251, TMRC30252, TMRC30266, TMRC30267, TMRC30268, TMRC30286, TMRC30326, TMRC30316, TMRC30317, TMRC30322, TMRC30323, TMRC30328, TMRC30318, TMRC30319, TMRC30324, TMRC30325, TMRC30320, TMRC30321, TMRC30327, TMRC30312, TMRC30297, TMRC30298, TMRC30299, TMRC30300, TMRC30295, TMRC30296, TMRC30303, TMRC30304, TMRC30301, TMRC30302, TMRC30314, TMRC30315, TMRC30313.
## subset_expt(): There were 68, now there are 14 samples.
In the previous block, we used a new invocation of ensembl-derived annotation data, this time we can just use our existing parasite gene annotations.
lp_macrophage <- create_expt(sample_sheet,
file_column = "lpanamensisv36hisatfile",
gene_info = all_lp_annot,
savefile = glue("rda/lp_macrophage-v{ver}.rda"),
annotation = "org.Lpanamensis.MHOMCOL81L13.v46.eg.db") %>%
set_expt_conditions(fact = "macrophagezymodeme") %>%
set_expt_batches(fact = "macrophagetreatment")
## 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: 69 rows(samples) and 83 columns(metadata fields).
## Warning in create_expt(sample_sheet, file_column = "lpanamensisv36hisatfile", :
## Some samples were removed when cross referencing the samples against the count
## data.
## Matched 8778 annotations and counts.
## Bringing together the count matrix and gene information.
## The final expressionset has 8778 features and 66 samples.
## The numbers of samples by condition are:
##
## none z2.2 z2.3
## 8 29 29
## The number of samples by batch are:
##
## inf inf_sb uninf uninf_sb
## 29 29 4 4
unfilt_written <- write_expt(
lp_macrophage,
excel = glue("analyses/macrophage_de/{ver}/read_counts/lp_macrophage_reads_unfiltered-v{ver}.xlsx"))
## Writing the first sheet, containing a legend and some summary data.
## The following samples have less than 5705.7 genes.
## [1] "TMRC30066" "TMRC30117" "TMRC30244" "TMRC30246" "TMRC30249" "TMRC30266"
## [7] "TMRC30268" "TMRC30326" "TMRC30323" "TMRC30319" "TMRC30325" "TMRC30327"
## [13] "TMRC30312" "TMRC30300" "TMRC30304" "TMRC30302" "TMRC30313" "TMRC30309"
## [19] "TMRC30292" "TMRC30331" "TMRC30332" "TMRC30330"
## 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.
## 175550 entries are 0. We are on a log scale, adding 1 to the data.
##
## Changed 175550 zero count features.
##
## Naively calculating coefficient of variation/dispersion with respect to condition.
##
## Finished calculating dispersion estimates.
##
## `geom_smooth()` using formula = 'y ~ x'
## This expressionset does not support lmer with condition+batch
## Error in density.default(x, adjust = adj) : 'x' contains missing values
## Error in density.default(x, adjust = adj) : 'x' contains missing values
## `geom_smooth()` using formula = 'y ~ x'
lp_macrophage_filt <- subset_expt(lp_macrophage, nonzero = 2500) %>%
semantic_expt_filter(semantic = c("amastin", "gp63", "leishmanolysin"),
semantic_column = "annot_gene_product")
## The samples (and read coverage) removed when filtering 2500 non-zero genes are:
## TMRC30066 TMRC30117 TMRC30244 TMRC30246 TMRC30266 TMRC30268 TMRC30326 TMRC30323
## 3080 1147 1662 2834 822 3444 375 84
## TMRC30319 TMRC30325 TMRC30327 TMRC30312 TMRC30304 TMRC30313 TMRC30309 TMRC30330
## 374 356 129 76 289 96 188 181
## by number of genes.
## subset_expt(): There were 66, now there are 50 samples.
## semantic_expt_filter(): Removed 68 genes.
data_structures <- c(data_structures, "lp_macrophage", "lp_macrophage_filt")
filt_written <- write_expt(lp_macrophage_filt,
excel = glue("analyses/macrophage_de/{ver}/read_counts/lp_macrophage_reads_filtered-v{ver}.xlsx"))
## Writing the first sheet, containing a legend and some summary data.
## The following samples have less than 5661.5 genes.
## [1] "TMRC30249" "TMRC30300" "TMRC30302" "TMRC30292" "TMRC30331" "TMRC30332"
## 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.
## 44583 entries are 0. We are on a log scale, adding 1 to the data.
##
## Changed 44583 zero count features.
##
## Naively calculating coefficient of variation/dispersion with respect to condition.
##
## Finished calculating dispersion estimates.
##
## `geom_smooth()` using formula = 'y ~ x'
## Error in density.default(x, adjust = adj) : 'x' contains missing values
## Error in density.default(x, adjust = adj) : 'x' contains missing values
## `geom_smooth()` using formula = 'y ~ x'
lp_macrophage <- lp_macrophage_filt
lp_macrophage_nosb <- subset_expt(lp_macrophage, subset="batch!='inf_sb'")
## The samples excluded are: TMRC30051, TMRC30062, TMRC30065, TMRC30069, TMRC30248, TMRC30249, TMRC30251, TMRC30252, TMRC30317, TMRC30321, TMRC30298, TMRC30300, TMRC30296, TMRC30302, TMRC30315, TMRC30294, TMRC30292, TMRC30308, TMRC30331, TMRC30332, TMRC30306.
## subset_expt(): There were 50, now there are 29 samples.
lp_nosb_write <- write_expt(
lp_macrophage_nosb,
excel = glue("analyses/macrophage_de/{ver}/read_counts/lp_macrophage_nosb_reads-v{ver}.xlsx"))
## Writing the first sheet, containing a legend and some summary data.
## 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.6396 entries are 0. We are on a log scale, adding 1 to the data.
## Changed 6396 zero count features.
## Naively calculating coefficient of variation/dispersion with respect to condition.
## Finished calculating dispersion estimates.
## `geom_smooth()` using formula = 'y ~ x'The expressionset has a minimal or missing set of conditions/batches.
## `geom_smooth()` using formula = 'y ~ x'
found_idx <- data_structures %in% ls()
if (sum(!found_idx) > 0) {
not_found <- data_structures[!found_idx]
warning("Some datastructures were not generated: ", toString(not_found), ".")
data_structures <- data_structures[found_idx]
}
save(list = data_structures, file = glue("rda/tmrc2_data_structures-v{ver}.rda"))
R version 4.3.3 (2024-02-29)
Platform: x86_64-conda-linux-gnu (64-bit)
locale: C
attached base packages: stats4, stats, graphics, grDevices, utils, datasets, methods and base
other attached packages: ruv(v.0.9.7.1), BiocParallel(v.1.36.0), variancePartition(v.1.32.5), org.Lpanamensis.MHOMCOL81L13.v68.eg.db(v.2024.08), AnnotationDbi(v.1.64.1), futile.logger(v.1.4.3), EuPathDB(v.1.6.0), GenomeInfoDbData(v.1.2.11), dplyr(v.1.1.4), Heatplus(v.3.10.0), ggplot2(v.3.5.1), hpgltools(v.1.0), Matrix(v.1.6-5), glue(v.1.7.0), SummarizedExperiment(v.1.32.0), GenomicRanges(v.1.54.1), GenomeInfoDb(v.1.38.8), IRanges(v.2.36.0), S4Vectors(v.0.40.2), MatrixGenerics(v.1.14.0), matrixStats(v.1.3.0), Biobase(v.2.62.0) and BiocGenerics(v.0.48.1)
loaded via a namespace (and not attached): fs(v.1.6.4), bitops(v.1.0-8), doParallel(v.1.0.17), HDO.db(v.0.99.1), httr(v.1.4.7), RColorBrewer(v.1.1-3), insight(v.0.20.2), numDeriv(v.2016.8-1.1), backports(v.1.5.0), tools(v.4.3.3), utf8(v.1.2.4), R6(v.2.5.1), statsExpressions(v.1.5.5), mgcv(v.1.9-1), lazyeval(v.0.2.2), withr(v.3.0.1), gridExtra(v.2.3), prettyunits(v.1.2.0), cli(v.3.6.3), performance(v.0.12.2), formatR(v.1.14), AnnotationHubData(v.1.35.0), labeling(v.0.4.3), sass(v.0.4.9), prismatic(v.1.1.2), BWStest(v.0.2.3), mvtnorm(v.1.2-5), genefilter(v.1.84.0), readr(v.2.1.5), pbapply(v.1.7-2), Rsamtools(v.2.18.0), yulab.utils(v.0.1.5), DOSE(v.3.28.2), stringdist(v.0.9.12), AnnotationForge(v.1.44.0), limma(v.3.58.1), RSQLite(v.2.3.7), generics(v.0.1.3), BiocIO(v.1.12.0), gtools(v.3.9.5), vroom(v.1.6.5), zip(v.2.3.1), GO.db(v.3.18.0), fansi(v.1.0.6), abind(v.1.4-5), lifecycle(v.1.0.4), edgeR(v.4.0.16), yaml(v.2.3.10), gplots(v.3.1.3.1), Rtsne(v.0.17), biocViews(v.1.70.0), qvalue(v.2.34.0), SparseArray(v.1.2.4), BiocFileCache(v.2.10.2), paletteer(v.1.6.0), grid(v.4.3.3), blob(v.1.2.4), promises(v.1.3.0), crayon(v.1.5.3), lattice(v.0.22-6), cowplot(v.1.1.3), GenomicFeatures(v.1.54.4), annotate(v.1.80.0), KEGGREST(v.1.42.0), zeallot(v.0.1.0), pillar(v.1.9.0), knitr(v.1.48), varhandle(v.2.0.6), fgsea(v.1.28.0), boot(v.1.3-30), rjson(v.0.2.21), corpcor(v.1.6.10), kSamples(v.1.2-10), codetools(v.0.2-20), fastmatch(v.1.1-4), data.table(v.1.15.4), Rdpack(v.2.6), vctrs(v.0.6.5), png(v.0.1-8), testthat(v.3.2.1.1), gtable(v.0.3.5), rematch2(v.2.1.2), datawizard(v.0.12.2), cachem(v.1.1.0), xfun(v.0.46), openxlsx(v.4.2.6.1), rbibutils(v.2.2.16), S4Arrays(v.1.2.1), mime(v.0.12), correlation(v.0.8.5), coda(v.0.19-4.1), survival(v.3.7-0), iterators(v.1.0.14), statmod(v.1.5.0), gmp(v.0.7-4), directlabels(v.2024.1.21), interactiveDisplayBase(v.1.40.0), nlme(v.3.1-165), pbkrtest(v.0.5.3), bit64(v.4.0.5), EnvStats(v.2.8.1), progress(v.1.2.3), filelock(v.1.0.3), rprojroot(v.2.0.4), bslib(v.0.8.0), KernSmooth(v.2.23-24), colorspace(v.2.1-1), DBI(v.1.2.3), tidyselect(v.1.2.1), bit(v.4.0.5), compiler(v.4.3.3), curl(v.5.2.1), rvest(v.1.0.4), httr2(v.1.0.2), graph(v.1.80.0), BiocCheck(v.1.38.2), xml2(v.1.3.6), desc(v.1.4.3), DelayedArray(v.0.28.0), plotly(v.4.10.4), bayestestR(v.0.14.0), rtracklayer(v.1.62.0), caTools(v.1.18.2), scales(v.1.3.0), remaCor(v.0.0.18), quadprog(v.1.5-8), RBGL(v.1.78.0), multcompView(v.0.1-10), rappdirs(v.0.3.3), stringr(v.1.5.1), digest(v.0.6.36), minqa(v.1.2.7), ggsankey(v.0.0.99999), rmarkdown(v.2.27), aod(v.1.3.3), RhpcBLASctl(v.0.23-42), XVector(v.0.42.0), htmltools(v.0.5.8.1), pkgconfig(v.2.0.3), lme4(v.1.1-35.5), highr(v.0.11), dbplyr(v.2.5.0), fastmap(v.1.2.0), rlang(v.1.1.4), htmlwidgets(v.1.6.4), shiny(v.1.9.1), SuppDists(v.1.1-9.7), farver(v.2.1.2), jquerylib(v.0.1.4), jsonlite(v.1.8.8), GOSemSim(v.2.28.1), RCurl(v.1.98-1.16), magrittr(v.2.0.3), patchwork(v.1.2.0), parameters(v.0.22.1), munsell(v.0.5.1), Rcpp(v.1.0.13), stringi(v.1.8.4), brio(v.1.1.5), zlibbioc(v.1.48.2), MASS(v.7.3-60), AnnotationHub(v.3.10.1), plyr(v.1.8.9), parallel(v.4.3.3), ggrepel(v.0.9.5), Biostrings(v.2.70.3), PMCMRplus(v.1.9.10), splines(v.4.3.3), pander(v.0.6.5), hms(v.1.1.3), locfit(v.1.5-9.10), fastcluster(v.1.2.6), RUnit(v.0.4.33), ggsignif(v.0.6.4), effectsize(v.0.8.9), reshape2(v.1.4.4), biomaRt(v.2.58.2), pkgload(v.1.4.0), futile.options(v.1.0.1), BiocVersion(v.3.18.1), XML(v.3.99-0.17), evaluate(v.0.24.0), lambda.r(v.1.2.4), BiocManager(v.1.30.23), nloptr(v.2.1.1), tzdb(v.0.4.0), foreach(v.1.5.2), httpuv(v.1.6.15), MatrixModels(v.0.5-3), BayesFactor(v.0.9.12-4.7), tidyr(v.1.3.1), purrr(v.1.0.2), BiocBaseUtils(v.1.4.0), broom(v.1.0.6), xtable(v.1.8-4), restfulr(v.0.0.15), Rmpfr(v.0.9-5), fANCOVA(v.0.6-1), later(v.1.3.2), viridisLite(v.0.4.2), OrganismDbi(v.1.44.0), tibble(v.3.2.1), lmerTest(v.3.1-3), ggstatsplot(v.0.12.4), memoise(v.2.0.1), GenomicAlignments(v.1.38.2), sva(v.3.50.0) and GSEABase(v.1.64.0)
## If you wish to reproduce this exact build of hpgltools, invoke the following:
## > git clone http://github.com/abelew/hpgltools.git
## > git reset
## This is hpgltools commit:
## Saving to 01datasets.rda.xz