sample_sheet <- glue::glue("sample_sheets/tmrc2_samples_{ver}.xlsx")
This document is intended to provide a general overview of the TMRC2 samples which have thus far been sequenced. 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:
The analyses in this document use the matrices of counts/gene from #3 and variants/position from #4 in order to provide some images and metrics describing the samples we have sequenced so far.
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.
tt <- sm(library(EuPathDB))
orgdb <- "org.Lpanamensis.MHOMCOL81L13.v46.eg.db"
tt <- sm(library(orgdb, character.only=TRUE))
pan_db <- org.Lpanamensis.MHOMCOL81L13.v46.eg.db
all_fields <- columns(pan_db)
all_lp_annot <- sm(load_orgdb_annotations(
pan_db,
keytype = "gid",
fields = c("annot_gene_entrez_id", "annot_gene_name",
"annot_strand", "annot_chromosome", "annot_cds_length",
"annot_gene_product")))$genes
lp_go <- sm(load_orgdb_go(pan_db))
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(EuPathDB::extract_eupath_orthologs(db = pan_db))
meta <- sm(EuPathDB::download_eupath_metadata(webservice="tritrypdb"))
lp_entry <- EuPathDB::get_eupath_entry(species="Leishmania panamensis", metadata=meta)
## Found the following hits: Leishmania panamensis MHOM/COL/81/L13, Leishmania panamensis strain MHOM/PA/94/PSC-1, choosing the first.
## Using: Leishmania panamensis MHOM/COL/81/L13.
colnames(lp_entry)
## [1] "AnnotationVersion" "AnnotationSource" "BiocVersion"
## [4] "DataProvider" "Genome" "GenomeSource"
## [7] "GenomeVersion" "NumArrayGene" "NumChipChipGene"
## [10] "NumChromosome" "NumCodingGene" "NumCommunity"
## [13] "NumContig" "NumEC" "NumEST"
## [16] "NumGene" "NumGO" "NumOrtholog"
## [19] "NumOtherGene" "NumPopSet" "NumProteomics"
## [22] "NumPseudogene" "NumRNASeq" "NumRTPCR"
## [25] "NumSNP" "NumTFBS" "Organellar"
## [28] "ReferenceStrain" "MegaBP" "PrimaryKey"
## [31] "ProjectID" "RecordClassName" "SourceID"
## [34] "SourceVersion" "TaxonomyID" "TaxonomyName"
## [37] "URLGenome" "URLGFF" "URLProtein"
## [40] "Coordinate_1_based" "Maintainer" "SourceUrl"
## [43] "Tags" "BsgenomePkg" "GrangesPkg"
## [46] "OrganismdbiPkg" "OrgdbPkg" "TxdbPkg"
## [49] "Taxon" "Genus" "Species"
## [52] "Strain" "BsgenomeFile" "GrangesFile"
## [55] "OrganismdbiFile" "OrgdbFile" "TxdbFile"
## [58] "GenusSpecies" "TaxonUnmodified" "TaxonCanonical"
## [61] "TaxonXref"
testing_panamensis <- "BSGenome.Leishmania.panamensis.MHOMCOL81L13.v53"
## testing_panamensis <- EuPathDB::make_eupath_bsgenome(entry=lp_entry, eu_version="v46")
library(as.character(testing_panamensis), character.only=TRUE)
## Loading required package: BSgenome
## Loading required package: Biostrings
## Loading required package: XVector
##
## Attaching package: 'Biostrings'
## The following object is masked from 'package:base':
##
## strsplit
## Loading required package: rtracklayer
lp_genome <- get0(as.character(testing_panamensis))
data_structures <- c(data_structures, "lp_genome")
## Error in eval(expr, envir, enclos): object 'data_structures' not found
Resequence samples: TMRC20002, TMRC20006, TMRC20004 (maybe TMRC20008 and TMRC20029)
The process of sample estimation takes two primary inputs:
An expressionset 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 first lines of the following block create the Expressionset. All of the following lines perform various normalizations and generate plots from it.
The following samples are much lower coverage:
20210610: I made some manual changes to the sample sheet which I downloaded, filling in some zymodeme with ‘unknown’
data_structures <- c()
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"),
"cf" = list(
"cure" = "#006f00",
"fail" = "#9dffa0",
"unknown" = "#cbcbcb",
"notapplicable" = "#000000"),
"susceptibility" = list(
"resistant" = "#8563a7",
"sensitive" = "#8d0000",
"ambiguous" = "#cbcbcb",
"unknown" = "#555555"))
data_structures <- c(data_structures, "color_choices")
sanitize_columns <- c("passagenumber", "clinicalresponse", "clinicalcategorical",
"zymodemecategorical")
lp_expt <- create_expt(sample_sheet,
gene_info = all_lp_annot,
annotation_name = orgdb,
savefile = glue::glue("rda/tmrc2_lp_expt_all_raw-v{ver}.rda"),
id_column = "hpglidentifier",
file_column = "lpanamensisv36hisatfile") %>%
set_expt_conditions(fact = "zymodemecategorical") %>%
subset_expt(nonzero = 8550) %>%
set_expt_colors(color_choices[["strain"]]) %>%
subset_expt(coverage = 5000000) %>%
set_expt_colors(color_choices[["strain"]]) %>%
semantic_expt_filter(semantic = c("amastin", "gp63", "leishmanolysin"),
semantic_column = "annot_gene_product") %>%
sanitize_expt_metadata(columns = sanitize_columns) %>%
set_expt_factors(columns = sanitize_columns, class = "factor")
## Reading the sample metadata.
## Dropped 11 rows from the sample metadata because the sample ID is blank.
## 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: 110 rows(samples) and 65 columns(metadata fields).
## Warning in create_expt(sample_sheet, gene_info = all_lp_annot, annotation_name
## = orgdb, : 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.
## Some annotations were lost in merging, setting them to 'undefined'.
## The final expressionset has 8778 features and 105 samples.
## The samples (and read coverage) removed when filtering 8550 non-zero genes are:
## TMRC20002 TMRC20006
## 11681227 6670348
## subset_expt(): There were 105, now there are 103 samples.
## Warning in set_expt_colors(., color_choices[["strain"]]): Colors for the
## following categories are not being used: null.
## The samples removed (and read coverage) when filtering samples with less than 5e+06 reads are:
## TMRC20004 TMRC20029
## 564812 1658096
## subset_expt(): There were 103, now there are 101 samples.
## Warning in set_expt_colors(., color_choices[["strain"]]): Colors for the
## following categories are not being used: null.
## semantic_expt_filter(): Removed 68 genes.
data_structures <- c(data_structures, "lp_expt")
save(list = "lp_expt", file = glue::glue("rda/tmrc2_lp_expt_all_sanitized-v{ver}.rda"))
table(pData(lp_expt)[["zymodemecategorical"]])
##
## b2904 unknown z10 z15 z20 z21 z22 z23 z24 z30
## 1 2 1 1 1 7 43 41 2 1
## z32
## 1
table(pData(lp_expt)[["clinicalresponse"]])
##
## cure failure
## 38 38
## laboratory line laboratory line miltefosine resistant
## 1 1
## nd reference strain
## 19 4
Column ‘Q’ in the sample sheet, make a categorical version of it with these parameters:
The sanitize_percent() function seeks to make the percentage values recorded by excel more reliable. Unfortunately, sometimes excel displays the value ‘49%’ when the information recorded in the worksheet is any one of the following:
st <- pData(lp_expt)[["susceptibilityinfectionreduction32ugmlsbvcurrentdata"]]
starting <- sanitize_percent(st)
st
## [1] "3%" "47%" "90%"
## [4] "19%" NA NA
## [7] NA NA NA
## [10] "In process" "92%" "36%"
## [13] "28%" "48%" "45%"
## [16] "85%" "52%" "0.83"
## [19] "60%" "0.18" "0.99"
## [22] "0.97" "0.59" "0.95"
## [25] "0.91" "In process" "99%"
## [28] "53%" "99%" "0.53"
## [31] "0.93" "50%" "94%"
## [34] "92%" "0.92" "68%"
## [37] "94%" "95%" "In process"
## [40] "0.36" "0.94" "48%"
## [43] "7%" "0.87" "44%"
## [46] "0.72" "0.14000000000000001" "97%"
## [49] "0.37" "In process" "0.78"
## [52] "In process" "0.83" "0.97"
## [55] "92%" "In process" "In process"
## [58] "In process" "In process" "In process"
## [61] "In process" "In process" "In process"
## [64] "In process" "In process" "In process"
## [67] "In process" "In process" "In process"
## [70] "0" "In process" "In process"
## [73] "99%" "0.44" "In process"
## [76] "In process" "In process" "85%"
## [79] "100%" "43%" "In process"
## [82] "98%" "In process" "88%"
## [85] "94%" "74%" "In process"
## [88] "38%" "In process" "52%"
## [91] "41%" "In process" "In process"
## [94] "In process" "100%" "In process"
## [97] "93%" "49%" "99%"
## [100] "99%" "63%"
starting
## [1] 0.03 0.47 0.90 0.19 NA NA NA NA NA NA 0.92 0.36 0.28 0.48 0.45
## [16] 0.85 0.52 0.83 0.60 0.18 0.99 0.97 0.59 0.95 0.91 NA 0.99 0.53 0.99 0.53
## [31] 0.93 0.50 0.94 0.92 0.92 0.68 0.94 0.95 NA 0.36 0.94 0.48 0.07 0.87 0.44
## [46] 0.72 0.14 0.97 0.37 NA 0.78 NA 0.83 0.97 0.92 NA NA NA NA NA
## [61] NA NA NA NA NA NA NA NA NA 0.00 NA NA 0.99 0.44 NA
## [76] NA NA 0.85 1.00 0.43 NA 0.98 NA 0.88 0.94 0.74 NA 0.38 NA 0.52
## [91] 0.41 NA NA NA 1.00 NA 0.93 0.49 0.99 0.99 0.63
sus_categorical <- starting
na_idx <- is.na(starting)
sum(na_idx)
## [1] 37
sus_categorical[na_idx] <- "unknown"
resist_idx <- starting <= 0.35
sus_categorical[resist_idx] <- "resistant"
indeterminant_idx <- starting >= 0.36 & starting <= 0.48
sus_categorical[indeterminant_idx] <- "ambiguous"
susceptible_idx <- starting >= 0.49
sus_categorical[susceptible_idx] <- "sensitive"
sus_categorical <- as.factor(sus_categorical)
pData(lp_expt)[["sus_category_historical"]] <- sus_categorical
table(sus_categorical)
## sus_categorical
## ambiguous resistant sensitive unknown
## 12 7 45 37
starting_current <- sanitize_percent(pData(lp_expt)[["susceptibilityinfectionreduction32ugmlsbvcurrentdata"]])
sus_categorical_current <- starting_current
na_idx <- is.na(starting_current)
sum(na_idx)
## [1] 37
sus_categorical_current[na_idx] <- "unknown"
resist_idx <- starting_current <= 0.35
sus_categorical_current[resist_idx] <- "resistant"
indeterminant_idx <- starting_current >= 0.36 & starting_current <= 0.48
sus_categorical_current[indeterminant_idx] <- "ambiguous"
susceptible_idx <- starting_current >= 0.49
sus_categorical_current[susceptible_idx] <- "sensitive"
sus_categorical_current <- as.factor(sus_categorical_current)
pData(lp_expt)[["sus_category_current"]] <- sus_categorical_current
table(sus_categorical_current)
## sus_categorical_current
## ambiguous resistant sensitive unknown
## 12 7 45 37
lp_strain <- lp_expt %>%
set_expt_batches(fact = sus_categorical_current) %>%
set_expt_colors(color_choices[["strain"]])
## Warning in set_expt_colors(., color_choices[["strain"]]): Colors for the
## following categories are not being used: null.
table(pData(lp_strain)[["condition"]])
##
## b2904 unknown z1.0 z1.5 z2.0 z2.1 z2.2 z2.3 z2.4 z3.0
## 1 2 1 1 1 7 43 41 2 1
## z3.2
## 1
save(list = "lp_strain", file = glue::glue("rda/tmrc2_lp_strain-v{ver}.rda"))
data_structures <- c(data_structures, "lp_strain")
lp_two_strains <- subset_expt(lp_strain, subset = "condition=='z2.3'|condition=='z2.2'")
## subset_expt(): There were 101, now there are 84 samples.
save(list = "lp_two_strains",
file = glue::glue("rda/tmrc2_lp_two_strains-v{ver}.rda"))
data_structures <- c(data_structures, "lp_two_strains")
lp_cf <- set_expt_conditions(lp_expt, fact = "clinicalcategorical") %>%
set_expt_batches(fact = sus_categorical_current) %>%
set_expt_colors(color_choices[["cf"]])
## Warning in set_expt_colors(., color_choices[["cf"]]): Colors for the following
## categories are not being used: notapplicable.
table(pData(lp_cf)[["condition"]])
##
## cure fail unknown
## 38 38 25
data_structures <- c(data_structures, "lp_cf")
save(list = "lp_cf",
file = glue::glue("rda/tmrc2_lp_cf-v{ver}.rda"))
lp_cf_known <- subset_expt(lp_cf, subset="condition!='unknown'")
## subset_expt(): There were 101, now there are 76 samples.
data_structures <- c(data_structures, "lp_cf_known")
save(list = "lp_cf_known",
file = glue::glue("rda/tmrc2_lp_cf_known-v{ver}.rda"))
lp_susceptibility <- set_expt_conditions(lp_expt, fact = "sus_category_current") %>%
set_expt_batches(fact = "clinicalcategorical") %>%
set_expt_colors(colors = color_choices[["susceptibility"]])
save(list = "lp_susceptibility",
file = glue::glue("rda/tmrc2_lp_susceptibility-v{ver}.rda"))
data_structures <- c(data_structures, "lp_susceptibility")
TODO: Do this with and without sva and compare the results.
lp_zymo <- subset_expt(lp_expt, subset = "condition=='z2.2'|condition=='z2.3'")
## subset_expt(): There were 101, now there are 84 samples.
data_structures <- c(data_structures, "lp_zymo")
save(list = "lp_zymo",
file = glue::glue("rda/tmrc2_lp_zymo-v{ver}.rda"))
macrophage_sheet <- "sample_sheets/tmrc2_macrophage_samples_202203.xlsx"
lp_macrophage <- create_expt(macrophage_sheet,
file_column="lpanamensisv36hisatfile",
gene_info=all_lp_annot,
savefile = glue::glue("rda/lp_macrophage-v{ver}.rda"),
annotation="org.Lpanamensis.MHOMCOL81L13.v46.eg.db") %>%
set_expt_conditions(fact="macrophagezymodeme") %>%
set_expt_batches(fact="macrophagetreatment") %>%
subset_expt(nonzero = 8550) %>%
semantic_expt_filter(semantic = c("amastin", "gp63", "leishmanolysin"),
semantic_column = "annot_gene_product")
## Reading the sample metadata.
## The sample definitions comprises: 28 rows(samples) and 68 columns(metadata fields).
## Warning in create_expt(macrophage_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.
## Some annotations were lost in merging, setting them to 'undefined'.
## The final expressionset has 8778 features and 23 samples.
## The samples (and read coverage) removed when filtering 8550 non-zero genes are:
## TMRC30051 TMRC30061 TMRC30062 TMRC30065 TMRC30066 TMRC30069 TMRC30117 TMRC30244
## 190577 636331 38270 34155 3080 424761 1147 1662
## TMRC30246 TMRC30249 TMRC30251 TMRC30252
## 2834 14507 84934 634958
## subset_expt(): There were 23, now there are 11 samples.
## semantic_expt_filter(): Removed 68 genes.
data_structures <- c(data_structures, "lp_macrophage")
hs_annot <- load_biomart_annotations(year="2020")
## 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")]
hs_macrophage <- create_expt(
macrophage_sheet,
gene_info = hs_annot,
file_column = "hg38100hisatfile") %>%
set_expt_conditions(fact="macrophagezymodeme") %>%
set_expt_batches(fact="macrophagetreatment") %>%
subset_expt(nonzero = 8550)
## Reading the sample metadata.
## The sample definitions comprises: 28 rows(samples) and 68 columns(metadata fields).
## Matched 21447 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 28 samples.
## The samples (and read coverage) removed when filtering 8550 non-zero genes are:
## named numeric(0)
## subset_expt(): There were 28, now there are 28 samples.
table(pData(hs_macrophage)$condition)
##
## none z2.2 z2.3
## 4 12 12
data_structures <- c(data_structures, "hs_macrophage")
lp_previous <- create_expt("sample_sheets/tmrc2_samples_20191203.xlsx",
file_column = "tophat2file",
savefile = glue::glue("rda/lp_previous-v{ver}.rda"))
## Reading the sample metadata.
## Dropped 13 rows from the sample metadata because the sample ID is blank.
## The sample definitions comprises: 50 rows(samples) and 38 columns(metadata fields).
## Warning in create_expt("sample_sheets/tmrc2_samples_20191203.xlsx", file_column
## = "tophat2file", : Some samples were removed when cross referencing the samples
## against the count data.
## Matched 8841 annotations and counts.
## Bringing together the count matrix and gene information.
## The final expressionset has 8841 features and 33 samples.
data_structures <- c(data_structures, "lp_previous")
tt <- lp_previous$expressionset
rownames(tt) <- gsub(pattern = "^exon_", replacement = "", x = rownames(tt))
rownames(tt) <- gsub(pattern = "\\.1$", replacement = "", x = rownames(tt))
lp_previous$expressionset <- tt
rm(tt)
## The next line drops the samples which are missing the SNP pipeline.
lp_snp <- subset_expt(lp_expt, subset="!is.na(pData(lp_expt)[['freebayessummary']])")
## subset_expt(): There were 101, now there are 101 samples.
new_snps <- count_expt_snps(lp_snp, annot_column = "freebayessummary", snp_column="PAIRED")
## New names:
## • `DP` -> `DP...3`
## • `RO` -> `RO...8`
## • `AO` -> `AO...9`
## • `QR` -> `QR...12`
## • `QA` -> `QA...13`
## • `DP` -> `DP...42`
## • `RO` -> `RO...43`
## • `QR` -> `QR...44`
## • `AO` -> `AO...45`
## • `QA` -> `QA...46`
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## • `DP` -> `DP...3`
## • `RO` -> `RO...8`
## • `AO` -> `AO...9`
## • `QR` -> `QR...12`
## • `QA` -> `QA...13`
## • `DP` -> `DP...42`
## • `RO` -> `RO...43`
## • `QR` -> `QR...44`
## • `AO` -> `AO...45`
## • `QA` -> `QA...46`
## Warning: NAs introduced by coercion
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## • `DP` -> `DP...3`
## • `RO` -> `RO...8`
## • `AO` -> `AO...9`
## • `QR` -> `QR...12`
## • `QA` -> `QA...13`
## • `DP` -> `DP...42`
## • `RO` -> `RO...43`
## • `QR` -> `QR...44`
## • `AO` -> `AO...45`
## • `QA` -> `QA...46`
data_structures <- c(data_structures, "new_snps")
old_snps <- count_expt_snps(lp_previous, annot_column = "bcftable", snp_column = 2)
## The rownames are missing the chromosome identifier,
## they probably came from an older version of this method.
data_structures <- c(data_structures, "old_snps")
save(list = "lp_snp",
file = glue::glue("rda/lp_snp-v{ver}.rda"))
data_structures <- c(data_structures, "lp_snp")
save(list = "new_snps",
file = glue::glue("rda/new_snps-v{ver}.rda"))
data_structures <- c(data_structures, "new_snps")
save(list = "old_snps",
file = glue::glue("rda/old_snps-v{ver}.rda"))
data_structures <- c(data_structures, "old_snps")
nonzero_snps <- exprs(new_snps) != 0
colSums(nonzero_snps)
## tmrc20001 tmrc20065 tmrc20005 tmrc20007 tmrc20008 tmrc20027 tmrc20028 tmrc20032
## 0 93649 0 0 0 351343 338580 146302
## tmrc20040 tmrc20066 tmrc20039 tmrc20037 tmrc20038 tmrc20067 tmrc20068 tmrc20041
## 58753 93615 25115 98958 97676 93954 96583 53184
## tmrc20015 tmrc20009 tmrc20010 tmrc20016 tmrc20011 tmrc20012 tmrc20013 tmrc20017
## 96398 15890 93816 146124 13914 456 94766 48288
## tmrc20014 tmrc20018 tmrc20019 tmrc20070 tmrc20020 tmrc20021 tmrc20022 tmrc20025
## 17245 140438 14829 97336 15484 101127 18143 364240
## tmrc20024 tmrc20036 tmrc20069 tmrc20033 tmrc20026 tmrc20031 tmrc20076 tmrc20073
## 18471 60087 18792 33663 15074 19139 18385 96169
## tmrc20055 tmrc20079 tmrc20071 tmrc20078 tmrc20094 tmrc20042 tmrc20058 tmrc20072
## 22246 96224 94353 18836 87878 19734 94524 50292
## tmrc20059 tmrc20048 tmrc20057 tmrc20088 tmrc20056 tmrc20060 tmrc20077 tmrc20074
## 94091 97164 48944 15594 22683 21506 18773 22132
## tmrc20063 tmrc20053 tmrc20052 tmrc20064 tmrc20075 tmrc20051 tmrc20050 tmrc20049
## 28254 20181 100709 93173 97982 94125 17200 16168
## tmrc20062 tmrc20110 tmrc20080 tmrc20043 tmrc20083 tmrc20054 tmrc20085 tmrc20046
## 93677 16997 96528 95623 21167 93603 89765 48608
## tmrc20093 tmrc20089 tmrc20047 tmrc20090 tmrc20044 tmrc20045 tmrc20061 tmrc20105
## 48254 90421 92637 91564 14861 50403 116906 86758
## tmrc20108 tmrc20109 tmrc20098 tmrc20096 tmrc20097 tmrc20101 tmrc20092 tmrc20082
## 97005 17932 92927 17534 46863 17753 16578 108121
## tmrc20102 tmrc20099 tmrc20100 tmrc20091 tmrc20084 tmrc20087 tmrc20103 tmrc20104
## 92380 91383 94381 15059 46548 14947 49368 94237
## tmrc20086 tmrc20107 tmrc20081 tmrc20106 tmrc20095
## 15813 95370 19533 18830 81200
## My old_snps is using an older annotation incorrectly, so fix it here:
Biobase::annotation(old_snps$expressionset) <- Biobase::annotation(new_snps$expressionset)
both_snps <- combine_expts(new_snps, old_snps)
save(list = "both_snps",
file = glue::glue("rda/both_snps-v{ver}.rda"))
data_structures <- c(data_structures, "both_snps")
save(list = data_structures, file = glue::glue("rda/tmrc2_data_structures-v{ver}.rda"))
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: stats4, stats, graphics, grDevices, utils, datasets, methods and base
other attached packages: BSGenome.Leishmania.panamensis.MHOMCOL81L13.v53(v.2021.07), BSgenome(v.1.64.0), rtracklayer(v.1.56.1), Biostrings(v.2.64.0), XVector(v.0.36.0), futile.logger(v.1.4.3), org.Lpanamensis.MHOMCOL81L13.v46.eg.db(v.2020.07), AnnotationDbi(v.1.58.0), EuPathDB(v.1.6.0), GenomeInfoDbData(v.1.2.8), Heatplus(v.3.4.0), hpgltools(v.1.0), testthat(v.3.1.4), 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), AnnotationForge(v.1.38.0), 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), lambda.r(v.1.2.4), 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), tzdb(v.0.3.0), 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), stringdist(v.0.9.8), 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), vroom(v.1.5.7), AnnotationHubData(v.1.26.1), 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), 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), AnnotationHub(v.3.4.0), rprojroot(v.2.0.3), polyclip(v.1.10-0), graph(v.1.74.0), Matrix(v.1.4-1), 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), blob(v.1.2.3), stringr(v.1.4.0), qvalue(v.2.28.0), readr(v.2.1.2), 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), lme4(v.1.1-30), Rsamtools(v.2.12.0), cli(v.3.3.0), urlchecker(v.1.0.1), patchwork(v.1.1.1), ps(v.1.7.1), formatR(v.1.12), MASS(v.7.3-58), mgcv(v.1.8-40), tidyselect(v.1.1.2), stringi(v.1.7.8), yaml(v.2.3.5), GOSemSim(v.2.22.0), locfit(v.1.5-9.6), ggrepel(v.0.9.1), biocViews(v.1.64.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), ggraph(v.2.0.5), digest(v.0.6.29), BiocManager(v.1.30.18), shiny(v.1.7.2), Rcpp(v.1.0.9), broom(v.1.0.0), BiocVersion(v.3.15.2), later(v.1.3.0), OrganismDbi(v.1.38.1), httr(v.1.4.3), Rdpack(v.2.4), colorspace(v.2.0-3), rvest(v.1.0.2), brio(v.1.1.3), XML(v.3.99-0.10), fs(v.1.5.2), splines(v.4.2.1), BiocCheck(v.1.32.0), yulab.utils(v.0.0.5), RBGL(v.1.72.0), 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), futile.options(v.1.0.1), jsonlite(v.1.8.0), nloptr(v.2.0.3), ggtree(v.3.4.1), tidygraph(v.1.2.1), ggfun(v.0.0.6), R6(v.2.5.1), RUnit(v.0.4.32), 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), BiocParallel(v.1.30.3), interactiveDisplayBase(v.1.34.0), 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), openxlsx(v.4.2.5), GO.db(v.3.15.0), 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), iterators(v.1.0.14), variancePartition(v.1.26.0), 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 605cc89b5f1cadea6923b53ac71e234ba0181fe7
## This is hpgltools commit: Wed Aug 10 22:39:40 2022 -0400: 605cc89b5f1cadea6923b53ac71e234ba0181fe7
message(paste0("Saving to ", savefile))
## Saving to tmrc2_datasets_202207.rda.xz
tmp <- sm(saveme(filename = savefile))
tmp <- loadme(filename = savefile)