L. panamensis 20201022: Sample Estimation

1 Introduction

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:

1. Default trimming was performed.
2. Hisat2 was used to map the remaining reads against the Leishmania panamensis genome revision 36.
3. The alignments from hisat2 were used to count reads/gene against the revision 36 annotations with htseq.
4. These alignments were also passed to the pileup functionality of samtools and the vcf/bcf utilities in order to make a matrix of all observed differences between each sample with respect to the reference.

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.

2 Annotations

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 <- 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

## Extracted all gene ids.

## Attempting to select: ANNOT_GENE_PRODUCT, ANNOT_GENE_TYPE, ANNOT_LOCATION_TEXT, ANNOT_GENE_ENTREZ_ID, ANNOT_GENE_NAME, ANNOT_STRAND, ANNOT_CHROMOSOME, ANNOT_CDS_LENGTH, ANNOT_GENE_PRODUCT

## 'select()' returned 1:1 mapping between keys and columns

lp_go <- load_orgdb_go(pan_db)
lp_lengths <- all_lp_annot[, c("gid", "annot_cds_length")]
colnames(lp_lengths)  <- c("ID", "length")

orthos <- sm(EuPathDB::extract_eupath_orthologs(db=pan_db))

hisat_annot <- all_lp_annot
## rownames(hisat_annot) <- paste0("exon_", rownames(hisat_annot), ".E1")

3 Sample Estimation

The process of sample estimation takes two primary inputs:

1. The sample sheet, which contains all the metadata we currently have on hand, including filenames for the outputs of #3 and #4 above.
2. The gene annotations.

An expressionset is primary 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.

3.1 Generate expressionsets

The first lines of the following block create the Expressionset. All of the following lines perform various normalizations and generate plots from it.

sample_sheet <- glue::glue("sample_sheets/tmrc2_samples_{ver}.xlsx")
lp_expt <- create_expt(sample_sheet,
                       gene_info=hisat_annot,
                       id_column="hpglidentifier",
                       file_column="lpanamensisv36hisatfile")

## Reading the sample metadata.

## Dropped 33 rows from the sample metadata because they were 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: 18 rows(samples) and 56 columns(metadata fields).

## Reading count tables.

## Reading count files with read.table().

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc20001/outputs/hisat2_lpanamensis_v36/concatenated.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows.

## preprocessing/tmrc20002/outputs/hisat2_lpanamensis_v36/concatenated.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20004/outputs/hisat2_lpanamensis_v36/concatenated.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20005/outputs/hisat2_lpanamensis_v36/concatenated.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20006/outputs/hisat2_lpanamensis_v36/concatenated.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20007/outputs/hisat2_lpanamensis_v36/concatenated.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20008/outputs/hisat2_lpanamensis_v36/concatenated.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20015/outputs/hisat2_lpanamensis_v36/r1.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20009/outputs/hisat2_lpanamensis_v36/r1.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20010/outputs/hisat2_lpanamensis_v36/r1.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20016/outputs/hisat2_lpanamensis_v36/r1.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20011/outputs/hisat2_lpanamensis_v36/r1.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20012/outputs/hisat2_lpanamensis_v36/r1.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20013/outputs/hisat2_lpanamensis_v36/r1.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20017/outputs/hisat2_lpanamensis_v36/r1.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20014/outputs/hisat2_lpanamensis_v36/r1.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20018/outputs/hisat2_lpanamensis_v36/r1.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## Finished reading count data.

## Warning in create_expt(sample_sheet, gene_info = hisat_annot, id_column =
## "hpglidentifier", : 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'.

## Saving the expressionset to 'expt.rda'.

## The final expressionset has 8778 rows and 17 columns.

lp_expt <- set_expt_conditions(lp_expt, fact="zymodemecategorical")

libsizes <- plot_libsize(lp_expt)

## The scale difference between the smallest and largest
## libraries is > 10. Assuming a log10 scale is better, set scale=FALSE if not.

libsizes$plot

## I think samples 7,10 should be removed at minimum, probably also 9,11
nonzero <- plot_nonzero(lp_expt)
nonzero$plot

plot_boxplot(lp_expt)

## This data will benefit from being displayed on the log scale.

## If this is not desired, set scale='raw'

## Some entries are 0.  We are on log scale, adding 1 to the data.

## Changed 1920 zero count features.

3.2 Distribution Visualization

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’.

all_norm <- normalize_expt(lp_expt, norm="quant", transform="log2", convert="cpm",
                           batch=FALSE, filter=TRUE)

## This function will replace the expt$expressionset slot with:

## log2(cpm(quant(cbcb(data))))

## It will save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep libsizes in mind
##  when invoking limma.  The appropriate libsize is non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Not correcting the count-data for batch effects.  If batch is
##  included in EdgerR/limma's model, then this is probably wise; but in extreme
##  batch effects this is a good parameter to play with.

## Step 1: performing count filter with option: cbcb

## Removing 155 low-count genes (8623 remaining).

## Step 2: normalizing the data with quant.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## transform_counts: Found 17 values equal to 0, adding 1 to the matrix.

## The method is: raw.

## Step 5: not doing batch correction.

plot_pca(all_norm, plot_title="PCA of parasite expression values")$plot

## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

