gbrowse 2.54+dfsg-2build1 / usr / bin / bed2gff3

#!/usr/bin/perl 

eval 'exec /usr/bin/perl  -S $0 ${1+"$@"}'
    if 0; # not running under some shell

# convert UCSC gene files into GFF3 data

use strict;
use File::Basename 'basename';
use Getopt::Long;

my $executable = basename($0);

my ($SRC,$ORIGIN);
GetOptions('src:s'    => \$SRC,
	   'origin:i' => \$ORIGIN,
	   ) or die <<USAGE;
Usage: $0 [options] ucsc_file1 ucsc_file2...

Convert UCSC Genome Browser BED-format gene files into GFF3 version files.
Only the gene IDs and their locations come through.  You have to get
the comments and aliases some other way.

Options:

    -src    <string>   Choose a source for the gene, default "UCSC"
    -origin <integer>  Choose a relative position to number from, default is "1"

The resulting file is in GFF3 format and should be loaded into a
Bio::DB::GFF database using the following command:

 bp_bulk_load_gff.pl -c -d db1 --maxfeature 1000000000 --gff3_munge file.gff

USAGE

$SRC    ||= 'UCSC';
$ORIGIN ||= 1;

print "##gff-version 3\n";

# automatically uncompress varous compression formats
foreach (@ARGV) {
  $_ = "gunzip     -c $_ |" if /\.gz$/;
  $_ = "uncompress -c $_ |" if /\.Z$/;
  $_ = "bunzip2    -c $_ |" if /\.bz2$/;
}

while (<>) {
  chomp;
  next if /^\#/;;
  next if /random/;  ## added line

  my ($chrom,$txStart,$txEnd,$id,$score,$strand,$cdsStart,$cdsEnd,
      $itemRGB,$exonCount,$exonSizes,$exonStarts) = split /\t/;
  my ($utr5_start,$utr5_end,$utr3_start,$utr3_end);

  # adjust for Jim's 0-based coordinates
  $txStart++;
  $cdsStart++;

  $txStart  -= $ORIGIN;
  $txEnd    -= $ORIGIN;
  $cdsStart -= $ORIGIN;
  $cdsEnd   -= $ORIGIN;

  # print the transcript
  print join("\t",$chrom,$SRC,'mRNA',$txStart,$txEnd,'.',$strand,'.',"ID=$id;Name=$id"),"\n";

  # now handle the CDS entries -- the tricky part is the need to keep
  # track of phase
  my $phase = 0;
  my @exon_starts = map {$_-$ORIGIN+$txStart} split ',',$exonStarts;
  my @exon_sizes  = map {$_-$ORIGIN} split ',',$exonSizes;
  my @exon_ends   = map {$exon_starts[$_]+$exon_sizes[$_]+1} (0..@exon_sizes);

  if ($strand eq '+') {
    for (my $i=0;$i<@exon_starts;$i++) {  # for each exon start
      my $exon_start = $exon_starts[$i] + 1;
      my $exon_end   = $exon_ends[$i];
      my (@utr_start,@utr_end,$cds_start,$cds_end);

      if ($exon_start < $cdsStart) { # in a 5' UTR
	push (@utr_start, $exon_start);
      } elsif ($exon_start > $cdsEnd) { 
	push (@utr_start, $exon_start);
      } else {
	$cds_start = $exon_start;
      }

      if ($exon_end < $cdsStart) {
	push (@utr_end, $exon_end);
      } elsif ($exon_end > $cdsEnd) {
	push (@utr_end, $exon_end);
      } else {
	$cds_end = $exon_end;
      }

      if ($utr_start[0] && !$utr_end[0]) { # half in half out on 5' end
	$utr_end[0]= $cdsStart - 1;
	$cds_start = $cdsStart;
	$cds_end   = $exon_end;
      }

      if ($utr_end[0] && !$utr_start[0]) { # half in half out on 3' end
	$utr_start[0]= $cdsEnd + 1;
	$cds_end     = $cdsEnd;
	$cds_start   = $exon_start;
      }

      # If the CDS is within the exon
      if (defined $utr_start[0] == defined $utr_end[0] && 
	  $utr_start[0] < $cdsStart && $utr_end[0] > $cdsEnd) {
	$utr_end[0]= $cdsStart - 1;
	$cds_start = $cdsStart;
	$cds_end   = $cdsEnd;
	
	push (@utr_start, $cdsEnd + 1);
	push (@utr_end, $exon_end);
      }


      die "programmer error, not an even number of utr_starts and
utr_ends"
	unless $#utr_start == $#utr_end;
      die "programmer error, cds_start and no cds_end" 
	unless defined $cds_start == defined $cds_end;

      for (my $i=0;$i<@utr_start;$i++) {  # for each utr start
	if (defined $utr_start[$i] && $utr_start[$i] <= $utr_end[$i] &&
$utr_start[$i] < $cdsStart) {
	  print join
("\t",$chrom,$SRC,"five_prime_UTR",$utr_start[$i],$utr_end[$i],'.',$strand,'.',"Parent=$id"),"\n"	
	} # end of if	    
      } # end of foreach

      if (defined $cds_start && $cds_start <= $cds_end) {
	print join
("\t",$chrom,$SRC,'CDS',$cds_start,$cds_end,'.',$strand,$phase,"Parent=$id"),"\n";
	$phase = (($cds_end-$cds_start+1-$phase)) % 3;
      }

      for (my $i=0;$i<@utr_start;$i++) {  # for each utr start
	if (defined $utr_start[$i] && $utr_start[$i] <= $utr_end[$i] &&
$utr_start[$i] > $cdsEnd) {
	  print join ("\t",$chrom,$SRC,"three_prime_UTR",,$utr_start[$i],
$utr_end[$i],'.',$strand,'.',"Parent=$id"),"\n"	
	}
      }
    } # end of for each exon
  } # matches if strand = +


  if ($strand eq '-') {
    my @lines;
    for (my $i=@exon_starts-1; $i>=0; $i--) { # count backwards
      my $exon_start = $exon_starts[$i] + 1;
      my $exon_end   = $exon_ends[$i];
      my (@utr_start,@utr_end,$cds_start,$cds_end);

      if ($exon_end > $cdsEnd) { # in a 5' UTR
	push (@utr_end,  $exon_end);
      } elsif ($exon_end < $cdsStart) {
	push (@utr_end,  $exon_end);
      } else {
	$cds_end = $exon_end;
      }

      if ($exon_start > $cdsEnd) {
	push (@utr_start, $exon_start);
      } elsif ($exon_start < $cdsStart) {
	push (@utr_start, $exon_start);
      } else {
	$cds_start = $exon_start;
      }

      if ($utr_start[0] && !$utr_end[0]) { # half in half out on 3' end
	$utr_end[0]   = $cdsStart - 1;
	$cds_start = $cdsStart;
	$cds_end   = $exon_end;
      }

      if ($utr_end[0] && !$utr_start[0]) { # half in half out on 5' end
	$utr_start[0] = $cdsEnd + 1;
	$cds_end   = $cdsEnd;
	$cds_start = $exon_start;
      }

      # If the CDS is within the exon  
      if (defined $utr_start[0] == defined $utr_end[0] && 
	  $utr_start[0] < $cdsStart && $utr_end[0] > $cdsEnd) {
	$utr_end[0]= $cdsStart - 1;
	$cds_start = $cdsStart;
	$cds_end   = $cdsEnd;
	
	push (@utr_start, $cdsEnd + 1);
	push (@utr_end, $exon_end);
      }

      die "programmer error, not an even number of utr_starts and
utr_ends"
	unless $#utr_start == $#utr_end;

      die "programmer error, cds_start and no cds_end" unless defined
$cds_start == defined $cds_end;

      for (my $i=0;$i<@utr_start;$i++) {  # for each utr start
	if (defined $utr_start[$i] && $utr_start[$i] <= $utr_end[$i] &&
$utr_start[$i] > $cdsEnd) {
	  unshift @lines,join
("\t",$chrom,$SRC,"five_prime_UTR",,$utr_start[$i],$utr_end[$i],'.',$strand,'.',"Parent=$id"),"\n"	
	}
      } # end of for

      if (defined $cds_start && $cds_start <= $cds_end) {
	unshift @lines,join
("\t",$chrom,$SRC,'CDS',$cds_start,$cds_end,'.',$strand,$phase,"Parent=$id"),"\n";
	$phase = (($cds_end-$cds_start+1-$phase)) % 3;
      }

      for (my $i=0;$i<@utr_start;$i++) {  # for each utr start
	if (defined $utr_start[$i] && $utr_start[$i] <= $utr_end[$i] &&
$utr_end[$i] < $cdsStart) {
	  unshift @lines,join
("\t",$chrom,$SRC,"three_prime_UTR",$utr_start[$i],$utr_end[$i],'.',$strand,'.',"Parent=$id"),"\n"	
	}
      } # end for
    }
    print @lines;
  }
} # end while <>

