"genes_with_missing_exons" not present in input gtf file

[Lil-Psilocybe]

Hello!

Thanks for a great tool! I successfully installed geneext and ran the test_data example without issue, however when working with my own data I got the same error as in issue #10 and got a lot of my genes without exons.

The genes in question came from a trinity assembly and are in the genes_with_missing_exons.txt like so:
head $geneext/genes_with_missing_exons.txt
TRINITY_DN26089_c0_g1_1
TRINITY_DN1846_c6_g1_1
TRINITY_DN11332_c0_g1_1
TRINITY_DN4372_c1_g1_1
TRINITY_DN9737_c0_g1_1
TRINITY_DN263_c7_g1_1
TRINITY_DN36841_c1_g1_1
TRINITY_DN4199_c0_g1_1
TRINITY_DN14827_c0_g1_1

When screening these genes in my gtf file, they do not return a result. However when I screen them without the final numbered suffix (i.e. TRINITY_DN26089_c0_g1, TRINITY_DN1846_c6_g1, TRINITY_DN11332_c0_g, etc), there is returned result with each of these genes having multiple exons.

I generated the input BAM file with the same gtf I am using for geneext, so I'm unsure where the problem could be.

I am working with SPLiT-Seq data from Parse BioSciences if that helps and the all_genes.csv file from the STAR alignment output does not contain the genes with this suffix either. Any input is greatly appreciated thanks!

[zolotarovgl]

Dear @Lil-Psilocybe,
thank you for raising this issue.
Most if not all issues come from a weirdly formatted input .gtf file - sometimes, there is no corresponding "gene" feature or missing exon features etc.
The reason why the GeneExt returns the names of those genes with the suffix _1 is probably due to an internal handling of IDs ( this may happen when the transcripts and genes have the same ids - i.e. both the gene and the transcript are called gene_1...)
What usually works out of the box, is trying to fix the .gtf with AGAT tool (https://agat.readthedocs.io/en/latest/agat_how_does_it_work.html)
Please, try running agat_convert_sp_gxf2gxf.pl on your input annotation before using GeneExt.

Also, I will greatly appreciate if you could paste some example genes "without exons" here from your gtf file to see what may be causing this error.

[Lil-Psilocybe]

Hi @zolotarovgl! Thank you so much for the quick reply over the holidays, i really appreciate it.

The gtf i'm using was actually generated from a gff using agat's agat_convert_sp_gff2gtf.pl and after checking out #10, I used agat_convert_sp_gxf2gxf.pl on both the original gff and non-working gtf, but no avail.

Interestingly, it appears to be the same genes that are genes_with_missing_exons.txt from running agat on either the gff and gtf, which tells us something useful.

Here's some output from 3 grep/less searches in the updated annotation:
some_gene_annots_with_missing_exons.txt.
It looks like I have some extra (even excessive) tags in the features column that may be contributing to this. When looking in my gtf for genes NOT in genes_with_missing_exons.txt they did not look noticeably different:
some_good_genes_with_exons.txt

We're updating gene annotations for a nonmodel species and need these tags for various troubleshooting.

Please let me know what you think and thank you for your time.

[Lil-Psilocybe]

Hello again! I wanted to update on this.

I was rerunning the data after using agat's agat_convert_sp_gxf2gxf.pl and then kept getting another error:
[E::hts_idx_push] Chromosome blocks not continuous [E::sam_index] Read '01_01_01__R__1_1_1__CATCATCC_AACGTGAT_AACGTGAT__ACCGTAGGAT__251210Cs__OH_@A01204:346:H3GFJDRX7:2:2137:17815:20290' with ref_name='Lachesis_group0__68_contigs__length_203855924', ref_length=203855924, flags=0, pos=128287724 cannot be indexed samtools index: failed to create index for "combined_bam_files.bam"

After doing some digging about [E::hts_idx_push] Chromosome blocks not continuous errors, it seemed the bam file I was working with was not sorted by coordinate, but rather by quername.

I ran samtools sort on my bam file and now geneext seems to be running fine! Fingers crossed, but it has been running for about 21 hours now so we'll see!

Will update when run finishes. Thank you!

[zolotarovgl]

Yes, @Lil-Psilocybe. I have never seen unsorted input bam files - usually both STARsolo and cellranger output the ones sorted by the position.
Thanks for spotting this, I will add this to the user manual.
Regarding "genes without exons" - I also got this error while working with the other genomes - it seems to be the error of gffutils - it infers the wrong gene IDs somehow by appending the _1 prefix (perhaps to make them different from the transcript IDs?)
Some genomes have nique features associated with genes and transcripts (geneA - transcriptA), while the others have the same id for both (gene_id "A", transcript_id "A") - this may be the reason.

For now, I have fixed this issue by doing some additional checks.
I would still appreciate if you could send me an example from the original .gtf file (the first chromosome is fine), so I can figure out why it has happened in your case.

Thank you for staying patient and not giving up on the tool!
Let me know if you spot any errors so I can adrress them.

Cheers,
Grisha

[Lil-Psilocybe]

Hello hello!

The run finished a couple days ago on the 17th. The BAM file I'm working with is really big and so the run took longer than I would have thought but it ran fine, which is great to see!

I wonder if the naming of isoforms may contributing to this discrepancy too (geneA - transcriptA, geneA - transcriptB, etc). My data is a bit funky in that it is "frankensteined" from combined sequences of a BRAKER run and also a trinity alignment; we are hoping to better capture some real genes that are not in the current assembly, but that we know are there.

Please find below ch1 from the original gtf I used that sparked the initial error.

Thanks for a necessary tool for sc/sn work! The community appreciates it.

ch1.gtf.txt