plot_tsne(all_norm, plot_title="TSNE of parasite expression values")$plot

plot_corheat(all_norm, title="Correlation heatmap of parasite expression values
(Same legend as above)")$plot

plot_sm(all_norm)$plot

## Performing correlation.

plot_variance_coefficients(all_norm)$plot

## Naively calculating coefficient of variation/dispersion with respect to condition.

## Finished calculating dispersion estimates.

plot_sample_cvheatmap(all_norm)$plot

## This function will replace the expt$expressionset slot with:

## cv(data)

## It will save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep libsizes in mind
##  when invoking limma.  The appropriate libsize is non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Leaving the data in its current base format, keep in mind that
##  some metrics are easier to see when the data is log2 transformed, but
##  EdgeR/DESeq do not accept transformed data.

## Leaving the data unconverted.  It is often advisable to cpm/rpkm
##  the data to normalize for sampling differences, keep in mind though that rpkm
##  has some annoying biases, and voom() by default does a cpm (though hpgl_voom()
##  will try to detect this).

## Leaving the data unnormalized.  This is necessary for DESeq, but
##  EdgeR/limma might benefit from normalization.  Good choices include quantile,
##  size-factor, tmm, etc.

## Not correcting the count-data for batch effects.  If batch is
##  included in EdgerR/limma's model, then this is probably wise; but in extreme
##  batch effects this is a good parameter to play with.

## Step 1: performing count filter with option: cv

## Removing 8603 low-count genes (20 remaining).

## Step 2: not normalizing the data.

## Step 3: not converting the data.

## Step 4: not transforming the data.

## The method is: raw.

## Step 5: not doing batch correction.

## The factor unknown has 3 rows.

## The factor z2.2 has 8 rows.

## The factor z2.3 has 6 rows.

## NULL

3.2.1 Notes

The following samples are much lower coverage:

• TMRC20002
• TMRC20006
• TMRC20007
• TMRC20008

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’.

4 Zymodeme analyses

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.

4.1 Differential expression

4.1.1 With respect to zymodeme attribution

zy_expt <- subset_expt(lp_expt, subset="condition=='z2.2'|condition=='z2.3'")

## Using a subset expression.

## There were 17, now there are 14 samples.

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")))

4.2 With respect to cure/failure

cf_expt <- set_expt_conditions(lp_expt, fact="clinicalcategorical")
cf_norm <- normalize_expt(cf_expt, filter=TRUE, norm="quant", transform="log2", convert="cpm", batch="sva")

## This function will replace the expt$expressionset slot with:

## log2(sva(cpm(quant(cbcb(data)))))

## It will save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep libsizes in mind
##  when invoking limma.  The appropriate libsize is non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Warning in normalize_expt(cf_expt, filter = TRUE, norm = "quant", transform =
## "log2", : Quantile normalization and sva do not always play well together.

## Step 1: performing count filter with option: cbcb

## Removing 155 low-count genes (8623 remaining).

## Step 2: normalizing the data with quant.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## transform_counts: Found 17 values equal to 0, adding 1 to the matrix.

## The method is: sva.

## Step 5: doing batch correction with sva.

## Using the current state of normalization.

## Passing the data to all_adjusters using the sva estimate type.

## batch_counts: Before batch/surrogate estimation, 145680 entries are x>1: 99%.

## batch_counts: Before batch/surrogate estimation, 17 entries are x==0: 0%.

## batch_counts: Before batch/surrogate estimation, 894 entries are 0<x<1: 1%.

## The be method chose 4 surrogate variables.

## Estimate type 'sva' is shorthand for 'sva_unsupervised'.

## Other sva options include: sva_supervised and svaseq.

## Attempting sva unsupervised surrogate estimation with 4 surrogates.

## There are 25 (0%) elements which are < 0 after batch correction.

## Setting low elements to zero.

plot_pca(cf_norm)$plot

cf_de <- all_pairwise(cf_expt, filter=TRUE, model_batch="svaseq")

## batch_counts: Before batch/surrogate estimation, 145827 entries are x>1: 99%.

## batch_counts: Before batch/surrogate estimation, 618 entries are x==0: 0%.

## The be method chose 1 surrogate variable.

## Attempting svaseq estimation with 1 surrogate.

## Plotting a PCA before surrogate/batch inclusion.

## Not putting labels on the PC plot.

## Using svaseq to visualize before/after batch inclusion.

## Performing a test normalization with: raw

## This function will replace the expt$expressionset slot with:

## log2(svaseq(cpm(cbcb(data))))

## It will save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep libsizes in mind
##  when invoking limma.  The appropriate libsize is non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Leaving the data unnormalized.  This is necessary for DESeq, but
##  EdgeR/limma might benefit from normalization.  Good choices include quantile,
##  size-factor, tmm, etc.

## Step 1: performing count filter with option: cbcb

## Removing 0 low-count genes (8623 remaining).

## Step 2: not normalizing the data.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## transform_counts: Found 618 values equal to 0, adding 1 to the matrix.

## The method is: svaseq.

## Step 5: doing batch correction with svaseq.

## Using the current state of normalization.

## Passing the data to all_adjusters using the svaseq estimate type.

## batch_counts: Before batch/surrogate estimation, 145262 entries are x>1: 99%.

## batch_counts: Before batch/surrogate estimation, 618 entries are x==0: 0%.

## batch_counts: Before batch/surrogate estimation, 711 entries are 0<x<1: 0%.

## The be method chose 4 surrogate variables.

## Attempting svaseq estimation with 4 surrogates.

## There are 61 (0%) elements which are < 0 after batch correction.

## Setting low elements to zero.

## Not putting labels on the PC plot.

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

cf_table <- combine_de_tables(cf_de, excel=glue::glue("excel/cf_tables-v{ver}.xlsx"))

## Deleting the file excel/cf_tables-v20201022.xlsx before writing the tables.

## Writing a legend of columns.

## Printing a pca plot before/after surrogates/batch estimation.

## Working on table 1/3: fail_vs_cure

## Working on table 2/3: unknown_vs_cure

## Working on table 3/3: unknown_vs_fail

## Adding venn plots for fail_vs_cure.

## Limma expression coefficients for fail_vs_cure; R^2: 0.984; equation: y = 0.967x + 0.2

## Deseq expression coefficients for fail_vs_cure; R^2: 0.985; equation: y = 0.972x + 0.298

## Edger expression coefficients for fail_vs_cure; R^2: 0.985; equation: y = 0.973x + 0.187

## Adding venn plots for unknown_vs_cure.

## Limma expression coefficients for unknown_vs_cure; R^2: 0.937; equation: y = 0.948x + 0.314

## Deseq expression coefficients for unknown_vs_cure; R^2: 0.922; equation: y = 0.983x + 0.187

## Edger expression coefficients for unknown_vs_cure; R^2: 0.922; equation: y = 0.983x + 0.166

## Adding venn plots for unknown_vs_fail.

## Limma expression coefficients for unknown_vs_fail; R^2: 0.921; equation: y = 0.961x + 0.246

## Deseq expression coefficients for unknown_vs_fail; R^2: 0.897; equation: y = 0.993x + 0.0751

## Edger expression coefficients for unknown_vs_fail; R^2: 0.897; equation: y = 0.992x + 0.071

## Writing summary information, compare_plot is: TRUE.

## Performing save of excel/cf_tables-v20201022.xlsx.

cf_sig <- extract_significant_genes(cf_table, excel=glue::glue("excel/cf_sig-v{ver}.xlsx"))

## Writing a legend of columns.

## Printing significant genes to the file: excel/cf_sig-v20201022.xlsx

## The up table fail_vs_cure is empty.

## The down table fail_vs_cure is empty.

## 2/3: Creating significant table up_limma_unknown_vs_cure

## 3/3: Creating significant table up_limma_unknown_vs_fail

## Printing significant genes to the file: excel/cf_sig-v20201022.xlsx

## The up table fail_vs_cure is empty.

## The down table fail_vs_cure is empty.

## 2/3: Creating significant table up_edger_unknown_vs_cure

## 3/3: Creating significant table up_edger_unknown_vs_fail

## Printing significant genes to the file: excel/cf_sig-v20201022.xlsx

## The up table fail_vs_cure is empty.

## The down table fail_vs_cure is empty.

## 2/3: Creating significant table up_deseq_unknown_vs_cure

## 3/3: Creating significant table up_deseq_unknown_vs_fail

## Printing significant genes to the file: excel/cf_sig-v20201022.xlsx

## The up table fail_vs_cure is empty.

## The down table fail_vs_cure is empty.

## 2/3: Creating significant table up_ebseq_unknown_vs_cure

## 3/3: Creating significant table up_ebseq_unknown_vs_fail

## Printing significant genes to the file: excel/cf_sig-v20201022.xlsx

## The up table fail_vs_cure is empty.

## The down table fail_vs_cure is empty.

## The up table unknown_vs_cure is empty.

## The down table unknown_vs_cure is empty.

## 3/3: Creating significant table up_basic_unknown_vs_fail

## Adding significance bar plots.

4.3 Ontology searches

colnames(lp_go)

## [1] "GID"      "GO"       "EVIDENCE" "ONTOLOGY"

lp_go <- lp_go[, c("GID", "GOALL")]

## Error in `[.data.frame`(lp_go, , c("GID", "GOALL")): undefined columns selected

colnames(lp_go) <- c("ID", "GO")

## 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))

4.3.1 A couple plots from the differential expression

4.3.1.1 Number of genes in agreement among DE methods, 2.3 more than 2.2

zy_table[["venns"]][[1]][["p_lfc1"]][["up_noweight"]]

4.3.1.2 Number of genes in agreement among DE methods, 2.2 more than 2.3

zy_table[["venns"]][[1]][["p_lfc1"]][["down_noweight"]]

4.3.1.3 MA plot of the differential expression between the zymodemes.

zy_table$plots[[1]][["deseq_ma_plots"]][["plot"]]

4.3.1.4 goseq ontology plots of groups of genes, 2.3 more than 2.2

zy_go_up$pvalue_plots$bpp_plot_over

4.3.1.5 goseq ontology plots of groups of genes, 2.2 more than 2.3

zy_go_down$pvalue_plots$bpp_plot_over

4.4 Zymodeme enzyme gene IDs

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:

1. ALAT: LPAL13_120010900 – alanine aminotransferase
2. ASAT: LPAL13_340013000 – aspartate aminotransferase
3. G6PD: LPAL13_000054100 – glucase-6-phosphate 1-dehydrogenase
4. NH: LPAL13_14006100, LPAL13_180018500 – inosine-guanine nucleoside hydrolase
5. MPI: LPAL13_320022300 (maybe) – mannose phosphate isomerase (I chose phosphomannose isomerase)

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.

4.4.1 Expression levels of zymodeme genes

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(all_norm, ids=my_genes, method="keep")

## Before removal, there were 8623 entries.

## Now there are 6 entries.

## Percent kept: 0.086, 0.081, 0.084, 0.084, 0.083, 0.086, 0.083, 0.084, 0.083, 0.084, 0.083, 0.083, 0.085, 0.085, 0.083, 0.083, 0.083

## Percent removed: 99.914, 99.919, 99.916, 99.916, 99.917, 99.914, 99.917, 99.916, 99.917, 99.916, 99.917, 99.917, 99.915, 99.915, 99.917, 99.917, 99.917

test <- plot_sample_heatmap(zymo_expt, row_label=my_names)

4.5 Empirically observed Zymodeme genes from differential expression analysis

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_000038400" "LPAL13_310031300" "LPAL13_000012000"
##  [5] "LPAL13_050005000" "LPAL13_210015500" "LPAL13_310039200" "LPAL13_340039600"
##  [9] "LPAL13_000038500" "LPAL13_000041000" "LPAL13_200013000" "LPAL13_330021800"
## [13] "LPAL13_340039700" "LPAL13_170015400" "LPAL13_240009700" "LPAL13_140019300"
## [17] "LPAL13_140019100" "LPAL13_320038700" "LPAL13_210005000" "LPAL13_280037900"
## [21] "LPAL13_230011200" "LPAL13_100010600" "LPAL13_140019200" "LPAL13_310028500"
## [25] "LPAL13_230011500" "LPAL13_000010600" "LPAL13_230011300" "LPAL13_350073200"
## [29] "LPAL13_250025700" "LPAL13_040007800" "LPAL13_200012900" "LPAL13_020006900"

upshared_expt <- exclude_genes_expt(all_norm, ids=rownames(up_shared), method="keep")

## Before removal, there were 8623 entries.

## Now there are 32 entries.

## Percent kept: 0.358, 0.288, 0.263, 0.255, 0.278, 0.264, 0.254, 0.372, 0.260, 0.366, 0.358, 0.270, 0.254, 0.367, 0.267, 0.267, 0.366

## Percent removed: 99.642, 99.712, 99.737, 99.745, 99.722, 99.736, 99.746, 99.628, 99.740, 99.634, 99.642, 99.730, 99.746, 99.633, 99.733, 99.733, 99.634

4.5.1 Heatmap of zymodeme gene expression increased in 2.3 vs. 2.2

test <- plot_sample_heatmap(upshared_expt, row_label=rownames(up_shared))

4.5.2 Heatmap of zymodeme gene expression increased in 2.2 vs. 2.3

down_shared <- shared_zymo[["downs"]][[1]][["data"]][["all"]]
downshared_expt <- exclude_genes_expt(all_norm, ids=rownames(down_shared), method="keep")

## Before removal, there were 8623 entries.

## Now there are 69 entries.

## Percent kept: 0.461, 0.816, 0.742, 0.757, 0.787, 0.751, 0.769, 0.423, 0.790, 0.406, 0.435, 0.791, 0.749, 0.433, 0.779, 0.796, 0.438

## Percent removed: 99.539, 99.184, 99.258, 99.243, 99.213, 99.249, 99.231, 99.577, 99.210, 99.594, 99.565, 99.209, 99.251, 99.567, 99.221, 99.204, 99.562

test <- plot_sample_heatmap(downshared_expt, row_label=rownames(down_shared))

5 SNP profiles

In this block, I am combining our previous samples and our new samples in the hopes of finding variant positions which help elucidate aspects of either the new or old samples. 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. We may be able 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).

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

