2
$\begingroup$

I am predicting genes using GeneMark-ES, I managed to generate a gtf file, but when I tried to extract the nt/protein sequence using get_sequence_from_GTF.pl, I am running into an error in the format of the gtf file, which I find puzzing as it was the gtf file was generated using gmes_petap.pl (an script in the GeneMark-ES package).

Here is the command I used and the error message which also prints the first line of the genemark.gtf file:

$ get_sequence_from_GTF.pl genemark.gtf DT_CAJQTU010000037_1.fasta

error in file format GTF: CAJQTU010000037.1 Desoria tigrina genome assembly, contig: utg000040l, whole genome shotgun sequence     GeneMark.hmm3   gene    2099    4039    .       +       .       gene_id "1_g"; transcript_id "1_t";

I think I have found the line in get_sequence_from_GTF.pl that is causing the error, but not being a perl monk, I am having troubles to read figure out what is wrong with my gtf file and more importantly figure out how to fix it.

...
if ( /^(\S+)\t\S+\t(\S+)\t(\d+)\t(\d+)\t.\t([\+-])\t([012.])\tgene_id \"(\S+)\"/ )
    {
    ... (some stuff)
    }
else
{
  if ( /^\s+$/ ) { next };
  if ( /^#/ ) { next };
  die( "error in file format GTF: $_\n" );
}
... 
$\endgroup$
0
1
$\begingroup$

I'm afraid this is a clear bug in the tool you are using. The regular expression requires that all fields in the GTF file contain no whitespace. The \S+ means "one or more non-whitespace characters" but your "source" field, the second field in the GTF file, contains spaces: Desoria tigrina genome assembly, contig: utg000040l, whole genome shotgun sequence.

The thing is, that is a perfectly valid field for a GTF file! GTF fields can contain spaces, they just can't contain tabs. I tried to find a repository for the software to suggest a fix, but found only an empty one. In any case, try changing the regex you show to this:

#if ( /^(\S+)\t\S+\t(\S+)\t(\d+)\t(\d+)\t.\t([\+-])\t([012.])\tgene_id \"(\S+)\"/ )
if ( /^([^\t]+)\t[^\t]+\t([^\t]+)\t(\d+)\t(\d+)\t.\t([\+-])\t([012.])\tgene_id \"(\S+)\"/ )
    {
    ... (some stuff)
    }

I can't guarantee that this will work as expected, please test! The details will depend on what they do with the captured groups (the some stuff), but at least that regex will match and I think I have kept the numbering of the capture groups.

$\endgroup$
0
$\begingroup$

You might use agat_sp_extract_sequences.pl from AGAT instead.

To extract and translate the coding regions:
agat_sp_extract_sequences.pl -g infile.gxf -f infile.fasta -p

see here for all possibilities with that script: https://agat.readthedocs.io/en/latest/tools/agat_sp_extract_sequences.html

Otherwise with the modification suggested by @terdon, your script should work fine.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.