Sequencing 3 Mycobacterium Tuberculosis strains.

1 Introduction

This project seeks to understand the effects of deleting/complementing pknF in Mtb. I cannot honestly say I know anything about the biology involved except the bits from the review Volker kindly sent me.

pknF is one of the surprisingly large number of serine threonine protein kinases in Mtb. Most bacteria have only a couple/few of these, but Mtb has 11, and uses them with similar frequency to the two component systems in order to respond to its changing environment during infection (I assume, I have not gotten that far in the review yet). In previous experiments performed by people in Volker’s lab, a deletion of this particular STPK causes a much more virulent phenotype during mouse infection.

2 Preprocessing

I wrote a quick sample sheet for these 9 samples and will use my cyoa toy to perform the various tasks to process the data.

2.1 What annotations do I have?

I am assuming the most appropriate annotations come from microbesonline? Let us check real quick. I should probably go back through my old logs to answer these questions.

annot <- load_microbesonline_annotations("CDC1551")

## Found 1 entry.

## Mycobacterium tuberculosis CDC1551Actinobacteriayes2007-05-08yes10429383331

## The species being downloaded is: Mycobacterium tuberculosis CDC1551

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

annot <- as.data.frame(annot)
rownames(annot) <- make.names(annot[["sysName"]], unique = TRUE)
go_db <- load_microbesonline_go(id = 83331, id_column = "sysName")

## The species being downloaded is: Mycobacterium tuberculosis CDC1551 and is being downloaded as 83331.tab.

h37_annot <- load_microbesonline_annotations("H37Rv")

## Found 1 entry.

## Mycobacterium tuberculosis H37RvActinobacteriayes2007-05-08yes10404783332

## The species being downloaded is: Mycobacterium tuberculosis H37Rv

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

h37_go <- load_microbesonline_go(id = 83332, id_column = "sysName")

## The species being downloaded is: Mycobacterium tuberculosis H37Rv and is being downloaded as 83332.tab.

That looks pretty reasonable to me, it even has COG groups and all the fun stuff I want.

2.2 Which genome am I using?

Before I start, I should recall which is my current reference genome and see if there is a better alternative since last I looked.

My current genome is strain H37Rv last downloaded in 2022, it is derived from accession NC_000962.3 / assembly GCF_000195955.2. Glancing at the gff file I have, I am guessing that the most appropriate tags to index against are either locus_tag(Rv0003) or gene(recF).

cd preprocessing
start=$(pwd)
for i in $(/bin/ls -d * | grep -v Undeter); do
    cd $i
    cyoa --method prnaseq --species mtuberculosis_h37rv \
         --gff_type gene --gff_tag locus_tag \
         --introns 0 --freebayes 1 \
         --input $(/bin/ls *.gz | tr '\n' ':' | sed 's/:$//g')
    cd $start
done

2.3 The correct Strain is CDC1551

Ok, so let us redo this accordingly.

cd preprocessing
start=$(pwd)
for i in $(/bin/ls -d * | grep -v Undeter); do
    cd $i
    cyoa --method prnaseq --species mycobacterium_tuberculosis_cdc1551 \
         --gff_type gene --gff_tag locus_tag \
         --introns 0 --freebayes 1 \
         --input $(/bin/ls *.gz | tr '\n' ':' | sed 's/:$//g')
    cd $start
done

3 Ortholog search of CDC1551/H37Rv

A few minutes ago I was making a couple of notes to myself, one of which was to explicitly note where PknF is in this strain. Sadly, when I looked into the annotations I downloaded it was unnamed. I was able to make a good guess given that I knew the location in H37Rv and saw a group of genes in CDC1551 in basically the same place in a similar orientation. Later, when the DE finished I was able to prove it by showing that MT1788 (position 1962996) is the most decreased in the Δ strain and most increased in the complement.

With that in mind, I was a little sad to not have the full set of H37Rv annotations available and so am spinning up my favorite ortholog search tool ‘orthofinder’. The following block is where I am going to do that.

I first need to extract the amino acid sequences from H37Rv (I have them in place already for CDC1551).

