S. agalactiae 202105: Poking at some TNSeq.

1 Grab annotation data for Streptococcus agalactiae.

1.1 Strain cjb111

I am going to use the gff file which I extracted from the NCBI genbank file which was used with htseq. In addition, microbesonline has some annotations for CJB111.

1.1.1 Microbesonline

Download the microbesonline annotations. In a moment we should be able to merge these with the gff/genbank data.

cjb_microbes <- load_microbesonline_annotations(species="CJB111")

## Found 1 entry.

## Streptococcus agalactiae CJB111Firmicutesyes2007-05-08noNANA2208342617

## The species being downloaded is: Streptococcus agalactiae CJB111

## Downloading: http://www.microbesonline.org/cgi-bin/genomeInfo.cgi?tId=342617;export=tab

cjb_microbes <- as.data.frame(cjb_microbes)
cjb_microbes[["sysName"]] <- make.names(cjb_microbes[["sysName"]], unique=TRUE)

1.1.2 GFF/Genbank

I downloaded Lindsey’s annotations and genbank data, here are the resulting annotations.

cjb_gff <- load_gff_annotations("reference/sagalactiae_cjb111.gff")

## Trying attempt: rtracklayer::import.gff3(gff, sequenceRegionsAsSeqinfo = TRUE)

## Had a successful gff import with rtracklayer::import.gff3(gff, sequenceRegionsAsSeqinfo = TRUE)

## Returning a df with 42 columns and 4347 rows.

rownames(cjb_gff) <- make.names(cjb_gff[["locus_tag"]], unique=TRUE)
## I am going to only pay attention to the first annotation for each locus tag from microbesonline.
## I should ask Lindsey if they would like me to merge these
cjb_annot <- cjb_gff
rownames(cjb_annot) <- paste0("gene-", make.names(cjb_annot[["Name"]], unique=TRUE))
##cjb_annot <- merge(cjb_gff, cjb_microbes, by.x="locus_tag", by.y="sysName")
## rownames(cjb_annot) <- make.names(cjb_annot[["locus_tag"]], unique=TRUE)

cjb_gb <- load_genbank_annotations(accession="CP063198")

## Loading required namespace: rentrez

## Done Parsing raw GenBank file text. [ 15.997 seconds ]

## 2021-06-02 13:25:00 Starting creation of gene GRanges

## 2021-06-02 13:25:01 Starting creation of CDS GRanges

## Translation product seems to be missing for 51 of 2010 CDS annotations. Setting to ''

## 2021-06-02 13:25:05 Starting creation of exon GRanges

## No exons read from genbank file. Assuming sections of CDS are full exons

## 2021-06-02 13:25:05 Starting creation of variant VRanges

## 2021-06-02 13:25:05 Starting creation of transcript GRanges

## No transcript features (mRNA) found, using spans of CDSs

## 2021-06-02 13:25:05 Starting creation of misc feature GRanges

## Warning in fill_stack_df(feats[!typs %in% c("gene", "exon", "CDS",
## "variation", : Got unexpected multi-value field(s) [ inference ]. The resulting
## column(s) will be of class CharacterList, rather than vector(s). Please contact
## the maintainer if multi-valuedness is expected/meaningful for the listed
## field(s).

## 2021-06-02 13:25:05 - Done creating GenBankRecord object [ 5.66 seconds ]

cjb_gb <- as.data.frame(cjb_gb[["genes"]])
rownames(cjb_gb) <- paste0("gene-", cjb_gb[["gene_id"]])

1.1.3 Pull them together

I used ggsearch to compare the microbesonline transcriptome to the genbank. The resulting tables should provide a bridge for using the microbesonline annotations.

I set a relatively strict cutoff and asked for only the best match.

ggsearch_result <- readr::read_tsv("reference/compare_cjb_annotations/search_parsed.txt")

## 
## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## cols(
##   QueryName = col_character(),
##   QueryLength = col_double(),
##   LibID = col_character(),
##   LibHitLength = col_double(),
##   Score = col_double(),
##   E = col_double(),
##   Identity = col_double(),
##   HitMatches = col_double(),
##   QueryStart = col_double(),
##   QueryEnd = col_double(),
##   LibStart = col_double(),
##   LibEnd = col_double(),
##   LStrand = col_double()
## )

genbank_ids <- readr::read_table2("reference/compare_cjb_annotations/gb_ids.txt", col_names = c("id", "sam"))

## 
## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## cols(
##   id = col_character(),
##   sam = col_character()
## )

genbank_bridge <- ggsearch_result[, c("QueryName", "LibID")]
colnames(genbank_bridge) <- c("id", "genbank")
genbank_bridge <- merge(genbank_bridge, genbank_ids, by = "id")
all_annot <- merge(cjb_gff, genbank_bridge, by.x = "row.names", by.y = "genbank")
all_annot <- merge(all_annot, cjb_microbes, by.x = "sam", by.y = "name")
rownames(all_annot) <- make.names(all_annot[["ID"]], unique = TRUE)
all_annot[["Row.names"]] <- NULL
rownames(all_annot) <- gsub(pattern = "\\.", replacement = "-", x = rownames(all_annot))

