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

Annotation of /FigKernelScripts/find_gaps.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.9 - (view) (download) (as text)

1 : overbeek 1.1 # -*- perl -*-
2 : gdpusch 1.5 ########################################################################
3 : overbeek 1.1 # Copyright (c) 2003-2006 University of Chicago and Fellowship
4 :     # for Interpretations of Genomes. All Rights Reserved.
5 :     #
6 :     # This file is part of the SEED Toolkit.
7 :     #
8 :     # The SEED Toolkit is free software. You can redistribute
9 :     # it and/or modify it under the terms of the SEED Toolkit
10 :     # Public License.
11 :     #
12 :     # You should have received a copy of the SEED Toolkit Public License
13 :     # along with this program; if not write to the University of Chicago
14 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
15 :     # Genomes at veronika@thefig.info or download a copy from
16 :     # http://www.theseed.org/LICENSE.TXT.
17 : gdpusch 1.5 ########################################################################
18 : overbeek 1.1
19 :     use FIG;
20 :     $fig = new FIG;
21 :    
22 :     use constant LEFT => 0;
23 :     use constant RIGHT => 1;
24 : gdpusch 1.7 use constant LEN => 2;
25 :     use constant FID => 3;
26 : overbeek 1.1
27 :    
28 :     $0 =~ m/([^\/]+)$/; $this_tool = $1;
29 : gdpusch 1.6 $usage = "$this_tool [-extend_5p=num] org_id contigs tbl_1 tbl_2 ... > gap_tbl";
30 : overbeek 1.2
31 : gdpusch 1.6 if (not @ARGV) {
32 :     die "\n\n\tusage: $usage\n\n";
33 : overbeek 1.2 }
34 : overbeek 1.1
35 :     $trouble = 0;
36 : gdpusch 1.6 $extend_5prime_end = 0;
37 : overbeek 1.1
38 : gdpusch 1.6 #...Process any flags and switches...
39 :     while ($ARGV[0] =~ m/^-/o) {
40 :     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 : overbeek 1.1 $trouble = 1;
77 : gdpusch 1.6 warn "\nERROR: File $ARGV[$i] does not exist\n";
78 : overbeek 1.1 }
79 :     }
80 : overbeek 1.4 die "\n\nAborting due to invalid args\n\n usage: $usage\n\n" if $trouble;
81 : overbeek 1.1
82 :     (($contigs_file = shift @ARGV) && (-s $contigs_file))
83 :     || die "Contigs file $contigs_file has zero size\n\n\tusage: $usage\n\n";
84 :    
85 :     ((@tbls = @ARGV) > 0) || die "\n\tusage: $usage\n\n";
86 :    
87 :     $len_of = &load_contig_lens($contigs_file);
88 : gdpusch 1.6 $regions = &load_regions($len_of, @tbls);
89 : overbeek 1.1 # die Dumper($features);
90 :    
91 : overbeek 1.3 $gap_num = 0;
92 : gdpusch 1.6 foreach $contig (sort keys %$len_of) {
93 : gdpusch 1.7 print STDERR "Processing $contig ...\n" if $ENV{VERBOSE};
94 : overbeek 1.4
95 : overbeek 1.1 $x = $regions->{$contig};
96 : overbeek 1.4 $x = defined($x) ? $x : [];
97 :    
98 :     # 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}]; }
100 : overbeek 1.1
101 : overbeek 1.4 unshift @$x, [0, 0];
102 :     push @$x, [1+$len_of->{$contig}, 1+$len_of->{$contig}];
103 : overbeek 1.1
104 : gdpusch 1.6 for ($i=1; $i < @$x; ++$i) {
105 : overbeek 1.1 $gap_beg = $x->[$i-1]->[RIGHT] + 1;
106 :     $gap_end = $x->[$i]->[LEFT] - 1;
107 : gdpusch 1.6 # $gap_len = $gap_end - $gap_beg + 1;
108 : overbeek 1.1 $gap_loc = "$contig\_$gap_beg\_$gap_end";
109 :    
110 : gdpusch 1.9 if (($gap_beg < 1) || ($gap_end > $len_of->{$contig}) || ($gap_beg > $gap_end)) {
111 :     print STDERR "Skipping out-of-bounds location: $gap_loc\n" if $ENV{VERBOSE};
112 : gdpusch 1.8 next;
113 :     }
114 :    
115 : gdpusch 1.9 print STDOUT (qq(fig|$org_id.gap.), (++$gap_num), qq(\t$gap_loc\n));
116 : overbeek 1.1 }
117 :     }
118 : overbeek 1.4 print STDERR "\n$this_tool done\n\n" if $ENV{VERBOSE};
119 : overbeek 1.1 exit(0);
120 :    
121 :    
122 :     sub load_contig_lens
123 :     {
124 :     my ($contigs_file) = @_;
125 :    
126 : overbeek 1.4 my $num_contigs = 0;
127 :     my $num_chars = 0;
128 :    
129 : overbeek 1.1 my $len_of = {};
130 :    
131 : overbeek 1.4 print STDERR "Loading $contigs_file ...\n" if $ENV{VERBOSE};
132 : overbeek 1.1 open (CONTIGS, "<$contigs_file") or die "could not open $contigs_file to read";
133 : gdpusch 1.6 while (($id, $seqP) = &FIG::read_fasta_record(\*CONTIGS)) {
134 : overbeek 1.4 ++$num_contigs;
135 :     $num_chars += $len_of->{$id} = length($$seqP);
136 : overbeek 1.1 }
137 :     close(CONTIGS) or die "could not close $contigs_file";
138 : gdpusch 1.5 print STDERR "Loaded $num_contigs contig lengths (total $num_chars chars) from $contigs_file\n\n" if $ENV{VERBOSE};
139 : overbeek 1.4
140 : overbeek 1.1 return $len_of;
141 :     }
142 :    
143 :     sub load_regions
144 :     {
145 : gdpusch 1.6 my ($len_of, @tbl_files) = @_;
146 : overbeek 1.1 my ($fid, $org, $locus, $contig, $beg, $end, $len);
147 :    
148 :     my $regions = {};
149 : gdpusch 1.6 foreach my $tbl_file (@tbl_files) {
150 : overbeek 1.1 my $num_fids = 0;
151 :     open(TBL,"<$tbl_file") || die "could not open $tbl_file";
152 :     print STDERR "Loading $tbl_file ...\n" if $ENV{VERBOSE};
153 :    
154 :     while (defined($entry = <TBL>))
155 :     {
156 :     ++$num_fids;
157 :     chomp $entry;
158 : gdpusch 1.6 if ($entry =~ m/^(\S+)\s+(\S+)/o) {
159 : overbeek 1.1 ($fid, $locus) = ($1, $2);
160 : gdpusch 1.6
161 :     $fid_type = "";
162 :     if ($fid =~ m/^fig\|\d+\.\d+\.(rna|peg|orf)\.\d+$/) {
163 :     $fid_type = $1;
164 :     }
165 :     else {
166 :     die "Unrecognized FID type: $fid";
167 : overbeek 1.1 }
168 :    
169 :     ($contig, $beg, $end) = $fig->boundaries_of($locus);
170 : gdpusch 1.6
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 : gdpusch 1.8 $beg = &FIG::max(1, &FIG::min($beg, $len_of->{$contig}));
175 : gdpusch 1.6 print STDERR "$beg\n" if $ENV{VERBOSE};
176 :     }
177 : overbeek 1.1 $len = 1 + abs($end-$beg);
178 :    
179 :     unless (defined($regions->{$contig})) { $regions->{$contig} = []; }
180 : gdpusch 1.7 push(@ { $regions->{$contig} }
181 :     , [ $left = &FIG::min($beg, $end), $right = &FIG::max($beg, $end), (1 + $right - $left), $fid ]
182 :     );
183 : overbeek 1.1 }
184 :     else
185 :     {
186 :     print STDERR "Skipping invalid entry $entry\n";
187 :     }
188 :     }
189 :     print STDERR "Loaded $num_fids features from $tbl_file\n\n" if $ENV{VERBOSE};
190 :     close(TBL);
191 :     }
192 :    
193 :     foreach my $contig (keys(%$regions))
194 :     {
195 : gdpusch 1.7 my $x = [ sort { ($a->[LEFT] <=> $b->[LEFT]) || ($b->[RIGHT] <=> $a->[RIGHT])
196 :     } @ { $regions->{$contig} }
197 :     ];
198 :    
199 : overbeek 1.1 for (my $i=1; $i < @$x; ++$i)
200 :     {
201 :     while (($i < @$x) && ($x->[$i-1]->[RIGHT] >= $x->[$i]->[LEFT]))
202 :     {
203 : gdpusch 1.7 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));
205 :    
206 :     $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 : overbeek 1.1 if ($ENV{VERBOSE} && ($ENV{VERBOSE} > 1));
211 : gdpusch 1.7
212 :     $x->[$i-1]->[FID] .= ",$x->[$i]->[FID]";
213 :     $x->[$i-1]->[LEN] = $new_len;
214 :    
215 : overbeek 1.1 splice @$x, $i, 1;
216 :     }
217 :     }
218 :    
219 :     $regions->{$contig} = $x;
220 : gdpusch 1.7 print STDERR "\nAfter mergers, $contig has ", (scalar @$x), " regions\n\n" if $ENV{VERBOSE};
221 : overbeek 1.1 }
222 :    
223 :     return $regions;
224 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3