mkdir orthologs
cd orthologs
cp ~/libraries/genome/mtub*.gb .
cp ~/libraries/genome/myco*.faa .
cyoa --method gb2gff --input *.gb
ls *.faa
cyoa --method orthofinder --input mtuberculosis_h37rv.faa:mycobacterium_tuberculosis_cdc1551.faa

For the curious, that last command generates the following on the computing cluster:

orthofinder -f outputs/50orthofinder/input \
  -o outputs/50orthofinder/output \
  2>>outputs/50orthofinder/orthofinder.stderr 1>>outputs/50orthofinder/orthofinder.stdout
mv outputs/50orthofinder/output/Results_Dec06/* outputs/50orthofinder/output
rmdir outputs/50orthofinder/output/Results_Dec06

It writes the input file for me so that I cannot introduce any annoying typeographical errrors[sic]. It also creates a little job that parses the quite confusing outputs from orthofinder into an (hopefully) easy to read table.

My hope is to shortly be able to append all the H37Rv IDs to my annotation set for CDC1551. We shall see, orthofinder takes ~ 10 minutes / 1000 genes.

4 Append information to the sample sheet

We have a relatively sparse sample sheet to start with, lets add some summary stats from cyoa to it.

There is one important piece of information that Daniel alerted me to which is not included in the gathered metadata (I think): these samples appear to have an extraorindarily large amount of potentially other bacterial sequence in them; and it appears to be dominated by Pseudomonas aeruginosa – I am assuming PA01 and will try mapping the reads against that strain to see what is going on.

Until those mappings are complete, I will continue on as normal, but keep in mind that the following will almost certain need to change.

spec <- make_rnaseq_spec()
written <- gather_preprocessing_metadata("sample_sheets/all_samples.xlsx",
                                         species = "mycobacterium_tuberculosis_cdc1551",
                                         specification = spec)

## Using provided specification

## Example filename: preprocessing/WT_1/outputs/*trimomatic/*-trimomatic.stderr.
## Example filename: preprocessing/WT_1/outputs/*trimomatic/*-trimomatic.stderr.

## The numerator column is: trimomatic_output.

## The denominator column is: trimomatic_input.

## Example filename: preprocessing/WT_1/outputs/*fastqc/*_fastqc/fastqc_data.txt.

## Skipping for now

## Not including new entries for: fastqc_most_overrepresented, it is empty.

## Example filename: preprocessing/WT_1/outputs/*kraken_*/kraken_report_matrix.tsv.

## Example filename: preprocessing/WT_1/outputs/*hisat2_mycobacterium_tuberculosis_cdc1551/hisat2_*rRNA*.stderr.

## Not including new entries for: hisat_rrna_single_concordant, it is empty.

## Example filename: preprocessing/WT_1/outputs/*hisat2_mycobacterium_tuberculosis_cdc1551/hisat2_*rRNA*.stderr.

## Not including new entries for: hisat_rrna_multi_concordant, it is empty.

## Missing data to calculate the ratio between: hisat_rrna_multi_concordant and trimomatic_output.

## Not including new entries for: hisat_rrna_percent, it is empty.

## Example filename: preprocessing/WT_1/outputs/*hisat2_mycobacterium_tuberculosis_cdc1551/hisat2_*genome*.stderr.
## Example filename: preprocessing/WT_1/outputs/*hisat2_mycobacterium_tuberculosis_cdc1551/hisat2_*genome*.stderr.
## Example filename: preprocessing/WT_1/outputs/*hisat2_mycobacterium_tuberculosis_cdc1551/hisat2_*genome*.stderr.
## Example filename: preprocessing/WT_1/outputs/*hisat2_mycobacterium_tuberculosis_cdc1551/hisat2_*genome*.stderr.
## Example filename: preprocessing/WT_1/outputs/*hisat2_mycobacterium_tuberculosis_cdc1551/hisat2_*genome*.stderr.

## The numerator column is: hisat_genome_single_concordant.

## The denominator column is: trimomatic_output.