## Clean up the table a little
all_annot[["start.y"]] <- NULL
all_annot[["strand.y"]] <- NULL
start_idx <- colnames(all_annot) == "start.x"
colnames(all_annot)[start_idx] <- "start"
strand_idx <- colnames(all_annot) == "strand.x"
colnames(all_annot)[strand_idx] <- "strand"

2 Create Expressionsets

The following block merges the various counts, annotations, and experimental metadata.

cjb_expt <- create_expt(metadata="sample_sheets/all_samples_cjb111.xlsx",
                        gene_info=all_annot, file_column="btv1m1sagalactiaecjb111")

## Reading the sample metadata.

## The sample definitions comprises: 25 rows(samples) and 12 columns(metadata fields).

## Matched 1842 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 2064 rows and 25 columns.

cjb_written <- write_expt(cjb_expt, batch="raw",
                           excel=glue::glue("excel/{rundate}-cjb_counts-v{ver}.xlsx"))

## Writing the first sheet, containing a legend and some summary data.

## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

## Warning: ggrepel: 21 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 21 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 12 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 12 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

## The following object is masked from 'package:S4Vectors':
## 
##     expand

## 
## Total:13 s

## Warning: ggrepel: 12 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 12 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'

## Warning: ggrepel: 12 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: ggrepel: 12 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## 
## Total:12 s

3 A Few diagnostic plots

cjb_written[["legend_plot"]]

cjb_written[["raw_libsize"]]

cjb_written[["raw_density"]]

## awesome

cjb_written[["norm_disheat"]]

cjb_written[["norm_corheat"]]

plot_topn(cjb_expt, num=40)$plot

## `geom_smooth()` using formula 'y ~ x'

## So, sample d3s3 has ~ 5 genes which dominate it.

cjb_n <- normalize_expt(cjb_expt, filter=TRUE, norm="quant", convert="cpm", transform="log2")

## Removing 117 low-count genes (1947 remaining).

pp(file = "images/all_samples_pca_nobatch.png", image = plot_pca(cjb_n)$plot)

cjb_nb <- normalize_expt(cjb_expt, filter=TRUE, convert="cpm",
                         transform="log2", batch="svaseq")

## Removing 117 low-count genes (1947 remaining).

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

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

## Setting 3683 low elements to zero.

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

plot_pca(cjb_nb)$plot

## Warning: ggrepel: 1 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

cjb_filt <- normalize_expt(cjb_expt, filter="simple")

## Removing 73 low-count genes (1991 remaining).

4 Check tnseq saturation

The essentiality programs result in a series of wig/tsv files which we may use to get a sense of how well they Bayesian methods worked with the data. This is especially important for the in-vivo data because they are very different from replicate to replicate and it is more difficult to get saturation.

These plots may be generated for either the individual samples, the combined samples, or both.

4.0.1 Saturation of the in vivo control

saturation_tn <- tnseq_saturation(
    "preprocessing/combined_tn/outputs/essentiality_sagalactiae_cjb111/consolidated-v1M1.wig")
saturation_tn[["plot"]]

## Warning: Removed 59310 rows containing non-finite values (stat_bin).

## Warning: Removed 59310 rows containing non-finite values (stat_density).

## Warning: Removed 2 rows containing missing values (geom_bar).

ess_tn <- plot_essentiality(
  "preprocessing/combined_tn/outputs/essentiality_sagalactiae_cjb111/mh_ess-consolidated-v1M1_gene_tas_m1.csv")
ess_tn[["zbar"]]

4.0.2 Saturation of the in vitro control

Note, that because the coverage is so much higher, I set -m to 8 here. This may be too high.

saturation_vitro_control <- tnseq_saturation(
    "preprocessing/combined_control/outputs/essentiality_sagalactiae_cjb111/consolidated-v1M1.wig")
saturation_vitro_control[["plot"]]

## Warning: Removed 41392 rows containing non-finite values (stat_bin).

## Warning: Removed 41392 rows containing non-finite values (stat_density).

## Warning: Removed 2 rows containing missing values (geom_bar).

ess_vitro_control <- plot_essentiality(
  "preprocessing/combined_control/outputs/essentiality_sagalactiae_cjb111/mh_ess-consolidated-v1M1_gene_tas_m8.csv")
ess_vitro_control[["zbar"]]

4.0.3 Saturation of the in vitro low calprotectin

saturation_vitro_low <- tnseq_saturation(
    "preprocessing/combined_low/outputs/essentiality_sagalactiae_cjb111/consolidated-v1M1.wig")
saturation_vitro_low[["plot"]]

## Warning: Removed 48470 rows containing non-finite values (stat_bin).

## Warning: Removed 48470 rows containing non-finite values (stat_density).

## Warning: Removed 2 rows containing missing values (geom_bar).

ess_vitro_low <- plot_essentiality(
  "preprocessing/combined_low/outputs/essentiality_sagalactiae_cjb111/mh_ess-consolidated-v1M1_gene_tas_m8.csv")
