Unexpected behavior: Nextclade ignores CDS unlinked to gene if gene is in gff3

[anna-parker]

I'm not sure if this is really a bug or just nextclade not behaving as I would expect from the documentation.

@chaoran-chen created gene maps for the H5N1 segment M as follows:

.	.	CDS	26	784	.	+	.	gene="M1"
.	.	gene	26	1007	.	+	.	gene=M2;ID=gene-M2
.	.	CDS	26	51	.	+	.	gene=M2;ID=cds-M2;Parent=gene-M2
.	.	CDS	740	1007	.	+	.	gene=M2;ID=cds-M2;Parent=gene-M2

From the documentation this looked accurate to me. However, I now realized that this leads nextclade to ignore the M1 gene:

Removing the second line (with gene) resolved the issue:

.	.	CDS	26	784	.	+	.	gene="M1"
.	.	CDS	26	51	.	+	.	gene=M2;ID=cds-M2
.	.	CDS	740	1007	.	+	.	gene=M2;ID=cds-M2

My naive assumption is that nextclade ignores the CDS that is not linked to a gene if a gene is present - however I would not expect this behavior from the docs.

This can be tested with examples here: nextstrain/nextclade_data#217

[ivan-aksamentov]

@anna-parker What would be expected behavior in your opinion?

To clarify, what I tried to code is that if there's only genes present, we add corresponding CDSes, if there's only CDSes, we add corresponding genes. With the the half-genes-half-CDS thing as in the first example - I don't know what to do. Perhaps, we should iterate over CDSes and see if there's a linked gene, and add a gene if there isn't? And then symmetrically, see if there's lonely genes without CDS? (i.e. make the decision based on individual gene or CDS level, rather than globally for all entries).

The same GFF file can often be interpreted in many ways. If you have a clear spec for the behavior here, let me know and I'll try to tweak tings according to it.

I think we might also consider to drop this complexity entirely and to enforce correct GFFs - containing all required genes and all CDSes. But it may become a hassle for users, because lots of stuff in databases is non-compliant and they will have to fix it manually to make it work. We are also partially guilty ourselves here, because the data we used in Nextclade v1 was a mess.

But in any case, at certain point, Nextclade should not attempt to salvage the broken data. There should be a line we define where we say "it's too broken, we are not going to fix this".

[anna-parker]

Thanks! Yes I agree, now I've also read the GFF3 docs this isn't really a spec-compliant GFF3 - from our side it would have been better to use the genbank original GFF3 in this case.

I think I was confused by the docs and the statement: The fundamental unit for Nextclade is a single CDS. - so I thought that all CDS would be translated regardless of whether they are in the same region as the listed gene. But yeah I agree this was not a spec-compliant GFF3 and I should have checked it better before using it.

It probably would have been helpful if I had seen a warning from nextclade that a CDS in my GFF3 was being ignored or that my GFF3 wasn't fully correct and then hopefully I would have spotted this.

[ivan-aksamentov]

Yep. Some sanity checks would be nice to add sometimes. I transferred issue to nextclade software repo for better visibility.