Did the stuff on this morning’s TODO which came out of this morning’s meeting: do a PCA without the oddball strains (already done in the worksheet), highlight reference strains, and add L.major IDs and Descriptions (done by appending a collapsed version of the ortholog data to the all_lp_annot data).
Fixed human IDs for the macrophage data.
Changed input metadata sheets: primarily because I only remembered yesterday to finish the SL search for samples >TMRC20095. They are running now and will be added momentarily (I will have to redownload the sheet).
Setting up to make a hclust/phylogenetic tree of strains, use these are reference: 2168(2.3), 2272(2.2), for other 2.x choose arbitrarily (lower numbers are better).
Added another sanitize columns call for Antimony vs. antimony and None vs. none in the TMRC2 macrophage samples.
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:
Our shared online sample sheet is nearly static at the time of this writing (202209), I expect at this point the only likely updates will be to annotate some strains as more or less susceptible to drug treatment.
sample_sheet <- "sample_sheets/ClinicalStrains_TMRC2_Frozen\ 21062023.xlsx"
macrophage_sheet <- glue("sample_sheets/tmrc2_macrophage_samples_202304_modified.xlsx")
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.
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
testing <- load_orgdb_annotations(pan_db, keytype = "gid", fields = "all")
## Selecting the following fields, this might be too many:
## ANNOT_BFD3_CDS, ANNOT_BFD3_MODEL, ANNOT_BFD6_CDS, ANNOT_BFD6_MODEL, ANNOT_CDS, ANNOT_CDS_LENGTH, ANNOT_CHROMOSOME, ANNOT_DIF_CDS, ANNOT_DIF_MODEL, ANNOT_EC_NUMBERS, ANNOT_EC_NUMBERS_DERIVED, ANNOT_EXON_COUNT, ANNOT_FC_BFD3_CDS, ANNOT_FC_BFD3_MODEL, ANNOT_FC_BFD6_CDS, ANNOT_FC_BFD6_MODEL, ANNOT_FC_DIF_CDS, ANNOT_FC_DIF_MODEL, ANNOT_FC_PF_CDS, ANNOT_FC_PF_MODEL, ANNOT_FIVE_PRIME_UTR_LENGTH, ANNOT_GENE_ENTREZ_ID, ANNOT_GENE_EXON_COUNT, ANNOT_GENE_HTS_NONCODING_SNPS, ANNOT_GENE_HTS_NONSYN_SYN_RATIO, ANNOT_GENE_HTS_NONSYNONYMOUS_SNPS, ANNOT_GENE_HTS_STOP_CODON_SNPS, ANNOT_GENE_HTS_SYNONYMOUS_SNPS, ANNOT_GENE_LOCATION_TEXT, ANNOT_GENE_NAME, ANNOT_GENE_ORTHOLOG_NUMBER, ANNOT_GENE_ORTHOMCL_NAME, ANNOT_GENE_PARALOG_NUMBER, ANNOT_GENE_PREVIOUS_IDS, ANNOT_GENE_PRODUCT, ANNOT_GENE_SOURCE_ID, ANNOT_GENE_TOTAL_HTS_SNPS, ANNOT_GENE_TRANSCRIPT_COUNT, ANNOT_GENE_TYPE, ANNOT_GO_COMPONENT, ANNOT_GO_FUNCTION, ANNOT_GO_ID_COMPONENT, ANNOT_GO_ID_FUNCTION, ANNOT_GO_ID_PROCESS, ANNOT_GO_PROCESS, ANNOT_HAS_MISSING_TRANSCRIPTS, ANNOT_INTERPRO_DESCRIPTION, ANNOT_INTERPRO_ID, ANNOT_IS_PSEUDO, ANNOT_ISOELECTRIC_POINT, ANNOT_LOCATION_TEXT, ANNOT_MATCHED_RESULT, ANNOT_MOLECULAR_WEIGHT, ANNOT_NO_TET_CDS, ANNOT_NO_TET_MODEL, ANNOT_ORGANISM, ANNOT_PF_CDS, ANNOT_PF_MODEL, ANNOT_PFAM_DESCRIPTION, ANNOT_PFAM_ID, ANNOT_PIRSF_DESCRIPTION, ANNOT_PIRSF_ID, ANNOT_PREDICTED_GO_COMPONENT, ANNOT_PREDICTED_GO_FUNCTION, ANNOT_PREDICTED_GO_ID_COMPONENT, ANNOT_PREDICTED_GO_ID_FUNCTION, ANNOT_PREDICTED_GO_ID_PROCESS, ANNOT_PREDICTED_GO_PROCESS, ANNOT_PROJECT_ID, ANNOT_PROSITEPROFILES_DESCRIPTION, ANNOT_PROSITEPROFILES_ID, ANNOT_PROTEIN_LENGTH, ANNOT_PROTEIN_SEQUENCE, ANNOT_SEQUENCE_ID, ANNOT_SIGNALP_PEPTIDE, ANNOT_SIGNALP_SCORES, ANNOT_SMART_DESCRIPTION, ANNOT_SMART_ID, ANNOT_SOURCE_ID, ANNOT_STRAND, ANNOT_SUPERFAMILY_DESCRIPTION, ANNOT_SUPERFAMILY_ID, ANNOT_THREE_PRIME_UTR_LENGTH, ANNOT_TIGRFAM_DESCRIPTION, ANNOT_TIGRFAM_ID, ANNOT_TM_COUNT, ANNOT_TRANS_FOUND_PER_GENE_INTERNAL, ANNOT_TRANSCRIPT_INDEX_PER_GENE, ANNOT_TRANSCRIPT_LENGTH, ANNOT_TRANSCRIPT_LINK, ANNOT_TRANSCRIPT_PRODUCT, ANNOT_TRANSCRIPT_SEQUENCE, ANNOT_TRANSCRIPTS_FOUND_PER_GENE, ANNOT_UNIPROT_ID, ANNOT_URI, ANNOT_WDK_WEIGHT
## Extracted all gene ids.
## Attempting to select: ANNOT_BFD3_CDS, ANNOT_BFD3_MODEL, ANNOT_BFD6_CDS, ANNOT_BFD6_MODEL, ANNOT_CDS, ANNOT_CDS_LENGTH, ANNOT_CHROMOSOME, ANNOT_DIF_CDS, ANNOT_DIF_MODEL, ANNOT_EC_NUMBERS, ANNOT_EC_NUMBERS_DERIVED, ANNOT_EXON_COUNT, ANNOT_FC_BFD3_CDS, ANNOT_FC_BFD3_MODEL, ANNOT_FC_BFD6_CDS, ANNOT_FC_BFD6_MODEL, ANNOT_FC_DIF_CDS, ANNOT_FC_DIF_MODEL, ANNOT_FC_PF_CDS, ANNOT_FC_PF_MODEL, ANNOT_FIVE_PRIME_UTR_LENGTH, ANNOT_GENE_ENTREZ_ID, ANNOT_GENE_EXON_COUNT, ANNOT_GENE_HTS_NONCODING_SNPS, ANNOT_GENE_HTS_NONSYN_SYN_RATIO, ANNOT_GENE_HTS_NONSYNONYMOUS_SNPS, ANNOT_GENE_HTS_STOP_CODON_SNPS, ANNOT_GENE_HTS_SYNONYMOUS_SNPS, ANNOT_GENE_LOCATION_TEXT, ANNOT_GENE_NAME, ANNOT_GENE_ORTHOLOG_NUMBER, ANNOT_GENE_ORTHOMCL_NAME, ANNOT_GENE_PARALOG_NUMBER, ANNOT_GENE_PREVIOUS_IDS, ANNOT_GENE_PRODUCT, ANNOT_GENE_SOURCE_ID, ANNOT_GENE_TOTAL_HTS_SNPS, ANNOT_GENE_TRANSCRIPT_COUNT, ANNOT_GENE_TYPE, ANNOT_GO_COMPONENT, ANNOT_GO_FUNCTION, ANNOT_GO_ID_COMPONENT, ANNOT_GO_ID_FUNCTION, ANNOT_GO_ID_PROCESS, ANNOT_GO_PROCESS, ANNOT_HAS_MISSING_TRANSCRIPTS, ANNOT_INTERPRO_DESCRIPTION, ANNOT_INTERPRO_ID, ANNOT_IS_PSEUDO, ANNOT_ISOELECTRIC_POINT, ANNOT_LOCATION_TEXT, ANNOT_MATCHED_RESULT, ANNOT_MOLECULAR_WEIGHT, ANNOT_NO_TET_CDS, ANNOT_NO_TET_MODEL, ANNOT_ORGANISM, ANNOT_PF_CDS, ANNOT_PF_MODEL, ANNOT_PFAM_DESCRIPTION, ANNOT_PFAM_ID, ANNOT_PIRSF_DESCRIPTION, ANNOT_PIRSF_ID, ANNOT_PREDICTED_GO_COMPONENT, ANNOT_PREDICTED_GO_FUNCTION, ANNOT_PREDICTED_GO_ID_COMPONENT, ANNOT_PREDICTED_GO_ID_FUNCTION, ANNOT_PREDICTED_GO_ID_PROCESS, ANNOT_PREDICTED_GO_PROCESS, ANNOT_PROJECT_ID, ANNOT_PROSITEPROFILES_DESCRIPTION, ANNOT_PROSITEPROFILES_ID, ANNOT_PROTEIN_LENGTH, ANNOT_PROTEIN_SEQUENCE, ANNOT_SEQUENCE_ID, ANNOT_SIGNALP_PEPTIDE, ANNOT_SIGNALP_SCORES, ANNOT_SMART_DESCRIPTION, ANNOT_SMART_ID, ANNOT_SOURCE_ID, ANNOT_STRAND, ANNOT_SUPERFAMILY_DESCRIPTION, ANNOT_SUPERFAMILY_ID, ANNOT_THREE_PRIME_UTR_LENGTH, ANNOT_TIGRFAM_DESCRIPTION, ANNOT_TIGRFAM_ID, ANNOT_TM_COUNT, ANNOT_TRANS_FOUND_PER_GENE_INTERNAL, ANNOT_TRANSCRIPT_INDEX_PER_GENE, ANNOT_TRANSCRIPT_LENGTH, ANNOT_TRANSCRIPT_LINK, ANNOT_TRANSCRIPT_PRODUCT, ANNOT_TRANSCRIPT_SEQUENCE, ANNOT_TRANSCRIPTS_FOUND_PER_GENE, ANNOT_UNIPROT_ID, ANNOT_URI, ANNOT_WDK_WEIGHT
## 'select()' returned 1:1 mapping between keys and columns
lp_go <- load_orgdb_go(pan_db)
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(EuPathDB::extract_eupath_orthologs(db = pan_db))
Recently there was a request to include the Leishmania major gene IDs and descriptions. Thus I will extract them along with the orthologs and append that to the annotations used.
Having spent the time to run the following code, I realized that the orthologs data structure above actually already has the gene IDs and descriptions.
Thus I will leave my query in place to extract the major annotations, but follow it up with a collapse of the major orthologs and appending of that to the panamensis annotations.
orgdb <- "org.Lmajor.Friedlin.v49.eg.db"
tt <- sm(library(orgdb, character.only = TRUE))
major_db <- org.Lmajor.Friedlin.v49.eg.db
all_fields <- columns(pan_db)
all_lm_annot <- sm(load_orgdb_annotations(
major_db,
keytype = "gid",
fields = c("annot_gene_entrez_id", "annot_gene_name",
"annot_strand", "annot_chromosome", "annot_cds_length",
"annot_gene_product")))$genes
wanted_orthos_idx <- orthos[["ORTHOLOGS_SPECIES"]] == "Leishmania major strain Friedlin"
sum(wanted_orthos_idx)
## [1] 10796
wanted_orthos <- orthos[wanted_orthos_idx, ]
wanted_orthos <- wanted_orthos[, c("GID", "ORTHOLOGS_ID", "ORTHOLOGS_NAME")]
collapsed_orthos <- wanted_orthos %>%
group_by(GID) %>%
summarise(collapsed_id = stringr::str_c(ORTHOLOGS_ID, collapse=" ; "),
collapsed_name = stringr::str_c(ORTHOLOGS_NAME, collapse=" ; "))
all_lp_annot <- merge(all_lp_annot, collapsed_orthos, by.x = "row.names",
by.y = "GID", all.x = TRUE)
rownames(all_lp_annot) <- all_lp_annot[["Row.names"]]
all_lp_annot[["Row.names"]] <- NULL
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.
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" "DataProvider" "Genome" "GenomeSource"
## [7] "GenomeVersion" "NumArrayGene" "NumChipChipGene" "NumChromosome" "NumCodingGene" "NumCommunity"
## [13] "NumContig" "NumEC" "NumEST" "NumGene" "NumGO" "NumOrtholog"
## [19] "NumOtherGene" "NumPopSet" "NumProteomics" "NumPseudogene" "NumRNASeq" "NumRTPCR"
## [25] "NumSNP" "NumTFBS" "Organellar" "ReferenceStrain" "MegaBP" "PrimaryKey"
## [31] "ProjectID" "RecordClassName" "SourceID" "SourceVersion" "TaxonomyID" "TaxonomyName"
## [37] "URLGenome" "URLGFF" "URLProtein" "Coordinate_1_based" "Maintainer" "SourceUrl"
## [43] "Tags" "BsgenomePkg" "GrangesPkg" "OrganismdbiPkg" "OrgdbPkg" "TxdbPkg"
## [49] "Taxon" "Genus" "Species" "Strain" "BsgenomeFile" "GrangesFile"
## [55] "OrganismdbiFile" "OrgdbFile" "TxdbFile" "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")
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 samples are much lower coverage:
There is a set of strains which acquired resistance in vitro. These are included in the dataset, but there are not likely enough of them to query that question explicitly.
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"),
"susceptibility" = list(
"resistant" = "#8563a7",
"sensitive" = "#8d0000",
"ambiguous" = "#cbcbcb",
"unknown" = "#555555"))
data_structures <- c(data_structures, "color_choices")
The data structure ‘lp_expt’ contains the data for all samples which have hisat2 count tables, and which pass a few initial quality tests (e.g. they must have more than 8550 genes with >0 counts and >5e6 reads which mapped to a gene); genes which are annotated with a few key redundant categories (leishmanolysin for example) are also culled.
There are a few metadata columns which we really want to make certain are standardized.
sanitize_columns <- c("passagenumber", "clinicalresponse", "clinicalcategorical",
"zymodemecategorical")
lp_expt <- create_expt(sample_sheet,
gene_info = all_lp_annot,
annotation_name = orgdb,
savefile = glue("rda/tmrc2_lp_expt_all_raw-v{ver}.rda"),
id_column = "hpglidentifier",
annotation = "org.Lpanamensis.MHOMCOL81L13.v46.eg.db",
file_column = "lpanamensisv36hisatfile") %>%
set_expt_conditions(fact = "zymodemecategorical", colors = color_choices[["strain"]]) %>%
subset_expt(nonzero = 8550) %>%
subset_expt(coverage = 5000000) %>%
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 72 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 numbers of samples by condition are:
##
## b2904 unknown z1.0 z1.5 z2.0 z2.1 z2.2 z2.3 z2.4 z3.0 z3.2
## 1 4 1 1 1 7 45 41 2 1 1
## The samples (and read coverage) removed when filtering 8550 non-zero genes are:
## subset_expt(): There were 105, now there are 103 samples.
## 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.
## semantic_expt_filter(): Removed 68 genes.
data_structures <- c(data_structures, "lp_expt")
save(list = "lp_expt", file = 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 z32
## 1 2 1 1 1 7 43 41 2 1 1
table(pData(lp_expt)[["clinicalresponse"]])
##
## cure failure failure miltefosine
## 40 36 1
## laboratory line laboratory line miltefosine resistant nd
## 1 1 18
## reference strain
## 4
ncol(exprs(lp_expt))
## [1] 101
All the following data will derive from this starting point.
Column ‘Q’ in the sample sheet, make a categorical version of it with these parameters:
Note that these cutoffs are only valid for the historical data. The newer susceptibility data uses a cutoff of 0.78 for sensitive. I will set ambiguous to 0.5 to 0.78?
max_resist_historical <- 0.35
min_sensitive_historical <- 0.49
## 202305: Removed ambiguous category for the current set.g
max_resist_current <- 0.76
min_sensitive_current <- 0.77
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:
Thus, the following block will sanitize these percentage values into a single decimal number and make a categorical variable from it using pre-defined values for resistant/ambiguous/sensitive. This categorical variable will be stored in a new column: ‘sus_category_historical’.
st <- pData(lp_expt)[["susceptibilityinfectionreduction32ugmlsbvhistoricaldata"]]
starting <- sanitize_percent(st)
st
## [1] "0.45" "0.14000000000000001" "0.97" NA NA NA
## [7] NA NA NA "0" "0.97" "0"
## [13] "0" "0.46" "0.45" "0.97" "0.56000000000000005" "0.99"
## [19] "0.46" "0.7" "0.99" "0.99" "0.45" "0.98"
## [25] "0.99" "0.49" "No data" "No data" "0.99" "0.66"
## [31] "0.99" "No data" "0.99" "1" "1" "0.94"
## [37] "0.94" "No data" "No data" "No data" "No data" "No data"
## [43] "No data" "No data" "No data" "No data" "No data" "No data"
## [49] "No data" "No data" "No data" "0.99" "0.99" "No data"
## [55] "0.98" "0.97" "0.96" "0.96" "0" "0"
## [61] "0" "0.06" "0.94" "0.94" "0.03" "0.94"
## [67] "0" "0.25" "0.95" "0.27" "No data" "No data"
## [73] "No data" "No data" "No data" "No data" "No data" "No data"
## [79] "No data" "No data" "No data" "No data" "No data" "No data"
## [85] "No data" "No data" "No data" "No data" "No data" "No data"
## [91] "No data" "No data" "No data" "No data" "No data" "No data"
## [97] "No data" "No data" "No data" "No data" "No data"
starting
## [1] 0.45 0.14 0.97 NA NA NA NA NA NA 0.00 0.97 0.00 0.00 0.46 0.45 0.97 0.56 0.99 0.46 0.70 0.99 0.99 0.45 0.98 0.99 0.49 NA NA
## [29] 0.99 0.66 0.99 NA 0.99 1.00 1.00 0.94 0.94 NA NA NA NA NA NA NA NA NA NA NA NA NA NA 0.99 0.99 NA 0.98 0.97
## [57] 0.96 0.96 0.00 0.00 0.00 0.06 0.94 0.94 0.03 0.94 0.00 0.25 0.95 0.27 NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [85] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
sus_categorical <- starting
na_idx <- is.na(starting)
sum(na_idx)
## [1] 55
sus_categorical[na_idx] <- "unknown"
resist_idx <- starting <= max_resist_historical
sus_categorical[resist_idx] <- "resistant"
indeterminant_idx <- starting > max_resist_historical &
starting < min_sensitive_historical
sus_categorical[indeterminant_idx] <- "ambiguous"
susceptible_idx <- starting >= min_sensitive_historical
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
## 5 12 29 55
The same process will be repeated for the current iteration of the sensitivity assay and stored in the ‘sus_category_current’ column.
starting_current <- sanitize_percent(pData(lp_expt)[["susceptibilityinfectionreduction32ugmlsbvcurrentdata"]])
sus_categorical_current <- starting_current
na_idx <- is.na(starting_current)
sum(na_idx)
## [1] 5
sus_categorical_current[na_idx] <- "unknown"
resist_idx <- starting_current <= max_resist_current
sus_categorical_current[resist_idx] <- "resistant"
indeterminant_idx <- starting_current > max_resist_current &
starting_current < min_sensitive_current
sus_categorical_current[indeterminant_idx] <- "ambiguous"
susceptible_idx <- starting_current >= min_sensitive_current
sus_categorical_current[susceptible_idx] <- "sensitive"
sus_categorical_current <- as.factor(sus_categorical_current)
pData(lp_expt)[["sus_category_current"]] <- sus_categorical_current
pData(lp_expt)[["susceptibility"]] <- sus_categorical_current
table(sus_categorical_current)
## sus_categorical_current
## resistant sensitive unknown
## 47 49 5
lp_sankey <- plot_meta_sankey(
lp_expt, factors = c("zymodemecategorical", "clinicalcategorical", "susceptibility"),
drill_down = TRUE, color_choices = color_choices, html = "meta_sankey.html")
## Warning: attributes are not identical across measure variables; they will be dropped
lp_sankey
## A sankey plot describing the metadata of 101 samples,
## including 46 out of 143 nodes and traversing metadata factors:
## zymodemecategorical, clinicalcategorical, susceptibility.
Here is a table of my current classifier’s interpretation of the strains.
table(pData(lp_expt)[["knnv2classification"]])
##
## z10 z21 z22 z23 z24 z32
## 3 6 47 41 2 2
merged_zymo <- lp_expt
pData(merged_zymo)[["zymodeme"]] <- as.character(pData(merged_zymo)[["zymodemecategorical"]])
z21_idx <- pData(merged_zymo)[["zymodeme"]] == "z21"
pData(merged_zymo)[z21_idx, "zymodeme"] <- "z22"
z24_idx <- pData(merged_zymo)[["zymodeme"]] == "z24"
pData(merged_zymo)[z24_idx, "zymodeme"] <- "z23"
keepers <- pData(merged_zymo)[["zymodeme"]] == "z22" |
pData(merged_zymo)[["zymodeme"]] == "z23"
merged_zymo <- merged_zymo[, keepers] %>%
set_expt_conditions(fact = "zymodeme", colors = color_choices[["zymo"]])
## Subsetting on samples.
## subset_expt(): There were 101, now there are 93 samples.
## The numbers of samples by condition are:
##
## z22 z23
## 50 43
two_sankey <- plot_meta_sankey(
merged_zymo, factors = c("zymodeme", "clinicalcategorical", "susceptibility"),
drill_down = TRUE, color_choices = color_choices, html = "meta_sankey.html")
two_sankey
## A sankey plot describing the metadata of 93 samples,
## including 16 out of 20 nodes and traversing metadata factors:
## zymodeme, clinicalcategorical, susceptibility.
In many queries, we will seek to compare only the two primary strains, zymodeme 2.2 and 2.3. The following block will extract only those samples.
lp_strain <- lp_expt %>%
set_expt_batches(fact = sus_categorical_current) %>%
set_expt_colors(color_choices[["strain"]])
## The number of samples by batch are:
##
## resistant sensitive unknown
## 47 49 5
table(pData(lp_strain)[["condition"]])
##
## b2904 unknown z1.0 z1.5 z2.0 z2.1 z2.2 z2.3 z2.4 z3.0 z3.2
## 1 2 1 1 1 7 43 41 2 1 1
save(list = "lp_strain", file = glue("rda/tmrc2_lp_strain-v{ver}.rda"))
data_structures <- c(data_structures, "lp_strain")
lp_two_strains <- merged_zymo
save(list = "lp_two_strains",
file = glue("rda/tmrc2_lp_two_strains-v{ver}.rda"))
data_structures <- c(data_structures, "lp_two_strains")
Clinical outcome is by far the most problematic comparison in this data, but here is the recategorization of the data using it:
lp_cf <- set_expt_conditions(lp_expt, fact = "clinicalcategorical",
colors = color_choices[["cf"]]) %>%
set_expt_batches(fact = sus_categorical_current)
## The numbers of samples by condition are:
##
## cure fail unknown
## 40 37 24
## Warning in set_expt_colors(new_expt, colors = colors): Colors for the following categories are not being used: notapplicable.
## The number of samples by batch are:
##
## resistant sensitive unknown
## 47 49 5
table(pData(lp_cf)[["condition"]])
##
## cure fail unknown
## 40 37 24
data_structures <- c(data_structures, "lp_cf")
save(list = "lp_cf",
file = 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 77 samples.
data_structures <- c(data_structures, "lp_cf_known")
save(list = "lp_cf_known",
file = glue("rda/tmrc2_lp_cf_known-v{ver}.rda"))
Use the factorized version of susceptibility to categorize the samples by the historical data.
lp_susceptibility_historical <- set_expt_conditions(
lp_expt, fact = "sus_category_historical", colors = color_choices[["susceptibility"]]) %>%
set_expt_batches(fact = "clinicalcategorical")
## The numbers of samples by condition are:
##
## ambiguous resistant sensitive unknown
## 5 12 29 55
## The number of samples by batch are:
##
## cure fail unknown
## 40 37 24
save(list = "lp_susceptibility_historical",
file = glue("rda/tmrc2_lp_susceptibility_historical-v{ver}.rda"))
data_structures <- c(data_structures, "lp_susceptibility_historical")
Use the factorized version of susceptibility to categorize the samples by the historical data.
This will likely be our canonical susceptibility dataset, so I will remove the suffix and just call it ‘lp_susceptibility’.
lp_susceptibility <- set_expt_conditions(
lp_expt, fact = "sus_category_current", colors = color_choices[["susceptibility"]]) %>%
set_expt_batches(fact = "clinicalcategorical")
## The numbers of samples by condition are:
##
## resistant sensitive unknown
## 47 49 5
## Warning in set_expt_colors(new_expt, colors = colors): Colors for the following categories are not being used: ambiguous.
## The number of samples by batch are:
##
## cure fail unknown
## 40 37 24
save(list = "lp_susceptibility",
file = glue("rda/tmrc2_lp_susceptibility-v{ver}.rda"))
data_structures <- c(data_structures, "lp_susceptibility")
I think this is redundant with a previous block, but I am leaving it until I am certain that it is not required in a following document.
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("rda/tmrc2_lp_zymo-v{ver}.rda"))
The following section will create some initial data structures of the observed variants in the parasite samples. This will include some of our 2016 samples for some classification queries.
I changed and improved the mapping and variant detection methods from what we used for the 2016 data. So some small changes will be required to merge them.
lp_previous <- create_expt("sample_sheets/tmrc2_samples_20191203.xlsx",
file_column = "tophat2file",
savefile = 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.
tt <- lp_previous$expressionset
rownames(tt) <- gsub(pattern = "^exon_", replacement = "", x = rownames(tt))
rownames(tt) <- gsub(pattern = "\\.1$", replacement = "", x = rownames(tt))
rownames(tt) <- gsub(pattern = "\\-1$", replacement = "", x = rownames(tt))
lp_previous$expressionset <- tt
rm(tt)
data_structures <- c(data_structures, "lp_previous")
The count_expt_snps() function uses our expressionset data and a metadata column in order to extract the mpileup or freebayes-based variant calls and create matrices of the likelihood that each position-per-sample is in fact a variant.
There is an important caveat here which changed on 202301: I was interpreting using the PAIRED tag, which is only used for, unsurprisingly, paired-end samples. A couple samples are not paired and so were failing silently. The QA tag looks like it is more appropriate and should work across both types. One way to find out, I am setting it here and will look to see if the results make more sense for my test samples (TMRC2001, TMRC2005, TMRC2007).
## 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 = "QA")
## 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:
## 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`
## Lets see if we get numbers which make sense.
summary(exprs(new_snps)[["tmrc20001"]]) ## My weirdo sample
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0 0.0 0.0 13.9 0.0 2217.0
summary(exprs(new_snps)[["tmrc20072"]]) ## Another sample chosen at random
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 64 0 247568
summary(exprs(new_snps)[["tmrc20021"]]) ## Another sample chosen at random
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 686 0 1708458
## Now that we are reasonably confident that things make more sense, lets save and move on...
data_structures <- c(data_structures, "new_snps")
tt <- normalize_expt(new_snps, transform = "log2")
## transform_counts: Found 146144951 values equal to 0, adding 1 to the matrix.
plot_boxplot(tt)
Now let us pull in the 2016 data.
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("rda/lp_snp-v{ver}.rda"))
data_structures <- c(data_structures, "lp_snp")
save(list = "new_snps",
file = glue("rda/new_snps-v{ver}.rda"))
data_structures <- c(data_structures, "new_snps")
save(list = "old_snps",
file = 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 tmrc20040 tmrc20066 tmrc20039 tmrc20037 tmrc20038 tmrc20067
## 74798 93649 12389 8675 975 351343 338580 146302 58753 93615 25115 98958 97676 93954
## tmrc20068 tmrc20041 tmrc20015 tmrc20009 tmrc20010 tmrc20016 tmrc20011 tmrc20012 tmrc20013 tmrc20017 tmrc20014 tmrc20018 tmrc20019 tmrc20070
## 96583 53184 96398 16031 94079 146124 14055 456 95040 48288 17245 140438 14829 97336
## tmrc20020 tmrc20021 tmrc20022 tmrc20025 tmrc20024 tmrc20036 tmrc20069 tmrc20033 tmrc20026 tmrc20031 tmrc20076 tmrc20073 tmrc20055 tmrc20079
## 15484 101127 18143 364240 18471 60087 18792 33663 15074 19139 17559 96169 22246 96224
## tmrc20071 tmrc20078 tmrc20094 tmrc20042 tmrc20058 tmrc20072 tmrc20059 tmrc20048 tmrc20057 tmrc20088 tmrc20056 tmrc20060 tmrc20077 tmrc20074
## 94353 18836 87878 19734 94524 50292 94091 97164 48944 15594 22683 21506 17925 21015
## tmrc20063 tmrc20053 tmrc20052 tmrc20064 tmrc20075 tmrc20051 tmrc20050 tmrc20049 tmrc20062 tmrc20110 tmrc20080 tmrc20043 tmrc20083 tmrc20054
## 28254 20181 100709 93173 96625 94125 17200 16168 93677 15487 95775 95623 21167 93603
## tmrc20085 tmrc20046 tmrc20093 tmrc20089 tmrc20047 tmrc20090 tmrc20044 tmrc20045 tmrc20061 tmrc20105 tmrc20108 tmrc20109 tmrc20098 tmrc20096
## 89765 48608 48254 90421 92637 91564 14861 50403 116906 86758 96831 17709 92927 17534
## tmrc20097 tmrc20101 tmrc20092 tmrc20082 tmrc20102 tmrc20099 tmrc20100 tmrc20091 tmrc20084 tmrc20087 tmrc20103 tmrc20104 tmrc20086 tmrc20107
## 46863 17753 16578 109909 92380 91383 94381 15059 46548 14947 49091 93979 15813 95144
## tmrc20081 tmrc20106 tmrc20095
## 19533 18545 81200
As far as I can tell, freebayes and mpileup are reasonably similar in their sensitivity/specificity; so combining the two datasets like this is expected to work with minimal problems. The most likely problem is that my mpileup-based pipeline is unable to handle indels.
## My old_snps is using an older annotation incorrectly, so fix it here:
#annotation(old_snps) <- annotation(new_snps)
both_snps <- combine_expts(new_snps, old_snps)
## Warning in combine_expts(new_snps, old_snps): There are many gene IDs which are not shared among the two datasets, this may fail.
## There are many gene IDs which are not shared among the two datasets.
## Here are some from the first: chr_LPAL13-SCAF000001_pos_1019_ref_G_alt_A, chr_LPAL13-SCAF000001_pos_1023_ref_C_alt_G, chr_LPAL13-SCAF000001_pos_1038_ref_C_alt_T, chr_LPAL13-SCAF000001_pos_1066_ref_T_alt_C, chr_LPAL13-SCAF000001_pos_106_ref_A_alt_G, chr_LPAL13-SCAF000001_pos_1092_ref_A_alt_G
## Here are some from the second: chr_LPAL13-SCAF000001_pos_1019_ref_G_alt_A, chr_LPAL13-SCAF000001_pos_106_ref_A_alt_G, chr_LPAL13-SCAF000001_pos_1092_ref_A_alt_G, chr_LPAL13-SCAF000001_pos_1138_ref_C_alt_A, chr_LPAL13-SCAF000001_pos_1149_ref_G_alt_A, chr_LPAL13-SCAF000001_pos_1183_ref_C_alt_T
## The numbers of samples by condition are:
##
## z2.3 z2.2 unknown z1.0 b2904 z3.0 z2.0 z1.5 z2.1 z2.4 z3.2 undef sh chr inf
## 41 43 2 1 1 1 1 1 7 2 1 0 13 14 6
save(list = "both_snps",
file = glue("rda/both_snps-v{ver}.rda"))
data_structures <- c(data_structures, "both_snps")
I am taking a heatmap from our variant data and manually identifying sample groups.
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.
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[["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(
macrophage_sheet,
gene_info = hs_annot,
file_column = "hg38100hisatfile") %>%
set_expt_conditions(fact = "macrophagetreatment") %>%
set_expt_batches(fact = "macrophagezymodeme") %>%
sanitize_expt_metadata(columns = sanitize_columns) %>%
subset_expt(nonzero = 12000)
## Reading the sample metadata.
## Dropped 1 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: 69 rows(samples) and 80 columns(metadata fields).
## Warning in create_expt(macrophage_sheet, gene_info = hs_annot, file_column = "hg38100hisatfile"): Even after changing the rownames in gene info,
## they do not match the count table.
## Even after changing the rownames in gene info, they do not match the count table.
## Here are the first few rownames from the count tables:
## ENSG00000000003, ENSG00000000005, ENSG00000000419, ENSG00000000457, ENSG00000000460, ENSG00000000938
## Here are the first few rownames from the gene information table:
## gene:ENSG00000004059, gene:ENSG00000003056, gene:ENSG00000173153, gene:ENSG00000004478, gene:ENSG00000003137, gene:ENSG00000003509
## 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:
## 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"
data_structures <- c(data_structures, "hs_macrophage")
Finally, split off the U937 samples.
hs_u937 <- subset_expt(hs_macrophage, subset = "typeofcells!='Macrophages'")
## subset_expt(): There were 68, now there are 14 samples.
data_structures <- c(data_structures, "hs_u937")
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(macrophage_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.
## Dropped 1 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: 69 rows(samples) and 80 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 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"))
## Deleting the file analyses/macrophage_de/202305/read_counts/lp_macrophage_reads_unfiltered-v202305.xlsx before writing the tables.
## 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" "TMRC30268" "TMRC30326" "TMRC30323" "TMRC30319" "TMRC30325"
## [12] "TMRC30327" "TMRC30312" "TMRC30300" "TMRC30304" "TMRC30302" "TMRC30313" "TMRC30309" "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.
## Not putting labels on the plot.
##
## 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.
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
## `geom_smooth()` using formula = 'y ~ x'
## Error in checkModelStatus(fit, showWarnings = showWarnings, colinearityCutoff = colinearityCutoff, :
## The variables specified in this model are redundant,
## so the design matrix is not full rank
## Retrying with only condition in the model.
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following object is masked from 'package:S4Vectors':
##
## expand
##
##
## Total:84 s
## 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'
## Error in checkModelStatus(fit, showWarnings = showWarnings, colinearityCutoff = colinearityCutoff, :
## The variables specified in this model are redundant,
## so the design matrix is not full rank
## Retrying with only condition in the model.
##
## Total:71 s
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:
## 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"))
## Deleting the file analyses/macrophage_de/202305/read_counts/lp_macrophage_reads_filtered-v202305.xlsx before writing the tables.
## 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.
## Not putting labels on the plot.
##
## 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'
##
## Total:106 s
## 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'
##
## Total:95 s
lp_macrophage <- lp_macrophage_filt
lp_macrophage_nosb <- subset_expt(lp_macrophage, subset="batch!='inf_sb'")
## 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"))
## Deleting the file analyses/macrophage_de/202305/read_counts/lp_macrophage_nosb_reads-v202305.xlsx before writing the tables.
## 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.Not putting labels on the plot.
## 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.
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
## `geom_smooth()` using formula = 'y ~ x'
## varpart sees only 1 batch, adjusting the model accordingly.
## Error in .fitExtractVarPartModel(exprObj, formula, data, REML = REML, :
## Variable in formula not found in data: condition
## Retrying with only condition in the model.
##
## Total:77 s
## `geom_smooth()` using formula = 'y ~ x'varpart sees only 1 batch, adjusting the model accordingly.
## Error in .fitExtractVarPartModel(exprObj, formula, data, REML = REML, :
## Variable in formula not found in data: condition
## Retrying with only condition in the model.
##
## Total:98 s
data_structures <- c(data_structures, "lp_macrophage_nosb")
spec <- make_rnaseq_spec()
test <- gather_preprocessing_metadata(macrophage_sheet, specification = spec)
## Using provided specification
## Dropped 1 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.
## Example filename: preprocessing/TMRC30051/outputs/*trimomatic/*-trimomatic.stderr.
## Example filename: preprocessing/TMRC30051/outputs/*trimomatic/*-trimomatic.stderr.
## The numerator column is: trimomatic_output.
## The denominator column is: trimomatic_input.
## Example filename: preprocessing/TMRC30051/outputs/*fastqc/*_fastqc/fastqc_data.txt.
## Skipping for now
## Not including new entries for: fastqc_most_overrepresented, it is empty.
## Example filename: preprocessing/TMRC30051/outputs/*hisat2_*/hisat2_*rRNA*.stderr.
## Not including new entries for: hisat_rrna_single_concordant, it is empty.
## Example filename: preprocessing/TMRC30051/outputs/*hisat2_*/hisat2_*rRNA*.stderr.
## Not including new entries for: hisat_rrna_multi_concordant, it is empty.
## Missing data to calculate the ratio between: hisat_rrna_multi_concordant and trimomatic_output.
## Not including new entries for: hisat_rrna_percent, it is empty.
## Example filename: preprocessing/TMRC30051/outputs/*hisat2_*/hisat2_*genome*.stderr.
## Example filename: preprocessing/TMRC30051/outputs/*hisat2_*/hisat2_*genome*.stderr.
## Example filename: preprocessing/TMRC30051/outputs/*hisat2_*/hisat2_*genome*.stderr.
## Example filename: preprocessing/TMRC30051/outputs/*hisat2_*/hisat2_*genome*.stderr.
## Example filename: preprocessing/TMRC30051/outputs/*hisat2_*/hisat2_*genome*.stderr.
## The numerator column is: hisat_genome_single_concordant.
## The denominator column is: trimomatic_output.
## Example filename: preprocessing/TMRC30051/outputs/*hisat2_*/*_genome*.count.xz.
## Example filename: preprocessing/TMRC30051/outputs/*hisat2_*/*_genome*.count.xz.
## Example filename: preprocessing/TMRC30051/outputs/*salmon_*/salmon_*.stderr.
## Example filename: preprocessing/TMRC30051/outputs/*salmon_*/quant.sf.
## Example filename: preprocessing/TMRC30051/outputs/*salmon_*/salmon_*.stderr.
## Example filename: preprocessing/TMRC30051/outputs/*salmon_*/salmon_*.stderr.
## Writing new metadata to: sample_sheets/tmrc2_macrophage_samples_202304_modified_modified.xlsx
## Deleting the file sample_sheets/tmrc2_macrophage_samples_202304_modified_modified.xlsx before writing the tables.
lp_meta <- pData(lp_macrophage)
lp_meta[["slvsreads_log"]] <- log10(lp_meta[["slvsreads"]])
inf_values <- is.infinite(lp_meta[["slvsreads_log"]])
lp_meta[inf_values, "slvsreads_log"] <- -10
color_vector <- as.character(color_choices[["strain"]])
names(color_vector) <- names(color_choices[["strain"]])
color_vector <- color_vector[c("z2.2", "z2.3", "unknown")]
names(color_vector) <- c("z2.2", "z2.3", "none")
sl_violin <- ggplot(lp_meta,
aes(x = .data[["condition"]], y = .data[["slvsreads_log"]],
fill = .data[["condition"]])) +
geom_violin() +
geom_point() +
scale_fill_manual(values = color_vector)
sl_violin
ggstatsplot::ggbetweenstats(lp_meta, x = "condition", y = "slvsreads_log")
save(list = data_structures, file = glue("rda/tmrc2_data_structures-v{ver}.rda"))
pander::pander(sessionInfo())
R version 4.2.0 (2022-04-22)
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-33), Matrix(v.1.5-4), BiocParallel(v.1.32.6), variancePartition(v.1.28.9), BSGenome.Leishmania.panamensis.MHOMCOL81L13.v53(v.2021.07), BSgenome(v.1.66.3), rtracklayer(v.1.58.0), Biostrings(v.2.66.0), XVector(v.0.38.0), futile.logger(v.1.4.3), org.Lmajor.Friedlin.v49.eg.db(v.2020.11), org.Lpanamensis.MHOMCOL81L13.v46.eg.db(v.2020.07), AnnotationDbi(v.1.60.2), EuPathDB(v.1.6.0), GenomeInfoDbData(v.1.2.9), hpgltools(v.1.0), testthat(v.3.1.8), Heatplus(v.3.6.0), ggplot2(v.3.4.2), glue(v.1.6.2), SummarizedExperiment(v.1.28.0), GenomicRanges(v.1.50.2), GenomeInfoDb(v.1.34.9), IRanges(v.2.32.0), S4Vectors(v.0.36.2), MatrixGenerics(v.1.10.0), matrixStats(v.0.63.0), Biobase(v.2.58.0) and BiocGenerics(v.0.44.0)
loaded via a namespace (and not attached): stringdist(v.0.9.10), corpcor(v.1.6.10), ps(v.1.7.5), Rsamtools(v.2.14.0), foreach(v.1.5.2), rprojroot(v.2.0.3), crayon(v.1.5.2), rbibutils(v.2.2.13), MASS(v.7.3-60), nlme(v.3.1-162), backports(v.1.4.1), sva(v.3.46.0), GOSemSim(v.2.24.0), rlang(v.1.1.1), HDO.db(v.0.99.1), nloptr(v.2.0.3), callr(v.3.7.3), limma(v.3.54.2), filelock(v.1.0.2), rjson(v.0.2.21), bit64(v.4.0.5), pbkrtest(v.0.5.2), parallel(v.4.2.0), processx(v.3.8.1), ggstatsplot(v.0.11.1), DOSE(v.3.24.2), tidyselect(v.1.2.0), usethis(v.2.1.6), BiocCheck(v.1.34.3), XML(v.3.99-0.14), tidyr(v.1.3.0), zoo(v.1.8-12), GenomicAlignments(v.1.34.1), MatrixModels(v.0.5-1), xtable(v.1.8-4), magrittr(v.2.0.3), evaluate(v.0.21), Rdpack(v.2.4), cli(v.3.6.1), zlibbioc(v.1.44.0), rstudioapi(v.0.14), miniUI(v.0.1.1.1), bslib(v.0.4.2), fastmatch(v.1.1-3), aod(v.1.3.2), lambda.r(v.1.2.4), treeio(v.1.22.0), shiny(v.1.7.4), xfun(v.0.39), parameters(v.0.21.0), pkgbuild(v.1.4.0), gson(v.0.1.0), caTools(v.1.18.2), tidygraph(v.1.2.3), AnnotationHubData(v.1.28.0), KEGGREST(v.1.38.0), clusterGeneration(v.1.3.7), tibble(v.3.2.1), interactiveDisplayBase(v.1.36.0), ggrepel(v.0.9.3), ape(v.5.7-1), png(v.0.1-8), zeallot(v.0.1.0), withr(v.2.5.0), bitops(v.1.0-7), ggforce(v.0.4.1), RBGL(v.1.74.0), plyr(v.1.8.8), GSEABase(v.1.60.0), coda(v.0.19-4), pillar(v.1.9.0), biocViews(v.1.66.3), gplots(v.3.1.3), cachem(v.1.0.8), GenomicFeatures(v.1.50.4), multcomp(v.1.4-23), fs(v.1.6.2), paletteer(v.1.5.0), clusterProfiler(v.4.6.2), RUnit(v.0.4.32), vctrs(v.0.6.2), ellipsis(v.0.3.2), generics(v.0.1.3), devtools(v.2.4.5), tools(v.4.2.0), remaCor(v.0.0.11), munsell(v.0.5.0), tweenr(v.2.0.2), fgsea(v.1.24.0), networkD3(v.0.4), emmeans(v.1.8.5), DelayedArray(v.0.24.0), fastmap(v.1.1.1), compiler(v.4.2.0), pkgload(v.1.3.2), httpuv(v.1.6.10), sessioninfo(v.1.2.2), plotly(v.4.10.1), gridExtra(v.2.3), edgeR(v.3.40.2), lattice(v.0.21-8), AnnotationForge(v.1.40.2), utf8(v.1.2.3), later(v.1.3.1), dplyr(v.1.1.2), prismatic(v.1.1.1), BiocFileCache(v.2.6.1), jsonlite(v.1.8.4), scales(v.1.2.1), graph(v.1.76.0), pbapply(v.1.7-0), tidytree(v.0.4.2), estimability(v.1.4.1), genefilter(v.1.80.3), lazyeval(v.0.2.2), promises(v.1.2.0.1), doParallel(v.1.0.17), effectsize(v.0.8.3), sandwich(v.3.0-2), rmarkdown(v.2.21), openxlsx(v.4.2.5.2), cowplot(v.1.1.1), Rtsne(v.0.16), pander(v.0.6.5), downloader(v.0.4), igraph(v.1.4.2), survival(v.3.5-5), yaml(v.2.3.7), htmltools(v.0.5.5), memoise(v.2.0.1), profvis(v.0.3.8), BiocIO(v.1.8.0), locfit(v.1.5-9.7), graphlayouts(v.1.0.0), quadprog(v.1.5-8), viridisLite(v.0.4.2), digest(v.0.6.31), RhpcBLASctl(v.0.23-42), mime(v.0.12), rappdirs(v.0.3.3), futile.options(v.1.0.1), bayestestR(v.0.13.1), RSQLite(v.2.3.1), yulab.utils(v.0.0.6), remotes(v.2.4.2), data.table(v.1.14.8), urlchecker(v.1.0.1), blob(v.1.2.4), labeling(v.0.4.2), rematch2(v.2.1.2), AnnotationHub(v.3.6.0), OrganismDbi(v.1.40.0), ggsankey(v.0.0.99999), RCurl(v.1.98-1.12), broom(v.1.0.4), hms(v.1.1.3), colorspace(v.2.1-0), BiocManager(v.1.30.20), aplot(v.0.1.10), sass(v.0.4.6), Rcpp(v.1.0.10), mvtnorm(v.1.1-3), enrichplot(v.1.18.4), fansi(v.1.0.4), tzdb(v.0.3.0), brio(v.1.1.3), R6(v.2.5.1), grid(v.4.2.0), lifecycle(v.1.0.3), formatR(v.1.14), statsExpressions(v.1.5.0), BayesFactor(v.0.9.12-4.4), zip(v.2.3.0), datawizard(v.0.7.1), curl(v.5.0.0), minqa(v.1.2.5), jquerylib(v.0.1.4), fastcluster(v.1.2.3), PROPER(v.1.30.0), qvalue(v.2.30.0), TH.data(v.1.1-2), desc(v.1.4.2), RColorBrewer(v.1.1-3), iterators(v.1.0.14), stringr(v.1.5.0), directlabels(v.2021.1.13), htmlwidgets(v.1.6.2), polyclip(v.1.10-4), biomaRt(v.2.54.1), purrr(v.1.0.1), shadowtext(v.0.1.2), gridGraphics(v.0.5-1), rvest(v.1.0.3), mgcv(v.1.8-42), insight(v.0.19.1), patchwork(v.1.1.2), codetools(v.0.2-19), GO.db(v.3.16.0), gtools(v.3.9.4), prettyunits(v.1.1.1), dbplyr(v.2.3.2), correlation(v.0.8.4), gtable(v.0.3.3), DBI(v.1.1.3), ggfun(v.0.0.9), httr(v.1.4.6), highr(v.0.10), KernSmooth(v.2.23-21), stringi(v.1.7.12), vroom(v.1.6.3), progress(v.1.2.2), reshape2(v.1.4.4), farver(v.2.1.1), annotate(v.1.76.0), viridis(v.0.6.3), ggtree(v.3.6.2), xml2(v.1.3.4), boot(v.1.3-28.1), restfulr(v.0.0.15), readr(v.2.1.4), ggplotify(v.0.1.0), BiocVersion(v.3.16.0), bit(v.4.0.5), scatterpie(v.0.1.9), ggraph(v.2.1.0), pkgconfig(v.2.0.3) and knitr(v.1.42)
message("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 d4aecb6c9f50378b78f8d3c152430d455ccf17a9
## This is hpgltools commit: Tue Jul 11 17:12:15 2023 -0400: d4aecb6c9f50378b78f8d3c152430d455ccf17a9
message("Saving to ", savefile)
## Saving to tmrc2_datasets_202305.rda.xz
# tmp <- sm(saveme(filename = savefile))
tmp <- loadme(filename = savefile)