search

gpertea / gffread

GFF/GTF utility providing format conversions, region filtering, FASTA sequence extraction and more
MIT License
353 stars 39 forks source link

double free or corruption (!prev) and gff corruption with ribosomal slippage #121

Open alephreish opened 9 months ago

alephreish commented 9 months ago

I do not entirely understand what is happening, but in some cases CDSs with ribosomal slippage from NCBI cause double free or corruption (!prev):

$ gffread --version
0.12.7
$ curl https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/900/114/235/GCF_900114235.1_IMG-taxon_2636415975_annotated_assembly/GCF_900114235.1_IMG-taxon_2636415975_annotated_assembly_genomic.fna.gz | gzip -cd > genome.fna
$ curl https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/900/114/235/GCF_900114235.1_IMG-taxon_2636415975_annotated_assembly/GCF_900114235.1_IMG-taxon_2636415975_annotated_assembly_genomic.gff.gz | gzip -cd > genome.fna
$ gffread -y- -g genome.fna genome.gff | seqkit stat
double free or corruption (!prev)
file  format  type     num_seqs  sum_len  min_len  avg_len  max_len
-     FASTA   Protein     2,302  680,383       30    295.6    1,443
$
$ # This is caused by locus BM212_RS11370:
$
$ grep BM212_RS11370 genome.gff
NZ_FOSJ01000070.1   RefSeq  gene    18  1144    .   +   .   ID=gene-BM212_RS11370;Name=BM212_RS11370;gbkey=Gene;gene_biotype=protein_coding;locus_tag=BM212_RS11370;old_locus_tag=SAMN04488569_10702
NZ_FOSJ01000070.1   Protein Homology    CDS 18  303 .   +   0   ID=cds-WP_091898643.1;Parent=gene-BM212_RS11370;Dbxref=GenBank:WP_091898643.1;Name=WP_091898643.1;Note=programmed frameshift;Ontology_term=GO:0004803;exception=ribosomal slippage;gbkey=CDS;go_function=transposase activity|0004803||IEA;inference=COORDINATES: similar to AA sequence:RefSeq:WP_076611996.1;locus_tag=BM212_RS11370;product=IS3 family transposase;protein_id=WP_091898643.1;transl_table=11
NZ_FOSJ01000070.1   Protein Homology    CDS 303 1144    .   +   2   ID=cds-WP_091898643.1;Parent=gene-BM212_RS11370;Dbxref=GenBank:WP_091898643.1;Name=WP_091898643.1;Note=programmed frameshift;Ontology_term=GO:0004803;exception=ribosomal slippage;gbkey=CDS;go_function=transposase activity|0004803||IEA;inference=COORDINATES: similar to AA sequence:RefSeq:WP_076611996.1;locus_tag=BM212_RS11370;product=IS3 family transposase;protein_id=WP_091898643.1;transl_table=11
$ grep -v BM212_RS11370 genome.gff | gffread -y- -g genome.fna | seqkit stat
file  format  type     num_seqs  sum_len  min_len  avg_len  max_len
-     FASTA   Protein     2,436  718,882       30    295.1    1,511

Notwithstanding the crash, the CDS is translated correctly:

$ grep BM212_RS11370 genome.gff | gffread -y- -g genome.fna
>gene-BM212_RS11370
MSTRRPRRTYTEEFKKQIVDLHKAGKSRKEIIEEYDLTGSAFDKWVRQHNQTGSFKEKDNLTPEQKELKE
LRKANTQLQMENDILKQAALIFGRKFEVIKRNKHKYSISAMCRALEISRGSYYYEVIKKESDAKLEQAII
EEFAKSKNNYGTRKLKKRLKKRAFTVSRRRIGRIMNKFHLVSNYDKRSYKPQQSGVNQAKTENALNREFD
QEQPMKAIVTDLTYVKVANKWFYVCFILDLFNREIIGYSAGPNKTADLVLQALATVKGDLHTVKVFHTDR
GKEFDNHILDELLGTFDIVRSLSRKGNPYDNAVAESTYKSFKFEFVYGHTFHTLYELQVQLMDYVHWWNP
VRPHGSLDYETPIDYREGWEQEQSE
double free or corruption (!prev)
Aborted (core dumped)

If translation is not requested, there is no crash but gffread produces weird CDS features:

$ grep BM212_RS11370 genome.gff | gffread
NZ_FOSJ01000070.1   RefSeq  gene    18  1144    .   +   .   ID=gene-BM212_RS11370;
NZ_FOSJ01000070.1   RefSeq  CDS 18  1144    .   +   .   Parent=gene-BM212_RS11370
NZ_FOSJ01000070.1   RefSeq  CDS 18  303 .   +   0   Parent=gene-BM212_RS11370
NZ_FOSJ01000070.1   RefSeq  CDS 303 1144    .   +   2   Parent=gene-BM212_RS11370
$ grep BM212_RS11370 genome.gff | gffread | gffread
NZ_FOSJ01000070.1   RefSeq  gene    18  1144    .   +   .   ID=gene-BM212_RS11370;
NZ_FOSJ01000070.1   RefSeq  CDS 18  1144    .   +   .   Parent=gene-BM212_RS11370
NZ_FOSJ01000070.1   RefSeq  CDS 18  303 .   +   0   Parent=gene-BM212_RS11370
NZ_FOSJ01000070.1   RefSeq  CDS 18  1144    .   +   2   Parent=gene-BM212_RS11370
NZ_FOSJ01000070.1   RefSeq  CDS 303 1144    .   +   0   Parent=gene-BM212_RS11370
$ grep BM212_RS11370 genome.gff | gffread | gffread | gffread
NZ_FOSJ01000070.1   RefSeq  gene    18  1144    .   +   .   ID=gene-BM212_RS11370;
NZ_FOSJ01000070.1   RefSeq  CDS 18  1144    .   +   .   Parent=gene-BM212_RS11370
NZ_FOSJ01000070.1   RefSeq  CDS 18  303 .   +   0   Parent=gene-BM212_RS11370
NZ_FOSJ01000070.1   RefSeq  CDS 18  1144    .   +   2   Parent=gene-BM212_RS11370
NZ_FOSJ01000070.1   RefSeq  CDS 18  1144    .   +   0   Parent=gene-BM212_RS11370
NZ_FOSJ01000070.1   RefSeq  CDS 303 1144    .   +   1   Parent=gene-BM212_RS11370