__END__

=head1 NAME

ucsc_genes2gff.pl - Convert UCSC Genome Browser-format gene files into GFF
files suitable for loading into gbrowse

=head1 SYNOPSIS

  % uscsc_genes2gff.pl [options] ucsc_file1 ucsc_file2...

Options:

    -src    <string>   Choose a source for the gene, default "UCSC"
    -origin <integer>  Choose a relative position to number from, default
is "1"

=head1 DESCRIPTION

This script massages the gene files available from the "tables" link
of the UCSC genome browser (genome.ucsc.edu) into a form suitable for
loading of gbrowse.  Warning: it only works with the gene tables.
Other tables, such as EST alignments, contours and repeats, have their
own formats which will require other scripts to parse.

To use this script, get one or more UCSC tables, either from the
"Tables" link on the browser, or from the UCSC Genome Browser FTP
site.  Give the table file as the argument to this script.  You may
want to provide an alternative "source" field.  Otherwise this script
defaults to "UCSC".

  % pucsc_genes2gff.pl -src RefSeq refseq_data.ucsc > refseq.gff

The resulting GFF file can then be loaded into a Bio::DB::GFF database
using the following command:

  % bulk_load_gff.pl -d <databasename> refseq.gff

=head1 SEE ALSO

L<Bio::DB::GFF>, L<bulk_load_gff.pl>, L<load_gff.pl>

=head1 AUTHOR

Lincoln Stein <lstein@cshl.org>.

Copyright (c) 2003 Cold Spring Harbor Laboratory

This library is free software; you can redistribute it and/or modify
it under the same terms as Perl itself.  See DISCLAIMER.txt for
disclaimers of warranty.

=cut