ess_vitro_low[["zbar"]]

4.0.4 Saturation of the in vitro high

saturation_vitro_high <- tnseq_saturation(
    "preprocessing/combined_high/outputs/essentiality_sagalactiae_cjb111/consolidated-v1M1.wig")
saturation_vitro_high[["plot"]]

## Warning: Removed 41526 rows containing non-finite values (stat_bin).

## Warning: Removed 41526 rows containing non-finite values (stat_density).

## Warning: Removed 2 rows containing missing values (geom_bar).

ess_vitro_high <- plot_essentiality(
  "preprocessing/combined_high/outputs/essentiality_sagalactiae_cjb111/mh_ess-consolidated-v1M1_gene_tas_m8.csv")

## Warning in lmrob.S(x, y, control = control): S refinements did not converge (to
## refine.tol=1e-07) in 200 (= k.max) steps

## Warning in lmrob.S(x, y, control = control): S refinements did not converge (to
## refine.tol=1e-07) in 200 (= k.max) steps

ess_vitro_high[["zbar"]]

4.0.5 D1 saturation

saturation_d1 <- tnseq_saturation(
    "preprocessing/combined_d1/outputs/essentiality_sagalactiae_cjb111/consolidated-v1M1.wig")
saturation_d1[["plot"]]

## Warning: Removed 18720 rows containing non-finite values (stat_bin).

## Warning: Removed 18720 rows containing non-finite values (stat_density).

## Warning: Removed 2 rows containing missing values (geom_bar).

ess_d1 <- plot_essentiality(
  "preprocessing/combined_d1/outputs/essentiality_sagalactiae_cjb111/mh_ess-consolidated-v1M1_gene_tas_m1.csv")

## Warning in lmrob.S(x, y, control = control): S refinements did not converge (to
## refine.tol=1e-07) in 200 (= k.max) steps

## Warning in lmrob.S(x, y, control = control): S refinements did not converge (to
## refine.tol=1e-07) in 200 (= k.max) steps

ess_d1[["zbar"]]

4.0.6 D3 saturation

saturation_d3 <- tnseq_saturation(
    "preprocessing/combined_d3/outputs/essentiality_sagalactiae_cjb111/consolidated-v1M1.wig")
saturation_d3[["plot"]]

## Warning: Removed 7013 rows containing non-finite values (stat_bin).

## Warning: Removed 7013 rows containing non-finite values (stat_density).

## Warning: Removed 2 rows containing missing values (geom_bar).

ess_d3 <- plot_essentiality(
  "preprocessing/combined_d3/outputs/essentiality_sagalactiae_cjb111/mh_ess-consolidated-v1M1_gene_tas_m1.csv")
ess_d3[["zbar"]]

4.0.7 CX saturation

saturation_cx <- tnseq_saturation(
    "preprocessing/combined_cx/outputs/essentiality_sagalactiae_cjb111/consolidated-v1M1.wig")
saturation_cx[["plot"]]

## Warning: Removed 7045 rows containing non-finite values (stat_bin).

## Warning: Removed 7045 rows containing non-finite values (stat_density).

## Warning: Removed 2 rows containing missing values (geom_bar).

ess_cx <- plot_essentiality(
  "preprocessing/combined_cx/outputs/essentiality_sagalactiae_cjb111/mh_ess-consolidated-v1M1_gene_tas_m1.csv")
ess_cx[["zbar"]]

saturation_ut <- tnseq_saturation(
    "preprocessing/combined_ut/outputs/essentiality_sagalactiae_cjb111/consolidated-v1M1.wig")
saturation_ut[["plot"]]

## Warning: Removed 2710 rows containing non-finite values (stat_bin).

## Warning: Removed 2710 rows containing non-finite values (stat_density).

## Warning: Removed 2 rows containing missing values (geom_bar).

ess_ut <- plot_essentiality(
  "preprocessing/combined_ut/outputs/essentiality_sagalactiae_cjb111/mh_ess-consolidated-v1M1_gene_tas_m1.csv")
ess_ut[["zbar"]]

saturation_vg <- tnseq_saturation(
    "preprocessing/combined_vg/outputs/essentiality_sagalactiae_cjb111/consolidated-v1M1.wig")
saturation_vg[["plot"]]

## Warning: Removed 6713 rows containing non-finite values (stat_bin).

## Warning: Removed 6713 rows containing non-finite values (stat_density).

## Warning: Removed 2 rows containing missing values (geom_bar).

ess_vg <- plot_essentiality(
  "preprocessing/combined_vg/outputs/essentiality_sagalactiae_cjb111/mh_ess-consolidated-v1M1_gene_tas_m1.csv")
ess_vg[["zbar"]]

These plots/estimations of essentiality are suspiciously similar. Something to me feels wrong, but I am not sure what it might be.

5 Changed genes