new_snps <- sm(count_expt_snps(lp_expt, 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_norm <- set_expt_conditions(both_norm, fact="strain")

5.1 Plot of SNP profiles for zymodemes

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).

tt <- plot_disheat(both_norm)

snp_sets <- get_snp_sets(both_snps, factor="condition")

## The factor z2.3 has 6 rows.

## The factor z2.2 has 8 rows.

## The factor unknown has 3 rows.

## The factor sh has 13 rows.

## The factor chr has 14 rows.

## The factor inf has 6 rows.

## Iterating over 722 elements.

both_expt <- combine_expts(lp_expt, old_expt)
snp_genes <- sm(snps_vs_genes(both_expt, snp_sets, expt_name_col="chromosome"))

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)))

6 Clinical response for new samples

clinical_sets <- get_snp_sets(new_snps, factor="clinicalresponse")

## The factor Cure has 6 rows.

## The factor Failure has 7 rows.

## The factor Laboratory line has 3 rows.

## The factor Laboratory line antimony resistant  has only 1 row.

## The factor Laboratory line miltefosine resistant  has only 1 row.

## The factor ND has only 1 row.

## Iterating over 620 elements.

clinical_genes <- sm(snps_vs_genes(lp_expt, clinical_sets, expt_name_col="chromosome"))
clinical_snps <- snps_intersections(lp_expt, clinical_sets, chr_column="chromosome")
head(as.data.frame(clinical_snps$inters[["Failure"]]))

