[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.7, Tue Oct 23 18:13:57 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";          print "fig|$org_id.gap.".(++$gap_num) . "\t$gap_loc\n";
111      }      }
112  }  }
113  print STDERR "$this_tool done\n\n";  print STDERR "\n$this_tool done\n\n" if $ENV{VERBOSE};
114  exit(0);  exit(0);
115    
116    
# Line 81  Line 118 
118  {  {
119      my ($contigs_file) = @_;      my ($contigs_file) = @_;
120    
121        my $num_contigs = 0;
122        my $num_chars   = 0;
123    
124      my $len_of = {};      my $len_of = {};
125    
126        print STDERR "Loading $contigs_file ...\n" if $ENV{VERBOSE};
127      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";
128      while (($id, $seqP) = &FIG::read_fasta_record(\*CONTIGS))      while (($id, $seqP) = &FIG::read_fasta_record(\*CONTIGS)) {
129      {          ++$num_contigs;
130          $len_of->{$id} = length($$seqP);          $num_chars += $len_of->{$id} = length($$seqP);
131      }      }
132      close(CONTIGS) or die "could not close $contigs_file";      close(CONTIGS) or die "could not close $contigs_file";
133        print STDERR "Loaded $num_contigs contig lengths (total $num_chars chars) from $contigs_file\n\n" if $ENV{VERBOSE};
134    
135      return $len_of;      return $len_of;
136  }  }
137    
138  sub load_regions  sub load_regions
139  {  {
140      my (@tbl_files) = @_;      my ($len_of, @tbl_files) = @_;
141      my ($fid, $org, $locus, $contig, $beg, $end, $len);      my ($fid, $org, $locus, $contig, $beg, $end, $len);
142    
143      my $regions  = {};      my $regions  = {};
144      foreach my $tbl_file (@tbl_files)      foreach my $tbl_file (@tbl_files) {
     {  
145          my $num_fids = 0;          my $num_fids = 0;
146          open(TBL,"<$tbl_file") || die "could not open $tbl_file";          open(TBL,"<$tbl_file") || die "could not open $tbl_file";
147          print STDERR "Loading $tbl_file ...\n" if $ENV{VERBOSE};          print STDERR "Loading $tbl_file ...\n" if $ENV{VERBOSE};
# Line 109  Line 150 
150          {          {
151              ++$num_fids;              ++$num_fids;
152              chomp $entry;              chomp $entry;
153              if ($entry =~ /^(\S+)\s+(\S+)/)              if ($entry =~ m/^(\S+)\s+(\S+)/o) {
             {  
154                  ($fid, $locus) = ($1, $2);                  ($fid, $locus) = ($1, $2);
155                  if ($fid !~ m/^fig\|(\d+\.\d+)\.(rna|peg|orf)\.(\d+)$/)  
156                  {                  $fid_type = "";
157                      die "Invalid FID $fid";                  if ($fid =~ m/^fig\|\d+\.\d+\.(rna|peg|orf)\.\d+$/) {
158                        $fid_type = $1;
159                    }
160                    else {
161                        die "Unrecognized FID type: $fid";
162                  }                  }
163    
164                  ($contig, $beg, $end) = $fig->boundaries_of($locus);                  ($contig, $beg, $end) = $fig->boundaries_of($locus);
165    
166                    if ($extend_5prime_end && ($fid_type eq qq(peg))) {
167                        print STDERR "$fid:\t$locus\tbeg = $beg --> " if $ENV{VERBOSE};
168                        $beg -= ($beg < $end) ? +$extend_5prime_end : -$extend_5prime_end;
169                        $beg = &FIG::max(0, &FIG::min($beg, $len_of->{$contig}));
170                        print STDERR "$beg\n" if $ENV{VERBOSE};
171                    }
172                  $len = 1 + abs($end-$beg);                  $len = 1 + abs($end-$beg);
173    
174                  unless (defined($regions->{$contig})) { $regions->{$contig} = []; }                  unless (defined($regions->{$contig})) { $regions->{$contig} = []; }
175                  push(@ { $regions->{$contig} }, [ &FIG::min($beg, $end), &FIG::max($beg, $end), $fid ] );                  push(@ { $regions->{$contig} }
176                         , [ $left = &FIG::min($beg, $end), $right = &FIG::max($beg, $end), (1 + $right - $left), $fid ]
177                         );
178              }              }
179              else              else
180              {              {
# Line 134  Line 187 
187    
188      foreach my $contig (keys(%$regions))      foreach my $contig (keys(%$regions))
189      {      {
190          my $x = [ sort { $a->[LEFT] <=> $b->[LEFT] }  @ { $regions->{$contig} } ];          my $x = [ sort { ($a->[LEFT] <=> $b->[LEFT]) || ($b->[RIGHT] <=> $a->[RIGHT])
191                             }  @ { $regions->{$contig} }
192                      ];
193    
194          for (my $i=1; $i < @$x; ++$i)          for (my $i=1; $i < @$x; ++$i)
195          {          {
196              while (($i < @$x) && ($x->[$i-1]->[RIGHT] >= $x->[$i]->[LEFT]))              while (($i < @$x) && ($x->[$i-1]->[RIGHT] >= $x->[$i]->[LEFT]))
197              {              {
198                  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]"
199                      if ($ENV{VERBOSE} && ($ENV{VERBOSE} > 1));                      if ($ENV{VERBOSE} && ($ENV{VERBOSE} > 1));
200    
201                  $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] );
202                    my $new_len = 1 + ($x->[$i-1]->[RIGHT] - $x->[$i-1]->[LEFT]);
203    
204                    print STDERR "\t--- new_len = $new_len\n"
205                        if ($ENV{VERBOSE} && ($ENV{VERBOSE} > 1));
206    
207                    $x->[$i-1]->[FID] .= ",$x->[$i]->[FID]";
208                    $x->[$i-1]->[LEN]  =   $new_len;
209    
210                  splice @$x, $i, 1;                  splice @$x, $i, 1;
211              }              }
212          }          }
213    
214          $regions->{$contig} = $x;          $regions->{$contig} = $x;
215          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};
216      }      }
217    
218      return $regions;      return $regions;

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3