For differential expression, I am going to assume until I hear otherwise, that my batch assignments are not correct and that the 1,2,3 assignments of the sample names do not actually delineate separate batches. Though, if they do delineate separate batches, it might be taken as a (very)small degree of evidence that 04 and 09 were switched.

5.1 Strain cjb111

combined_essentiality <- list(
  "vitro_control" = "preprocessing/combined_control/outputs/essentiality_sagalactiae_cjb111/mh_ess-consolidated-v1M1_gene_tas_m8.csv",
  "vitro_low" = "preprocessing/combined_low/outputs/essentiality_sagalactiae_cjb111/mh_ess-consolidated-v1M1_gene_tas_m8.csv",
  "vitro_high" = "preprocessing/combined_high/outputs/essentiality_sagalactiae_cjb111/mh_ess-consolidated-v1M1_gene_tas_m8.csv",
  "vivo_cx" = "preprocessing/combined_cx/outputs/essentiality_sagalactiae_cjb111/mh_ess-consolidated-v1M1_gene_tas_m2.csv",
  "vivo_d1" = "preprocessing/combined_d1/outputs/essentiality_sagalactiae_cjb111/mh_ess-consolidated-v1M1_gene_tas_m2.csv",
  "vivo_d3" = "preprocessing/combined_d3/outputs/essentiality_sagalactiae_cjb111/mh_ess-consolidated-v1M1_gene_tas_m2.csv",
  "vivo_ut" = "preprocessing/combined_ut/outputs/essentiality_sagalactiae_cjb111/mh_ess-consolidated-v1M1_gene_tas_m2.csv",
  "vivo_vg" = "preprocessing/combined_vg/outputs/essentiality_sagalactiae_cjb111/mh_ess-consolidated-v1M1_gene_tas_m2.csv",  
  "vivo_tn" = "preprocessing/combined_tn/outputs/essentiality_sagalactiae_cjb111/mh_ess-consolidated-v1M1_gene_tas_m2.csv")
ess_table <- data.frame()
for (f in 1:length(combined_essentiality)) {
  name <- names(combined_essentiality)[f]
  column_names <- c("orf", "k", "n", "r", "s", "zbar", "call")
  names <- paste0(name, "_", column_names)
  r <- readr::read_tsv(combined_essentiality[[f]], comment="#", col_names=names)
  colnames(r)[1] <- "orf"
  if (f == 1) {
    ess_table <- r
  } else {
    ess_table <- merge(ess_table, r, by="orf")
  }
}

## 
## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## cols(
##   vitro_control_orf = col_character(),
##   vitro_control_k = col_double(),
##   vitro_control_n = col_double(),
##   vitro_control_r = col_double(),
##   vitro_control_s = col_double(),
##   vitro_control_zbar = col_double(),
##   vitro_control_call = col_character()
## )

## 
## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## cols(
##   vitro_low_orf = col_character(),
##   vitro_low_k = col_double(),
##   vitro_low_n = col_double(),
##   vitro_low_r = col_double(),
##   vitro_low_s = col_double(),
##   vitro_low_zbar = col_double(),
##   vitro_low_call = col_character()
## )

## 
## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## cols(
##   vitro_high_orf = col_character(),
##   vitro_high_k = col_double(),
##   vitro_high_n = col_double(),
##   vitro_high_r = col_double(),
##   vitro_high_s = col_double(),
##   vitro_high_zbar = col_double(),
##   vitro_high_call = col_character()
## )

## 
## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## cols(
##   vivo_cx_orf = col_character(),
##   vivo_cx_k = col_double(),
##   vivo_cx_n = col_double(),
##   vivo_cx_r = col_double(),
##   vivo_cx_s = col_double(),
##   vivo_cx_zbar = col_double(),
##   vivo_cx_call = col_character()
## )

## 
## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## cols(
##   vivo_d1_orf = col_character(),
##   vivo_d1_k = col_double(),
##   vivo_d1_n = col_double(),
##   vivo_d1_r = col_double(),
##   vivo_d1_s = col_double(),
##   vivo_d1_zbar = col_double(),
##   vivo_d1_call = col_character()
## )

## 
## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## cols(
##   vivo_d3_orf = col_character(),
##   vivo_d3_k = col_double(),
##   vivo_d3_n = col_double(),
##   vivo_d3_r = col_double(),
##   vivo_d3_s = col_double(),
##   vivo_d3_zbar = col_double(),
##   vivo_d3_call = col_character()
## )

## 
## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## cols(
##   vivo_ut_orf = col_character(),
##   vivo_ut_k = col_double(),
##   vivo_ut_n = col_double(),
##   vivo_ut_r = col_double(),
##   vivo_ut_s = col_double(),
##   vivo_ut_zbar = col_double(),
##   vivo_ut_call = col_character()
## )

## 
## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## cols(
##   vivo_vg_orf = col_character(),
##   vivo_vg_k = col_double(),
##   vivo_vg_n = col_double(),
##   vivo_vg_r = col_double(),
##   vivo_vg_s = col_double(),
##   vivo_vg_zbar = col_double(),
##   vivo_vg_call = col_character()
## )