##                                       seqnames  start    end width strand
## chr_LpaL13-01_pos_100120_ref_A_alt_C LpaL13-01 100120 100121     2      +
## chr_LpaL13-01_pos_100402_ref_A_alt_C LpaL13-01 100402 100403     2      +
## chr_LpaL13-01_pos_10056_ref_T_alt_A  LpaL13-01  10056  10057     2      +
## chr_LpaL13-01_pos_10158_ref_A_alt_G  LpaL13-01  10158  10159     2      +
## chr_LpaL13-01_pos_104367_ref_T_alt_C LpaL13-01 104367 104368     2      +
## chr_LpaL13-01_pos_112913_ref_A_alt_G LpaL13-01 112913 112914     2      +

head(as.data.frame(clinical_snps$inters[["Cure"]]))

##                                       seqnames  start    end width strand
## chr_LpaL13-02_pos_58689_ref_A_alt_G  LpaL13-02  58689  58690     2      +
## chr_LpaL13-02_pos_58964_ref_G_alt_C  LpaL13-02  58964  58965     2      +
## chr_LpaL13-03_pos_159681_ref_T_alt_C LpaL13-03 159681 159682     2      +
## chr_LpaL13-03_pos_189630_ref_C_alt_A LpaL13-03 189630 189631     2      +
## chr_LpaL13-04_pos_51133_ref_A_alt_G  LpaL13-04  51133  51134     2      +
## chr_LpaL13-07_pos_353845_ref_A_alt_G LpaL13-07 353845 353846     2      +

