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.

3.0.1 Read in orthogroup information

The orthologs produced by orthofinder appear to pick up single-groups for 3343 genes, which seems pretty good, I am assuming the remainder of genes are multi-gene families which will also be nicely identified.

single_orthologs <- as.data.frame(readr::read_tsv("orthologs/outputs/50orthofinder/orthogroups_single_named.tsv"))

## New names:
## • `` -> `...4`

## Warning: One or more parsing issues, call `problems()` on your data frame for details, e.g.:
##   dat <- vroom(...)
##   problems(dat)

## Rows: 3343 Columns: 4
## ── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (3): Orthogroup, mtuberculosis_h37rv, ...4
## lgl (1): mycobacterium_tuberculosis_cdc1551
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

rownames(single_orthologs) <- single_orthologs[["Orthogroup"]]
single_orthologs <- single_orthologs[, c(2, 4)]
colnames(single_orthologs) <- c("h37rv", "cdc1551")

## There are a bunch of 1:many mappings and many:many mappings it turns out!
## and I need to reformat the text file containing them, supido orthofinder
## added some random characters in the lists of genes which make it weird to parse.
all_orthogroups <- as.data.frame(readr::read_tsv("orthologs/outputs/50orthofinder/orthogroups_all_named.tsv"))

## New names:
## • `` -> `...2`

## Warning: One or more parsing issues, call `problems()` on your data frame for details, e.g.:
##   dat <- vroom(...)
##   problems(dat)

## Rows: 3551 Columns: 2
## ── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (2): Orthogroup, ...2
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Ok, so I need to do a little work to clean up the (1|many):(1|many) mappings; until then let us just merge the 1:1 mappings into the CDC1551 annotations.

I think I will need to get the NIH IDs from the genbank or gff file before merging to the microbesonline annotations. An important caveat I did not consider before: the ortholog searches were explicitly on CDS/proteins and my general annotations are not.

h37_cds <- load_gff_annotations("orthologs/mtuberculosis_h37rv_all.gff", type = "CDS")

## Returning a df with 30 columns and 3906 rows.

h37_cds_singles <- merge(h37_cds, single_orthologs, by.x = "protein_id",
                         by.y = "h37rv", all.x = TRUE)

cdc_cds <- load_gff_annotations("orthologs/mycobacterium_tuberculosis_cdc1551_cds.gff", type = "CDS")

## Returning a df with 18 columns and 4189 rows.

h37_ncbi_cdc <- merge(h37_cds_singles, cdc_cds, by.x = "cdc1551", by.y = "protein_id")

single_gene_mapping <- h37_ncbi_cdc[, c("locus_tag.x", "locus_tag.y", "gene.x",
                                        "product.x", "gene.y", "product.y")]
colnames(single_gene_mapping) <- c("h37rv_tag", "cdc1551_tag", "h37_gene",
                                   "h37_product", "cdc_gene", "cdc_product")

3.1 Add the single gene maps to the annotations

Oh, the way I did the above, this will almost certainly end up with redundant columns; but since those columns are coming from different sources it will be a good control to ensure that everything worked properly (e.g. if they prove not to be redundant then we know there were shenanigans in my ortholog mapping).

annot <- merge(annot, single_gene_mapping, by.x = "row.names",
               by.y = "cdc1551_tag", all.x = TRUE)
rownames(annot) <- annot[["Row.names"]]
annot[["Row.names"]] <- NULL

4 A slightly random aside, rRNA and 6S

Given the pecularities we observed vis a vis the Pseudomonas reads, there is explicit interest from April to know how well the ribosomal depletion went; so I did a little poking and found that the ribosomal genes are named rrf(5S), rrs(16S), and rrl(23S). When I make a bacterial ribosomal database I usually like to include a few tRNAs and 6S, even though 5S/6S/tRNAs should be size-selected out of the libraries pretty consistently. It turns out the biology of RNAs which interact with σ in Mtb is kind of wild and there is no known 6s, but a somewhat similar RNA named ‘Ms1’.

Anyhow, I made a rRNA database for Mtb and here is a query against it to get a sense of likely rRNA reads in the samples.

