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" );
}
...
Improve this question
$\endgroup$
0

2 Answers 2

Reset to default
2
$\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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ I agree with this 100%. I read the regex and didn't spot the reply. There are quite a few problems, the use of a bare "." just means 'any character' rather '\.' which in this example would fit (albeit it might be intentional). The final section '; transcript_id "1_t";' isn't captured. However, the key problem is \S which will trip over the space characters in "CAJQTU010000037.1 Desoria tigrina genome assembly, contig: utg000040l, whole genome shotgun sequence". The use of '[^\t]', i.e. not a tab character will work. Personally I would split on the tab character. $\endgroup$
    – M__ ♦
    Mar 29, 2022 at 3:32
  • $\begingroup$ Just to clarify this "\t.\t" means any single character between two tabs. However, this "[012.]" might mean a \. i.e. just a full-stop, I can't remember. If it doesn't [.] is exactly the same as [012.] or even [a-z|0-9] and beyond. Anyway this is the accepted answer. $\endgroup$
    – M__ ♦
    Mar 29, 2022 at 3:39
1
$\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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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