head(clinical_snps$gene_summaries$Failure)

## LPAL13_110010300 LPAL13_340060500 LPAL13_000027300 LPAL13_270015100 
##               18               17               15               15 
## LPAL13_310006100 LPAL13_310013800 
##               14               14

head(clinical_snps$gene_summaries$Cure)

## LPAL13_200020000 LPAL13_000018700 LPAL13_080009800 LPAL13_110015800 
##                7                5                4                4 
## LPAL13_000030400 LPAL13_230026400 
##                3                3

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"]]), 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

7 Zymodeme for new samples

The heatmap produced here should show the variants only for the zymodeme genes.

new_sets <- get_snp_sets(new_snps, factor="phenotypiccharacteristics")

## The factor 2.2 has 8 rows.

## The factor 2.3 has 6 rows.

## The factor Laboratory line has 3 rows.

## The factor Laboratory line antimony resistant  has only 1 row.

## The factor Laboratory line miltefosine resistant  has only 1 row.

## Iterating over 620 elements.

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)

## This function will replace the expt$expressionset slot with:

## log2(cpm(quant(cbcb(data))))

## It will save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep libsizes in mind
##  when invoking limma.  The appropriate libsize is non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Not correcting the count-data for batch effects.  If batch is
##  included in EdgerR/limma's model, then this is probably wise; but in extreme
##  batch effects this is a good parameter to play with.

## Step 1: performing count filter with option: cbcb

## Removing 0 low-count genes (133649 remaining).

## Step 2: normalizing the data with quant.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## transform_counts: Found 209094 values equal to 0, adding 1 to the matrix.

## The method is: raw.

## Step 5: not doing batch correction.

new_zymo_norm <- set_expt_conditions(new_zymo_norm, fact="phenotypiccharacteristics")
zymo_heat <- plot_disheat(new_zymo_norm)

zymo_subset <- snp_subset_genes(lp_expt, new_snps,
                                genes=c("LPAL13_120010900", "LPAL13_340013000", "LPAL13_000054100",
                                        "LPAL13_140006100", "LPAL13_180018500", "LPAL13_320022300"))