## Example filename: preprocessing/WT_1/outputs/*hisat2_mycobacterium_tuberculosis_cdc1551/mycobacterium_tuberculosis_cdc1551_genome*.count.xz.

## Example filename: preprocessing/WT_1/outputs/*hisat2_mycobacterium_tuberculosis_cdc1551/mycobacterium_tuberculosis_cdc1551_genome*.count.xz.

## Example filename: preprocessing/WT_1/outputs/*salmon_*/salmon_*.stderr.

## Not including new entries for: salmon_stranded, it is empty.

## Example filename: preprocessing/WT_1/outputs/*salmon_mycobacterium_tuberculosis_cdc1551/quant.sf.

## Warning in dispatch_filename_search(meta, input_file_spec, verbose = verbose, : The input file is NA for:
## preprocessing/WT_1/outputs/*salmon_mycobacterium_tuberculosis_cdc1551/quant.sf.

## Warning in dispatch_filename_search(meta, input_file_spec, verbose = verbose, : The input file is NA for:
## preprocessing/WT_2/outputs/*salmon_mycobacterium_tuberculosis_cdc1551/quant.sf.

## Warning in dispatch_filename_search(meta, input_file_spec, verbose = verbose, : The input file is NA for:
## preprocessing/WT_3/outputs/*salmon_mycobacterium_tuberculosis_cdc1551/quant.sf.

## Warning in dispatch_filename_search(meta, input_file_spec, verbose = verbose, : The input file is NA for:
## preprocessing/Mutant_pknF_1/outputs/*salmon_mycobacterium_tuberculosis_cdc1551/quant.sf.

## Warning in dispatch_filename_search(meta, input_file_spec, verbose = verbose, : The input file is NA for:
## preprocessing/Mutant_pknF_2/outputs/*salmon_mycobacterium_tuberculosis_cdc1551/quant.sf.

## Warning in dispatch_filename_search(meta, input_file_spec, verbose = verbose, : The input file is NA for:
## preprocessing/Mutant_pknF_3/outputs/*salmon_mycobacterium_tuberculosis_cdc1551/quant.sf.

## Warning in dispatch_filename_search(meta, input_file_spec, verbose = verbose, : The input file is NA for:
## preprocessing/Comp_pknF_1/outputs/*salmon_mycobacterium_tuberculosis_cdc1551/quant.sf.

## Warning in dispatch_filename_search(meta, input_file_spec, verbose = verbose, : The input file is NA for:
## preprocessing/Comp_pknF_2/outputs/*salmon_mycobacterium_tuberculosis_cdc1551/quant.sf.

## Warning in dispatch_filename_search(meta, input_file_spec, verbose = verbose, : The input file is NA for:
## preprocessing/Comp_pknF_3/outputs/*salmon_mycobacterium_tuberculosis_cdc1551/quant.sf.

## Not including new entries for: salmon_count_table, it is empty.

## Example filename: preprocessing/WT_1/outputs/*salmon_*/salmon_*.stderr.

## Not including new entries for: salmon_mapped, it is empty.

## Example filename: preprocessing/WT_1/outputs/*salmon_*/salmon_*.stderr.

## Not including new entries for: salmon_percent, it is empty.

## Writing new metadata to: sample_sheets/all_samples_modified.xlsx

## Deleting the file sample_sheets/all_samples_modified.xlsx before writing the tables.

5 Create an expressionset

expt <- create_expt(written[["new_file"]], gene_info = annot, file_column = "hisatcounttable")

## Reading the sample metadata.

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

## Matched 4245 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 4245 features and 9 samples.

6 Poke at it

plot_libsize(expt)

## Library sizes of 9 samples, 
## ranging from 3,338,710 to 13,076,335.

plot_nonzero(expt)

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
## A non-zero genes plot of 9 samples.
## These samples have an average 8.867 CPM coverage and 4235 genes observed, ranging from 4227 to
## 4242.

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

## Removing 91 low-count genes (4154 remaining).

plot_corheat(norm)

## A heatmap of pairwise sample correlations ranging from: 
## 0.924732991192705 to 0.987981556516807.

plot_disheat(norm)

