#!/usr/bin/env perl

use strict;
use warnings;
use Getopt::Long;
use List::Util qw(min max);
use feature "say";

my $usage = <<EOS;
  Synopsis: cat GFF_FILE.gff3 | marker_gff_to_bed.pl -fai GENOME.fna.fai [options]
  
  Read marker coordinates in GFF format and a genome fai index file, and return
  a four-column BED file with genomic coordinates flanking the marker coordinates.
  Adjust from GFF 1-based, closed [start, end] to BED 0-based, half-open [start-1, end)
  Unless the "-merge" flag is set, two sequences are returned for each marker: 
  one before (>ID.UP), one after (>ID.DN). 
  The fai file (or a file with seqid in column 1 and molecule length in column 2) 
  is used to identify markers near the ends of molecules.
  The optional -pad flag indicates the distance up- and down-stream from the marker
  to be extracted (default 1000). If the marker position is less than (pad) units from 
  the start of a molecule, then the padded start will be set to 0; and if the marker
  position is less than (pad) units from the end of the molecule, then the padded end
  will be set to the length of that molecule.
  
  Required:
    stream of GFF data on STDIN
    -fai   Index file generated by faidx
       
  Options:
    -pad   The flanking distance up- and down-stream from the marker. Default 1000.
    -merge (boolean; default false). If set, then one sequence will be returned for each marker;
             otherwise, return one beUP (>ID.UP), one DNer (>ID.DN).
    -out   File to write to; otherwise, to stdout.
    -help  (boolean) This message.
EOS

my ($fai_file, $out_file, $merge, $help);
my $pad = 1000;

GetOptions (
  "fai=s" =>  \$fai_file,   
  "pad:i" =>  \$pad,  
  "merge" =>  \$merge,
  "out:s" =>  \$out_file,   
  "help" =>   \$help,
);

die "$usage" unless $fai_file;
die "$usage" if $help;

open my $FAI_FH, "<", $fai_file or die "Can't open in $fai_file: $!\n";

my $OUT_FH;
if ( $out_file ){
  open ( $OUT_FH, ">", $out_file ) or die "Can't open out $out_file: $!\n";
}
else {
  open ( $OUT_FH, ">&", \*STDOUT) or die;
}

# Read in fai file and get molecule lengths
my %seqid_len;
while (<$FAI_FH>) {
  chomp;
  my @fields = split(/\t/, $_);
  my ($seqid, $len) = ($fields[0], $fields[1]);
  $seqid_len{$seqid} = $len;
  #say "$seqid\t$len";
}

# Read in the GFF;
while (<>) {
  s/\r?\n\z//; # CRLF to LF
  chomp;
  my $line = $_;
  if ( $line =~ /^#/ ) { next } # Skip comment and spacer lines
  else { # body of the GFF
    my @parts = split(/\t/, $line);
    if (scalar(@parts)<9){ next }
    my $col9 = $parts[8];
    my @col9_attrs = split(/;/, $col9);
    my ($seqid, $mrk_start, $mrk_end) = ($parts[0], $parts[3], $parts[4]);
    my $name;
    for my $attr (@col9_attrs){
      if ($attr =~ /Name=(.+)/){
        $name = $1;
        $name =~ s/([^;]+);.+/$1/;
      }
    }
    my $pad_start = max(0, $mrk_start-$pad-1);
    my $pad_end = min($mrk_end+1+$pad, $seqid_len{$seqid});
    if ($merge){ # Return one sequence per marker, including flanking upstream and down
      say $OUT_FH join("\t", $seqid, $pad_start, $pad_end, $name);
    }
    else { # Return two sequences per marker: one upstream, one down
      say $OUT_FH join("\t", $seqid, $pad_start, $mrk_end-1, "$name.UP");
      say $OUT_FH join("\t", $seqid, $mrk_end+1, $pad_end, "$name.DN");
    }
  }
}

__END__

Steven Cannon
2023-03-09 Start 

