1
$\begingroup$

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

$\endgroup$
2
  • $\begingroup$ Are you going for the short format or long format? Can you provide a spec for what you'd like the output to look like? ZFF is not a common file format. $\endgroup$
    – conchoecia
    Nov 8 '21 at 6:35
  • $\begingroup$ I would go the long format. The specs are described here. $\endgroup$
    – user977828
    Nov 8 '21 at 7:11
0
$\begingroup$

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...

$\endgroup$
3
  • $\begingroup$ Originally, I had a GTF but I failed with gffread --gene2exon -E traingenes.gtf -o traingenes.gff3. What did I miss? $\endgroup$
    – user977828
    Nov 8 '21 at 23:36
  • $\begingroup$ What do you mean 'failed'? $\endgroup$
    – Dan Bolser
    Nov 9 '21 at 13:30
  • $\begingroup$ Sorry, I meant it still kept CDS rather than changing them to Exon. What did I miss? $\endgroup$
    – user977828
    Nov 10 '21 at 0:34
0
$\begingroup$

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.

$\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.