## Warning in .Seqinfo.mergexy(x, y): Each of the 2 combined objects has sequence levels not in the other:
##   - in 'x': LPAL13-SCAF000002, LPAL13-SCAF000003, LPAL13-SCAF000005, LPAL13-SCAF000009, LPAL13-SCAF000014, LPAL13-SCAF000015, LPAL13-SCAF000018, LPAL13-SCAF000019, LPAL13-SCAF000020, LPAL13-SCAF000022, LPAL13-SCAF000023, LPAL13-SCAF000026, LPAL13-SCAF000029, LPAL13-SCAF000032, LPAL13-SCAF000035, LPAL13-SCAF000036, LPAL13-SCAF000037, LPAL13-SCAF000038, LPAL13-SCAF000042, LPAL13-SCAF000043, LPAL13-SCAF000045, LPAL13-SCAF000047, LPAL13-SCAF000049, LPAL13-SCAF000050, LPAL13-SCAF000054, LPAL13-SCAF000056, LPAL13-SCAF000057, LPAL13-SCAF000058, LPAL13-SCAF000066, LPAL13-SCAF000067, LPAL13-SCAF000070, LPAL13-SCAF000073, LPAL13-SCAF000082, LPAL13-SCAF000083, LPAL13-SCAF000085, LPAL13-SCAF000086, LPAL13-SCAF000088, LPAL13-SCAF000090, LPAL13-SCAF000091, LPAL13-SCAF000092, LPAL13-SCAF000095, LPAL13-SCAF000098, LPAL13-SCAF000101, LPAL13-SCAF000103, LPAL13-SCAF000106, LPAL13-SCAF000109, LPAL13-SCAF000111, LPAL13-SCAF000112, LPAL13-SCAF000113, LPAL13-SCAF000118, LPAL13-SCAF000125, LPAL13-SCAF000126, LPAL13-SCAF000138, LPAL13-SCAF000139, LPAL13-SCAF000140, LPAL13-SCAF000141, LPAL13-SCAF000144, LPAL13-SCAF000145, LPAL13-SCAF000147, LPAL13-SCAF000148, LPAL13-SCAF000150, LPAL13-SCAF000151, LPAL13-SCAF000152, LPAL13-SCAF000154, LPAL13-SCAF000155, LPAL13-SCAF000156, LPAL13-SCAF000158, LPAL13-SCAF000159, LPAL13-SCAF000160, LPAL13-SCAF000161, LPAL13-SCAF000163, LPAL13-SCAF000164, LPAL13-SCAF000167, LPAL13-SCAF000168, LPAL13-SCAF000169, LPAL13-SCAF000170, LPAL13-SCAF000175, LPAL13-SCAF000177, LPAL13-SCAF000178, LPAL13-SCAF000179, LPAL13-SCAF000180, LPAL13-SCAF000184, LPAL13-SCAF000185, LPAL13-SCAF000189, LPAL13-SCAF000190, LPAL13-SCAF000192, LPAL13-SCAF000195, LPAL13-SCAF000196, LPAL13-SCAF000198, LPAL13-SCAF000199, LPAL13-SCAF000204, LPAL13-SCAF000207, LPAL13-SCAF000210, LPAL13-SCAF000212, LPAL13-SCAF000214, LPAL13-SCAF000215, LPAL13-SCAF000216, LPAL13-SCAF000218, LPAL13-SCAF000219, LPAL13-SCAF000221, LPAL13-SCAF000223, LPAL13-SCAF000224, LPAL13-SCAF000226, LPAL13-SCAF000228, LPAL13-SCAF000234, LPAL13-SCAF000236, LPAL13-SCAF000240, LPAL13-SCAF000241, LPAL13-SCAF000242, LPAL13-SCAF000243, LPAL13-SCAF000244, LPAL13-SCAF000246, LPAL13-SCAF000247, LPAL13-SCAF000251, LPAL13-SCAF000252, LPAL13-SCAF000254, LPAL13-SCAF000255, LPAL13-SCAF000257, LPAL13-SCAF000258, LPAL13-SCAF000260, LPAL13-SCAF000262, LPAL13-SCAF000263, LPAL13-SCAF000269, LPAL13-SCAF000270, LPAL13-SCAF000272, LPAL13-SCAF000274, LPAL13-SCAF000275, LPAL13-SCAF000276, LPAL13-SCAF000277, LPAL13-SCAF000278, LPAL13-SCAF000279, LPAL13-SCAF000280, LPAL13-SCAF000282, LPAL13-SCAF000283, LPAL13-SCAF000284, LPAL13-SCAF000289, LPAL13-SCAF000290, LPAL13-SCAF000293, LPAL13-SCAF000294, LPAL13-SCAF000297, LPAL13-SCAF000298, LPAL13-SCAF000299, LPAL13-SCAF000304, LPAL13-SCAF000306, LPAL13-SCAF000307, LPAL13-SCAF000311, LPAL13-SCAF000312, LPAL13-SCAF000315, LPAL13-SCAF000323, LPAL13-SCAF000325, LPAL13-SCAF000327, LPAL13-SCAF000331, LPAL13-SCAF000332, LPAL13-SCAF000333, LPAL13-SCAF000334, LPAL13-SCAF000336, LPAL13-SCAF000341, LPAL13-SCAF000342, LPAL13-SCAF000344, LPAL13-SCAF000348, LPAL13-SCAF000349, LPAL13-SCAF000350, LPAL13-SCAF000351, LPAL13-SCAF000352, LPAL13-SCAF000355, LPAL13-SCAF000357, LPAL13-SCAF000359, LPAL13-SCAF000360, LPAL13-SCAF000361, LPAL13-SCAF000362, LPAL13-SCAF000365, LPAL13-SCAF000366, LPAL13-SCAF000371, LPAL13-SCAF000375, LPAL13-SCAF000376, LPAL13-SCAF000377, LPAL13-SCAF000378, LPAL13-SCAF000379, LPAL13-SCAF000380, LPAL13-SCAF000381, LPAL13-SCAF000382, LPAL13-SCAF000383, LPAL13-SCAF000385, LPAL13-SCAF000386, LPAL13-SCAF000387, LPAL13-SCAF000390, LPAL13-SCAF000392, LPAL13-SCAF000393, LPAL13-SCAF000394, LPAL13-SCAF000395, LPAL13-SCAF000396, LPAL13-SCAF000398, LPAL13-SCAF000399, LPAL13-SCAF000406, LPAL13-SCAF000407, LPAL13-SCAF000408, LPAL13-SCAF000409, LPAL13-SCAF000411, LPAL13-SCAF000412, LPAL13-SCAF000413, LPAL13-SCAF000416, LPAL13-SCAF000418, LPAL13-SCAF000422, LPAL13-SCAF000423, LPAL13-SCAF000425, LPAL13-SCAF000427, LPAL13-SCAF000431, LPAL13-SCAF000435, LPAL13-SCAF000438, LPAL13-SCAF000439, LPAL13-SCAF000441, LPAL13-SCAF000442, LPAL13-SCAF000443, LPAL13-SCAF000444, LPAL13-SCAF000445, LPAL13-SCAF000449, LPAL13-SCAF000450, LPAL13-SCAF000451, LPAL13-SCAF000452, LPAL13-SCAF000454, LPAL13-SCAF000455, LPAL13-SCAF000457, LPAL13-SCAF000458, LPAL13-SCAF000462, LPAL13-SCAF000464, LPAL13-SCAF000466, LPAL13-SCAF000467, LPAL13-SCAF000472, LPAL13-SCAF000473, LPAL13-SCAF000474, LPAL13-SCAF000475, LPAL13-SCAF000476, LPAL13-SCAF000479, LPAL13-SCAF000481, LPAL13-SCAF000482, LPAL13-SCAF000485, LPAL13-SCAF000489, LPAL13-SCAF000494, LPAL13-SCAF000497, LPAL13-SCAF000498, LPAL13-SCAF000499, LPAL13-SCAF000504, LPAL13-SCAF000506, LPAL13-SCAF000510, LPAL13-SCAF000513, LPAL13-SCAF000514, LPAL13-SCAF000516, LPAL13-SCAF000517, LPAL13-SCAF000518, LPAL13-SCAF000519, LPAL13-SCAF000520, LPAL13-SCAF000521, LPAL13-SCAF000524, LPAL13-SCAF000525, LPAL13-SCAF000526, LPAL13-SCAF000530, LPAL13-SCAF000531, LPAL13-SCAF000545, LPAL13-SCAF000546, LPAL13-SCAF000550, LPAL13-SCAF000557, LPAL13-SCAF000565, LPAL13-SCAF000571, LPAL13-SCAF000581, LPAL13-SCAF000584, LPAL13-SCAF000589, LPAL13-SCAF000592, LPAL13-SCAF000595, LPAL13-SCAF000596, LPAL13-SCAF000597, LPAL13-SCAF000602, LPAL13-SCAF000604, LPAL13-SCAF000606, LPAL13-SCAF000608, LPAL13-SCAF000609, LPAL13-SCAF000612, LPAL13-SCAF000613, LPAL13-SCAF000615, LPAL13-SCAF000620, LPAL13-SCAF000621, LPAL13-SCAF000623, LPAL13-SCAF000624, LPAL13-SCAF000629, LPAL13-SCAF000631, LPAL13-SCAF000632, LPAL13-SCAF000633, LPAL13-SCAF000634, LPAL13-SCAF000635, LPAL13-SCAF000640, LPAL13-SCAF000642, LPAL13-SCAF000647, LPAL13-SCAF000648, LPAL13-SCAF000657, LPAL13-SCAF000658, LPAL13-SCAF000660, LPAL13-SCAF000662, LPAL13-SCAF000663, LPAL13-SCAF000664, LPAL13-SCAF000665, LPAL13-SCAF000667, LPAL13-SCAF000669, LPAL13-SCAF000670, LPAL13-SCAF000675, LPAL13-SCAF000676, LPAL13-SCAF000677, LPAL13-SCAF000678, LPAL13-SCAF000683, LPAL13-SCAF000684, LPAL13-SCAF000685, LPAL13-SCAF000686, LPAL13-SCAF000687, LPAL13-SCAF000689, LPAL13-SCAF000690, LPAL13-SCAF000691, LPAL13-SCAF000693, LPAL13-SCAF000694, LPAL13-SCAF000699, LPAL13-SCAF000701, LPAL13-SCAF000702, LPAL13-SCAF000703, LPAL13-SCAF000705, LPAL13-SCAF000706, LPAL13-SCAF000708, LPAL13-SCAF000709, LPAL13-SCAF000710, LPAL13-SCAF000712, LPAL13-SCAF000715, LPAL13-SCAF000718, LPAL13-SCAF000721, LPAL13-SCAF000725, LPAL13-SCAF000728, LPAL13-SCAF000729, LPAL13-SCAF000730, LPAL13-SCAF000731, LPAL13-SCAF000733, LPAL13-SCAF000736, LPAL13-SCAF000739, LPAL13-SCAF000740, LPAL13-SCAF000741, LPAL13-SCAF000742, LPAL13-SCAF000743, LPAL13-SCAF000745, LPAL13-SCAF000746, LPAL13-SCAF000747, LPAL13-SCAF000749, LPAL13-SCAF000750, LPAL13-SCAF000751, LPAL13-SCAF000752, LPAL13-SCAF000753, LPAL13-SCAF000754, LPAL13-SCAF000755, LPAL13-SCAF000756, LPAL13-SCAF000757, LPAL13-SCAF000758, LPAL13-SCAF000759, LPAL13-SCAF000763, LPAL13-SCAF000764, LPAL13-SCAF000765, LPAL13-SCAF000766, LPAL13-SCAF000767, LPAL13-SCAF000768, LPAL13-SCAF000769, LPAL13-SCAF000770, LPAL13-SCAF000771, LPAL13-SCAF000773, LPAL13-SCAF000774, LPAL13-SCAF000775, LPAL13-SCAF000776, LPAL13-SCAF000777, LPAL13-SCAF000778, LPAL13-SCAF000784, LPAL13-SCAF000785, LPAL13-SCAF000786, LPAL13-SCAF000787, LPAL13-SCAF000790, LPAL13-SCAF000794, LPAL13-SCAF000795, LPAL13-SCAF000796, LPAL13-SCAF000803, LPAL13-SCAF000807, LPAL13-SCAF000809, LPAL13-SCAF000810, LPAL13-SCAF000811, LPAL13-SCAF000814, LPAL13-SCAF000817, LPAL13-SCAF000819, LPAL13-SCAF000820
##   - in 'y': LPAL13-SCAF000074, LPAL13-SCAF000136, LPAL13-SCAF000230, LPAL13-SCAF000286, LPAL13-SCAF000328, LPAL13-SCAF000539, LPAL13-SCAF000653
##   Make sure to always combine/compare objects based on the same reference
##   genome (use suppressWarnings() to suppress this warning).