## A heatmap of pairwise sample distances ranging from: 
## 23.2338031782721 to 58.1252671359335.

plot_pca(norm)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by wt, delta, complement
## Shapes are defined by b1, b2, b3.

nb <- normalize_expt(expt, transform = "log2", convert = "cpm", batch = "limma", filter = TRUE)

## Removing 91 low-count genes (4154 remaining).
## If you receive a warning: 'NANs produced', one potential reason is that the data was quantile normalized.
## Setting 2 low elements to zero.
## transform_counts: Found 2 values equal to 0, adding 1 to the matrix.

plot_pca(nb)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by wt, delta, complement
## Shapes are defined by b1, b2, b3.

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

## Removing 91 low-count genes (4154 remaining).
## Setting 2 low elements to zero.
## transform_counts: Found 2 values equal to 0, adding 1 to the matrix.

plot_pca(ns)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by wt, delta, complement
## Shapes are defined by b1, b2, b3.

Well, that is not great. I arbitrarily tried a batch in model and sva adjustment, but am not including them without further information.

7 varpart

varpart <- simple_varpart(expt)

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

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

## 
## Total:34 s

varpart

## The result of using variancePartition with the model:
## ~ condition + batch

8 Kraken matrix

genus_expt <- create_expt("sample_sheets/all_samples.xlsx",
                          file_column = "krakentsv", file_type = "table")

## Reading the sample metadata.

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

## Error in create_expt("sample_sheets/all_samples.xlsx", file_column = "krakentsv", : This requires either a count dataframe/matrix or column containing the filenames.

genus_norm <- normalize_expt(genus_expt, convert = "cpm")

## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'state': object 'genus_expt' not found

plot_disheat(genus_norm)

## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt_data' in selecting a method for function 'plot_heatmap': object 'genus_norm' not found

genus_normv2 <- normalize_expt(genus_expt, convert = "cpm", transform = "log2")

## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'state': object 'genus_expt' not found

plot_pca(genus_normv2)

## Error in eval(expr, envir, enclos): object 'genus_normv2' not found

plot_libsize(genus_expt)

## Error in h(simpleError(msg, call)): error in evaluating the argument 'data' in selecting a method for function 'plot_libsize': object 'genus_expt' not found

head(exprs(genus_expt))

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'head': error in evaluating the argument 'object' in selecting a method for function 'exprs': object 'genus_expt' not found

exprs(genus_expt)["Pseudomonas", ]

## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'exprs': object 'genus_expt' not found

9 Differential Expression

I am going to assume for the moment that we can learn interesting things from the data, but that it will require more intervention than usual for a bacterial dataset.

keepers <- list(
  "delta_wt" = c("delta", "wt"),
  "comp_wt" = c("complement", "wt"),
  "delta_comp" = c("delta", "complement"))
de <- all_pairwise(expt, keepers = keepers, filter = TRUE, model_batch = "svaseq")

## 
##         wt      delta complement 
##          3          3          3

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

## Setting 2 low elements to zero.

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

tables <- combine_de_tables(
  de, keepers = keepers, excel = glue("excel/de_tables-v{ver}.xlsx"))

## Checking limma for name delta_wt:wt_vs_delta

## Checking deseq for name delta_wt:wt_vs_delta

## Checking edger for name delta_wt:wt_vs_delta

## Checking ebseq for name delta_wt:wt_vs_delta

## Checking noiseq for name delta_wt:wt_vs_delta

## Checking basic for name delta_wt:wt_vs_delta

## Checking limma for name comp_wt:wt_vs_complement

## Checking deseq for name comp_wt:wt_vs_complement

## Checking edger for name comp_wt:wt_vs_complement

## Checking ebseq for name comp_wt:wt_vs_complement

## Checking noiseq for name comp_wt:wt_vs_complement

## Checking basic for name comp_wt:wt_vs_complement

## Checking limma for name delta_comp:complement_vs_delta

## Checking deseq for name delta_comp:complement_vs_delta

## Checking edger for name delta_comp:complement_vs_delta

