sanitize-gtf-for-featureCounts

#! /usr/bin/env perl

# ------------------------------------------------------------------------
# Script to remove and flag GTF entries that will cause featureCounts
# to fail.
# ------------------------------------------------------------------------

use strict;
use warnings FATAL => 'all';
#use Carp::Always;

# use FindBin;
# use lib "$FindBin::Bin";
# use Xyzzy;

use constant { TRUE => 1, FALSE => 0 };

# ------------------------------------------------------------------------
# Process the command line
# ------------------------------------------------------------------------

use File::Basename;
use Getopt::Std;

our $opt_a;
our $opt_f;
our $opt_h;

sub usage {
  my $progname = basename($0);
  print STDERR "Usage: $progname [options] < annotations.gtf\n";
  print STDERR "-a ATTRIBUTE[,ATTRIBUTE,...] - attributes to retain\n";
  print STDERR "-f FEATURE[,FEATURE,...] - features to retain\n";
  print STDERR "-h - print help\n";
  exit(@_);
}

my $stat = getopts('a:f:h');
if (!$stat) {
  usage(1);
}
if ($opt_h) {
  usage();
}
if (scalar(@ARGV) != 0) {
  usage(1);
}

# ------------------------------------------------------------------------

my %target_attributes;
if ( $opt_a ) {
  foreach my $attribute ( split(/,/,$opt_a) ) {
    $target_attributes{$attribute} = TRUE;
  }
}

my %target_features;
if ( $opt_f ) {
  foreach my $feature ( split(/,/,$opt_f) ) {
    $target_features{$feature} = TRUE;
  }
}

# ------------------------------------------------------------------------

use constant { NO_VALUE => ";no-value;" };

sub parse_gtf_attributes {
  my ($raw_attributes) = @_;
  my $attributes = {};
  foreach my $key_val (split(/; */,$raw_attributes)) {
    my ($key,$val);
    if ( $key_val =~ /^([^ ]+) (.*)/ ) {
      ($key,$val) = ($1,$2);
      $val =~ s/^"//;
      $val =~ s/"$//;
      $val =~ s/" +"/,/g;
    } else {
      ($key,$val) = ($key_val, NO_VALUE);
    }
    $attributes->{$key} = $val;
  }
  return $attributes;
}

sub unparse_gtf_attributes {
  my ($attributes) = @_;
  my @l;
  foreach my $key (sort (keys %$attributes) ) {
    my $val = $attributes->{$key};
    if ( $val eq NO_VALUE ) {
      push @l, $key;
    } else {
      if ( $val =~ /[ "'\\]/ ) {
	$val = '"'.$val.'"';
      }
      push @l, $key." ".$val;
    }
  }
  return join(";",@l);
}

# ------------------------------------------------------------------------

while (<STDIN>) {
  chomp;
  if ( /^#/ ) {
    next;
  }
  my ($accession,$source,$feature,$start,$end,
      $score,$strand,$frame,$raw_attributes) = split(/\t/,$_);
  if ($opt_f and !$target_features{$feature}) {
    next;
  }
  my $in_attributes = parse_gtf_attributes($raw_attributes);
  if ($start >= $end) {
    print STDERR "Skipping non-positive length feature: ", $in_attributes->{gene_id},"\n";
    next;
  }
  my $out_attributes;
  if (!$opt_a) {
    $out_attributes = $in_attributes;
  } else {
    foreach my $k ( keys($in_attributes->%*) ) {
      if ($target_attributes{$k}) {
	$out_attributes->{$k} = $in_attributes->{$k};
      }
    }
  }
  print join("\t",$accession,$source,$feature,$start,$end,
	     $score,$strand,$frame,unparse_gtf_attributes($out_attributes)),"\n";
}