cd preprocessing
start=$(pwd)
for i in $(/bin/ls -d * | grep -v Undetermined); do
    cd $i
    cyoa --method hisat --libtype rRNA --species mtuberculosis_h37rv \
         --input $(/bin/ls *-trimmed.fastq.xz | tr '\n' ':' | sed 's/:$//g')
    cd $start
done

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

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

7 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.1 Check the rRNA content

The following does not work because I fell into the assumption that these are H37Rv samples again! Doofus.

rrna_loci <- c("Rvnr01", "Rvnr02", "Rvnr03")
rrna_expt <- expt[rrna_loci, ]

## Subsetting on genes.

## remove_genes_expt(), before removal, there were 4245 genes, now there are 0.

## There are 9 samples which kept less than 90 percent counts.

##          WT_1          WT_2          WT_3 Mutant_pknF_1 Mutant_pknF_2 Mutant_pknF_3   Comp_pknF_1   Comp_pknF_2   Comp_pknF_3 
##             0             0             0             0             0             0             0             0             0

While I was typing the above incorrectly, the alignments started finishing the query against my Mtb rRNA database. It looks like ~ 3% of the reads map to rRNA.

8 varpart

varpart <- simple_varpart(expt)

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

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

## 
## Total:42 s

varpart

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

9 Kraken matrix

genus_expt <- create_expt(written[["new_file"]],
                          file_column = "krakenmatrix", file_type = "table")

## Reading the sample metadata.

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

## Warning in first_rownames != current_rownames: longer object length is not a multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, : The file:
## /home/trey/sshfs/scratch/atb/rnaseq/mycobacterium_tuberculosis_2023/preprocessing/WT_2/outputs/05kraken_bacteria/kraken_report_matrix.tsv has mismatched
## rownames.

## Warning in first_rownames != current_rownames: longer object length is not a multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, : The file:
## /home/trey/sshfs/scratch/atb/rnaseq/mycobacterium_tuberculosis_2023/preprocessing/WT_3/outputs/05kraken_bacteria/kraken_report_matrix.tsv has mismatched
## rownames.

## Warning in first_rownames != current_rownames: longer object length is not a multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, : The file:
## /home/trey/sshfs/scratch/atb/rnaseq/mycobacterium_tuberculosis_2023/preprocessing/Mutant_pknF_1/outputs/05kraken_bacteria/kraken_report_matrix.tsv has
## mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, : The file:
## /home/trey/sshfs/scratch/atb/rnaseq/mycobacterium_tuberculosis_2023/preprocessing/Mutant_pknF_2/outputs/05kraken_bacteria/kraken_report_matrix.tsv has
## mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, : The file:
## /home/trey/sshfs/scratch/atb/rnaseq/mycobacterium_tuberculosis_2023/preprocessing/Mutant_pknF_3/outputs/05kraken_bacteria/kraken_report_matrix.tsv has
## mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, : The file:
## /home/trey/sshfs/scratch/atb/rnaseq/mycobacterium_tuberculosis_2023/preprocessing/Comp_pknF_1/outputs/05kraken_bacteria/kraken_report_matrix.tsv has
## mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, : The file:
## /home/trey/sshfs/scratch/atb/rnaseq/mycobacterium_tuberculosis_2023/preprocessing/Comp_pknF_2/outputs/05kraken_bacteria/kraken_report_matrix.tsv has
## mismatched rownames.

## Warning in first_rownames != current_rownames: longer object length is not a multiple of shorter object length

## Warning in read_counts_expt(sample_ids, filenames, countdir = countdir, : The file:
## /home/trey/sshfs/scratch/atb/rnaseq/mycobacterium_tuberculosis_2023/preprocessing/Comp_pknF_3/outputs/05kraken_bacteria/kraken_report_matrix.tsv has
## mismatched rownames.

## Warning in create_expt(written[["new_file"]], file_column = "krakenmatrix", : There are some NAs in this data, the 'handle_nas' parameter may be required.

## Matched 191 annotations and counts.

## Bringing together the count matrix and gene information.

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

## The final expressionset has 191 features and 9 samples.

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

## A heatmap of pairwise sample distances ranging from: 
## 4945.06863345337 to 102810.646015848.

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

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

plot_libsize(genus_expt)

## Library sizes of 9 samples, 
## ranging from 371,301 to 2,637,312.

head(exprs(genus_expt))