## Checking ebseq for name delta_comp:delta_vs_complement

## Checking noiseq for name delta_comp:complement_vs_delta

## Checking basic for name delta_comp:complement_vs_delta

## About to start combine_mapped_table()

## Inverting column: basic_logfc.

## Inverting column: deseq_logfc.

## Inverting column: ebseq_logfc.

## Inverting column: edger_logfc.

## Inverting column: limma_logfc.

## Inverting column: noiseq_logfc.

## Finished combine_mapped_table()

## About to start combine_mapped_table()

## Inverting column: basic_logfc.

## Inverting column: deseq_logfc.

## Inverting column: ebseq_logfc.

## Inverting column: edger_logfc.

## Inverting column: limma_logfc.

## Inverting column: noiseq_logfc.

## Finished combine_mapped_table()

## About to start combine_mapped_table()

## Inverting column: basic_logfc.

## Inverting column: deseq_logfc.

## Inverting column: edger_logfc.

## Inverting column: limma_logfc.

## Inverting column: noiseq_logfc.

## Finished combine_mapped_table()

sig <- extract_significant_genes(
  tables, according_to = "deseq", excel = glue("excel/de_sig-v{ver}.xlsx"))

9.1 The intersection of Δ/WT and Δ/comp

library(UpSetR)
first_up <- sig[["deseq"]][["ups"]][["delta_wt"]]
second_up <- sig[["deseq"]][["ups"]][["delta_comp"]]

first_down <- sig[["deseq"]][["downs"]][["delta_wt"]]
second_down <- sig[["deseq"]][["downs"]][["delta_comp"]]

comparison <- list(
  "d_w_up" = rownames(first_up),
  "d_c_up" = rownames(second_up),
  "d_w_down" = rownames(first_down),
  "d_c_down" = rownames(second_down))
query <- UpSetR::fromList(comparison)
sets <- UpSetR::upset(query)
sets

## So, there are 2 shared down and 12 up.  Also 1 weirdo (which given
## its proximity to the operon containing PknF is likely interesting)

overlap_ids <- overlap_groups(comparison)
interesting_up <- overlap_geneids(overlap_ids, "d_w_up:d_c_up")
interesting_up

##  d_w_up3  d_w_up6 d_w_up18 d_w_up20 d_w_up29 d_w_up30 d_w_up37 d_w_up45 d_w_up57 d_w_up71 d_w_up90 d_w_up95 
## "MT3214" "MT0995" "MT2020" "MT0423" "MT1255" "MT0424" "MT2579" "MT2574" "MT3712" "MT2576" "MT1256" "MT0595"

interesting_down <- overlap_geneids(overlap_ids, "d_w_down:d_c_down")
interesting_down

## d_c_up25 d_c_up26 
## "MT1702" "MT0995"

weirdo <- overlap_geneids(overlap_ids, "d_w_up:d_c_down")
weirdo

##  d_w_up2 
## "MT1790"

10 Ontology enrichment

Given the relative sparsity of the annotations for CDC1551, I might either go back to my H37Rv mapping data (which was pretty decent), or use the orthologs and their annotations; but since I already have the data in its current form loaded, let us just try it and see.

10.1 Volker’s explicit query first

There are only 15 genes in the set of interesting intersections; so I am going to group them together.

In my annotation block above I downloaded the H37Rv and CDC1551 annotations/GO categories. There appears to be a bug in the GO loader, it needed the actual ID number instead of a unique string like the annotations; happily the annotation load prints out the ID number so I did not have to think very hard.

annot[["length"]] <- abs(annot[["start"]] - annot[["stop"]])
length_db <- as.data.frame(annot[, c("name", "length")])
colnames(length_db) <- c("ID", "length")
query_genes <- unique(c(interesting_up, interesting_down, weirdo))
goseq_query <- simple_goseq(query_genes, go_db = go_db,
                            length_db = length_db, min_xref = 4)

## Found 9 go_db genes and 11 length_db genes out of 14.

## Testing that go categories are defined.

## Removing undefined categories.

## Gathering synonyms.

## Gathering category definitions.

goseq_query

