Hi! What does rnaseq use to annotate the reference transcriptome? I want to quantify prokaryotic data with Salmon (which requires a reference transcriptome) but according to this: https://nf-co.re/rnaseq/latest/docs/usage/ , the genome is annotated using featureCounts.

From the doc page you are quoting:

Since v3.0 of the pipeline, featureCounts is no longer used to perform gene/transcript quantification, however it is still used to generate QC metrics based on biotype information available within GFF/GTF genome annotation files.

(Emphasis mine.) The default quantification is done by Salmon. You don't need to supply a transcriptome file if you have a genome + annotation (but you may; otherwise it will be autogenerated from genome and GTF/GFF). So far, so good. The tricky part, though, with bacteria is that the GTF often does not contain the

exon

entries and

transcript_id

attributes that the pipeline expects by default (often there are only

CDS

entries). While you could use

CDS

instead of

exon

(I think there are parameters), you still need the transcripts annotated (they are basically equivalent to genes; the concept of operons is not something standard tools know about, and you they would not be annotated either way), Your best bet is to either find a transcript that conforms to the expectations (Ensembl may be your best bet, if they have your species/strain) or make one yourself with the use of some scripting. 🤞

If I want to run Salmon is mapping mode instead of selective alignment mode (which requires the SAM/BAM alignment files + genome to act as the decoy) I would need a transcriptome. But I agree that salmon cannot account for operons in bacteria

If your GTF has transcript information and you have a compatible genome FASTA, the pipeline should automatically generate a transcriptome file, which would also be used for Salmon in pseudo-aligner mode. I prefer to map to the genome, because the QC information is much better that way, but with bacterial genomes (which have no intronic and almost no intergenic regions) the difference may not be that large.

From what I understand, the issue for bacteria is not whether a transcriptome can be generated, but that standard bacterial annotations (if a GTF/GFF3 file is available) lack expected transcript structure and biologically true transcript annotations (e.g. operons, UTRs, various RNAs). As a result, the autogenerated “transcriptome” is effectively CDS-as-transcript, which Salmon can quantify but that does not resolve operons mapping issue (dense, overlapping transcripts). Unfortunately, I believe the problem comes down to the available annotated transcripts being not fully representative. I understand the issue a bit better now, I'll try to work around it!

As far as I know, there is no actual annotation for operons in any GTF, at least I've never seen one. In many cases it is not even known or certain if certain genes form an operon. But you can do quantification of RNAseq results at the gene level (you just call each gene also a transcript, and if you don't know the UTRs, the CDS boundaries can serve as the gene/transcript boundaries -- the difference is typically pretty small). If you have a bunch of neighboring genes with similar expression, you are likely looking at an operon, but that knowledge is not really necessary to interpret the results when comparing gene expression between sample groups. There may be reads that start in one gene and go into the next (if both are in an operon), but in the grand scheme of things those are not that many and Salmon will probably assign them to one or the other.