## 
## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## cols(
##   vivo_tn_orf = col_character(),
##   vivo_tn_k = col_double(),
##   vivo_tn_n = col_double(),
##   vivo_tn_r = col_double(),
##   vivo_tn_s = col_double(),
##   vivo_tn_zbar = col_double(),
##   vivo_tn_call = col_character()
## )

rownames(ess_table) <- gsub(x=ess_table[["orf"]], pattern="^cds_", replacement="")
ess_table[["orf"]] <- NULL

written <- write_xlsx(data = ess_table, excel = "excel/essentiality_results_one_table.xlsx")

6 Perform DEish analysis and include the essentiality results

tn_de <- all_pairwise(cjb_filt, model_batch="svaseq")

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

## Plotting a PCA before surrogate/batch inclusion.

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

## Performing a test normalization with: raw

## Removing 44 low-count genes (1947 remaining).

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

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

## Setting 3683 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

keepers <- list(
  "ut" = c("ut", "control"),
  "d1" = c("d1", "control"),
  "d3" = c("d3", "control"),
  "cx" = c("cx", "control"),
  "vg" = c("vg", "control"),
  "low" = c("cal_low", "control"),
  "high" = c("cal_high", "control"))
tn_tables <- combine_de_tables(tn_de, keepers=keepers, extra_annot = ess_table,
                               excel=glue::glue("excel/vivo_fitness-v{ver}.xlsx"))

## Deleting the file excel/vivo_fitness-v202105.xlsx before writing the tables.

7 Circos

This worksheet was revived primarily because of an email query from Lindsey about some circos plots:

For the circos plot, I’m thinking nonessential (black), unknown (gray), and essential (red). For samples I think we would need to include mRPMI, high calprotectin, and low calprotectin from the in vitro data and each of these samples from the in vivo (D1 swabs, D3 swabs, vagina, cervix, uterus)? It sounds like it could be very complicated to me, do you think it will be possible to see anything clearly?

The following block will attempt to answer this question… I am going to use our in-vitro circos plot as a template and attempt to add to it the above…

combined_expt <- create_expt("sample_sheets/sagalactiae_combined_samples.xlsx",
                             gene_info=all_annot)

## Reading the sample metadata.

## The sample definitions comprises: 9 rows(samples) and 11 columns(metadata fields).

## Matched 1842 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 2064 rows and 9 columns.

combined_written <- write_expt(combined_expt, excel = "excel/combined_expt.xlsx")

## Deleting the file excel/combined_expt.xlsx before writing the tables.

## Writing the first sheet, containing a legend and some summary data.

## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 9 rows containing missing values (geom_text).

## Warning: Removed 9 rows containing missing values (geom_text).

## Warning in doTryCatch(return(expr), name, parentenv, handler): display list
## redraw incomplete

## varpart sees only 1 batch, adjusting the model accordingly.

## Error in if (ncol(exprObj) != nrow(data)) { : argument is of length zero

## Retrying with only condition in the model.

## Error in checkModelStatus(fit, showWarnings = showWarnings, colinearityCutoff = colinearityCutoff) : 
##   Colinear score = 1 > 0.999 
## Covariates in the formula are so strongly correlated that the
## parameter estimates from this model are not meaningful.
## Dropping one or more of the covariates will fix this problem

## Attempting again with only condition failed.

## Error in simple_varpart(expt, predictor = NULL, factors = c("condition",  : 
##   
## Error in dstat0[, i] : subscript out of bounds

## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

## Warning in doTryCatch(return(expr), name, parentenv, handler): display list
## redraw incomplete

## varpart sees only 1 batch, adjusting the model accordingly.

## Error in if (ncol(exprObj) != nrow(data)) { : argument is of length zero

## Retrying with only condition in the model.

## Error in checkModelStatus(fit, showWarnings = showWarnings, colinearityCutoff = colinearityCutoff) : 
##   Colinear score = 1 > 0.999 
## Covariates in the formula are so strongly correlated that the
## parameter estimates from this model are not meaningful.
## Dropping one or more of the covariates will fix this problem

## Attempting again with only condition failed.

## Error in simple_varpart(norm_data, predictor = NULL, factors = c("condition",  : 
## 

combined_norm <- normalize_expt(combined_expt, convert="rpkm",
                                transform="log2", na_to_zero=TRUE)

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

combined_exprs <- exprs(combined_norm)

ess_circos <- ess_table[, c("vitro_control_call", "vitro_low_call", "vitro_high_call",
                            "vivo_cx_call", "vivo_d1_call", "vivo_d3_call", "vivo_tn_call",
                            "vivo_ut_call", "vivo_vg_call")]
rownames(ess_circos) <- rownames(ess_table)

colors <- c("990000", "000000", "000099", "777777")
names(colors) <- c("E", "NE", "S", "U")

## A Reminder of the contrasts in tn_tables[["data"]]
##  "ut" = c("ut", "control"),
##  "d1" = c("d1", "control"),
##  "d3" = c("d3", "control"),
##  "cx" = c("cx", "control"),
##  "vg" = c("vg", "control"))