## A set of ontologies produced by gprofiler using 14
## with significance cutoff 0.1.
## There are 4 MF hits, 1, BP hits, and 1 CC hits.

## Category mfp_plot_over is the most populated with 3 hits.

Given that there are only a few genes in the above query, let us also do the three primary contrasts separately as I think they will prove more interesting.

goseq_up_delta_wt <- simple_goseq(sig[["deseq"]][["ups"]][["delta_wt"]], go_db = go_db,
                                  length_db = length_db)

## Found 67 go_db genes and 83 length_db genes out of 108.

## Testing that go categories are defined.

## Removing undefined categories.

## Gathering synonyms.

## Gathering category definitions.

goseq_up_delta_wt

## A set of ontologies produced by gprofiler using 108
## with significance cutoff 0.1.
## There are 27 MF hits, 12, BP hits, and 2 CC hits.

## Category mfp_plot_over is the most populated with 11 hits.

goseq_down_delta_wt <- simple_goseq(sig[["deseq"]][["downs"]][["delta_wt"]], go_db = go_db,
                                    length_db = length_db)

## Found 3 go_db genes and 17 length_db genes out of 17.

goseq_down_delta_wt

## NULL

goseq_up_comp_wt <- simple_goseq(sig[["deseq"]][["ups"]][["comp_wt"]], go_db = go_db,
                                  length_db = length_db, min_xref = 20)

## Found 27 go_db genes and 51 length_db genes out of 54.

## Testing that go categories are defined.

## Removing undefined categories.

## Gathering synonyms.

## Gathering category definitions.

goseq_up_comp_wt

## A set of ontologies produced by gprofiler using 54
## with significance cutoff 0.1.
## There are 20 MF hits, 11, BP hits, and 0 CC hits.

## Category mfp_plot_over is the most populated with 8 hits.

goseq_down_comp_wt <- simple_goseq(sig[["deseq"]][["downs"]][["comp_wt"]], go_db = go_db,
                                    length_db = length_db, min_xref = 5)

## Found 6 go_db genes and 24 length_db genes out of 25.

## Testing that go categories are defined.

## Removing undefined categories.

## Gathering synonyms.

## Gathering category definitions.

goseq_down_comp_wt

## A set of ontologies produced by gprofiler using 25
## with significance cutoff 0.1.
## There are 4 MF hits, 3, BP hits, and 0 CC hits.

## Category mfp_plot_over is the most populated with 3 hits.

goseq_up_delta_comp <- simple_goseq(sig[["deseq"]][["ups"]][["delta_comp"]], go_db = go_db,
                                    length_db = length_db, min_xref = 10)

## Found 20 go_db genes and 32 length_db genes out of 36.

## Testing that go categories are defined.

## Removing undefined categories.

## Gathering synonyms.

## Gathering category definitions.

goseq_up_delta_comp

## A set of ontologies produced by gprofiler using 36
## with significance cutoff 0.1.
## There are 8 MF hits, 3, BP hits, and 1 CC hits.

## Category mfp_plot_over is the most populated with 7 hits.

goseq_down_delta_comp <- simple_goseq(sig[["deseq"]][["downs"]][["delta_comp"]], go_db = go_db,
                                    length_db = length_db, min_xref = 5)

## Found 7 go_db genes and 31 length_db genes out of 32.

## Testing that go categories are defined.

## Removing undefined categories.

## Gathering synonyms.

## Gathering category definitions.

goseq_down_delta_comp

## A set of ontologies produced by gprofiler using 32
## with significance cutoff 0.1.
## There are 4 MF hits, 2, BP hits, and 1 CC hits.

## Category mfp_plot_over is the most populated with 4 hits.

11 Circos

I think it would be nice to have the wt mean expression values as a baseline

means <- mean_by_factor(norm)[["medians"]]

## The factor wt has 3 rows.

## The factor delta has 3 rows.

## The factor complement has 3 rows.