Interestingly, there is a different transposase CDS in this example with ribosome slippage that does not cause the crash:

$ grep BM212_RS06860 genome.gff
NZ_FOSJ01000022.1   RefSeq  gene    35670   37183   .   +   .   ID=gene-BM212_RS06860;Name=BM212_RS06860;gbkey=Gene;gene_biotype=protein_coding;locus_tag=BM212_RS06860
NZ_FOSJ01000022.1   Protein Homology    CDS 35670   36288   .   +   0   ID=cds-WP_091897549.1;Parent=gene-BM212_RS06860;Dbxref=GenBank:WP_091897549.1;Name=WP_091897549.1;Note=programmed frameshift;Ontology_term=GO:0004803;exception=ribosomal slippage;gbkey=CDS;go_function=transposase activity|0004803||IEA;inference=COORDINATES: similar to AA sequence:RefSeq:WP_016179228.1;locus_tag=BM212_RS06860;product=IS3 family transposase;protein_id=WP_091897549.1;transl_table=11
NZ_FOSJ01000022.1   Protein Homology    CDS 36288   37183   .   +   2   ID=cds-WP_091897549.1;Parent=gene-BM212_RS06860;Dbxref=GenBank:WP_091897549.1;Name=WP_091897549.1;Note=programmed frameshift;Ontology_term=GO:0004803;exception=ribosomal slippage;gbkey=CDS;go_function=transposase activity|0004803||IEA;inference=COORDINATES: similar to AA sequence:RefSeq:WP_016179228.1;locus_tag=BM212_RS06860;product=IS3 family transposase;protein_id=WP_091897549.1;transl_table=11
$ grep BM212_RS06860 genome.gff | gffread -y- -g genome.fna
>gene-BM212_RS06860
MYSYEERMKAVKLYIKYEHSLASVIHELGYPSSKALYKWFKEYEENGDLSKTSKRKREFTQEEKQFAVNH
YIEYGRNYSRTVRMLGYPDRNILSDWCKELAPEARKIRRSKLNYSQEQKKEAVIHLVSRETSVQKIAERL
QVTRKTLYDWKEELIGEELPPNMTNMPDSPQLEGLKNEVDVLKQEVYRLRMEKDILEKSAELVKKNGGIN
PKHLSNKEKTRVIDALKTFYPLALLLENLDLVKSSYFYHRSQAKLPDKYADLRVMLKDIFIESRRTYGYR
RMHDVLKKQNIVVSEKVVRRLMKEEKLIVHSVTKKKYSSYAGEITPATPNLIQRNFKAEAPDKKWLTDIT
EFRIPAGKVYLSPIIDCFDGALVSWTISTQPNAELVNSMFDQALDTLADNSKPIIHSDRGAHYRWPGWIE
RIKTNGLTQSMSKKGCSPDNSACEGFFGRLKNEFFYGFKWDNISIDAFIELLDDYLHWYNNKRIKLSLGG
LSPIDYRKKLGLIA

But gffread still corrupts the gff in this case as well:

$ grep BM212_RS06860 genome.gff | gffread
NZ_FOSJ01000022.1   RefSeq  gene    35670   37183   .   +   .   ID=gene-BM212_RS06860;
NZ_FOSJ01000022.1   RefSeq  CDS 35670   37183   .   +   .   Parent=gene-BM212_RS06860
NZ_FOSJ01000022.1   RefSeq  CDS 35670   36288   .   +   0   Parent=gene-BM212_RS06860
NZ_FOSJ01000022.1   RefSeq  CDS 36288   37183   .   +   2   Parent=gene-BM212_RS06860
$ grep BM212_RS06860 genome.gff | gffread | gffread -y- -g genome.fna
malloc(): corrupted top size
Aborted (core dumped)
alephreish commented 9 months ago

agat does not crash but the translation is incorrect:

$ agat_sp_extract_sequences.pl -g genome.gff -f genome.fna -t cds -p -o genome.faa
$ seqkit grep -rp BM212_RS11370 genome.faa 
>gene-BM212_RS11370 gene=nbis-gene-2141 seq_id=NZ_FOSJ01000070.1 type=cds
MSTRRPRRTYTEEFKKQIVDLHKAGKSRKEIIEEYDLTGSAFDKWVRQHNQTGSFKEKDN
LTPEQKELKELRKANTQLQMENDILKQAALIFGRKSK*SNVININILYQRCAVPLRSVED
LIIMK**RKKVTRNLNKRLLKSSLKAKTIMARVN*KND*RNVRLLCPVGGLGAL*TSFIW
YPTMINGHTNRNRVASIKRRLKTH*TESSIKSSR*KQSSRT*PMSKLPTSGSMSVLF*TY
LTARLSGIPQGQIRQLTWFCRLSLQLRETYIRSRYSILTGERNSTTTS*MNYWVPSILCA
R*VEKEILTIMP*RNLRINHLSLNLSTVTHSIHSMNCKSN*WTMFIGGILFAHMDPWTTK
RLSIIVKVGNRNSLN