low_df <- tn_tables[["data"]][["low_cal_vs_control"]]
high_df <- tn_tables[["data"]][["high_cal_vs_control"]]
ut_df <- tn_tables[["data"]][["ut_vs_control"]]
vg_df <- tn_tables[["data"]][["vg_vs_control"]]
cx_df <- tn_tables[["data"]][["cx_vs_control"]]
d1_df <- tn_tables[["data"]][["d1_vs_control"]]
d3_df <- tn_tables[["data"]][["d3_vs_control"]]

circos_cfg <- circos_prefix(all_annot, name="cjb", stop_column = "end", cog_column = "COGFun")

## This assumes you have a colors.conf in circos/colors/ and fonts.conf in circos/fonts/

## It also assumes you have conf/ideogram.conf, conf/ticks.conf, and conf/housekeeping.conf

## It will write circos/conf/cjb.conf with a reasonable first approximation config file.

## Wrote karyotype to circos/conf/ideograms/cjb.conf

## This should match the karyotype= line in cjb.conf

## Wrote ticks to circos/conf/ticks_cjb.conf

cjb_kary <- circos_karyotype(circos_cfg, fasta = "reference/sagalactiae_cjb111.fasta")

## Wrote karyotype to circos/conf/karyotypes/cjb.conf

## This should match the karyotype= line in cjb.conf

cjb_plus_minus <- circos_plus_minus(circos_cfg, width = 0.06, thickness = 40,
                                    url_string = "http://www.microbesonline.org/cgi-bin/fetchLocus.cgi?locus=[id]")

## Writing data file: circos/data/cjb_plus_go.txt with the + strand GO data.

## Writing data file: circos/data/cjb_minus_go.txt with the - strand GO data.

## Wrote the +/- config files.  Appending their inclusion to the master file.

## Returning the inner width: 0.88.  Use it as the outer for the next ring.

##cjb_control <- circos_hist(circos_cfg, combined_exprs, colname="combined_control",
##                           basename="control", outer=cjb_plus_minus,
##                           fill_color="vvdpblue", width=0.06)
##cjb_low_rpkm <- circos_hist(circos_cfg, combined_exprs, colname="combined_low",
##                            basename="low", outer=cjb_control_rpkm,
##                            fill_color="vdpblue", width=0.06)
##cjb_high_rpkm <- circos_hist(circos_cfg, combined_exprs, colname="combined_high",
##                             basename="high", outer=cjb_low_rpkm,
##                             fill_color="dpblue", width=0.06)
##cjb_cx <- circos_hist(circos_cfg, combined_exprs, colname="combined_cx",
##                      basename="cx", outer=cjb_high_rpkm,
##                      fill_color="dpblue", width=0.06)
##cjb_d1 <- circos_hist(circos_cfg, combined_exprs, colname="combined_d1",
##                      basename="d1", outer=cjb_cx,
##                      fill_color="dpblue", width=0.06)
##cjb_d3 <- circos_hist(circos_cfg, combined_exprs, colname="combined_d3",
##                      basename="d3", outer=cjb_d1,
##                      fill_color="dpblue", width=0.06)
##cjb_ut <- circos_hist(circos_cfg, combined_exprs, colname="combined_ut",
##                      basename="ut", outer=cjb_d3,
##                      fill_color="dpblue", width=0.06)
##cjb_vg <- circos_hist(circos_cfg, combined_exprs, colname="combined_vg",
##                      fill_color="dpblue", width=0.06)
cjb_control_tile <- circos_tile(circos_cfg, ess_circos, colname="vitro_control_call",
                                colors=colors, basename="control_tile",
                                outer=cjb_plus_minus, width=0.06)

## Writing data file: circos/data/cjbvitro_control_call_tile.txt with the control_tilevitro_control_call column.

## Returning the inner width: 0.82.  Use it as the outer for the next ring.

cjb_high_tile <- circos_tile(circos_cfg, ess_circos, colname="vitro_high_call",
                                colors=colors, basename="high_tile",
                             outer=cjb_control_tile, width=0.06)

## Writing data file: circos/data/cjbvitro_high_call_tile.txt with the high_tilevitro_high_call column.

## Returning the inner width: 0.76.  Use it as the outer for the next ring.

cjb_low_tile <- circos_tile(circos_cfg, ess_circos, colname="vitro_low_call",
                                colors=colors, basename="low_tile",
                                outer=cjb_high_tile, width=0.06)

## Writing data file: circos/data/cjbvitro_low_call_tile.txt with the low_tilevitro_low_call column.

## Returning the inner width: 0.7.  Use it as the outer for the next ring.

cjb_tn_tile <- circos_tile(circos_cfg, ess_circos, colname="vivo_tn_call",
                           colors=colors, basename="tn_tile",
                           outer=cjb_low_tile, width=0.06)

## Writing data file: circos/data/cjbvivo_tn_call_tile.txt with the tn_tilevivo_tn_call column.

## Returning the inner width: 0.64.  Use it as the outer for the next ring.

