This repository uses conda to manage software installations. You can find operating system-specific instructions for installing miniconda here.
conda env create --name microberna --file environment.yml
conda activate microberna
Once the base environment is installed, the analysis is automated with snakemake.
snakemake -s best_genome_reference.snakefile -j 1 --use-conda --rerun-incomplete
Snakemake can parallelize job submission and modulate resource usage (RAM, CPUs). We used the command below in a slurm cluster, but other cluster engines are also supported.
snakemake -s best_genome_reference.snakefile -j 16 --use-conda --rerun-incomplete --latency-wait 15 -
-resources mem_mb=200000 --cluster "sbatch -t 720 -J comp -p bmm -n 1 -N 1 -c {threads} --mem={resources.mem_mb}" -k
The goal of this repository is to generate an automated, lightweight, and generalized pipeline to create a compendium of isolate bacterial and archaeal RNA-seq data. In model organisms like human and mouse, compendia like recount2 (and recount3), GTEx, and The Cancer Genome Atlas (TCGA) have been rich community resources that have allowed researchers to interface with and draw insights from gene expression data. Standard pre-processing enables integrated analysis across data sets derived from different studies, tissues, or organisms. This repository aims to produce a compendium that integrates the thousands of publicly available isolate bacterial and archaeal RNA-seq data via uniform pre-processing.
While the exact output formats will become more clear as the pre-processing pipeline is finalized, the pipeline aims to produce annotation-dependent and annotation-independent products.
Annotation-dependent products
sourmash gather
to determine the taxonomic profile of each isolate RNA-seq data set using genomes in the GTDB database (rs202) and the human transcriptome (GCF_000001405.39_GRCh38.p13) as references.tximport
to produce ortholog-summarized counts for each RNA-seq sample. Eggnog provides annotations for orthologs in COG, GO, KEGG, EC, CAZy, BiGG Reaction, and PFAM databases.Annotation-independent products
As of January 2022, there were 54,445[^samples] publicly available bacterial and archaeal RNA seq samples from 5,751[^experiments] experiments on the Sequence Read Archive (SRA). These samples were attributable to 1,722 distinct user-supplied organism names. The median number of times an organism name was observed was 7.
Using text-based matching, we converted organism names in SRA metadata to species names in the GTDB taxonomy. Through text-based matching alone, 3,274 RNA-seq samples derived from 344 organism names did not have a reference genome in the GTDB database. 2,524 GTDB species were represented in the 51,171 samples with matches, with a median of 29 observations per species. The majority of samples were Bacterial (50,182 Bacterial, 989 Archaeal).
A major challenge associated with building this compendium is the selection of reference sequences and the integration of annotations between reference sequences. Here, we focus on generating a reference transcriptome, as we have elected to perform transcript quantification using Salmon because it is lightweight and scales to thousands of samples.
The reference transcriptome is an important component of transcriptome quantification. If the reference does not contain all of the sequences that are contained within the sample, quantification will be incomplete (see here and here). This is nicely illustrated in this figure by Lisa K. Johnson (original notebook here). Johnson demonstrates that has Jaccard similarity decreases between an RNA-seq sample and the reference used to quantify it, mapping rates decrease. As such, it is important to have both a complete (all genes present) and similar (high percent identity) reference transcriptome with regards to each sample.
Given the importance of having a representative and accurate transcriptome, one goal of this repository is to determine the best approach for reference transcriptome generation for any publicly available bacterial or archaeal RNA-seq sample. Thus far, we have identified two approaches for reference transcriptome generation that we plan to test: a single pangenome approach and a best genome approach.
The single pangenome approach will generate a single pangenome using all "reference" genomes for a species. We have selected GTDB rs202 as our reference genome database, as it encompasses some unculturable organisms (unlike RefSeq), but is still quality controlled more extensively than e.g. GenBank. Annotated features within this pangenome (CDS, tRNA, rRNA, tmRNA, ncRNA) will be used to create a reference transcriptome that will be used to quantify all RNA-seq libraries from this species. We suspect that this approach might be successful at capturing microbes that have mixes of genes (e.g. accessory elements) which have not previously be observed in conjunction in a single genome before; a similar method was successful in the analysis of RNA-seq samples from Staphylococcus aureus strains without reference genomes
The best genome approach will select the best genome and use it as the reference transcriptome for each RNA-seq library.
To achieve this, we will use sourmash gather
to identify the reference genome in GTDB that best covers the k-mers in a given RNA-seq library.
Here, "covers" refers to a minimum set cover (see here, or the minimum number of reference genomes in a database needed to cover the k-mers in the RNA-seq data set.
Using this information, we will download this genome and identify, index, and annotate its open reading frames.
To harmonize annotations across data sets, we will rely on ortholog annotations produced by e.g. eggnog.
One benefit of this approach is the direct use of model organism transcriptomes (e.g. Escherichia coli k12, Pseudomonas aeruginosa PAO1 or PA14), which could be exported as sub-compendia and would better integrate with existing and well-known community annotation resources.
Prior to evaluation, we can foresee different challenges associated with each approach.
Challenges with single pangenome approach:
Challenges with best genome approach:
Three studies representing 29 samples have been submitted to the SRA for F. prausnitzii A2-165. While some of the samples are contaminated with human or other organisms, only strain A2-165 was detected among F. prausnitzii reads. These RNA-seq samples are a good test case to determine whether the pangenome approach works as well as the best genome approach. We selected three samples, each from a different study, and used these to investigate this question.
SRA identifier | Annotated strain |
---|---|
ERX4307280 | A2-165 |
SRX10245671 | A2-165 |
SRX3847835 | A2-165 |
See results here
Five studies representing 21 samples have been submitted to the SRA and labelled as Bacillus pumilus. We selected one sample from each of these studies to test. Four of the samples were annotated by the submitters as B. pumilus, illustrating the necessity of determining the reference transcriptome by sequence comparison against a database. The difference in both the total number of Bacillus strains detected, as well as the identity of those strains (not shown), suggests that each transcriptome captures a slightly different strain, and that none of those strains are not represented in a single reference genome. The mix of strains present makes these samples a good test case to determine whether the pangenome approach allows more reads to be recovered over the best genome approach.
SRA identifier | Top annotated strain | Total Bacillus strains annotated |
---|---|---|
SRX4378884 | Bacillus altitudinis strain=BA06 | 5 |
SRX5179263 | Bacillus altitudinis strain=BA06 | 6 |
SRX5678458 | Bacillus altitudinis strain=BA06 | 8 |
SRX7088415 | Bacillus altitudinis strain=BA06 | 5 |
Recently, compendia of P. aeruginosa were developed by mapping RNA-seq samples against PAO1 or PA14 references using salmon. We identified 6 samples analyzed by this study to use as test cases for this compendium.
SRA identifier | Annotated strain |
---|---|
ERX3558803 | PAO1 |
SRX540112 | PAO1 |
SRX589549 | PA14 |
SRX4624095 | PA14 |
SRX396878 | clinical isolate |
SRX5123759 | clinical isolate |
P. aeruginosa is a good test case for this pipeline because:
Salmon in part uses k-mers in RNA-seq reads and transcripts in the reference transcriptome to quantify transcript abundance in a sample. One user-controllable parameter is the minimum k-size for an acceptable match. We need to test how this parameter behaves in relation to the read size and the similarity between the RNA-seq sample and the selected reference transcriptome.
The salmon documentation states:
This will build the mapping-based index, using an auxiliary k-mer hash over k-mers of length 31. While the mapping algorithms will make used of arbitrarily long matches between the query and reference, the k size selected here will act as the minimum acceptable length for a valid match. Thus, a smaller value of k may slightly improve sensitivty. We find that a k of 31 seems to work well for reads of 75bp or longer, but you might consider a smaller k if you plan to deal with shorter reads. Also, a shoter value of k may improve sensitivity even more when using selective alignment (enabled via the –validateMappings flag). So, if you are seeing a smaller mapping rate than you might expect, consider building the index with a slightly smaller k.
Read length: Estimating average read length by dividing SRA metadata Total Bases by Total Spots,
It may be prudent to require read length > 31 base pairs after trimming for inclusion in the compendia. In either case, we need to test k size across read lengths to see if there are biases.
ANI between genome in the RNA-seq sample and the reference transcriptome: Decreased relatedness between the organism represented by the RNA-seq sample and that represented by the reference transcriptome should decrease mapping rates, but decreasing the minimum acceptable k-mer length could improve mapping rates in this case. We need to test whether decreasing the k size increases mapping rate without increasing off target mapping, for both RNA-seq samples that are closely and distantly related to the selected reference.
It would be great if there is a universal best k-size to use.
I will probably evaluate this after I have gather results for a few species. An ideal data set would have varied jaccard similarity relative to the reference pangenome (or other reference transcriptome used).
General ideas
Single pangenome ideas
[^samples]: Search criteria on the SRA ("Bacteria"[Organism] OR "Archaea"[Organism]) AND ("biomol rna"[Properties] AND "platform illumina"[Properties] AND "filetype fastq"[Properties]), and additionally filtered on Library Source "Transcriptomic", Library Strategy "RNA-Seq" or "OTHER", and Library Selection "cDNA", "other", "RANDOM", "unspecified", "RT-PCR", "RANDOM PCR", "PCR"
[^experiments]: Measured as distinct study accessions