Hi, I plan to run rnaseq pipeline with bacterial RNA-Seq data. The documentation says: Some bacterial GTFs from RefSeq contain gene entries with empty

transcript_id

fields, which can cause Salmon errors. If you encounter issues, check your annotation file and consider filtering problematic entries. The annotation GFF I have seemed to be generated using Prodigal and there is no

transcript_id

at all:

##gff-version 3
##sequence-region 1 1 2842373
1	Prodigal:002006	CDS	1	1410	.	+	0	ID=327701E_Td35405_00001;Name=dnaA;db_xref=COG:COG0593;gene=dnaA;inference=ab initio prediction:Prodigal:002006,similar to AA sequence:UniProtKB:O83047;locus_tag=327701E_Td35405_00001;product=Chromosomal replication initiator protein DnaA
1	Prodigal:002006	CDS	1415	1870	.	-	0	ID=327701E_Td35405_00002;inference=ab initio prediction:Prodigal:002006;locus_tag=327701E_Td35405_00002;product=hypothetical protein
1	Prodigal:002006	CDS	2015	5464	.	+	0	ID=327701E_Td35405_00003;eC_number=2.7.7.7;Name=dnaE;db_xref=COG:COG0587;gene=dnaE;inference=ab initio prediction:Prodigal:002006,similar to AA sequence:UniProtKB:Q9XDH5;locus_tag=327701E_Td35405_00003;product=DNA polymerase III subunit alpha
1	Prodigal:002006	CDS	5505	6086	.	+	0	ID=327701E_Td35405_00004;inference=ab initio prediction:Prodigal:002006;locus_tag=327701E_Td35405_00004;product=hypothetical protein

Could I simply change

ID

to

transcript_ID

?

No. You have a GFF3, which is different from a G*T*F in both the structure and the identifier names it uses. While a GTF has all of the information (e.g, for an exon) in one line (you can see which gene and transcript it belongs to via the

gene_id

and

transcript_id

), the GFF3 is hierarchical (more like a tree) -- each entity (exon, transcript, gene, etc.) has an

ID

, and most have also a

Parent

which contains the

ID

of whatever structure contains the entity. E.g., you have several `exon`s (labeled as such in column 3), and their

Parent

is an

mRNA

(labeled in column 3, this is similar to a transcript in the GTF, but GFF3s often distinguish between multiple types of transcripts). The

Parent

of the

mRNA

, in turn, is a

gene

, which is typically the top of the hierarchy, in which case it does not have a

Parent

. If you randomly replace

ID

with in such a file with something else, you will completely mess it up. You are not showing the whole file -- do you have lines with something like

mRNA

,

ncRNA

,

tRNA

,

rRNA

,

lnc_RNA

,

miRNA

, etc., or even

transcript

? If so, and if all of these lines have a non-empty

ID

, you don't have to worry about missing transcript IDs. I would run the pipeline with this GFF3 (use the

--gff

parameter) and see what happens. If it works, great, if you get an error, let us know what the message says and we'll take it from there.

Ah, sorry I messed it up. I meant to say GFF. And I just noticed that the

transcript_id

requirement is for GTF. I do have other features besides CDS:

2578 CDS
   1 repeat_region
   6 rRNA
   1 tmRNA
   44 tRNA

Thanks a lot for your detailed explanation! Sounds like I can use this GFF directly.

Yes, you should (ideally). The

rRNA

,

tmRNA

, and

tRNA

are your "transcripts" (maybe there are more types?).

repeat_region

is just an annotation, which can be ignored, and

CDS

is part of the coding sequence, i.e. similar to an

exon

, and I don't think these would be used for RNAseq at all. If these are really all the types you have, I see a problem, though, because I think you need "exons", "transcripts", and "genes". I'm not sure if or how the latest changes in the pipelines have changed these requirements. All you have is

CDS

, but these don't have

Parents

. In bacteria, an exon is the same as a transcript or a gene, because there is no splicing, and the CDS is just the open reading frame within this exon/transcript/gene, and since the exact transcription start and end are typically unknown (and very close by), it is usually considered the equivalent of the exon/transcript/gene. A worst case scenario is that you'll have to make an

exon

out of every

CDS

and create a

Parent

for it in form of an

mRNA

, as well as a

Parent

for that in form of a

gene

. But if the pipeline does what I hope it does, you don't have worry and it will just use the

CDS

for everything automatically. 🤞

There are no "parents".

CDS

is basically all it has. Would you recommend I run some extra annotation to generate a proper GFF that is more compatible with the rnaseq pipeline?

The older versions of the pipeline would require that, but I think that the newest version (when run with

-profile prokaryotic

) will know what to do with those seemingly "orphaned"

CDS

entries, so ideally it will be fine.

OK. I'll give it a try then. Thanks!

Yep, we've tried to improve prokaryotic operation, but it hasn't been widely tested and prokaryotic annotation files can be quite variable. The way we hope it should work is documented here: https://nf-co.re/rnaseq/3.23.0/docs/usage/#prokaryotic-genome-annotations

Hi @Jonathan Manning @Thomas Danhorn I see that in addition to the CDS entries in the gff, other features, such as tmRNA, tRNA, rRNA, are also captured and converted into transcripts.fasta. This is good, but I just wonder if this sentence in the documentation is not that stringent anymore:

• CDS features required: The annotation must contain CDS (coding sequence) features. The pipeline extracts transcripts from these.

Good point: https://github.com/nf-core/rnaseq/pull/1796