##                  WT_1 WT_2 WT_3 Mutant_pknF_1 Mutant_pknF_2 Mutant_pknF_3 Comp_pknF_1 Comp_pknF_2 Comp_pknF_3
## Acetobacter        96   29   66           116            81            53          27         146         143
## Achromobacter      26   20   15            34            28            21           6          49          59
## Acidiferrobacter   26    8   12            18            21            10           2          23          31
## Acidihalobacter    14    5   11            23            38            11           4          38          40
## Acinetobacter     819  314  448           770           691           477         147        1130         962
## Actinobacillus     26   20   26            36            36            10           4          28          32

exprs(genus_expt)["Pseudomonas", ]

##          WT_1          WT_2          WT_3 Mutant_pknF_1 Mutant_pknF_2 Mutant_pknF_3   Comp_pknF_1   Comp_pknF_2   Comp_pknF_3 
##       1509718        546307       1020882       1635517       1502338        902653        341293       2477879       2388934

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

## Deleting the file excel/de_tables-v20231207.xlsx before writing the tables.

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

## Deleting the file excel/de_sig-v20231207.xlsx before writing the tables.

10.0.1 Test batch-in-model

Given the PCA above, it seems at least possible that including batch in model might provide more reliable results than sva. So, the following block does that. Spoiler alert: I don’t think it helped.

de_batch <- all_pairwise(expt, keepers = keepers, filter = TRUE, model_batch = TRUE)

## 
##         wt      delta complement 
##          3          3          3 
## 
## b1 b2 b3 
##  3  3  3

de_batch

## A pairwise differential expression with results from: basic, deseq, ebseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: batch in model/limma.

## The primary analysis performed 15 comparisons.

tables_batch <- combine_de_tables(de_batch, keepers = keepers,
                                  excel = glue("excel/de_tables_batch-v{ver}.xlsx"))

## Deleting the file excel/de_tables_batch-v20231207.xlsx before writing the tables.

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

tables_batch

## A set of combined differential expression results.

##                          table deseq_sigup deseq_sigdown edger_sigup edger_sigdown limma_sigup limma_sigdown
## 1         wt_vs_delta-inverted          31             1          16             3           2             1
## 2    wt_vs_complement-inverted          17             7          14             6           2             1
## 3 complement_vs_delta-inverted          12             9          11             8           1             1

sig_batch <- extract_significant_genes(tables_batch, according_to = "deseq",
                                       excel = glue("excel/de_sig_batch-v{ver}.xlsx"))

## Deleting the file excel/de_sig_batch-v20231207.xlsx before writing the tables.

sig_batch

## A set of genes deemed significant according to deseq.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##            deseq_up deseq_down
## delta_wt         31          1
## comp_wt          17          7
## delta_comp       12          9

10.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"

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

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

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

12.1 Collect the STPKs

stpk_operons <- c("pknA", "pknB", ## ~15k, perhaps ending a 7 gene operon (fhaA, fhaB, pstP, rodA, rbpA, pknA, pknB)
                 ## No pknC
                 "pknD", ## ~ 1038500, perhaps 2 gene operon with pstS2
                 "pknE", ## ~ 2000000, perhaps with vapB34, vapC34, Rv1742, pknE
                 "pknF", ## ~ 2000000, close to E, pknF, Rv1747, Rv1748
                 "pknG", ## ~ 2000000, Rv0412c, hlnH, pknG
                 "pknH", ## ~ 1415000, maybe 8 genes, Rv1273c, Rv1272c, Rv1271c, lprA, Rv1269c, Rv1268c, embR, pknH
                 "pknI", ## ~ 3222000, ffh, Rv2915c, pknI, Rv2913c, Rv2912c
                 "pknJ", ## ~ 2345000, Rv2082, Rv2083, Rv2084, Rv2085, Rv2086, Rv2087, pknJ
                 "pknK", ## ~ 3442500, pknK, Rv3079c
                 "pknL",
                 "lppH" ## when searching for pknM, lpphH, Rv3577, arsB2 4015000
                 )
mapped_stpk <- single_gene_mapping[["h37_gene"]] %in% stpk_operons
mapped_ids <- single_gene_mapping[mapped_stpk, ]
mapped_ids

