sample_sheet <- glue::glue("sample_sheets/tmrc2_samples_20210817.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))
tt <- sm(library(org.Lpanamensis.MHOMCOL81L13.v46.eg.db))
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))
hisat_annot <- all_lp_annot
## rownames(hisat_annot) <- paste0("exon_", rownames(hisat_annot), ".E1")
meta <- EuPathDB::download_eupath_metadata(webservice="tritrypdb")
## Appending to an existing file: EuPathDB/metadata/biocv3.13_tritrypdbv54_metadata.csv
## Appending to an existing file: EuPathDB/metadata/GRanges_biocv3.13_tritrypdbv54_metadata.csv
## Appending to an existing file: EuPathDB/metadata/OrgDb_biocv3.13_tritrypdbv54_metadata.csv
## Appending to an existing file: EuPathDB/metadata/TxDb_biocv3.13_tritrypdbv54_metadata.csv
## Appending to an existing file: EuPathDB/metadata/OrganismDbi_biocv3.13_tritrypdbv54_metadata.csv
## Appending to an existing file: EuPathDB/metadata/BSgenome_biocv3.13_tritrypdbv54_metadata.csv
## Appending to an existing file: EuPathDB/metadata/biocv3.13_tritrypdbv54_invalid_metadata.csv
## Appending to an existing file: EuPathDB/metadata/GRanges_biocv3.13_tritrypdbv54_invalid_metadata.csv
## Appending to an existing file: EuPathDB/metadata/OrgDb_biocv3.13_tritrypdbv54_invalid_metadata.csv
## Appending to an existing file: EuPathDB/metadata/TxDb_biocv3.13_tritrypdbv54_invalid_metadata.csv
## Appending to an existing file: EuPathDB/metadata/OrganismDbi_biocv3.13_tritrypdbv54_invalid_metadata.csv
## Appending to an existing file: EuPathDB/metadata/BSgenome_biocv3.13_tritrypdbv54_invalid_metadata.csv
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.
lp_entry
## AnnotationVersion AnnotationSource BiocVersion DataProvider
## 6 Dec 16, 2014 GeneDB 3.13 TriTrypDB
## Genome GenomeSource GenomeVersion
## 6 TriTrypDB-54_LpanamensisMHOMCOL81L13 GenBank GCA_000340495.1
## NumArrayGene NumChipChipGene NumChromosome NumCodingGene NumCommunity
## 6 <NA> <NA> 36 8665 12
## NumContig NumEC NumEST NumGene NumGO NumOrtholog NumOtherGene NumPopSet
## 6 0 62 <NA> 8778 3930 8665 113 37300
## NumProteomics NumPseudogene NumRNASeq NumRTPCR NumSNP NumTFBS Organellar
## 6 <NA> 0 <NA> NA 5337 NA no
## ReferenceStrain MegaBP PrimaryKey ProjectID
## 6 yes 31.26 Leishmania panamensis MHOM/COL/81/L13 TriTrypDB
## RecordClassName SourceID SourceVersion
## 6 OrganismRecordClasses.OrganismRecordClass NCBITAXON_1055687 54
## TaxonomyID TaxonomyName
## 6 1295824 Leishmania panamensis MHOM/COL/81/L13
## URLGenome
## 6 http://TriTrypDB.org/common/downloads/release-54/LpanamensisMHOMCOL81L13/fasta/data/TriTrypDB-54_LpanamensisMHOMCOL81L13_Genome.fasta
## URLGFF
## 6 http://TriTrypDB.org/common/downloads/release-54/LpanamensisMHOMCOL81L13/gff/data/TriTrypDB-54_LpanamensisMHOMCOL81L13.gff
## URLProtein
## 6 http://TriTrypDB.org/common/downloads/release-54/LpanamensisMHOMCOL81L13/fasta/data/TriTrypDB-54_LpanamensisMHOMCOL81L13_AnnotatedProteins.fasta
## Coordinate_1_based Maintainer
## 6 TRUE Keith Hughitt <khughitt@umd.edu>
## SourceUrl
## 6 http://TriTrypDB.org/common/downloads/release-54/LpanamensisMHOMCOL81L13/gff/data/TriTrypDB-54_LpanamensisMHOMCOL81L13.gff
## Tags
## 6 Annotation:EuPathDB:Eukaryote:Pathogen:Parasite:Trypanosome:Kinetoplastid:Leishmania
## BsgenomePkg
## 6 BSGenome.Leishmania.panamensis.MHOMCOL81L13.v54
## GrangesPkg
## 6 GRanges.Leishmania.panamensis.MHOMCOL81L13.v54.rda
## OrganismdbiPkg
## 6 tritrypdb.Leishmania.panamensis.MHOMCOL81L13.v54
## OrgdbPkg
## 6 org.Lpanamensis.MHOMCOL81L13.v54.eg.db
## TxdbPkg
## 6 TxDb.Leishmania.panamensis.MHOMCOL81L13.TriTrypDB.v54
## Taxon Genus Species Strain
## 6 Leishmania panamensis MHOMCOL81L13 Leishmania panamensis MHOMCOL81L13
## BsgenomeFile
## 6 EuPathDB/BSgenome/3.13/BSGenome.Leishmania.panamensis.MHOMCOL81L13.v54/single_sequences.2bit
## GrangesFile
## 6 EuPathDB/GRanges/3.13/GRanges.Leishmania.panamensis.MHOMCOL81L13.v54.rda
## OrganismdbiFile
## 6 EuPathDB/OrganismDbi/3.13/tritrypdb.Leishmania.panamensis.MHOMCOL81L13.v54/graphInfo.rda
## OrgdbFile
## 6 EuPathDB/OrgDb/3.13/org.Lpanamensis.MHOMCOL81L13.v54.eg.sqlite
## TxdbFile
## 6 EuPathDB/TxDb/3.13/TxDb.Leishmania.panamensis.MHOMCOL81L13.TriTrypDB.v54.sqlite
## GenusSpecies TaxonUnmodified
## 6 Leishmania panamensis Leishmania panamensis MHOM/COL/81/L13
## TaxonCanonical TaxonXref
## 6 Leishmania panamensis MHOM/COL/81/L13 Leishmania panamensis MHOM/COL/81/L13
## testing_panamensis <- EuPathDB::make_eupath_bsgenome(entry=lp_entry)
library(as.character(testing_panamensis), character.only=TRUE)
## Error in library(as.character(testing_panamensis), character.only = TRUE): object 'testing_panamensis' not found
genome <- get0(as.character(testing_panamensis))
## Error in get0(as.character(testing_panamensis)): object 'testing_panamensis' 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’
sanitize_columns <- c("passagenumber", "clinicalresponse", "clinicalcategorical",
"zymodemecategorical", "phenotypiccharacteristics")
lp_expt <- sm(create_expt(sample_sheet,
gene_info = hisat_annot,
id_column = "hpglidentifier",
file_column = "lpanamensisv36hisatfile")) %>%
set_expt_conditions(fact = "zymodemecategorical") %>%
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")
## The samples (and read coverage) removed when filtering 8550 non-zero genes are:
## TMRC20002 TMRC20006
## 11681227 6670348
## subset_expt(): There were 74, now there are 72 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 72, now there are 70 samples.
## semantic_expt_filter(): Removed 68 genes.
libsizes <- plot_libsize(lp_expt)
pp(file = "images/lp_expt_libsizes.png", image = libsizes$plot, width = 14, height = 9)
## I think samples 7,10 should be removed at minimum, probably also 9,11
nonzero <- plot_nonzero(lp_expt)
pp(file = "images/lp_nonzero.png", image = nonzero$plot, width = 9, height = 9)
## Warning: ggrepel: 49 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
## Warning: ggrepel: 53 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
lp_box <- plot_boxplot(lp_expt)
## 5252 entries are 0. We are on a log scale, adding 1 to the data.
pp(file = "images/lp_expt_boxplot.png", image = lp_box, width = 12, height = 9)
filter_plot <- plot_libsize_prepost(lp_expt)
filter_plot$lowgene_plot
## Warning: Using alpha for a discrete variable is not advised.
filter_plot$count_plot
Najib’s favorite plots are of course the PCA/TNSE. These are nice to look at in order to get a sense of the relationships between samples. They also provide a good opportunity to see what happens when one applies different normalizations, surrogate analyses, filters, etc. In addition, one may set different experimental factors as the primary ‘condition’ (usually the color of plots) and surrogate ‘batches’.
Column ‘Q’ in the sample sheet, make a categorical version of it with these parameters:
starting <- as.numeric(pData(lp_expt)[["susceptibilityinfectionreduction32ugmlsbvhistoricaldata"]])
sus_categorical <- starting
na_idx <- is.na(starting)
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"
pData(lp_expt$expressionset)[["sus_category"]] <- sus_categorical
clinical_colors <- list(
"z2.3" = "#874400",
"z2.2" = "#df7000",
"unknown" = "#cbcbcb",
"null" = "#000000")
clinical_samples <- lp_expt %>%
set_expt_batches(fact = sus_categorical) %>%
set_expt_colors(clinical_colors)
clinical_norm <- sm(normalize_expt(clinical_samples, norm = "quant", transform = "log2",
convert = "cpm", batch = FALSE, filter = TRUE))
zymo_pca <- plot_pca(clinical_norm, plot_title = "PCA of parasite expression values",
plot_labels = FALSE)
pp(file = "images/zymo_pca_sus_shape.png", image = zymo_pca$plot)
zymo_3dpca <- plot_3d_pca(zymo_pca)
zymo_3dpca$plot
clinical_n <- sm(normalize_expt(clinical_samples, transform = "log2",
convert = "cpm", batch = FALSE, filter = TRUE))
zymo_tsne <- plot_tsne(clinical_n, plot_title = "TSNE of parasite expression values")
zymo_tsne$plot
clinical_nb <- normalize_expt(clinical_samples, convert = "cpm", transform = "log2",
filter = TRUE, batch = "svaseq")
## Removing 142 low-count genes (8568 remaining).
## batch_counts: Before batch/surrogate estimation, 989 entries are x==0: 0%.
## batch_counts: Before batch/surrogate estimation, 3326 entries are 0<x<1: 1%.
## Setting 335 low elements to zero.
## transform_counts: Found 335 values equal to 0, adding 1 to the matrix.
clinical_nb_pca <- plot_pca(clinical_nb, plot_title = "PCA of parasite expression values",
plot_labels = FALSE)
pp(file = "images/clinical_nb_pca_sus_shape.png", image = clinical_nb_pca$plot)
clinical_nb_tsne <- plot_tsne(clinical_nb, plot_title = "TSNE of parasite expression values")
clinical_nb_tsne$plot
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
## Warning: ggrepel: 70 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
corheat <- plot_corheat(clinical_norm, plot_title = "Correlation heatmap of parasite
expression values
")
corheat$plot
plot_sm(clinical_norm)$plot
## Performing correlation.
cf_colors <- list(
"cure" = "#006f00",
"fail" = "#9dffa0",
"unknown" = "#cbcbcb",
"notapplicable" = "#000000")
cf_expt <- set_expt_conditions(lp_expt, fact = "clinicalcategorical") %>%
set_expt_batches(fact = sus_categorical) %>%
set_expt_colors(cf_colors)
cf_norm <- normalize_expt(cf_expt, convert = "cpm", transform = "log2",
norm = "quant", filter = TRUE)
## Removing 142 low-count genes (8568 remaining).
## transform_counts: Found 2 values equal to 0, adding 1 to the matrix.
start_cf <- plot_pca(cf_norm, plot_title = "PCA of parasite expression values",
plot_labels = FALSE)
pp(file = "images/cf_sus_shape.png", image = start_cf$plot)
cf_nb <- normalize_expt(cf_expt, convert = "cpm", transform = "log2",
norm = "quant", filter = TRUE, batch = "svaseq")
## Warning in normalize_expt(cf_expt, convert = "cpm", transform = "log2", :
## Quantile normalization and sva do not always play well together.
## Removing 142 low-count genes (8568 remaining).
## batch_counts: Before batch/surrogate estimation, 2 entries are x==0: 0%.
## batch_counts: Before batch/surrogate estimation, 4074 entries are 0<x<1: 1%.
## Setting 181 low elements to zero.
## transform_counts: Found 181 values equal to 0, adding 1 to the matrix.
cf_nb_pca <- plot_pca(cf_nb, plot_title = "PCA of parasite expression values",
plot_labels = FALSE)
pp(file = "images/cf_sus_share_nb.png", image = cf_nb_pca$plot)
cf_norm <- normalize_expt(cf_expt, transform = "log2", convert = "cpm",
filter = TRUE, norm = "quant")
## Removing 142 low-count genes (8568 remaining).
## transform_counts: Found 2 values equal to 0, adding 1 to the matrix.
test <- pca_information(cf_norm,
expt_factors = c("clinicalcategorical", "zymodemecategorical",
"pathogenstrain", "passagenumber"),
num_components = 6, plot_pcas = TRUE)
test$anova_p
## PC1 PC2 PC3 PC4 PC5 PC6
## clinicalcategorical 0.2746798 0.0762058 0.3880 2.535e-04 0.34372 8.192e-01
## zymodemecategorical 0.0004562 0.8405440 0.4855 1.353e-01 0.33859 5.403e-02
## pathogenstrain 0.7985483 0.9035654 0.7089 6.322e-06 0.01760 5.689e-01
## passagenumber 0.6040298 0.0001338 0.9479 1.223e-02 0.07891 8.438e-05
test$cor_heatmap
sus_colors <- list(
"resistant" = "#8563a7",
"sensitive" = "#8d0000",
"ambiguous" = "#cbcbcb",
"unknown" = "#000000")
sus_expt <- set_expt_conditions(lp_expt, fact = "sus_category") %>%
set_expt_batches(fact = "zymodemecategorical") %>%
set_expt_colors(colors = sus_colors)
sus_norm <- normalize_expt(sus_expt, transform = "log2", convert = "cpm",
norm = "quant", filter = TRUE)
## Removing 142 low-count genes (8568 remaining).
## transform_counts: Found 2 values equal to 0, adding 1 to the matrix.
sus_pca <- plot_pca(sus_norm, plot_title = "PCA of parasite expression values",
plot_labels = FALSE)
pp(file = "images/sus_norm_pca.png", image = sus_pca[["plot"]])
sus_nb <- normalize_expt(sus_expt, transform = "log2", convert = "cpm",
batch = "svaseq", filter = TRUE)
## Removing 142 low-count genes (8568 remaining).
## batch_counts: Before batch/surrogate estimation, 989 entries are x==0: 0%.
## batch_counts: Before batch/surrogate estimation, 3326 entries are 0<x<1: 1%.
## Setting 222 low elements to zero.
## transform_counts: Found 222 values equal to 0, adding 1 to the matrix.
sus_nb_pca <- plot_pca(sus_nb, plot_title = "PCA of parasite expression values",
plot_labels = FALSE)
pp(file = "images/sus_nb_pca.png", image = sus_nb_pca[["plot"]])
At this time, we do not have very many samples, so the set of metrics/plots is fairly limited. There is really only one factor in the metadata which we can use for performing differential expression analyses, the ‘zymodeme’.
The following sections perform a series of analyses which seek to elucidate differences between the zymodemes 2.2 and 2.3 either through differential expression or variant profiles.
TODO: Do this with and without sva and compare the results.
zy_expt <- subset_expt(lp_expt, subset = "condition=='z2.2'|condition=='z2.3'")
## subset_expt(): There were 70, now there are 34 samples.
zy_norm <- normalize_expt(zy_expt, filter = TRUE, convert = "cpm", norm = "quant")
## Removing 166 low-count genes (8544 remaining).
zy_de_nobatch <- sm(all_pairwise(zy_expt, filter = TRUE, model_batch = "svaseq"))
zy_de <- sm(all_pairwise(zy_expt, filter = TRUE, model_batch = "svaseq"))
zy_table <- sm(combine_de_tables(zy_de, excel = glue::glue("excel/zy_tables-v{ver}.xlsx")))
zy_sig <- sm(extract_significant_genes(zy_table, excel = glue::glue("excel/zy_sig-v{ver}.xlsx")))
pp(file = "images/zymo_ma.png", image = zy_table[["plots"]][["z23_vs_z22"]][["deseq_ma_plots"]][["plot"]])
In contrast, we can search for genes which are differentially expressed with respect to cure/failure status.
cf_de <- sm(all_pairwise(cf_expt, filter = TRUE, model_batch = "svaseq"))
cf_table <- sm(combine_de_tables(cf_de, excel = glue::glue("excel/cf_tables-v{ver}.xlsx")))
cf_sig <- sm(extract_significant_genes(cf_table, excel = glue::glue("excel/cf_sig-v{ver}.xlsx")))
Finally, we can use our category of susceptibility and look for genes which change from sensitive to resistant. Keep in mind, though, that for the moment we have a lot of ambiguous and unknown strains.
sus_de <- sm(all_pairwise(sus_expt, filter = TRUE, model_batch = "svaseq"))
sus_table <- sm(combine_de_tables(sus_de, excel = glue::glue("excel/sus_tables-v{ver}.xlsx")))
sus_sig <- sm(extract_significant_genes(sus_table, excel = glue::glue("excel/sus_sig-v{ver}.xlsx")))
knitr::kable(head(sus_sig$deseq$ups$sensitive_vs_resistant, n = 20))
gid | annotgeneproduct | annotgenetype | chromosome | start | end | strand | annotgeneentrezid | annotgenename | annotstrand | annotchromosome | annotcdslength | length | deseq_logfc | deseq_adjp | edger_logfc | edger_adjp | limma_logfc | limma_adjp | basic_nummed | basic_denmed | basic_numvar | basic_denvar | basic_logfc | basic_t | basic_p | basic_adjp | deseq_basemean | deseq_lfcse | deseq_stat | deseq_p | ebseq_fc | ebseq_logfc | ebseq_c1mean | ebseq_c2mean | ebseq_mean | ebseq_var | ebseq_postfc | ebseq_ppee | ebseq_ppde | ebseq_adjp | edger_logcpm | edger_lr | edger_p | limma_ave | limma_t | limma_b | limma_p | limma_adjp_ihw | deseq_adjp_ihw | edger_adjp_ihw | ebseq_adjp_ihw | basic_adjp_ihw | lfc_meta | lfc_var | lfc_varbymed | p_meta | p_var | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
LPAL13_000044900 | LPAL13_000044900 | actin-related protein 2, putative | protein coding | LPAL13_SCAF000645 | 507 | 1685 | - | reverse | Not Assigned | 1179.0 | 1178 | 28.310 | 0e+00 | 13.450 | 0.0000 | 9.1070 | 0.2437 | 3.9690 | -4.1720 | 15.754 | 0.1283 | 8.141 | 10.160 | 0.0000 | 0.0000 | 858.70 | 1.1810 | 23.960 | 0 | 117202.52 | 16.839 | 0.0000 | 1172.015 | 813.899 | 5.036e+05 | 330.256 | 1.0000 | 0e+00 | 1.0000 | 5.0000 | 56.74 | 0e+00 | 1.4520 | 1.603 | -4.2930 | 0.1136 | 2.445e-01 | 5.558e-123 | 6.257e-11 | 0.000e+00 | 7.112e-08 | 14.150 | 8.959e+00 | 6.329e-01 | 3.787e-02 | 4.302e-03 | ||
LPAL13_000035800 | LPAL13_000035800 | hypothetical protein | protein coding | LPAL13_SCAF000500 | 737 | 1006 | - | reverse | Not Assigned | 270.0 | 269 | 14.870 | 0e+00 | 14.070 | 0.0000 | 10.4000 | 0.2983 | 5.3200 | -3.9370 | 15.228 | 0.3932 | 9.258 | 11.530 | 0.0000 | 0.0000 | 2938.00 | 1.1720 | 12.690 | 0 | 29230.09 | 14.835 | 0.1384 | 4339.122 | 3013.322 | 1.388e+07 | 1218.932 | 0.0000 | 0e+00 | 0.0000 | 6.7700 | 84.28 | 0e+00 | 2.3910 | 1.446 | -4.4250 | 0.1527 | 4.223e-01 | 2.867e-33 | 3.683e-16 | 0.000e+00 | 9.311e-09 | 15.970 | 7.884e+00 | 4.936e-01 | 5.090e-02 | 7.772e-03 | ||
LPAL13_320026300 | LPAL13_320026300 | hypothetical protein, conserved | protein coding | LpaL13_32 | 754268 | 755485 | - | reverse | 32 | 1218.0 | 1217 | 14.060 | 0e+00 | 13.250 | 0.0000 | 9.8700 | 0.3157 | 4.7870 | -4.0020 | 18.217 | 0.6715 | 8.789 | 9.890 | 0.0000 | 0.0000 | 1552.00 | 1.1420 | 12.310 | 0 | 16393.05 | 14.001 | 0.1294 | 2285.903 | 1587.472 | 2.058e+06 | 622.804 | 0.0000 | 0e+00 | 0.0000 | 5.8510 | 73.68 | 0e+00 | 2.1190 | 1.400 | -4.4620 | 0.1661 | 3.124e-01 | 2.373e-31 | 3.939e-14 | 0.000e+00 | 5.456e-08 | 12.340 | 1.164e-01 | 9.438e-03 | 5.537e-02 | 9.196e-03 | ||
LPAL13_000053200 | LPAL13_000053200 | hypothetical protein | protein coding | LPAL13_SCAF000804 | 5037 | 5249 | - | reverse | Not Assigned | 213.0 | 212 | 8.847 | 0e+00 | 10.250 | 0.0000 | 5.8640 | 0.0345 | 1.0040 | -4.1720 | 9.054 | 0.1283 | 5.176 | 8.466 | 0.0000 | 0.0000 | 76.55 | 1.1180 | 7.915 | 0 | 12683.81 | 13.631 | 0.0000 | 126.828 | 88.075 | 7.967e+03 | 36.435 | 1.0000 | 0e+00 | 1.0000 | 1.5600 | 49.10 | 0e+00 | -0.8649 | 2.777 | -2.9970 | 0.0071 | 3.422e-02 | 1.764e-12 | 1.145e-09 | 0.000e+00 | 1.092e-06 | 8.176 | 8.389e-02 | 1.026e-02 | 2.361e-03 | 1.672e-05 | ||
LPAL13_000051300 | LPAL13_000051300 | hypothetical protein, conserved | protein coding | LPAL13_SCAF000772 | 11 | 2344 | + | forward | Not Assigned | 2334.0 | 2333 | 8.358 | 0e+00 | 9.303 | 0.0000 | 3.8240 | 0.2550 | 0.2392 | -3.9740 | 10.240 | 0.5654 | 4.213 | 6.206 | 0.0000 | 0.0000 | 138.30 | 1.2740 | 6.560 | 0 | 1623.37 | 10.665 | 0.0951 | 170.535 | 118.456 | 6.846e+04 | 53.119 | 0.0000 | 0e+00 | 0.0000 | 2.4210 | 29.89 | 0e+00 | -1.0680 | 1.570 | -4.3230 | 0.1211 | 2.551e-01 | 1.216e-08 | 3.080e-06 | 0.000e+00 | 4.597e-05 | 7.011 | 2.071e+00 | 2.954e-01 | 4.037e-02 | 4.888e-03 | ||
LPAL13_000017600 | LPAL13_000017600 | hypothetical protein, conserved | protein coding | LPAL13_SCAF000146 | 359 | 586 | + | forward | Not Assigned | 228.0 | 227 | 6.566 | 0e+00 | 6.547 | 0.0000 | 5.9330 | 0.0833 | 4.4480 | -1.2190 | 4.355 | 2.8119 | 5.666 | 8.643 | 0.0000 | 0.0000 | 628.20 | 0.6955 | 9.441 | 0 | 80.11 | 6.324 | 12.0678 | 967.574 | 675.614 | 3.961e+05 | 63.578 | 0.0000 | 1e+00 | 0.0000 | 4.5480 | 54.14 | 0e+00 | 2.3880 | 2.306 | -3.3600 | 0.0242 | 1.191e-01 | 1.019e-17 | 1.364e-10 | 9.176e-01 | 1.268e-06 | 6.766 | 2.103e+00 | 3.109e-01 | 8.050e-03 | 1.944e-04 | ||
LPAL13_000040700 | LPAL13_000040700 | hypothetical protein, conserved | protein coding | LPAL13_SCAF000598 | 54 | 1067 | + | forward | Not Assigned | 1014.0 | 1013 | 6.558 | 0e+00 | 7.957 | 0.0000 | 3.0190 | 0.0829 | -1.2480 | -4.1720 | 6.491 | 0.1283 | 2.923 | 5.613 | 0.0000 | 0.0002 | 20.55 | 1.1800 | 5.557 | 0 | 2601.36 | 11.345 | 0.0000 | 26.004 | 18.058 | 4.051e+02 | 7.897 | 1.0000 | 0e+00 | 1.0000 | -0.1262 | 28.95 | 0e+00 | -2.3590 | 2.309 | -3.6060 | 0.0240 | 8.201e-02 | 1.917e-06 | 4.657e-06 | 0.000e+00 | 2.261e-04 | 5.728 | 1.196e+00 | 2.089e-01 | 8.003e-03 | 1.922e-04 | ||
LPAL13_300029400 | LPAL13_300029400 | hypothetical protein, conserved | protein coding | LpaL13_30 | 853953 | 854150 | - | reverse | 30 | 198.0 | 197 | 6.169 | 0e+00 | 6.102 | 0.0000 | 4.8680 | 0.0042 | 1.7290 | -2.3910 | 1.488 | 1.7754 | 4.120 | 8.766 | 0.0000 | 0.0000 | 90.87 | 0.7525 | 8.198 | 0 | 59.50 | 5.895 | 2.0963 | 125.320 | 87.668 | 9.233e+03 | 22.842 | 0.0000 | 1e+00 | 0.0000 | 1.7520 | 47.86 | 0e+00 | 0.0911 | 3.752 | -1.0820 | 0.0004 | 4.212e-03 | 3.350e-13 | 1.785e-09 | 8.838e-01 | 6.533e-06 | 5.860 | 4.150e-01 | 7.082e-02 | 1.221e-04 | 4.470e-08 | ||
LPAL13_000011700 | LPAL13_000011700 | hypothetical protein | protein coding | LPAL13_SCAF000076 | 101 | 364 | - | reverse | Not Assigned | 264.0 | 263 | 5.905 | 1e-04 | 7.455 | 0.0000 | 2.6870 | 0.0504 | -1.3930 | -4.1720 | 6.742 | 0.1283 | 2.779 | 5.239 | 0.0000 | 0.0004 | 14.11 | 1.2890 | 4.581 | 0 | 2444.21 | 11.255 | 0.0000 | 24.432 | 16.967 | 4.690e+02 | 7.590 | 1.0000 | 0e+00 | 1.0000 | -0.5868 | 23.32 | 0e+00 | -2.9030 | 2.580 | -3.2820 | 0.0121 | 5.072e-02 | 1.116e-04 | 4.581e-05 | 0.000e+00 | 4.319e-04 | 5.160 | 1.655e+00 | 3.208e-01 | 4.025e-03 | 4.854e-05 | ||
LPAL13_080010600 | LPAL13_080010600 | hypothetical protein, conserved | protein coding | LpaL13_08 | 195555 | 195749 | - | reverse | 8 | 195.0 | 194 | 5.774 | 0e+00 | 7.291 | 0.0000 | 2.3560 | 0.0586 | -1.8800 | -4.1720 | 4.901 | 0.1283 | 2.292 | 5.029 | 0.0000 | 0.0005 | 11.19 | 1.1660 | 4.954 | 0 | 1847.09 | 10.851 | 0.0000 | 18.461 | 12.820 | 4.800e+02 | 6.130 | 1.0000 | 0e+00 | 1.0000 | -0.9234 | 24.49 | 0e+00 | -3.1040 | 2.502 | -3.4150 | 0.0148 | 5.863e-02 | 3.267e-05 | 3.031e-05 | 0.000e+00 | 6.681e-04 | 4.688 | 2.810e+00 | 5.995e-01 | 4.920e-03 | 7.261e-05 | ||
LPAL13_040019400 | LPAL13_040019400 | hypothetical protein | protein coding | LpaL13_04 | 440768 | 441127 | - | reverse | 4 | 360.0 | 359 | 5.665 | 0e+00 | 5.547 | 0.0000 | 3.4780 | 0.0263 | -0.4351 | -3.4850 | 1.758 | 1.0862 | 3.050 | 7.418 | 0.0000 | 0.0000 | 34.01 | 0.9723 | 5.827 | 0 | 48.41 | 5.597 | 0.6939 | 34.070 | 23.872 | 1.800e+03 | 8.761 | 0.0000 | 0e+00 | 0.0000 | 0.3889 | 28.41 | 0e+00 | -1.6670 | 2.915 | -2.8310 | 0.0048 | 3.758e-02 | 7.294e-07 | 6.064e-06 | 0.000e+00 | 8.857e-06 | 4.958 | 3.756e-02 | 7.575e-03 | 1.605e-03 | 7.728e-06 | ||
LPAL13_200050100 | LPAL13_200050100 | hypothetical protein | protein coding | LpaL13_20.1 | 1627529 | 1627717 | + | forward | 20.1 | 189.0 | 188 | 5.513 | 0e+00 | 5.463 | 0.0000 | 4.7340 | 0.0024 | 2.4580 | -1.9700 | 1.008 | 2.4447 | 4.428 | 8.642 | 0.0000 | 0.0000 | 118.50 | 0.6171 | 8.934 | 0 | 26.30 | 4.717 | 7.3797 | 194.358 | 137.226 | 2.167e+04 | 18.302 | 0.0000 | 1e+00 | 0.0000 | 2.1670 | 54.24 | 0e+00 | 0.8180 | 4.011 | -0.3392 | 0.0002 | 2.416e-03 | 7.046e-16 | 1.364e-10 | 8.708e-01 | 3.310e-05 | 5.220 | 8.145e-01 | 1.560e-01 | 5.103e-05 | 7.813e-09 | ||
LPAL13_350011800 | LPAL13_350011800 | hypothetical protein, conserved | protein coding | LpaL13_35 | 171009 | 171242 | + | forward | 35 | 234.0 | 233 | 5.141 | 0e+00 | 5.126 | 0.0000 | 4.3400 | 0.0062 | 2.8760 | -0.8627 | 2.424 | 0.1736 | 3.739 | 11.140 | 0.0000 | 0.0000 | 178.20 | 0.5811 | 8.847 | 0 | 31.92 | 4.997 | 9.4333 | 301.448 | 212.221 | 5.943e+04 | 24.086 | 0.0000 | 1e+00 | 0.0000 | 2.7240 | 54.09 | 0e+00 | 1.2480 | 3.584 | -0.8512 | 0.0006 | 6.229e-03 | 1.282e-15 | 1.364e-10 | 9.176e-01 | 4.879e-09 | 4.952 | 7.332e-01 | 1.481e-01 | 2.108e-04 | 1.333e-07 | ||
LPAL13_080010800 | LPAL13_080010800 | hypothetical protein | protein coding | LpaL13_08 | 199409 | 199792 | - | reverse | 8 | 384.0 | 383 | 5.012 | 1e-04 | 6.449 | 0.0001 | 1.6610 | 0.2389 | -2.3550 | -4.1720 | 4.137 | 0.1283 | 1.817 | 4.316 | 0.0002 | 0.0023 | 10.46 | 1.0880 | 4.608 | 0 | 1048.69 | 10.034 | 0.0000 | 10.477 | 7.276 | 1.077e+02 | 3.599 | 1.0000 | 0e+00 | 1.0000 | -0.8685 | 22.29 | 0e+00 | -3.0780 | 1.618 | -4.2920 | 0.1102 | 2.393e-01 | 1.006e-04 | 6.986e-05 | 0.000e+00 | 2.524e-03 | 3.942 | 3.004e+00 | 7.622e-01 | 3.674e-02 | 4.048e-03 | ||
LPAL13_170014500 | LPAL13_170014500 | hypothetical protein, conserved | protein coding | LpaL13_17 | 361708 | 362040 | + | forward | 17 | 333.0 | 332 | 4.932 | 1e-04 | 4.711 | 0.0012 | 2.8030 | 0.0333 | -0.6524 | -3.1970 | 6.988 | 1.5170 | 2.544 | 3.938 | 0.0004 | 0.0038 | 21.46 | 1.0430 | 4.730 | 0 | 43.02 | 5.427 | 1.0218 | 44.379 | 31.131 | 1.643e+03 | 10.389 | 0.0000 | 0e+00 | 0.0000 | -0.2584 | 15.30 | 1e-04 | -2.3840 | 2.796 | -2.9090 | 0.0067 | 3.341e-02 | 6.295e-05 | 1.234e-03 | 0.000e+00 | 3.646e-03 | 4.132 | 5.989e-02 | 1.450e-02 | 2.273e-03 | 1.487e-05 | ||
LPAL13_000011800 | LPAL13_000011800 | hypothetical protein, conserved | protein coding | LPAL13_SCAF000076 | 446 | 640 | - | reverse | Not Assigned | 195.0 | 194 | 4.806 | 2e-04 | 5.382 | 0.0004 | 0.9836 | 0.4742 | -2.4970 | -3.9580 | 3.651 | 0.6042 | 1.461 | 3.260 | 0.0025 | 0.0153 | 11.25 | 1.0820 | 4.441 | 0 | 56.34 | 5.816 | 0.1523 | 9.134 | 6.389 | 7.337e+01 | 2.967 | 0.9995 | 5e-04 | 0.9995 | -0.8436 | 17.82 | 0e+00 | -3.0280 | 1.041 | -4.6760 | 0.3014 | 6.606e-01 | 1.913e-04 | 4.376e-04 | 1.114e-02 | 1.713e-02 | 3.267 | 3.339e+00 | 1.022e+00 | 1.005e-01 | 3.028e-02 | ||
LPAL13_000014000 | LPAL13_000014000 | hypothetical protein | protein coding | LPAL13_SCAF000119 | 655 | 942 | + | forward | Not Assigned | 288.0 | 287 | 4.413 | 0e+00 | 4.400 | 0.0000 | 4.0190 | 0.0151 | 2.4340 | -1.2390 | 1.518 | 1.8989 | 3.673 | 7.604 | 0.0000 | 0.0000 | 130.40 | 0.5404 | 8.166 | 0 | 19.75 | 4.304 | 9.7345 | 192.485 | 136.644 | 1.297e+04 | 14.826 | 0.0000 | 1e+00 | 0.0000 | 2.3000 | 49.02 | 0e+00 | 1.1140 | 3.185 | -1.7980 | 0.0022 | 1.518e-02 | 3.350e-13 | 1.145e-09 | 9.176e-01 | 3.693e-05 | 4.411 | 8.730e-01 | 1.979e-01 | 7.287e-04 | 1.593e-06 | ||
LPAL13_000026500 | LPAL13_000026500 | hypothetical protein | protein coding | LPAL13_SCAF000301 | 144 | 494 | - | reverse | Not Assigned | 351.0 | 350 | 4.293 | 0e+00 | 4.241 | 0.0001 | 2.5410 | 0.1555 | 0.1369 | -2.4930 | 5.421 | 2.0348 | 2.630 | 4.149 | 0.0003 | 0.0028 | 44.64 | 0.8060 | 5.326 | 0 | 21.34 | 4.415 | 2.5893 | 55.452 | 39.299 | 1.544e+03 | 9.519 | 0.0000 | 0e+00 | 0.0000 | 0.9069 | 22.12 | 0e+00 | -0.8699 | 1.921 | -3.9520 | 0.0590 | 1.557e-01 | 5.506e-06 | 7.424e-05 | 0.000e+00 | 2.670e-03 | 3.646 | 2.689e-02 | 7.375e-03 | 1.966e-02 | 1.159e-03 | ||
LPAL13_000035500 | LPAL13_000035500 | hypothetical protein, conserved | protein coding | LPAL13_SCAF000492 | 7045 | 7410 | + | forward | Not Assigned | 366.0 | 365 | 4.211 | 0e+00 | 4.212 | 0.0000 | 3.7780 | 0.0333 | 4.5240 | 0.8413 | 2.356 | 0.6018 | 3.683 | 9.543 | 0.0000 | 0.0000 | 518.90 | 0.6035 | 6.977 | 0 | 18.47 | 4.207 | 48.2560 | 891.445 | 633.804 | 3.430e+05 | 17.782 | 0.0000 | 1e+00 | 0.0000 | 4.2780 | 35.09 | 0e+00 | 2.8460 | 2.797 | -2.3710 | 0.0067 | 3.335e-02 | 9.969e-10 | 3.974e-07 | 9.176e-01 | 2.696e-08 | 4.168 | 7.708e-01 | 1.849e-01 | 2.234e-03 | 1.497e-05 | ||
LPAL13_220019500 | LPAL13_220019500 | hypothetical protein | protein coding | LpaL13_22 | 578260 | 578538 | + | forward | 22 | 279.0 | 278 | 3.915 | 0e+00 | 3.913 | 0.0000 | 3.0880 | 0.0423 | 3.5420 | 0.2509 | 2.124 | 0.7532 | 3.291 | 8.402 | 0.0000 | 0.0000 | 289.60 | 0.5056 | 7.744 | 0 | 15.50 | 3.955 | 27.6226 | 428.407 | 305.945 | 7.136e+04 | 14.131 | 0.0000 | 1e+00 | 0.0000 | 3.4440 | 45.27 | 0e+00 | 2.3790 | 2.676 | -2.6320 | 0.0093 | 4.198e-02 | 7.595e-12 | 5.076e-09 | 8.708e-01 | 3.683e-07 | 3.826 | 4.090e-01 | 1.069e-01 | 3.112e-03 | 2.906e-05 |
knitr::kable(head(sus_sig$deseq$downs$sensitive_vs_resistant, n = 20))
gid | annotgeneproduct | annotgenetype | chromosome | start | end | strand | annotgeneentrezid | annotgenename | annotstrand | annotchromosome | annotcdslength | length | deseq_logfc | deseq_adjp | edger_logfc | edger_adjp | limma_logfc | limma_adjp | basic_nummed | basic_denmed | basic_numvar | basic_denvar | basic_logfc | basic_t | basic_p | basic_adjp | deseq_basemean | deseq_lfcse | deseq_stat | deseq_p | ebseq_fc | ebseq_logfc | ebseq_c1mean | ebseq_c2mean | ebseq_mean | ebseq_var | ebseq_postfc | ebseq_ppee | ebseq_ppde | ebseq_adjp | edger_logcpm | edger_lr | edger_p | limma_ave | limma_t | limma_b | limma_p | limma_adjp_ihw | deseq_adjp_ihw | edger_adjp_ihw | ebseq_adjp_ihw | basic_adjp_ihw | lfc_meta | lfc_var | lfc_varbymed | p_meta | p_var | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
LPAL13_000033300 | LPAL13_000033300 | hypothetical protein, conserved | protein coding | LPAL13_SCAF000463 | 551 | 811 | + | forward | Not Assigned | 261.0 | 260 | -5.366 | 0.0006 | -5.252 | 0.0016 | -6.005 | 0.0000 | -3.7290 | 3.4740 | 10.5855 | 0.0633 | -7.202 | -10.990 | 0 | 0e+00 | 127.700 | 1.3070 | -4.106 | 0.0000 | 0.1238 | -3.013 | 311.93 | 38.622 | 122.132 | 2.245e+04 | 0.1327 | 0.0000 | 0.0000 | 0.0000 | 2.2460 | 14.750 | 0.0001 | -1.1070 | -5.720 | 5.0450 | 0e+00 | 3.109e-05 | 7.620e-04 | 1.557e-03 | 0.000e+00 | 2.310e-08 | -5.541 | 0.000e+00 | 0.000e+00 | 5.432e-05 | 3.885e-09 | ||
LPAL13_000038400 | LPAL13_000038400 | expression-site associated gene (esag3), putative | protein coding | LPAL13_SCAF000573 | 101 | 1360 | + | forward | Not Assigned | 1260.0 | 1259 | -2.987 | 0.0000 | -2.981 | 0.0000 | -3.311 | 0.0001 | 4.6320 | 8.2380 | 3.1197 | 0.0291 | -3.606 | -10.100 | 0 | 0e+00 | 3589.000 | 0.5604 | -5.329 | 0.0000 | 0.1737 | -2.526 | 8823.76 | 1532.438 | 3760.342 | 1.554e+07 | 0.1785 | 0.0000 | 0.0000 | 0.0000 | 7.0530 | 31.760 | 0.0000 | 5.7730 | -5.315 | 5.1220 | 0e+00 | 8.555e-05 | 5.448e-06 | 1.508e-06 | 0.000e+00 | 7.112e-08 | -3.046 | 2.268e-02 | -7.444e-03 | 4.680e-07 | 5.059e-13 | ||
LPAL13_350063000 | LPAL13_350063000 | hypothetical protein | protein coding | LpaL13_35 | 1964328 | 1964543 | - | reverse | 35 | 216.0 | 215 | -2.831 | 0.0000 | -2.804 | 0.0000 | -3.432 | 0.0000 | -2.3300 | 1.2080 | 2.1382 | 0.2275 | -3.538 | -10.860 | 0 | 0e+00 | 20.770 | 0.4869 | -5.813 | 0.0000 | 0.1384 | -2.853 | 55.15 | 7.623 | 22.146 | 6.307e+02 | 0.1595 | 0.0000 | 1.0000 | 0.0000 | -0.3699 | 32.430 | 0.0000 | -1.4710 | -6.977 | 8.0100 | 0e+00 | 1.358e-06 | 6.482e-07 | 1.165e-06 | 8.838e-01 | 4.879e-09 | -3.042 | 6.226e-03 | -2.047e-03 | 6.694e-09 | 2.932e-17 | ||
LPAL13_140019300 | LPAL13_140019300 | bt1 family, putative | protein coding | LpaL13_14 | 530784 | 531350 | + | forward | 14 | 567.0 | 566 | -2.670 | 0.0000 | -2.666 | 0.0000 | -2.508 | 0.0000 | 4.6450 | 7.0660 | 0.4797 | 1.1123 | -2.420 | -6.978 | 0 | 2e-04 | 1818.000 | 0.3821 | -6.989 | 0.0000 | 0.1702 | -2.555 | 4611.25 | 784.787 | 1953.984 | 5.164e+06 | 0.1758 | 0.0000 | 1.0000 | 0.0000 | 6.0720 | 56.420 | 0.0000 | 5.3600 | -6.856 | 11.0800 | 0e+00 | 1.496e-06 | 9.490e-10 | 6.257e-11 | 1.000e+00 | 1.499e-04 | -2.652 | 1.598e-01 | -6.026e-02 | 8.713e-10 | 2.270e-18 | ||
LPAL13_000012000 | LPAL13_000012000 | hypothetical protein | protein coding | LPAL13_SCAF000080 | 710 | 1159 | - | reverse | Not Assigned | 450.0 | 449 | -2.643 | 0.0007 | -2.634 | 0.0003 | -3.063 | 0.0036 | 0.0950 | 3.9560 | 7.6621 | 0.1804 | -3.861 | -6.796 | 0 | 0e+00 | 203.300 | 0.6482 | -4.077 | 0.0000 | 0.2237 | -2.160 | 451.41 | 100.978 | 208.054 | 4.893e+04 | 0.2351 | 0.1874 | 0.8126 | 0.1874 | 2.9190 | 18.660 | 0.0000 | 1.3190 | -3.821 | 0.2633 | 3e-04 | 3.601e-03 | 8.461e-04 | 3.059e-04 | 8.579e-01 | 1.862e-05 | -2.788 | 1.668e-02 | -5.981e-03 | 1.174e-04 | 2.285e-08 | ||
LPAL13_310035500 | LPAL13_310035500 | hypothetical protein | protein coding | LpaL13_31 | 1198439 | 1198957 | - | reverse | 31 | 519.0 | 518 | -2.628 | 0.0033 | -2.493 | 0.0023 | -3.214 | 0.0000 | -4.1830 | -0.4027 | 4.0925 | 0.4977 | -3.780 | -8.269 | 0 | 0e+00 | 7.018 | 0.7348 | -3.576 | 0.0003 | 0.2939 | -1.767 | 18.43 | 5.409 | 9.386 | 3.219e+02 | 0.3477 | 0.0000 | 0.0000 | 0.0000 | -1.9350 | 13.780 | 0.0002 | -3.1660 | -7.180 | 6.1600 | 0e+00 | 7.301e-07 | 3.277e-03 | 2.291e-03 | 0.000e+00 | 3.375e-07 | -2.778 | 0.000e+00 | 0.000e+00 | 1.846e-04 | 3.067e-08 | ||
LPAL13_310039200 | LPAL13_310039200 | hypothetical protein | protein coding | LpaL13_31 | 1301745 | 1301972 | - | reverse | 31 | 228.0 | 227 | -2.429 | 0.0000 | -2.425 | 0.0000 | -2.486 | 0.0000 | 1.2200 | 3.7790 | 1.3494 | 0.2183 | -2.559 | -9.417 | 0 | 0e+00 | 188.800 | 0.4143 | -5.864 | 0.0000 | 0.2457 | -2.025 | 396.35 | 97.358 | 188.717 | 3.434e+04 | 0.2570 | 0.4442 | 0.5558 | 0.4442 | 2.8160 | 38.010 | 0.0000 | 1.9640 | -5.821 | 6.9370 | 0e+00 | 3.275e-05 | 5.090e-07 | 1.181e-07 | 5.203e-01 | 2.310e-08 | -2.476 | 8.014e-02 | -3.237e-02 | 6.034e-08 | 1.000e-14 | ||
LPAL13_000012100 | LPAL13_000012100 | hypothetical protein | protein coding | LPAL13_SCAF000080 | 1637 | 1894 | - | reverse | Not Assigned | 258.0 | 257 | -2.286 | 0.0082 | -2.272 | 0.0075 | -3.355 | 0.0003 | -2.2140 | 1.2000 | 6.2602 | 0.6808 | -3.414 | -6.109 | 0 | 0e+00 | 30.200 | 0.7042 | -3.246 | 0.0012 | 0.2746 | -1.865 | 66.01 | 18.121 | 32.755 | 1.804e+03 | 0.3055 | 0.0537 | 0.9463 | 0.0537 | 0.2122 | 10.880 | 0.0010 | -1.4340 | -4.809 | 2.2850 | 0e+00 | 3.152e-04 | 8.203e-03 | 7.497e-03 | 8.335e-01 | 4.173e-05 | -2.689 | 7.115e-02 | -2.646e-02 | 7.163e-04 | 3.850e-07 | ||
LPAL13_340039600 | LPAL13_340039600 | hypothetical protein | protein coding | LpaL13_34 | 1247554 | 1247757 | - | reverse | 34 | 204.0 | 203 | -2.241 | 0.0002 | -2.237 | 0.0001 | -2.628 | 0.0015 | 1.2500 | 4.2050 | 3.6080 | 0.0559 | -2.955 | -7.645 | 0 | 0e+00 | 218.200 | 0.5128 | -4.369 | 0.0000 | 0.2254 | -2.149 | 518.32 | 116.821 | 239.500 | 4.536e+04 | 0.2307 | 0.0000 | 1.0000 | 0.0000 | 3.0090 | 21.460 | 0.0000 | 1.9570 | -4.190 | 1.3890 | 1e-04 | 1.508e-03 | 2.494e-04 | 9.710e-05 | 9.635e-01 | 4.535e-06 | -2.376 | 8.498e-03 | -3.576e-03 | 3.288e-05 | 1.872e-09 | ||
LPAL13_310031000 | LPAL13_310031000 | hypothetical protein, conserved | protein coding | LpaL13_31 | 1075172 | 1075459 | - | reverse | 31 | 288.0 | 287 | -2.221 | 0.0000 | -2.201 | 0.0000 | -2.810 | 0.0000 | -2.0060 | 0.9976 | 3.5024 | 0.5671 | -3.003 | -6.860 | 0 | 0e+00 | 26.100 | 0.4570 | -4.860 | 0.0000 | 0.2698 | -1.890 | 55.87 | 15.064 | 27.532 | 1.035e+03 | 0.2951 | 0.1254 | 0.8746 | 0.1254 | 0.0381 | 23.870 | 0.0000 | -1.1800 | -5.973 | 5.4530 | 0e+00 | 1.497e-05 | 3.843e-05 | 3.755e-05 | 8.682e-01 | 6.085e-06 | -2.400 | 6.650e-04 | -2.772e-04 | 7.662e-07 | 3.423e-13 | ||
LPAL13_000038500 | LPAL13_000038500 | hypothetical protein | protein coding | LPAL13_SCAF000575 | 39 | 251 | + | forward | Not Assigned | 213.0 | 212 | -2.098 | 0.0050 | -2.083 | 0.0077 | -3.311 | 0.0001 | -1.9360 | 1.4300 | 4.6504 | 0.6757 | -3.367 | -6.768 | 0 | 0e+00 | 31.310 | 0.6124 | -3.427 | 0.0006 | 0.2833 | -1.820 | 77.47 | 21.939 | 38.906 | 2.414e+03 | 0.3096 | 0.1593 | 0.8407 | 0.1593 | 0.2145 | 10.810 | 0.0010 | -1.3130 | -5.467 | 4.1160 | 0e+00 | 5.635e-05 | 4.994e-03 | 7.747e-03 | 8.406e-01 | 6.533e-06 | -2.471 | 1.211e-01 | -4.903e-02 | 5.407e-04 | 2.584e-07 | ||
LPAL13_310031300 | LPAL13_310031300 | hypothetical protein, conserved | protein coding | LpaL13_31 | 1084772 | 1085059 | - | reverse | 31 | 288.0 | 287 | -2.033 | 0.0047 | -2.029 | 0.0032 | -3.046 | 0.0004 | -1.0710 | 2.0840 | 3.9543 | 0.7623 | -3.155 | -6.616 | 0 | 0e+00 | 61.360 | 0.5900 | -3.446 | 0.0006 | 0.2465 | -2.020 | 130.45 | 32.149 | 62.187 | 5.901e+03 | 0.2666 | 0.0529 | 0.9471 | 0.0529 | 1.1930 | 12.960 | 0.0003 | -0.2414 | -4.723 | 2.5090 | 0e+00 | 5.553e-04 | 4.724e-03 | 3.223e-03 | 8.478e-01 | 8.857e-06 | -2.305 | 3.960e-02 | -1.718e-02 | 3.001e-04 | 7.796e-08 | ||
LPAL13_050005000 | LPAL13_050005000 | hypothetical protein | protein coding | LpaL13_05 | 3394 | 3612 | - | reverse | 5 | 219.0 | 218 | -1.983 | 0.0077 | -1.977 | 0.0046 | -2.732 | 0.0002 | 0.2033 | 2.6750 | 2.0316 | 0.1724 | -2.472 | -7.940 | 0 | 0e+00 | 88.260 | 0.6063 | -3.271 | 0.0011 | 0.2820 | -1.826 | 177.99 | 50.188 | 89.238 | 6.186e+03 | 0.2912 | 0.0006 | 0.9994 | 0.0006 | 1.7080 | 12.080 | 0.0005 | 0.4586 | -5.053 | 3.9660 | 0e+00 | 2.371e-04 | 9.310e-03 | 4.591e-03 | 8.708e-01 | 6.620e-07 | -2.227 | 1.857e-02 | -8.335e-03 | 5.284e-04 | 2.851e-07 | ||
LPAL13_140019100 | LPAL13_140019100 | bt1 family, putative | protein coding | LpaL13_14 | 525164 | 525514 | + | forward | 14 | 351.0 | 350 | -1.975 | 0.0000 | -1.972 | 0.0000 | -2.055 | 0.0000 | 3.9170 | 5.9980 | 0.3507 | 0.5490 | -2.081 | -8.232 | 0 | 0e+00 | 857.500 | 0.3048 | -6.479 | 0.0000 | 0.2333 | -2.100 | 1937.62 | 451.966 | 905.916 | 6.953e+05 | 0.2374 | 0.0000 | 1.0000 | 0.0000 | 4.9890 | 50.690 | 0.0000 | 4.5600 | -7.326 | 12.9600 | 0e+00 | 5.315e-07 | 1.926e-08 | 5.922e-10 | 9.176e-01 | 2.350e-05 | -1.982 | 9.112e-03 | -4.598e-03 | 1.551e-10 | 3.741e-20 | ||
LPAL13_340039700 | LPAL13_340039700 | snare domain containing protein, putative | protein coding | LpaL13_34 | 1248192 | 1248947 | - | reverse | 34 | 756.0 | 755 | -1.781 | 0.0000 | -1.778 | 0.0000 | -1.882 | 0.0000 | 4.6130 | 6.6510 | 0.6868 | 0.0521 | -2.038 | -11.360 | 0 | 0e+00 | 1351.000 | 0.3095 | -5.753 | 0.0000 | 0.2859 | -1.806 | 2810.80 | 803.693 | 1416.976 | 1.132e+06 | 0.2899 | 0.0000 | 1.0000 | 0.0000 | 5.6450 | 38.840 | 0.0000 | 5.2460 | -6.273 | 8.7710 | 0e+00 | 7.177e-06 | 8.329e-07 | 8.056e-08 | 9.573e-01 | 4.879e-09 | -1.773 | 2.901e-02 | -1.636e-02 | 1.256e-08 | 2.072e-16 | ||
LPAL13_170015400 | LPAL13_170015400 | hypothetical protein, conserved | protein coding | LpaL13_17 | 395975 | 396307 | + | forward | 17 | 333.0 | 332 | -1.674 | 0.0000 | -1.669 | 0.0000 | -1.669 | 0.0002 | 1.2440 | 3.2670 | 0.9459 | 0.1363 | -2.023 | -9.027 | 0 | 0e+00 | 148.400 | 0.2999 | -5.584 | 0.0000 | 0.3146 | -1.669 | 270.10 | 84.958 | 141.528 | 1.297e+04 | 0.3199 | 0.0000 | 1.0000 | 0.0000 | 2.4610 | 31.800 | 0.0000 | 1.9810 | -5.001 | 4.0120 | 0e+00 | 1.920e-04 | 1.774e-06 | 1.495e-06 | 8.708e-01 | 6.699e-08 | -1.685 | 3.478e-02 | -2.064e-02 | 1.442e-06 | 6.065e-12 | ||
LPAL13_140019200 | LPAL13_140019200 | inositol-3-phosphate synthase | protein coding | LpaL13_14 | 527711 | 529291 | + | INO1 | forward | 14 | 1581.0 | 1580 | -1.465 | 0.0000 | -1.461 | 0.0000 | -1.529 | 0.0000 | 8.8250 | 10.3900 | 0.1718 | 0.4796 | -1.566 | -6.969 | 0 | 2e-04 | 18910.000 | 0.2607 | -5.619 | 0.0000 | 0.3510 | -1.510 | 36839.52 | 12931.847 | 20236.970 | 1.798e+08 | 0.3540 | 0.0000 | 1.0000 | 0.0000 | 9.4520 | 45.990 | 0.0000 | 9.2330 | -6.630 | 10.2000 | 0e+00 | 3.648e-06 | 1.506e-06 | 3.767e-09 | 1.000e+00 | 2.363e-04 | -1.541 | 2.454e-02 | -1.592e-02 | 8.612e-09 | 9.489e-17 | |
LPAL13_320038700 | LPAL13_320038700 | hypothetical protein, conserved | protein coding | LpaL13_32 | 1175024 | 1175257 | + | forward | 32 | 234.0 | 233 | -1.426 | 0.0000 | -1.421 | 0.0000 | -1.458 | 0.0000 | 2.5500 | 3.9570 | 0.4259 | 0.1120 | -1.407 | -8.529 | 0 | 0e+00 | 255.900 | 0.2296 | -6.211 | 0.0000 | 0.4144 | -1.271 | 436.05 | 180.702 | 258.724 | 2.162e+04 | 0.4183 | 0.0011 | 0.9989 | 0.0011 | 3.2490 | 40.400 | 0.0000 | 3.0040 | -6.105 | 8.1100 | 0e+00 | 1.204e-05 | 8.363e-08 | 4.223e-08 | 8.838e-01 | 2.066e-07 | -1.482 | 2.663e-02 | -1.797e-02 | 1.899e-08 | 1.040e-15 | ||
LPAL13_040007800 | LPAL13_040007800 | hypothetical protein, conserved | protein coding | LpaL13_04 | 77524 | 78306 | + | forward | 4 | 783.0 | 782 | -1.387 | 0.0000 | -1.383 | 0.0000 | -1.273 | 0.0000 | 6.5520 | 7.7820 | 0.1432 | 0.3791 | -1.230 | -6.135 | 0 | 6e-04 | 4045.000 | 0.2359 | -5.881 | 0.0000 | 0.3806 | -1.394 | 7003.17 | 2665.491 | 3990.892 | 7.653e+06 | 0.3857 | 0.0001 | 0.9999 | 0.0001 | 7.2260 | 43.070 | 0.0000 | 7.0020 | -5.799 | 6.9070 | 0e+00 | 2.572e-05 | 4.795e-07 | 1.328e-08 | 1.000e+00 | 6.587e-04 | -1.374 | 6.214e-02 | -4.523e-02 | 6.531e-08 | 1.200e-14 | ||
LPAL13_350013200 | LPAL13_350013200 | hypothetical protein, conserved | protein coding | LpaL13_35 | 223837 | 224070 | + | forward | 35 | 234.0 | 233 | -1.377 | 0.0114 | -1.363 | 0.0201 | -2.035 | 0.0003 | -2.0610 | -0.1712 | 1.3327 | 0.7337 | -1.890 | -5.455 | 0 | 3e-04 | 10.700 | 0.4406 | -3.126 | 0.0018 | 0.3587 | -1.479 | 24.88 | 8.918 | 13.797 | 2.950e+02 | 0.3915 | 0.7783 | 0.2217 | 0.7783 | -1.1510 | 8.503 | 0.0035 | -2.0600 | -4.819 | 1.7080 | 0e+00 | 3.139e-04 | 1.141e-02 | 2.008e-02 | 2.557e-01 | 2.467e-04 | -1.577 | 3.478e-02 | -2.205e-02 | 1.776e-03 | 3.127e-06 |
sus_ma <- sus_table[["plots"]][["sensitive_vs_resistant"]][["deseq_ma_plots"]][["plot"]]
pp(file = "images/sus_ma.png", image = sus_ma)
## test <- ggplt(sus_ma)
Now let us look for ontology categories which are increased in the 2.3 samples followed by the 2.2 samples.
## Gene categories more represented in the 2.3 group.
zy_go_up <- sm(simple_goseq(sig_genes = zy_sig[["deseq"]][["ups"]][[1]],
go_db = lp_go, length_db = lp_lengths))
## Gene categories more represented in the 2.2 group.
zy_go_down <- sm(simple_goseq(sig_genes = zy_sig[["deseq"]][["downs"]][[1]],
go_db = lp_go, length_db = lp_lengths))
In the function ‘combined_de_tables()’ above, one of the tasks performed is to look at the agreement among DESeq2, limma, and edgeR. The following show a couple of these for the set of genes observed with a fold-change >= |2| and adjusted p-value <= 0.05.
zy_table[["venns"]][[1]][["p_lfc1"]][["up_noweight"]]
zy_table[["venns"]][[1]][["p_lfc1"]][["down_noweight"]]
zy_go_up$pvalue_plots$bpp_plot_over
zy_go_down$pvalue_plots$bpp_plot_over
Remind myself, the data structures are (zy|sus)_(de|table|sig).
zy_df <- zy_table[["data"]][["z23_vs_z22"]]
sus_df <- sus_table[["data"]][["sensitive_vs_resistant"]]
both_df <- merge(zy_df, sus_df, by = "row.names")
plot_df <- both_df[, c("deseq_logfc.x", "deseq_logfc.y")]
rownames(plot_df) <- both_df[["Row.names"]]
colnames(plot_df) <- c("z23_vs_z22", "sensitive_vs_resistant")
compare <- plot_linear_scatter(plot_df)
## Warning in plot_multihistogram(df): NAs introduced by coercion
pp(file = "images/compare_sus_zy.png", image = compare$scatter)
compare$cor
##
## Pearson's product-moment correlation
##
## data: df[, 1] and df[, 2]
## t = -188, df = 8542, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.9014 -0.8932
## sample estimates:
## cor
## -0.8974
Najib read me an email listing off the gene names associated with the zymodeme classification. I took those names and cross referenced them against the Leishmania panamensis gene annotations and found the following:
They are:
Given these 6 gene IDs (NH has two gene IDs associated with it), I can do some looking for specific differences among the various samples.
The following creates a colorspace (red to green) heatmap showing the observed expression of these genes in every sample.
my_genes <- c("LPAL13_120010900", "LPAL13_340013000", "LPAL13_000054100",
"LPAL13_140006100", "LPAL13_180018500", "LPAL13_320022300",
"other")
my_names <- c("ALAT", "ASAT", "G6PD", "NHv1", "NHv2", "MPI", "other")
zymo_expt <- exclude_genes_expt(zy_norm, ids = my_genes, method = "keep")
## Before removal, there were 8544 genes, now there are 6.
## There are 34 samples which kept less than 90 percent counts.
## TMRC20001 TMRC20065 TMRC20005 TMRC20066 TMRC20039 TMRC20037 TMRC20038 TMRC20067
## 0.1313 0.1250 0.1325 0.1059 0.1303 0.1102 0.1129 0.1165
## TMRC20068 TMRC20041 TMRC20015 TMRC20009 TMRC20010 TMRC20016 TMRC20011 TMRC20012
## 0.1155 0.1181 0.1147 0.1137 0.1098 0.1059 0.1103 0.1207
## TMRC20013 TMRC20017 TMRC20014 TMRC20018 TMRC20021 TMRC20022 TMRC20077 TMRC20074
## 0.1205 0.1064 0.1089 0.1147 0.1063 0.1310 0.1221 0.1209
## TMRC20063 TMRC20053 TMRC20052 TMRC20064 TMRC20075 TMRC20051 TMRC20050 TMRC20049
## 0.1169 0.1184 0.1105 0.1140 0.1111 0.1285 0.1155 0.1400
## TMRC20062 TMRC20054
## 0.1288 0.1279
zymo_heatmap <- plot_sample_heatmap(zymo_expt, row_label = my_names)
zymo_heatmap
In contrast, the following plots take the set of genes which are shared among all differential expression methods (|lfc| >= 1.0 and adjp <= 0.05) and use them to make categories of genes which are increased in 2.3 or 2.2.
shared_zymo <- intersect_significant(zy_table)
up_shared <- shared_zymo[["ups"]][[1]][["data"]][["all"]]
rownames(up_shared)
## [1] "LPAL13_000033300" "LPAL13_000012000" "LPAL13_310031300" "LPAL13_000038500"
## [5] "LPAL13_000038400" "LPAL13_000012100" "LPAL13_340039600" "LPAL13_310031000"
## [9] "LPAL13_310039200" "LPAL13_050005000" "LPAL13_350063000" "LPAL13_210015500"
## [13] "LPAL13_140019300" "LPAL13_180013900" "LPAL13_340039700" "LPAL13_170015400"
## [17] "LPAL13_270034100" "LPAL13_350013200" "LPAL13_250006300" "LPAL13_140019100"
## [21] "LPAL13_350012400" "LPAL13_350073400" "LPAL13_330021800" "LPAL13_240009700"
## [25] "LPAL13_000052700" "LPAL13_140019200" "LPAL13_250025700" "LPAL13_320038700"
## [29] "LPAL13_330021900" "LPAL13_210005000" "LPAL13_350073200" "LPAL13_310032500"
## [33] "LPAL13_230011200" "LPAL13_310028500" "LPAL13_230011400" "LPAL13_230011500"
## [37] "LPAL13_160014500" "LPAL13_050009600" "LPAL13_230011300" "LPAL13_040007800"
## [41] "LPAL13_160014100"
upshared_expt <- exclude_genes_expt(zy_norm, ids = rownames(up_shared), method = "keep")
## Before removal, there were 8544 genes, now there are 41.
## There are 34 samples which kept less than 90 percent counts.
## TMRC20001 TMRC20065 TMRC20005 TMRC20066 TMRC20039 TMRC20037 TMRC20038 TMRC20067
## 0.32583 0.40481 0.09692 0.35825 0.14393 0.38093 0.50866 0.29289
## TMRC20068 TMRC20041 TMRC20015 TMRC20009 TMRC20010 TMRC20016 TMRC20011 TMRC20012
## 0.34705 0.14461 0.39129 0.12494 0.37945 0.27684 0.13830 0.11794
## TMRC20013 TMRC20017 TMRC20014 TMRC20018 TMRC20021 TMRC20022 TMRC20077 TMRC20074
## 0.32652 0.13750 0.14642 0.30204 0.33880 0.11629 0.11344 0.14166
## TMRC20063 TMRC20053 TMRC20052 TMRC20064 TMRC20075 TMRC20051 TMRC20050 TMRC20049
## 0.12976 0.15457 0.40198 0.36990 0.30504 0.56640 0.12676 0.14631
## TMRC20062 TMRC20054
## 0.59300 0.49938
We can plot a quick heatmap to get a sense of the differences observed between the genes which are different between the two zymodemes.
high_23_heatmap <- plot_sample_heatmap(upshared_expt, row_label = rownames(up_shared))
high_23_heatmap
down_shared <- shared_zymo[["downs"]][[1]][["data"]][["all"]]
downshared_expt <- exclude_genes_expt(zy_norm, ids = rownames(down_shared), method = "keep")
## Before removal, there were 8544 genes, now there are 67.
## There are 34 samples which kept less than 90 percent counts.
## TMRC20001 TMRC20065 TMRC20005 TMRC20066 TMRC20039 TMRC20037 TMRC20038 TMRC20067
## 0.2732 0.2311 0.7574 0.2684 0.7461 0.2422 0.2348 0.2786
## TMRC20068 TMRC20041 TMRC20015 TMRC20009 TMRC20010 TMRC20016 TMRC20011 TMRC20012
## 0.2430 0.7511 0.2312 0.7425 0.2165 0.2537 0.6578 0.6401
## TMRC20013 TMRC20017 TMRC20014 TMRC20018 TMRC20021 TMRC20022 TMRC20077 TMRC20074
## 0.2153 0.7571 0.7355 0.2065 0.2035 0.7990 0.6559 0.7732
## TMRC20063 TMRC20053 TMRC20052 TMRC20064 TMRC20075 TMRC20051 TMRC20050 TMRC20049
## 0.7174 0.6859 0.2162 0.2381 0.2163 0.2301 0.7039 0.8006
## TMRC20062 TMRC20054
## 0.2268 0.2470
high_22_heatmap <- plot_sample_heatmap(downshared_expt, row_label = rownames(down_shared))
high_22_heatmap
Now I will combine our previous samples and our new samples in the hopes of finding variant positions which help elucidate currently unknown aspects of either group via their clustering to known samples from the other group. In other words, we do not know the zymodeme annotations for the old samples nor the strain identities (or the shortcut ‘chronic vs. self-healing’) for the new samples. I hope to make educated guesses given the variant profiles. There are some differences in how the previous and current data sets were analyzed (though I have since redone the old samples so it should be trivial to remove those differences now).
I added our 2016 data to a specific TMRC2 sample sheet, dated 20191203. Thus I will load the data here. That previous data was mapped using tophat, so I will also need to make some changes to the gene names to accomodate the two mappings.
old_expt <- sm(create_expt("sample_sheets/tmrc2_samples_20191203.xlsx",
file_column = "tophat2file"))
tt <- lp_expt[["expressionset"]]
rownames(tt) <- gsub(pattern = "^exon_", replacement = "", x = rownames(tt))
rownames(tt) <- gsub(pattern = "\\.E1$", replacement = "", x = rownames(tt))
lp_expt$expressionset <- tt
tt <- old_expt$expressionset
rownames(tt) <- gsub(pattern = "^exon_", replacement = "", x = rownames(tt))
rownames(tt) <- gsub(pattern = "\\.1$", replacement = "", x = rownames(tt))
old_expt$expressionset <- tt
rm(tt)
One other important caveat, we have a group of new samples which have not yet run through the variant search pipeline, so I need to remove them from consideration. Though it looks like they finished overnight…
## The next line drops the samples which are missing the SNP pipeline.
lp_snp <- subset_expt(lp_expt, subset="!is.na(pData(lp_expt)[['bcftable']])")
## subset_expt(): There were 70, now there are 66 samples.
new_snps <- sm(count_expt_snps(lp_snp, annot_column = "bcftable"))
old_snps <- sm(count_expt_snps(old_expt, annot_column = "bcftable", snp_column = 2))
both_snps <- combine_expts(new_snps, old_snps)
both_norm <- sm(normalize_expt(both_snps, transform = "log2", convert = "cpm", filter = TRUE))
## strains <- both_norm[["design"]][["strain"]]
both_strain <- set_expt_conditions(both_norm, fact = "strain")
The data structure ‘both_norm’ now contains our 2016 data along with the newer data collected since 2019.
The following plot shows the SNP profiles of all samples (old and new) where the colors at the top show either the 2.2 strains (orange), 2.3 strains (green), the previous samples (purple), or the various lab strains (pink etc).
old_new_variant_heatmap <- plot_disheat(both_norm)
pp(file = "images/raw_snp_disheat.png", image = old_new_variant_heatmap,
height = 12, width = 12)
The function get_snp_sets() takes the provided metadata factor (in this case ‘condition’) and looks for variants which are exclusive to each element in it. In this case, this is looking for differences between 2.2 and 2.3, as well as the set shared among them.
snp_sets <- get_snp_sets(both_snps, factor = "condition")
## The factor z2.3 has 18 rows.
## The factor z2.2 has 13 rows.
## The factor unknown has 16 rows.
## The factor null has 19 rows.
## The factor sh has 13 rows.
## The factor chr has 14 rows.
## The factor inf has 6 rows.
## Iterating over 727 elements.
both_expt <- combine_expts(lp_expt, old_expt)
snp_genes <- sm(snps_vs_genes(both_expt, snp_sets, expt_name_col = "chromosome"))
## I think we have some metrics here we can plot...
snp_subset <- sm(snp_subset_genes(
both_expt, both_snps,
genes = c("LPAL13_120010900", "LPAL13_340013000", "LPAL13_000054100",
"LPAL13_140006100", "LPAL13_180018500", "LPAL13_320022300")))
zymo_heat <- plot_sample_heatmap(snp_subset, row_label = rownames(exprs(snp_subset)))
zymo_heat
Didn’t I create a set of densities by chromosome? Oh I think they come in from get_snp_sets()
clinical_sets <- get_snp_sets(new_snps, factor = "clinicalresponse")
## The factor cure has 26 rows.
## The factor failure has 21 rows.
## The factor laboratory line has only 1 row.
## The factor laboratory line miltefosine resistant has only 1 row.
## The factor nd has 4 rows.
## The factor notapplicable has 9 rows.
## The factor reference strain has 4 rows.
## Iterating over 695 elements.
density_vec <- clinical_sets[["density"]]
chromosome_idx <- grep(pattern = "LpaL", x = names(density_vec))
density_df <- as.data.frame(density_vec[chromosome_idx])
density_df[["chr"]] <- rownames(density_df)
colnames(density_df) <- c("density_vec", "chr")
ggplot(density_df, aes_string(x = "chr", y = "density_vec")) +
ggplot2::geom_col() +
ggplot2::theme(axis.text = ggplot2::element_text(size = 10, colour = "black"),
axis.text.x = ggplot2::element_text(angle = 90, vjust = 0.5))
## clinical_written <- write_variants(new_snps)
clinical_genes <- sm(snps_vs_genes(lp_expt, clinical_sets, expt_name_col = "chromosome"))
snp_density <- merge(as.data.frame(clinical_genes[["summary_by_gene"]]),
as.data.frame(fData(lp_expt)),
by = "row.names")
snp_density <- snp_density[, c(1, 2, 4, 15)]
colnames(snp_density) <- c("name", "snps", "product", "length")
snp_density[["product"]] <- tolower(snp_density[["product"]])
snp_density[["length"]] <- as.numeric(snp_density[["length"]])
snp_density[["density"]] <- snp_density[["snps"]] / snp_density[["length"]]
snp_idx <- order(snp_density[["density"]], decreasing = TRUE)
snp_density <- snp_density[snp_idx, ]
removers <- c("amastin", "gp63", "leishmanolysin")
for (r in removers) {
drop_idx <- grepl(pattern = r, x = snp_density[["product"]])
snp_density <- snp_density[!drop_idx, ]
}
## Filter these for [A|a]mastin gp63 Leishmanolysin
clinical_snps <- snps_intersections(lp_expt, clinical_sets, chr_column = "chromosome")
fail_ref_snps <- as.data.frame(clinical_snps[["inters"]][["failure, reference strain"]])
cure_snps <- as.data.frame(clinical_snps[["inters"]][["cure"]])
head(fail_ref_snps)
## seqnames start end width strand
## chr_LpaL13-13_pos_167047_ref_G_alt_C LpaL13-13 167047 167048 2 +
## chr_LpaL13-15_pos_42885_ref_A_alt_G LpaL13-15 42885 42886 2 +
## chr_LpaL13-20.1_pos_111781_ref_T_alt_C LpaL13-20.1 111781 111782 2 +
## chr_LpaL13-25_pos_147406_ref_G_alt_C LpaL13-25 147406 147407 2 +
## chr_LpaL13-29_pos_649951_ref_T_alt_C LpaL13-29 649951 649952 2 +
## chr_LpaL13-30_pos_597834_ref_T_alt_G LpaL13-30 597834 597835 2 +
head(cure_snps)
## seqnames start end width strand
## chr_LpaL13-08_pos_184791_ref_T_alt_A LpaL13-08 184791 184792 2 +
## chr_LpaL13-20.1_pos_369935_ref_C_alt_T LpaL13-20.1 369935 369936 2 +
## chr_LpaL13-20.1_pos_370282_ref_C_alt_T LpaL13-20.1 370282 370283 2 +
## chr_LpaL13-20.1_pos_371356_ref_T_alt_C LpaL13-20.1 371356 371357 2 +
## chr_LpaL13-20.1_pos_380785_ref_A_alt_G LpaL13-20.1 380785 380786 2 +
## chr_LpaL13-20.1_pos_382801_ref_A_alt_C LpaL13-20.1 382801 382802 2 +
annot <- fData(lp_expt)
clinical_interest <- as.data.frame(clinical_snps[["gene_summaries"]][["cure"]])
clinical_interest <- merge(clinical_interest,
as.data.frame(clinical_snps[["gene_summaries"]][["failure, reference strain"]]),
by = "row.names")
rownames(clinical_interest) <- clinical_interest[["Row.names"]]
clinical_interest[["Row.names"]] <- NULL
colnames(clinical_interest) <- c("cure_snps","fail_snps")
annot <- merge(annot, clinical_interest, by = "row.names")
rownames(annot) <- annot[["Row.names"]]
annot[["Row.names"]] <- NULL
fData(lp_expt$expressionset) <- annot
The heatmap produced here should show the variants only for the zymodeme genes.
I am thinking that if we find clusters of locations which are variant, that might provide some PCR testing possibilities.
new_sets <- get_snp_sets(new_snps, factor = "phenotypiccharacteristics")
## The factor 22 has 13 rows.
## The factor 23 has 18 rows.
## The factor laboratory line has only 1 row.
## The factor laboratory line miltefosine resistant has only 1 row.
## The factor notapplicable has 29 rows.
## The factor reference strain has 4 rows.
## Iterating over 695 elements.
summary(new_sets)
## Length Class Mode
## medians 7 data.frame list
## possibilities 6 -none- character
## intersections 58 -none- list
## chr_data 695 -none- list
## set_names 64 -none- list
## invert_names 64 -none- list
## density 695 -none- numeric
## 1000000: 2.2
## 0100000: 2.3
summary(new_sets[["intersections"]][["10000"]])
## Length Class Mode
## 0 NULL NULL
summary(new_sets[["intersections"]][["01000"]])
## Length Class Mode
## 0 NULL NULL
Thus we see that there are 511 variants associated with 2.2 and 49,790 associated with 2.3.
The following function uses the positional data to look for sequential mismatches associated with zymodeme in the hopes that there will be some regions which would provide good potential targets for a PCR-based assay.
sequential_variants <- function(snp_sets, conditions = NULL, minimum = 3, maximum_separation = 3) {
if (is.null(conditions)) {
conditions <- 1
}
intersection_sets <- snp_sets[["intersections"]]
intersection_names <- snp_sets[["set_names"]]
chosen_intersection <- 1
if (is.numeric(conditions)) {
chosen_intersection <- conditions
} else {
intersection_idx <- intersection_names == conditions
chosen_intersection <- names(intersection_names)[intersection_idx]
}
possible_positions <- intersection_sets[[chosen_intersection]]
position_table <- data.frame(row.names = possible_positions)
pat <- "^chr_(.+)_pos_(.+)_ref_.*$"
position_table[["chr"]] <- gsub(pattern = pat, replacement = "\\1", x = rownames(position_table))
position_table[["pos"]] <- as.numeric(gsub(pattern = pat, replacement = "\\2", x = rownames(position_table)))
position_idx <- order(position_table[, "chr"], position_table[, "pos"])
position_table <- position_table[position_idx, ]
position_table[["dist"]] <- 0
last_chr <- ""
for (r in 1:nrow(position_table)) {
this_chr <- position_table[r, "chr"]
if (r == 1) {
position_table[r, "dist"] <- position_table[r, "pos"]
last_chr <- this_chr
next
}
if (this_chr == last_chr) {
position_table[r, "dist"] <- position_table[r, "pos"] - position_table[r - 1, "pos"]
} else {
position_table[r, "dist"] <- position_table[r, "pos"]
}
last_chr <- this_chr
}
## Working interactively here.
doubles <- position_table[["dist"]] == 1
doubles <- position_table[doubles, ]
write.csv(doubles, "doubles.csv")
one_away <- position_table[["dist"]] == 2
one_away <- position_table[one_away, ]
write.csv(one_away, "one_away.csv")
two_away <- position_table[["dist"]] == 3
two_away <- position_table[two_away, ]
write.csv(two_away, "two_away.csv")
combined <- rbind(doubles, one_away)
combined <- rbind(combined, two_away)
position_idx <- order(combined[, "chr"], combined[, "pos"])
combined <- combined[position_idx, ]
this_chr <- ""
for (r in 1:nrow(combined)) {
this_chr <- combined[r, "chr"]
if (r == 1) {
combined[r, "dist_pair"] <- combined[r, "pos"]
last_chr <- this_chr
next
}
if (this_chr == last_chr) {
combined[r, "dist_pair"] <- combined[r, "pos"] - combined[r - 1, "pos"]
} else {
combined[r, "dist_pair"] <- combined[r, "pos"]
}
last_chr <- this_chr
}
dist_pair_maximum <- 1000
dist_pair_minimum <- 200
dist_pair_idx <- combined[["dist_pair"]] <= dist_pair_maximum &
combined[["dist_pair"]] >= dist_pair_minimum
remaining <- combined[dist_pair_idx, ]
no_weak_idx <- grepl(pattern="ref_(G|C)", x=rownames(remaining))
remaining <- remaining[no_weak_idx, ]
print(head(table(position_table[["dist"]])))
sequentials <- position_table[["dist"]] <= maximum_separation
message("There are ", sum(sequentials), " candidate regions.")
## The following can tell me how many runs of each length occurred, that is not quite what I want.
## Now use run length encoding to find the set of sequential sequentials!
rle_result <- rle(sequentials)
rle_values <- rle_result[["values"]]
## The following line is equivalent to just leaving values alone:
## true_values <- rle_result[["values"]] == TRUE
rle_lengths <- rle_result[["lengths"]]
true_sequentials <- rle_lengths[rle_values]
rle_idx <- cumsum(rle_lengths)[which(rle_values)]
position_table[["last_sequential"]] <- 0
count <- 0
for (r in rle_idx) {
count <- count + 1
position_table[r, "last_sequential"] <- true_sequentials[count]
}
message("The maximum sequential set is: ", max(position_table[["last_sequential"]]), ".")
wanted_idx <- position_table[["last_sequential"]] >= minimum
wanted <- position_table[wanted_idx, c("chr", "pos")]
return(wanted)
}
zymo22_sequentials <- sequential_variants(new_sets, conditions = "22", minimum=1, maximum_separation=2)
dim(zymo22_sequentials)
## 7 candidate regions for zymodeme 2.2 -- thus I am betting that the reference strain is a 2.2
zymo23_sequentials <- sequential_variants(new_sets, conditions = "23",
minimum = 2, maximum_separation = 2)
dim(zymo23_sequentials)
## In contrast, there are lots (587) of interesting regions for 2.3!
The first 4 candidate regions from my set of remaining: * Chr Pos. Distance * LpaL13-15 238433 448 * LpaL13-18 142844 613 * LpaL13-29 830342 252 * LpaL13-33 1331507 843
Lets define a couple of terms: * Third: Each of the 4 above positions. * Second: Third - Distance * End: Third + PrimerLen * Start: Second - Primerlen
In each instance, these are the last positions, so we want to grab three things:
## * LpaL13-15 238433 448
first_candidate_chr <- genome[["LpaL13_15"]]
primer_length <- 22
amplicon_length <- 448
first_candidate_third <- 238433
first_candidate_second <- first_candidate_third - amplicon_length
first_candidate_start <- first_candidate_second - primer_length
first_candidate_end <- first_candidate_third + primer_length
first_candidate_region <- subseq(first_candidate_chr, first_candidate_start, first_candidate_end)
first_candidate_region
first_candidate_5p <- subseq(first_candidate_chr, first_candidate_start, first_candidate_second)
as.character(first_candidate_5p)
first_candidate_3p <- spgs::reverseComplement(subseq(first_candidate_chr, first_candidate_third, first_candidate_end))
first_candidate_3p
## * LpaL13-18 142844 613
second_candidate_chr <- genome[["LpaL13_18"]]
primer_length <- 22
amplicon_length <- 613
second_candidate_third <- 142844
second_candidate_second <- second_candidate_third - amplicon_length
second_candidate_start <- second_candidate_second - primer_length
second_candidate_end <- second_candidate_third + primer_length
second_candidate_region <- subseq(second_candidate_chr, second_candidate_start, second_candidate_end)
second_candidate_region
second_candidate_5p <- subseq(second_candidate_chr, second_candidate_start, second_candidate_second)
as.character(second_candidate_5p)
second_candidate_3p <- spgs::reverseComplement(subseq(second_candidate_chr, second_candidate_third, second_candidate_end))
second_candidate_3p
## * LpaL13-29 830342 252
third_candidate_chr <- genome[["LpaL13_29"]]
primer_length <- 22
amplicon_length <- 252
third_candidate_third <- 830342
third_candidate_second <- third_candidate_third - amplicon_length
third_candidate_start <- third_candidate_second - primer_length
third_candidate_end <- third_candidate_third + primer_length
third_candidate_region <- subseq(third_candidate_chr, third_candidate_start, third_candidate_end)
third_candidate_region
third_candidate_5p <- subseq(third_candidate_chr, third_candidate_start, third_candidate_second)
as.character(third_candidate_5p)
third_candidate_3p <- spgs::reverseComplement(subseq(third_candidate_chr, third_candidate_third, third_candidate_end))
third_candidate_3p
## You are a garbage polypyrimidine tract.
## Which is actually interesting if the mutations mess it up.
## * LpaL13-33 1331507 843
fourth_candidate_chr <- genome[["LpaL13_33"]]
primer_length <- 22
amplicon_length <- 843
fourth_candidate_third <- 1331507
fourth_candidate_second <- fourth_candidate_third - amplicon_length
fourth_candidate_start <- fourth_candidate_second - primer_length
fourth_candidate_end <- fourth_candidate_third + primer_length
fourth_candidate_region <- subseq(fourth_candidate_chr, fourth_candidate_start, fourth_candidate_end)
fourth_candidate_region
fourth_candidate_5p <- subseq(fourth_candidate_chr, fourth_candidate_start, fourth_candidate_second)
as.character(fourth_candidate_5p)
fourth_candidate_3p <- spgs::reverseComplement(subseq(fourth_candidate_chr, fourth_candidate_third, fourth_candidate_end))
fourth_candidate_3p
I made a fun little function which should find regions which have lots of variants associated with a given experimental factor.
pheno <- subset_expt(lp_expt, subset = "condition=='z2.2'|condition=='z2.3'")
pheno <- subset_expt(pheno, subset = "!is.na(pData(pheno)[['bcftable']])")
pheno_snps <- sm(count_expt_snps(pheno, annot_column = "bcftable"))
fun_stuff <- snp_density_primers(pheno_snps,
bsgenome="BSGenome.Leishmania.panamensis.MHOMCOL81L13.v53",
gff="reference/TriTrypDB-53_LpanamensisMHOMCOL81L13.gff")
drop_scaffolds <- grepl(x = rownames(fun_stuff$favorites), pattern = "SCAF")
favorite_primer_regions <- fun_stuff[["favorites"]][!drop_scaffolds, ]
favorite_primer_regions[["bin"]] <- rownames(favorite_primer_regions)
library(dplyr)
favorite_primer_regions <- favorite_primer_regions %>%
relocate(bin)
Here is my note from our meeting:
Cross reference primers to DE genes of 2.2/2.3 and/or resistance/suscpetible, add a column to the primer spreadsheet with the DE genes (in retrospect I am guessing this actually means to put the logFC as a column.
One nice thing, I did a semantic removal on the lp_expt, so the set of logFC/pvalues should not have any of the offending types; thus I should be able to automagically get rid of them in the merge.
logfc <- zy_table[["data"]][["z23_vs_z22"]]
logfc_columns <- logfc[, c("deseq_logfc", "deseq_adjp")]
colnames(logfc_columns) <- c("z23_logfc", "z23_adjp")
new_table <- merge(favorite_primer_regions, logfc_columns,
by.x = "closest_gene_before_id", by.y = "row.names")
sus <- sus_table[["data"]][["sensitive_vs_resistant"]]
sus_columns <- sus[, c("deseq_logfc", "deseq_adjp")]
colnames(sus_columns) <- c("sus_logfc", "sus_adjp")
new_table <- merge(new_table, sus_columns,
by.x = "closest_gene_before_id", by.y = "row.names") %>%
relocate(bin)
written <- write_xlsx(data=new_table,
excel="excel/favorite_primers_xref_zy_sus.
## Error: <text>:13:29: unexpected INCOMPLETE_STRING
## 12: written <- write_xlsx(data=new_table,
## 13: excel="excel/favorite_primers_xref_zy_sus.
## ^
We can cross reference the variants against the zymodeme status and plot a heatmap of the results and hopefully see how they separate.
snp_genes <- sm(snps_vs_genes(lp_expt, new_sets, expt_name_col = "chromosome"))
new_zymo_norm <- normalize_expt(new_snps, filter = TRUE, convert = "cpm", norm = "quant", transform = TRUE)
## Removing 0 low-count genes (568413 remaining).
## transform_counts: Found 18767635 values equal to 0, adding 1 to the matrix.
new_zymo_norm <- set_expt_conditions(new_zymo_norm, fact = "phenotypiccharacteristics")
zymo_heat <- plot_disheat(new_zymo_norm)
zymo_heat[["plot"]]
Now let us try to make a heatmap which includes some of the annotation data.
des <- both_norm[["design"]]
undef_idx <- is.na(des[["strain"]])
des[undef_idx, "strain"] <- "unknown"
##hmcols <- colorRampPalette(c("yellow","black","darkblue"))(256)
correlations <- hpgl_cor(exprs(both_norm))
zymo_missing_idx <- is.na(des[["phenotypiccharacteristics"]])
des[["phenotypiccharacteristics"]] <- as.character(des[["phenotypiccharacteristics"]])
des[["clinicalcategorical"]] <- as.character(des[["clinicalcategorical"]])
des[zymo_missing_idx, "phenotypiccharacteristics"] <- "unknown"
mydendro <- list(
"clustfun" = hclust,
"lwd" = 2.0)
col_data <- as.data.frame(des[, c("phenotypiccharacteristics", "clinicalcategorical")])
unknown_clinical <- is.na(col_data[["clinicalcategorical"]])
row_data <- as.data.frame(des[, c("strain")])
colnames(col_data) <- c("zymodeme", "outcome")
col_data[unknown_clinical, "outcome"] <- "undefined"
colnames(row_data) <- c("strain")
myannot <- list(
"Col" = list("data" = col_data),
"Row" = list("data" = row_data))
myclust <- list("cuth" = 1.0,
"col" = BrewerClusterCol)
mylabs <- list(
"Row" = list("nrow" = 4),
"Col" = list("nrow" = 4))
hmcols <- colorRampPalette(c("darkblue", "beige"))(240)
map1 <- annHeatmap2(
correlations,
dendrogram = mydendro,
annotation = myannot,
cluster = myclust,
labels = mylabs,
## The following controls if the picture is symmetric
scale = "none",
col = hmcols)
## Warning in breakColors(breaks, col): more colors than classes: ignoring 29 last
## colors
pp(file = "images/dendro_heatmap.png", image = map1, height = 20, width = 20)
## annotated Heatmap
##
## Rows: 'dendrogram' with 2 branches and 99 members total, at height 6.095
## 11 annotation variable(s)
## Cols: 'dendrogram' with 2 branches and 99 members total, at height 6.095
## 12 annotation variable(s)
Print the larger heatmap so that all the labels appear. Keep in mind that as we get more samples, this image needs to continue getting bigger.
big heatmap
xref_prop <- table(pheno_snps[["conditions"]])
## Error in eval(quote(list(...)), env): object 'pheno_snps' not found
pheno_snps$conditions
## Error in eval(expr, envir, enclos): object 'pheno_snps' not found
idx_tbl <- exprs(pheno_snps) > 5
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'exprs': object 'pheno_snps' not found
new_tbl <- data.frame(row.names = rownames(exprs(pheno_snps)))
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'rownames': error in evaluating the argument 'object' in selecting a method for function 'exprs': object 'pheno_snps' not found
for (n in names(xref_prop)) {
new_tbl[[n]] <- 0
idx_cols <- which(pheno_snps[["conditions"]] == n)
prop_col <- rowSums(idx_tbl[, idx_cols]) / xref_prop[n]
new_tbl[n] <- prop_col
}
## Error in eval(expr, envir, enclos): object 'xref_prop' not found
keepers <- grepl(x = rownames(new_tbl), pattern = "LpaL13")
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'grepl': error in evaluating the argument 'x' in selecting a method for function 'rownames': object 'new_tbl' not found
new_tbl <- new_tbl[keepers, ]
## Error in eval(expr, envir, enclos): object 'new_tbl' not found
new_tbl[["strong22"]] <- 1.001 - new_tbl[["z2.2"]]
## Error in eval(expr, envir, enclos): object 'new_tbl' not found
new_tbl[["strong23"]] <- 1.001 - new_tbl[["z2.3"]]
## Error in eval(expr, envir, enclos): object 'new_tbl' not found
s22_na <- new_tbl[["strong22"]] > 1
## Error in eval(expr, envir, enclos): object 'new_tbl' not found
new_tbl[s22_na, "strong22"] <- 1
## Error in new_tbl[s22_na, "strong22"] <- 1: object 'new_tbl' not found
s23_na <- new_tbl[["strong23"]] > 1
## Error in eval(expr, envir, enclos): object 'new_tbl' not found
new_tbl[s23_na, "strong23"] <- 1
## Error in new_tbl[s23_na, "strong23"] <- 1: object 'new_tbl' not found
new_tbl[["SNP"]] <- rownames(new_tbl)
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'rownames': object 'new_tbl' not found
new_tbl[["Chromosome"]] <- gsub(x = new_tbl[["SNP"]], pattern = "chr_(.*)_pos_.*", replacement = "\\1")
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'gsub': object 'new_tbl' not found
new_tbl[["Position"]] <- gsub(x = new_tbl[["SNP"]], pattern = ".*_pos_(\\d+)_.*", replacement = "\\1")
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'gsub': object 'new_tbl' not found
new_tbl <- new_tbl[, c("SNP", "Chromosome", "Position", "strong22", "strong23")]
## Error in eval(expr, envir, enclos): object 'new_tbl' not found
library(CMplot)
## Much appreciate for using CMplot.
## Full description, Bug report, Suggestion and the latest codes:
## https://github.com/YinLiLin/CMplot
simplify <- new_tbl
## Error in eval(expr, envir, enclos): object 'new_tbl' not found
simplify[["strong22"]] <- NULL
## Error in simplify[["strong22"]] <- NULL: object 'simplify' not found
CMplot(simplify, bin.size = 100000)
## Error in is.data.frame(x): object 'simplify' not found
CMplot(new_tbl, plot.type="m", multracks=TRUE, threshold = c(0.01, 0.05),
threshold.lwd=c(1,1), threshold.col=c("black","grey"),
amplify=TRUE, bin.size=10000,
chr.den.col=c("darkgreen", "yellow", "red"),
signal.col=c("red", "green", "blue"),
signal.cex=1, file="jpg", memo="", dpi=300, file.output=TRUE, verbose=TRUE)
## Error in is.data.frame(x): object 'new_tbl' not found
This tool looks a little opaque, but provides sample data with things that make sense to me and should be pretty easy to recapitulate in our data.
## For this, let us use the 'new_snps' data structure.
## Caveat here: these need to be coerced to numbers.
my_covariates <- pData(new_snps)[, c("phenotypiccharacteristics", "clinicalcategorical")]
for (col in colnames(my_covariates)) {
my_covariates[[col]] <- as.numeric(as.factor(my_covariates[[col]]))
}
my_covariates <- t(my_covariates)
my_geneloc <- fData(lp_expt)[, c("gid", "chromosome", "start", "end")]
colnames(my_geneloc) <- c("geneid", "chr", "left", "right")
my_ge <- exprs(normalize_expt(lp_expt, transform = "log2", filter = TRUE, convert = "cpm"))
used_samples <- tolower(colnames(my_ge)) %in% colnames(exprs(new_snps))
my_ge <- my_ge[, used_samples]
my_snpsloc <- data.frame(rownames = rownames(exprs(new_snps)))
## Oh, caveat here: Because of the way I stored the data,
## I could have duplicate rows which presumably will make matrixEQTL sad
my_snpsloc[["chr"]] <- gsub(pattern = "^chr_(.+)_pos(.+)_ref_.*$", replacement = "\\1",
x = rownames(my_snpsloc))
my_snpsloc[["pos"]] <- gsub(pattern = "^chr_(.+)_pos(.+)_ref_.*$", replacement = "\\2",
x = rownames(my_snpsloc))
test <- duplicated(my_snpsloc)
## Each duplicated row would be another variant at that position;
## so in theory we would do a rle to number them I am guessing
## However, I do not have different variants so I think I can ignore this for the moment
## but will need to make my matrix either 0 or 1.
if (sum(test) > 0) {
message("There are: ", sum(duplicated), " duplicated entries.")
keep_idx <- ! test
my_snpsloc <- my_snpsloc[keep_idx, ]
}
my_snps <- exprs(new_snps)
one_idx <- my_snps > 0
my_snps[one_idx] <- 1
## Ok, at this point I think I have all the pieces which this method wants...
## Oh, no I guess not; it actually wants the data as a set of filenames...
library(MatrixEQTL)
write.table(my_snps, "eqtl/snps.tsv", na = "NA", col.names = TRUE, row.names = TRUE, sep = "\t", quote = TRUE)
## readr::write_tsv(my_snps, "eqtl/snps.tsv", )
write.table(my_snpsloc, "eqtl/snpsloc.tsv", na = "NA", col.names = TRUE, row.names = TRUE, sep = "\t", quote = TRUE)
## readr::write_tsv(my_snpsloc, "eqtl/snpsloc.tsv")
write.table(as.data.frame(my_ge), "eqtl/ge.tsv", na = "NA", col.names = TRUE, row.names = TRUE, sep = "\t", quote = TRUE)
## readr::write_tsv(as.data.frame(my_ge), "eqtl/ge.tsv")
write.table(as.data.frame(my_geneloc), "eqtl/geneloc.tsv", na = "NA", col.names = TRUE, row.names = TRUE, sep = "\t", quote = TRUE)
## readr::write_tsv(as.data.frame(my_geneloc), "eqtl/geneloc.tsv")
write.table(as.data.frame(my_covariates), "eqtl/covariates.tsv", na = "NA", col.names = TRUE, row.names = TRUE, sep = "\t", quote = TRUE)
## readr::write_tsv(as.data.frame(my_covariates), "eqtl/covariates.tsv")
useModel = modelLINEAR # modelANOVA, modelLINEAR, or modelLINEAR_CROSS
# Genotype file name
SNP_file_name = "eqtl/snps.tsv"
snps_location_file_name = "eqtl/snpsloc.tsv"
expression_file_name = "eqtl/ge.tsv"
gene_location_file_name = "eqtl/geneloc.tsv"
covariates_file_name = "eqtl/covariates.tsv"
# Output file name
output_file_name_cis = tempfile()
output_file_name_tra = tempfile()
# Only associations significant at this level will be saved
pvOutputThreshold_cis = 0.1
pvOutputThreshold_tra = 0.1
# Error covariance matrix
# Set to numeric() for identity.
errorCovariance = numeric()
# errorCovariance = read.table("Sample_Data/errorCovariance.txt");
# Distance for local gene-SNP pairs
cisDist = 1e6
## Load genotype data
snps = SlicedData$new()
snps$fileDelimiter = "\t" # the TAB character
snps$fileOmitCharacters = "NA" # denote missing values;
snps$fileSkipRows = 1 # one row of column labels
snps$fileSkipColumns = 1 # one column of row labels
snps$fileSliceSize = 2000 # read file in slices of 2,000 rows
snps$LoadFile(SNP_file_name)
## Load gene expression data
gene = SlicedData$new()
gene$fileDelimiter = "\t" # the TAB character
gene$fileOmitCharacters = "NA" # denote missing values;
gene$fileSkipRows = 1 # one row of column labels
gene$fileSkipColumns = 1 # one column of row labels
gene$fileSliceSize = 2000 # read file in slices of 2,000 rows
gene$LoadFile(expression_file_name)
## Load covariates
cvrt = SlicedData$new()
cvrt$fileDelimiter = "\t" # the TAB character
cvrt$fileOmitCharacters = "NA" # denote missing values;
cvrt$fileSkipRows = 1 # one row of column labels
cvrt$fileSkipColumns = 1 # one column of row labels
if(length(covariates_file_name) > 0) {
cvrt$LoadFile(covariates_file_name)
}
## Run the analysis
snpspos = read.table(snps_location_file_name, header = TRUE, stringsAsFactors = FALSE)
genepos = read.table(gene_location_file_name, header = TRUE, stringsAsFactors = FALSE)
me = Matrix_eQTL_main(
snps = snps,
gene = gene,
cvrt = cvrt,
output_file_name = output_file_name_tra,
pvOutputThreshold = pvOutputThreshold_tra,
useModel = useModel,
errorCovariance = errorCovariance,
verbose = TRUE,
output_file_name.cis = output_file_name_cis,
pvOutputThreshold.cis = pvOutputThreshold_cis,
snpspos = snpspos,
genepos = genepos,
cisDist = cisDist,
pvalue.hist = "qqplot",
min.pv.by.genesnp = FALSE,
noFDRsaveMemory = FALSE);
if (!isTRUE(get0("skip_load"))) {
pander::pander(sessionInfo())
message(paste0("This is hpgltools commit: ", get_git_commit()))
message(paste0("Saving to ", savefile))
tmp <- sm(saveme(filename = savefile))
}
## If you wish to reproduce this exact build of hpgltools, invoke the following:
## > git clone http://github.com/abelew/hpgltools.git
## > git reset ce88fe653383489942b667e9b834f802690f27d1
## This is hpgltools commit: Wed Nov 3 13:30:35 2021 -0400: ce88fe653383489942b667e9b834f802690f27d1
## Saving to tmrc2_02sample_estimation_v202108.rda.xz
tmp <- loadme(filename = savefile)