## Before removal, there were 133649 entries.

## Now there are 17 entries.

## Percent kept: 0.037, 0.016, 0.000, 0.042, 0.081, 0.053, 0.046, 0.026, 0.000, 0.028, 0.026, 0.025, 0.000, 0.029, 0.020, 0.018, 0.027

## Percent removed: 99.963, 99.984, 100.000, 99.958, 99.919, 99.947, 99.954, 99.974, 100.000, 99.972, 99.974, 99.975, 100.000, 99.971, 99.980, 99.982, 99.973

zymo_subset <- set_expt_conditions(zymo_subset, fact="phenotypiccharacteristics")
## zymo_heat <- plot_sample_heatmap(zymo_subset, row_label=rownames(exprs(snp_subset)))

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[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"))(170)
map1 <- annHeatmap2(
  correlations,
  dendrogram=mydendro,
  annotation=myannot,
  cluster=myclust,
  labels=mylabs,
  col=hmcols)

## Warning in breakColors(breaks, col): more colors than classes: ignoring 3 last
## colors

plot(map1)

8 Using Variant profiles to make guesses about strains and chronic/self-healing

The following uses the same information to make some guesses about the strains used in the new samples.

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))

mydendro <- list(
  "clustfun" = hclust,
  "lwd" = 2.0)