##      h37rv_tag cdc1551_tag h37_gene                                        h37_product cdc_gene                                            cdc_product
## 14     Rv0014c      MT0017     pknB               serine/threonine-protein kinase PknB     <NA>                        serine/threonine protein kinase
## 15     Rv0015c      MT0018     pknA               serine/threonine-protein kinase PknA     <NA>                        serine/threonine protein kinase
## 346    Rv0410c      MT0423     pknG               serine/threonine-protein kinase PknG     <NA>                        serine/threonine protein kinase
## 778    Rv0931c      MT0958     pknD               serine/threonine-protein kinase PknD     <NA>                        serine/threonine protein kinase
## 1051   Rv1266c      MT1304     pknH               serine/threonine-protein kinase PknH     <NA>                        serine/threonine protein kinase
## 1462    Rv1743      MT1785     pknE               serine/threonine-protein kinase PknE     <NA>                        serine/threonine protein kinase
## 1465    Rv1746      MT1788     pknF               serine/threonine-protein kinase PknF     <NA>                        serine/threonine protein kinase
## 1755    Rv2088      MT2149     pknJ transmembrane serine/threonine-protein kinase PknJ     <NA>                        serine/threonine protein kinase
## 1832    Rv2176      MT2232     pknL               serine/threonine-protein kinase PknL     <NA>                        serine/threonine protein kinase
## 2474   Rv2914c      MT2982     pknI               serine/threonine-protein kinase PknI     <NA>                        serine/threonine protein kinase
## 2614   Rv3080c      MT3165     pknK               serine/threonine-protein kinase PknK     <NA> serine/threonine protein kinase / MalT-related protein
## 3040    Rv3576      MT3681     lppH                                   lipoprotein LppH     <NA>                                   hypothetical protein

mapped_ids[["cdc1551_tag"]]

##  [1] "MT0017" "MT0018" "MT0423" "MT0958" "MT1304" "MT1785" "MT1788" "MT2149" "MT2232" "MT2982" "MT3165" "MT3681"

## Oh no shit the reason I couldn't find them before was the difference between pknB and PknB, damnit.
## So that is annoying I likely did not need to do the silly ortholog map at all.

Let us add a column to the annotations noting the STPKs and potentially include them in the circos plot.

The current map is in order from outer->inner

• • strand genes colored by COGFun
• • strand genes
• WT gene expression on log scale, color scale chosen automatically as heatmap
• Δ/WT comparison on log scale, color scale chosen automatically as heatmap
• comp/WT comparison, Ibid.
• Δ/comp comparison, Ibid.
• ideally, tiles showing the location of the 11 STPKs and 1 rando which might also be a STPK.

annot[["chromosome"]] <- "AE000516"
annot[["stpk"]] <- "NULL"
annot[mapped_ids[["cdc1551_tag"]], "stpk"] <- "stpk"
stpk_rows <- annot[mapped_ids[["cdc1551_tag"]], ]

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.

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

## This should match the ideogram= line in pknf.conf

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

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

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

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

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

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

## Writing data file: circos/data/pknf_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.

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

## Assuming the input is a dataframe.

## Writing data file: circos/data/pknf_mean_wtwt_heatmap.txt with the mean_wtwt column.

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

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

## Assuming the input is a dataframe.

## Writing data file: circos/data/pknf_delta_wtdeseq_logfc_heatmap.txt with the delta_wtdeseq_logfc column.

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

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

## Assuming the input is a dataframe.

## Writing data file: circos/data/pknf_comp_wtdeseq_logfc_heatmap.txt with the comp_wtdeseq_logfc column.

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

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

## Assuming the input is a dataframe.

## Writing data file: circos/data/pknf_delta_compdeseq_logfc_heatmap.txt with the delta_compdeseq_logfc column.

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

tile_colors <- "0000AA"
names(tile_colors) <- c("stpk")
stpk_tiles <- circos_tile(pknf_cfg, stpk_rows, colname = "stpk",
                          colors = tile_colors, width = 0.125,
                          padding = 0, thickness = 125, margin = 0.00,
                          basename = "stpk", outer = delta_comp_hist)

## Writing data file: circos/data/pknfstpk_tile.txt with the stpkstpk column.

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

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

One of the primary reasons I usually do histograms instead of heatmaps: I do not really understand how to lay out the color range so that it is consistent across different ranges of logFC values (I suppose I could just normalize them to -1<=x<=1?

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