[Bio] / FigKernelScripts / find_gaps.pl Repository:
ViewVC logotype

Diff of /FigKernelScripts/find_gaps.pl

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.2, Thu Jun 8 17:43:12 2006 UTC revision 1.8, Thu Dec 20 07:53:36 2007 UTC
# Line 1  Line 1 
1  # -*- perl -*-  # -*- perl -*-
2  #  ########################################################################
3  # Copyright (c) 2003-2006 University of Chicago and Fellowship  # Copyright (c) 2003-2006 University of Chicago and Fellowship
4  # for Interpretations of Genomes. All Rights Reserved.  # for Interpretations of Genomes. All Rights Reserved.
5  #  #
# Line 14  Line 14 
14  # at info@ci.uchicago.edu or the Fellowship for Interpretation of  # at info@ci.uchicago.edu or the Fellowship for Interpretation of
15  # Genomes at veronika@thefig.info or download a copy from  # Genomes at veronika@thefig.info or download a copy from
16  # http://www.theseed.org/LICENSE.TXT.  # http://www.theseed.org/LICENSE.TXT.
17  #  ########################################################################
   
18    
19  use FIG;  use FIG;
20  $fig = new FIG;  $fig = new FIG;
21    
22  use constant LEFT   =>  0;  use constant LEFT   =>  0;
23  use constant RIGHT  =>  1;  use constant RIGHT  =>  1;
24  use constant FID    =>  2;  use constant LEN    =>  2;
25    use constant FID    =>  3;
26    
 use constant TRAN   =>  0;  
 use constant OTHERS =>  1;  