annot[["chromosome"]] <- "AE000516"
delta_wt_df <- tables[["data"]][["delta_wt"]][, c("deseq_logfc", "deseq_adjp")]
comp_wt_df <- tables[["data"]][["comp_wt"]][, c("deseq_logfc", "deseq_adjp")]
delta_comp_df <- tables[["data"]][["delta_comp"]][, c("deseq_logfc", "deseq_adjp")]

pknf_cfg <- circos_prefix(annot, name="pknf", cog_column = "COGFun",
                          start_column = "start", stop_column = "stop",
                          strand_column = "strand", chr_column = "chromosome")

## 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/pknf.conf with a reasonable first approximation config file.

## Error in path.package("hpgldata"): none of the packages are loaded

pknf_kary <- circos_karyotype(
  pknf_cfg, fasta = "~/libraries/genome/mycobacterium_tuberculosis_cdc1551.fasta")

## Error in eval(expr, envir, enclos): object 'pknf_cfg' not found

plus_minus <- circos_plus_minus(pknf_cfg, width = 0.06, thickness = 40)

## Error in eval(expr, envir, enclos): object 'pknf_cfg' not found

wt_exprs_heat <- circos_heatmap(pknf_cfg, means, colname = "wt", basename = "mean_wt",
                             outer = plus_minus)

## Error in eval(expr, envir, enclos): object 'pknf_cfg' not found

delta_wt_hist <- circos_heatmap(pknf_cfg, delta_wt_df, colname = "deseq_logfc",
                             basename = "delta_wt", outer = wt_exprs_heat)

## Error in eval(expr, envir, enclos): object 'pknf_cfg' not found

comp_wt_hist <- circos_heatmap(pknf_cfg, comp_wt_df, colname = "deseq_logfc",
                            basename = "comp_wt", outer = delta_wt_hist)

## Error in eval(expr, envir, enclos): object 'pknf_cfg' not found

delta_comp_hist <- circos_heatmap(pknf_cfg, delta_comp_df, colname = "deseq_logfc",
                               basename = "delta_comp", outer = comp_wt_hist)

## Error in eval(expr, envir, enclos): object 'pknf_cfg' not found

finish <- circos_suffix(pknf_cfg)

## Error in eval(expr, envir, enclos): object 'pknf_cfg' not found

made <- circos_make(pknf_cfg, target = "pknf")

## Error in eval(expr, envir, enclos): object 'pknf_cfg' not found

The following is pulled from a Pseudomonas tiling of specific genes of interest; I have a feeling that the various STPKs may want to be added in a similar fashion

other <- c("gene1652095", "gene1652097", "gene1652099", "gene1652101",
           "gene1652103", "gene1652105", "gene1652107", "gene1652109",
           "gene1652111", "gene1652113", "gene1652115", "gene1652117",
           "gene1652119", "gene1652121", "gene1652123", "gene1652125",
           "gene1652127", "gene1652129", "gene1652131", "gene1652133",
           "gene1652135", "gene1652137", "gene1652139", "gene1652141",
           "gene1652143", "gene1652145", "gene1652147", "gene1652149",
           "gene1652151", "gene1652153", "gene1652155", "gene1652157",
           "gene1652159", "gene1652161", "gene1652163", "gene1652165")
## also BLUE
pel <- c("gene1654787", "gene1654789", "gene1654793",
         "gene1654795", "gene1654797", "gene1654799")
extra[down_orn, "flag"] <- "down_orn"
extra[up_orn, "flag"] <- "up_orn"
extra[effectors, "flag"] <- "effectors"
extra[other, "flag"] <- "other"
extra[pel, "flag"] <- "pel"
keep_idx <- extra[["flag"]] != ""
extra <- extra[keep_idx, ]
summary(as.factor(extra$flag))
tile_colors <-        c("AA0000",   "0000AA",    "00AA00", "00AA00", "0000AA")
names(tile_colors) <- c("down_orn", "effectors", "other",  "pel",    "up_orn")


pa13_tiles <- circos_tile(extra, colname="flag", colors=tile_colors, width=0.125,
                          padding=0, thickness=125, margin=0.00,
                          cfgout=circos_pa13, basename="tile", outer=pa13_st_wtmt)