cjb_d1_tile <- circos_tile(circos_cfg, ess_circos, colname="vivo_d1_call",
                           colors=colors, basename="d1_tile",
                           outer=cjb_tn_tile, width=0.06)

## Writing data file: circos/data/cjbvivo_d1_call_tile.txt with the d1_tilevivo_d1_call column.

## Returning the inner width: 0.58.  Use it as the outer for the next ring.

cjb_d3_tile <- circos_tile(circos_cfg, ess_circos, colname="vivo_d3_call",
                           colors=colors, basename="d3_tile",
                           outer=cjb_d1_tile, width=0.06)

## Writing data file: circos/data/cjbvivo_d3_call_tile.txt with the d3_tilevivo_d3_call column.

## Returning the inner width: 0.52.  Use it as the outer for the next ring.

cjb_vg_tile <- circos_tile(circos_cfg, ess_circos, colname="vivo_vg_call",
                           colors=colors, basename="vg_tile",
                           outer=cjb_d3_tile, width=0.06)

## Writing data file: circos/data/cjbvivo_vg_call_tile.txt with the vg_tilevivo_vg_call column.

## Returning the inner width: 0.46.  Use it as the outer for the next ring.

cjb_cx_tile <- circos_tile(circos_cfg, ess_circos, colname="vivo_cx_call",
                           colors=colors, basename="cx_tile",
                           outer=cjb_vg_tile, width=0.06)

## Writing data file: circos/data/cjbvivo_cx_call_tile.txt with the cx_tilevivo_cx_call column.

## Returning the inner width: 0.4.  Use it as the outer for the next ring.

cjb_ut_tile <- circos_tile(circos_cfg, ess_circos, colname="vivo_ut_call",
                           colors=colors, basename="ut_tile",
                           outer=cjb_cx_tile, width=0.06)

## Writing data file: circos/data/cjbvivo_ut_call_tile.txt with the ut_tilevivo_ut_call column.

## Returning the inner width: 0.34.  Use it as the outer for the next ring.

cjb_suffix <- circos_suffix(circos_cfg)
##made <- circos_make(circos_cfg, target="cjb")

pander::pander(sessionInfo())

R version 4.0.3 (2020-10-10)

Platform: x86_64-pc-linux-gnu (64-bit)

locale: LC_CTYPE=en_US.UTF-8, LC_NUMERIC=C, LC_TIME=en_US.UTF-8, LC_COLLATE=en_US.UTF-8, LC_MONETARY=en_US.UTF-8, LC_MESSAGES=en_US.UTF-8, LC_PAPER=en_US.UTF-8, LC_NAME=C, LC_ADDRESS=C, LC_TELEPHONE=C, LC_MEASUREMENT=en_US.UTF-8 and LC_IDENTIFICATION=C

attached base packages: splines, stats4, parallel, stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: ruv(v.0.9.7.1), lme4(v.1.1-26), Matrix(v.1.3-2), BiocParallel(v.1.24.1), variancePartition(v.1.20.0), hpgltools(v.1.0), testthat(v.3.0.2), SummarizedExperiment(v.1.20.0), GenomicRanges(v.1.42.0), GenomeInfoDb(v.1.26.2), IRanges(v.2.24.1), S4Vectors(v.0.28.1), MatrixGenerics(v.1.2.1), matrixStats(v.0.58.0), Biobase(v.2.50.0) and BiocGenerics(v.0.36.0)