27    
28  $0 =~ m/([^\/]+)$/;   $this_tool = $1;  $0 =~ m/([^\/]+)$/;   $this_tool = $1;
29  $usage = "$this_tool  org_id  contigs  tbl_1 tbl_2 ... > gap_tbl";  $usage = "$this_tool  [-extend_5p=num] org_id  contigs  tbl_1 tbl_2 ... > gap_tbl";
30    
31  if ((not @ARGV) || ($ARGV[0] =~ m/-h(elp)?/)) {  if (not @ARGV) {
32      die "\n   usage: $usage\n\n";      die "\n\n\tusage: $usage\n\n";
33  }  }
34    
35  $trouble = 0;  $trouble = 0;
36  ($org_id = shift @ARGV) || (($trouble = 1) && (warn "No Org-ID given"));  $extend_5prime_end = 0;
37    
38  for ($i=0; $i < @ARGV; ++$i)  #...Process any flags and switches...
39  {  while ($ARGV[0] =~ m/^-/o) {
40      if (!-e $ARGV[$i])      if ($ARGV[$i] =~ m/-help/) {
41      {          print STDERR "\n\n\tusage: $usage\n\n";
42            exit(0);
43        }
44        elsif ($ARGV[$i] =~ m/-extend_5p=(\d+)/o) {
45            $extend_5prime_end = $1;
46            print STDERR "PEG 5-prime end will be extended $extend_5prime_end bp\n" if $ENV{VERBOSE};
47        }
48        else {
49            $trouble = 1;
50            warn "Unrecognized flag: $ARGV[$i]\n";
51        }
52    
53        shift @ARGV;
54    }
55    
56    
57    #...Process Org-ID argument...
58    if (not @ARGV) {
59        warn "\n\nNo argument list given\n\n";
60        die "\n\n\tusage: $usage\n\n";
61    }
62    else {
63        if ($ARGV[0] !~ m/^\d+\.\d+$/o) {
64            $trouble = 1;
65            warn "\nOrg-ID $ARGV[0] is malformed\n";
66        }
67        else {
68            $org_id = shift @ARGV;
69        }
70    }
71    
72    
73    #...Check existence of file arguments...
74    for ($i=0; $i < @ARGV; ++$i) {
75        if (!-e $ARGV[$i]) {
76          $trouble = 1;          $trouble = 1;
77          print STDERR "ERROR: File $ARGV[$i] does not exist\n";          warn "\nERROR: File $ARGV[$i] does not exist\n";
78      }      }
79  }  }
80  die "\n\nAborting due to invalid args\n\nusage: $usage\n\n" if $trouble;  die "\n\nAborting due to invalid args\n\nusage: $usage\n\n" if $trouble;
# Line 53  Line 85 
85  ((@tbls = @ARGV) > 0) || die "\n\tusage: $usage\n\n";  ((@tbls = @ARGV) > 0) || die "\n\tusage: $usage\n\n";
86    
87  $len_of  = &load_contig_lens($contigs_file);  $len_of  = &load_contig_lens($contigs_file);
88  $regions = &load_regions(@tbls);  $regions = &load_regions($len_of, @tbls);
89  # die Dumper($features);  # die Dumper($features);
90    
91  foreach $contig (keys %$regions)  $gap_num = 0;
92  {  foreach $contig (sort keys %$len_of) {
93        print STDERR "Processing $contig ...\n" if $ENV{VERBOSE};
94    
95      $x = $regions->{$contig};      $x = $regions->{$contig};
96        $x = defined($x) ? $x : [];
97    
98      if ($x->[0]  > 1)                  { unshift @$x, [0, 0]; }  #   if ($x->[0]  > 1)                  { unshift @$x, [0, 0]; }
99      if ($x->[-1] > $len_of->{$contig}) { push    @$x, [1+$len_of->{$contig}, 1+$len_of->{$contig}]; }  #   if ($x->[-1] > $len_of->{$contig}) { push    @$x, [1+$len_of->{$contig}, 1+$len_of->{$contig}]; }
100    
101      for ($i=1; $i < @$x; ++$i)      unshift @$x, [0, 0];
102      {      push    @$x, [1+$len_of->{$contig}, 1+$len_of->{$contig}];
103    
104        for ($i=1; $i < @$x; ++$i) {
105          $gap_beg = $x->[$i-1]->[RIGHT] + 1;          $gap_beg = $x->[$i-1]->[RIGHT] + 1;
106          $gap_end = $x->[$i]->[LEFT] - 1;          $gap_end = $x->[$i]->[LEFT] - 1;
107          $gap_len = $gap_end - $gap_beg + 1;  #       $gap_len = $gap_end - $gap_beg + 1;
108          $gap_loc = "$contig\_$gap_beg\_$gap_end";          $gap_loc = "$contig\_$gap_beg\_$gap_end";
109    
110          print "fig|$org_id.gap.".($i+1) . "\t$gap_loc\n";          if (($gap_beg < 1) || ($gap_end > $len_of->{$contig})) {
111                print STDERR "Skipping out-of-bounds location: $gap_loc\n";
112                next;
113            }
114    
115            print "fig|$org_id.gap.".(++$gap_num) . "\t$gap_loc\n";
116      }      }
117  }  }
118  print STDERR "$this_tool done\n\n";  print STDERR "\n$this_tool done\n\n" if $ENV{VERBOSE};
119  exit(0);  exit(0);
120    
121    
# Line 81  Line 123 
123  {  {
124      my ($contigs_file) = @_;      my ($contigs_file) = @_;
125    
126        my $num_contigs = 0;
127        my $num_chars   = 0;
128    
129      my $len_of = {};      my $len_of = {};
130    
131        print STDERR "Loading $contigs_file ...\n" if $ENV{VERBOSE};
132      open (CONTIGS, "<$contigs_file") or die "could not open $contigs_file to read";      open (CONTIGS, "<$contigs_file") or die "could not open $contigs_file to read";
133      while (($id, $seqP) = &FIG::read_fasta_record(\*CONTIGS))      while (($id, $seqP) = &FIG::read_fasta_record(\*CONTIGS)) {
134      {          ++$num_contigs;
135          $len_of->{$id} = length($$seqP);          $num_chars += $len_of->{$id} = length($$seqP);
136      }      }
137      close(CONTIGS) or die "could not close $contigs_file";      close(CONTIGS) or die "could not close $contigs_file";
138        print STDERR "Loaded $num_contigs contig lengths (total $num_chars chars) from $contigs_file\n\n" if $ENV{VERBOSE};
139    
140      return $len_of;      return $len_of;
141  }  }
142    
143  sub load_regions  sub load_regions
144  {  {
145      my (@tbl_files) = @_;      my ($len_of, @tbl_files) = @_;
146      my ($fid, $org, $locus, $contig, $beg, $end, $len);      my ($fid, $org, $locus, $contig, $beg, $end, $len);
147    
148      my $regions  = {};      my $regions  = {};
149      foreach my $tbl_file (@tbl_files)      foreach my $tbl_file (@tbl_files) {
     {  
150          my $num_fids = 0;          my $num_fids = 0;
151          open(TBL,"<$tbl_file") || die "could not open $tbl_file";          open(TBL,"<$tbl_file") || die "could not open $tbl_file";
152          print STDERR "Loading $tbl_file ...\n" if $ENV{VERBOSE};          print STDERR "Loading $tbl_file ...\n" if $ENV{VERBOSE};
# Line 109  Line 155 
155          {          {
156              ++$num_fids;              ++$num_fids;
157              chomp $entry;              chomp $entry;
158              if ($entry =~ /^(\S+)\s+(\S+)/)              if ($entry =~ m/^(\S+)\s+(\S+)/o) {
             {  
159                  ($fid, $locus) = ($1, $2);                  ($fid, $locus) = ($1, $2);
160                  if ($fid !~ m/^fig\|(\d+\.\d+)\.(rna|peg|orf)\.(\d+)$/)  
161                  {                  $fid_type = "";
162                      die "Invalid FID $fid";                  if ($fid =~ m/^fig\|\d+\.\d+\.(rna|peg|orf)\.\d+$/) {
163                        $fid_type = $1;
164                    }
165                    else {
166                        die "Unrecognized FID type: $fid";
167                  }                  }
168    
169                  ($contig, $beg, $end) = $fig->boundaries_of($locus);                  ($contig, $beg, $end) = $fig->boundaries_of($locus);
170    
171                    if ($extend_5prime_end && ($fid_type eq qq(peg))) {
172                        print STDERR "$fid:\t$locus\tbeg = $beg --> " if $ENV{VERBOSE};
173                        $beg -= ($beg < $end) ? +$extend_5prime_end : -$extend_5prime_end;
174                        $beg = &FIG::max(1, &FIG::min($beg, $len_of->{$contig}));
175                        print STDERR "$beg\n" if $ENV{VERBOSE};
176                    }
177                  $len = 1 + abs($end-$beg);                  $len = 1 + abs($end-$beg);
178    
179                  unless (defined($regions->{$contig})) { $regions->{$contig} = []; }                  unless (defined($regions->{$contig})) { $regions->{$contig} = []; }
180                  push(@ { $regions->{$contig} }, [ &FIG::min($beg, $end), &FIG::max($beg, $end), $fid ] );                  push(@ { $regions->{$contig} }
181                         , [ $left = &FIG::min($beg, $end), $right = &FIG::max($beg, $end), (1 + $right - $left), $fid ]
182                         );
183              }              }
184              else              else
185              {              {
# Line 134  Line 192 
192    
193      foreach my $contig (keys(%$regions))      foreach my $contig (keys(%$regions))
194      {      {
195          my $x = [ sort { $a->[LEFT] <=> $b->[LEFT] }  @ { $regions->{$contig} } ];          my $x = [ sort { ($a->[LEFT] <=> $b->[LEFT]) || ($b->[RIGHT] <=> $a->[RIGHT])
196                             }  @ { $regions->{$contig} }
197                      ];
198    
199          for (my $i=1; $i < @$x; ++$i)          for (my $i=1; $i < @$x; ++$i)
200          {          {
201              while (($i < @$x) && ($x->[$i-1]->[RIGHT] >= $x->[$i]->[LEFT]))              while (($i < @$x) && ($x->[$i-1]->[RIGHT] >= $x->[$i]->[LEFT]))
202              {              {
203                  print STDERR "Merging $i:\t$x->[$i-1]->[FID]\t$x->[$i]->[FID]\n"                  print STDERR "Merging $i:\t$x->[$i-1]->[FID]:$x->[$i-1]->[LEN]\t$x->[$i]->[FID]:$x->[$i]->[LEN]"
204                      if ($ENV{VERBOSE} && ($ENV{VERBOSE} > 1));                      if ($ENV{VERBOSE} && ($ENV{VERBOSE} > 1));
205    
206                  $x->[$i-1]->[RIGHT] = &FIG::max( $x->[$i-1]->[RIGHT], $x->[$i]->[RIGHT] );                  $x->[$i-1]->[RIGHT] = &FIG::max( $x->[$i-1]->[RIGHT], $x->[$i]->[RIGHT] );
207                    my $new_len = 1 + ($x->[$i-1]->[RIGHT] - $x->[$i-1]->[LEFT]);
208    
209                    print STDERR "\t--- new_len = $new_len\n"
210                        if ($ENV{VERBOSE} && ($ENV{VERBOSE} > 1));
211    
212                    $x->[$i-1]->[FID] .= ",$x->[$i]->[FID]";
213                    $x->[$i-1]->[LEN]  =   $new_len;
214    
215                  splice @$x, $i, 1;                  splice @$x, $i, 1;
216              }              }
217          }          }
218    
219          $regions->{$contig} = $x;          $regions->{$contig} = $x;
220          print STDERR "After merger, $contig has ", (scalar @$x), " regions\n\n" if $ENV{VERBOSE};          print STDERR "\nAfter mergers, $contig has ", (scalar @$x), " regions\n\n" if $ENV{VERBOSE};
221      }      }
222    
223      return $regions;      return $regions;

Legend:
Removed from v.1.2  
changed lines
  Added in v.1.8

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3