Perl convert GFF3 file

Question

I have the following GFF3 file which I would like to convert to ZFF file with the below script:

$ head -n 20 traingenes.gff3 
##gff-version 3
# gffread v0.12.7
# gffread -E traingenes.gtf -o traingenes.gff3
NbLab330C00 GeneMark.hmm3   mRNA    74501   76501   .   +   .   ID=1_t;geneID=1_g
NbLab330C00 GeneMark.hmm3   CDS 74501   74512   .   +   0   Parent=1_t
NbLab330C00 GeneMark.hmm3   CDS 75568   75999   .   +   0   Parent=1_t
NbLab330C00 GeneMark.hmm3   CDS 76151   76501   .   +   0   Parent=1_t
NbLab330C00 GeneMark.hmm3   mRNA    76535   77079   .   +   .   ID=2_t;geneID=2_g
NbLab330C00 GeneMark.hmm3   CDS 76535   76639   .   +   0   Parent=2_t
NbLab330C00 GeneMark.hmm3   CDS 76702   77079   .   +   0   Parent=2_t
NbLab330C00 GeneMark.hmm3   mRNA    93763   100703  .   -   .   ID=3_t;geneID=3_g
NbLab330C00 GeneMark.hmm3   CDS 93763   93837   .   -   0   Parent=3_t
NbLab330C00 GeneMark.hmm3   CDS 93915   94031   .   -   0   Parent=3_t
NbLab330C00 GeneMark.hmm3   CDS 94351   94430   .   -   2   Parent=3_t
NbLab330C00 GeneMark.hmm3   CDS 95483   95589   .   -   1   Parent=3_t
NbLab330C00 GeneMark.hmm3   CDS 95697   95746   .   -   0   Parent=3_t
NbLab330C00 GeneMark.hmm3   CDS 100464  100703  .   -   0   Parent=3_t

#!/usr/bin/env perl
use strict;
# 
# source: https://biowize.wordpress.com/2012/06/01/training-the-snap-ab-initio-gene-predictor/ 
#
my @exons;
my $gene_count = 0;
my $current_seq = "";
while (my $line = <STDIN>) {
  if ($line =~ m/^###/) {
    flush(\@exons);
    next;
  }
  my @fields = split(/\t/, $line);
  if ($fields[2] eq "mRNA") {
    flush(\@exons);
  } elsif ($fields[2] eq "exon") {
    if ($fields[0] ne $current_seq) {
      $current_seq = $fields[0];
      printf(">%s\n", $current_seq);
    }
    push(@exons, \@fields);
  }
}
flush();
 
sub flush
{
  my $num_exons = scalar(@exons);
  return if($num_exons == 0);
   
  my $group = sprintf("%s.%d", $exons[0]->[0], $gene_count);
  $gene_count++;
   
  if ($num_exons == 1) {
my($start, $end) = ($exons[0]->[3], $exons[0]->[4]);
if ($exons[0]->[6] eq "-") {
  ($start, $end) = ($exons[0]->[4], $exons[0]->[3]);
}
printf("Esngl\t%lu\t%lu\t%s\n", $start, $end, $group);
  } else {
@exons = reverse(@exons) if($exons[0]->[6] eq "-");
for (my $i = 0; $i < $num_exons; $i++) {
  my $exon_type = "Exon";
  if ($i == 0) {
    $exon_type = "Einit";
  } elsif ($i == $num_exons - 1) {
    $exon_type = "Eterm";
      }
       
      my($start, $end) = ($exons[$i]->[3], $exons[$i]->[4]);
      if ($exons[0]->[6] eq "-") {
    ($start, $end) = ($exons[$i]->[4], $exons[$i]->[3]);
      }
       
      printf("%s\t%lu\t%lu\t%s\n", $exon_type, $start, $end, $group);
    }
  }
  @exons = ();
}

However, when I run it like perl gff3_to_zff.pl < traingenes.gff3 >> traingenes.ann an empty file get generated.

What did I miss?

Thank you in advance,

Michal

Answer 1

Your file looks like it is CDS feature based, not exon feature based, which seems to be what the script expects.

There is a tool to convert CDS based GFF to exon based that I've used successfully in the past...

• Could be GFF utilities?
• Actually, I think it's GenomeTools

Answer 2

The easiest is to use agat_convert_sp_gff2zff.pl from AGAT
(Edit Actually you will have to first standardize your file because you do not have gene features, neither gene_id or Parent atttributes in your mRNA features. AGAT will be confused. So, in your case, prior the gff2zff conversion you will have to run agat_convert_sp_gxf2gxf.pl -c geneID --gff file.gff -o file_clean.gff)

Otherwise you can use your script but you need to add the exon as suggested by @Dan.
To do so you can just replace the CDS by exon using a sed command, it should work.
You can also use agat_convert_sp_gxf2gxf.pl (with parameter -c geneID) it will add the gene and exon features.