loaded via a namespace (and not attached): rappdirs(v.0.3.3), rtracklayer(v.1.50.0), R.methodsS3(v.1.8.1), tidyr(v.1.1.2), ggplot2(v.3.3.3), bit64(v.4.0.5), knitr(v.1.31), DelayedArray(v.0.16.1), R.utils(v.2.10.1), data.table(v.1.14.0), RCurl(v.1.98-1.2), doParallel(v.1.0.16), generics(v.0.1.0), GenomicFeatures(v.1.42.1), preprocessCore(v.1.52.1), callr(v.3.5.1), cowplot(v.1.1.1), usethis(v.2.0.1), RSQLite(v.2.2.3), shadowtext(v.0.0.7), bit(v.4.0.4), enrichplot(v.1.10.2), xml2(v.1.3.2), assertthat(v.0.2.1), viridis(v.0.5.1), xfun(v.0.21), hms(v.1.0.0), jquerylib(v.0.1.3), evaluate(v.0.14), IHW(v.1.18.0), DEoptimR(v.1.0-8), fansi(v.0.4.2), progress(v.1.2.2), caTools(v.1.18.1), dbplyr(v.2.1.0), igraph(v.1.2.6), DBI(v.1.1.1), geneplotter(v.1.68.0), purrr(v.0.3.4), ellipsis(v.0.3.1), selectr(v.0.4-2), dplyr(v.1.0.4), backports(v.1.2.1), annotate(v.1.68.0), biomaRt(v.2.46.3), vctrs(v.0.3.6), remotes(v.2.2.0), cachem(v.1.0.4), withr(v.2.4.1), ggforce(v.0.3.2), BSgenome(v.1.58.0), robustbase(v.0.93-7), GenomicAlignments(v.1.26.0), fdrtool(v.1.2.16), prettyunits(v.1.1.1), DOSE(v.3.16.0), crayon(v.1.4.1), genefilter(v.1.72.1), edgeR(v.3.32.1), pkgconfig(v.2.0.3), slam(v.0.1-48), labeling(v.0.4.2), tweenr(v.1.0.1), nlme(v.3.1-152), pkgload(v.1.2.0), devtools(v.2.3.2), rlang(v.0.4.10), lifecycle(v.1.0.0), downloader(v.0.4), BiocFileCache(v.1.14.0), directlabels(v.2021.1.13), rprojroot(v.2.0.2), polyclip(v.1.10-0), graph(v.1.68.0), lpsymphony(v.1.18.0), boot(v.1.3-27), processx(v.3.4.5), viridisLite(v.0.3.0), bitops(v.1.0-6), R.oo(v.1.24.0), KernSmooth(v.2.23-18), pander(v.0.6.3), Biostrings(v.2.58.0), blob(v.1.2.1), stringr(v.1.4.0), qvalue(v.2.22.0), readr(v.1.4.0), scales(v.1.1.1), memoise(v.2.0.0), magrittr(v.2.0.1), plyr(v.1.8.6), gplots(v.3.1.1), zlibbioc(v.1.36.0), compiler(v.4.0.3), scatterpie(v.0.1.5), RColorBrewer(v.1.1-2), DESeq2(v.1.30.1), Rsamtools(v.2.6.0), cli(v.2.3.1), XVector(v.0.30.0), ps(v.1.5.0), MASS(v.7.3-53.1), mgcv(v.1.8-34), tidyselect(v.1.1.0), stringi(v.1.5.3), highr(v.0.8), yaml(v.2.2.1), GOSemSim(v.2.16.1), askpass(v.1.1), locfit(v.1.5-9.4), ggrepel(v.0.9.1), grid(v.4.0.3), sass(v.0.3.1), VariantAnnotation(v.1.36.0), fastmatch(v.1.1-0), tools(v.4.0.3), rstudioapi(v.0.13), foreach(v.1.5.1), genbankr(v.1.18.0), gridExtra(v.2.3), farver(v.2.0.3), Rtsne(v.0.15), ggraph(v.2.0.5), digest(v.0.6.27), rvcheck(v.0.1.8), BiocManager(v.1.30.10), quadprog(v.1.5-8), Rcpp(v.1.0.6), broom(v.0.7.5), httr(v.1.4.2), AnnotationDbi(v.1.52.0), colorspace(v.2.0-0), rvest(v.0.3.6), XML(v.3.99-0.5), fs(v.1.5.0), RBGL(v.1.66.0), statmod(v.1.4.35), PROPER(v.1.22.0), graphlayouts(v.0.7.1), sessioninfo(v.1.1.1), xtable(v.1.8-4), jsonlite(v.1.7.2), nloptr(v.1.2.2.2), tidygraph(v.1.2.0), corpcor(v.1.6.9), R6(v.2.5.0), Vennerable(v.3.1.0.9000), pillar(v.1.5.0), htmltools(v.0.5.1.1), glue(v.1.4.2), fastmap(v.1.1.0), minqa(v.1.2.4), clusterProfiler(v.3.18.1), codetools(v.0.2-18), fgsea(v.1.16.0), pkgbuild(v.1.2.0), utf8(v.1.1.4), lattice(v.0.20-41), bslib(v.0.2.4), tibble(v.3.0.6), sva(v.3.38.0), pbkrtest(v.0.5-0.1), curl(v.4.3), rentrez(v.1.2.3), colorRamps(v.2.3), gtools(v.3.8.2), zip(v.2.1.1), GO.db(v.3.12.1), openxlsx(v.4.2.3), openssl(v.1.4.3), survival(v.3.2-7), limma(v.3.46.0), rmarkdown(v.2.7), desc(v.1.2.0), munsell(v.0.5.0), DO.db(v.2.9), fastcluster(v.1.1.25), GenomeInfoDbData(v.1.2.4), iterators(v.1.0.13), reshape2(v.1.4.4) and gtable(v.0.3.0)

message(paste0("This is hpgltools commit: ", get_git_commit()))

## If you wish to reproduce this exact build of hpgltools, invoke the following:

## > git clone http://github.com/abelew/hpgltools.git

## > git reset 72947fcc6afe09da22d71967059edd84e3063341

## This is hpgltools commit: Tue Jun 1 15:57:56 2021 -0400: 72947fcc6afe09da22d71967059edd84e3063341

this_save <- paste0(gsub(pattern="\\.Rmd", replace="", x=rmd_file), "-v", ver, ".rda.xz")
message(paste0("Saving to ", this_save))

## Saving to index_202105-v202105.rda.xz

tmp <- sm(saveme(filename=this_save))

loadme(filename=this_save)