col_data <- as.data.frame(des[, c("condition")])
row_data <- as.data.frame(des[, c("strain")])
colnames(col_data) <- c("condition")
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"))(170)
map1 <- annHeatmap2(
  correlations,
  dendrogram=mydendro,
  annotation=myannot,
  cluster=myclust,
  labels=mylabs,
  col=hmcols)

## Warning in breakColors(breaks, col): more colors than classes: ignoring 3 last
## colors

plot(map1)

pheno <- subset_expt(lp_expt, subset="condition=='z2.2'|condition=='z2.3'")

## Using a subset expression.

## There were 17, now there are 14 samples.

pheno_snps <- sm(count_expt_snps(pheno, annot_column="bcftable"))

xref_prop <- table(pheno_snps$conditions)
pheno_snps$conditions

##  [1] "z2.3" "z2.2" "z2.2" "z2.2" "z2.3" "z2.2" "z2.3" "z2.3" "z2.2" "z2.2"
## [11] "z2.3" "z2.2" "z2.2" "z2.3"

idx_tbl <- exprs(pheno_snps) > 5
new_tbl <- data.frame(row.names=rownames(exprs(pheno_snps)))
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
}
new_tbl[["ratio"]] <- (new_tbl[["z2.2"]] - new_tbl[["z2.3"]])
keepers <- grepl(x=rownames(new_tbl), pattern="LpaL13")
new_tbl <- new_tbl[keepers, ]
new_tbl[["SNP"]] <- rownames(new_tbl)
new_tbl[["Chromosome"]] <- gsub(x=new_tbl[["SNP"]], pattern="chr_(.*)_pos_.*", replacement="\\1")
new_tbl[["Position"]] <- gsub(x=new_tbl[["SNP"]], pattern=".*_pos_(\\d+)_.*", replacement="\\1")
new_tbl <- new_tbl[, c("SNP", "Chromosome", "Position", "ratio")]
library(CMplot)

## Much appreciate for using CMplot.

## Full description, Bug report, Suggestion and the latest codes:

## https://github.com/YinLiLin/CMplot

CMplot(new_tbl)

##  SNP-Density Plotting.

##  Circular-Manhattan Plotting ratio.
##  Rectangular-Manhattan Plotting ratio.

##  QQ Plotting ratio.

##  Plots are stored in: /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019

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 052640a0b091e9a740e487957e087f265d0c74b5

## This is hpgltools commit: Thu Dec 3 13:42:02 2020 -0500: 052640a0b091e9a740e487957e087f265d0c74b5

## Saving to 01_annotation_v20201022.rda.xz

tmp <- loadme(filename=savefile)