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

Annotation of /FigKernelScripts/find_gaps.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 # -*- perl -*-
2 :     #
3 :     # 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 :     #
18 :    
19 :    
20 :     use FIG;
21 :     $fig = new FIG;
22 :    
23 :     use constant LEFT => 0;
24 :     use constant RIGHT => 1;
25 :     use constant FID => 2;
26 :    
27 :     use constant TRAN => 0;
28 :     use constant OTHERS => 1;
29 :    
30 :     $0 =~ m/([^\/]+)$/; $this_tool = $1;
31 : overbeek 1.2 $usage = "$this_tool org_id contigs tbl_1 tbl_2 ... > gap_tbl";
32 :    
33 :     if ((not @ARGV) || ($ARGV[0] =~ m/-h(elp)?/)) {
34 :     die "\n usage: $usage\n\n";
35 :     }
36 : overbeek 1.1
37 :     $trouble = 0;
38 : overbeek 1.4 ($org_id = shift @ARGV) || (($trouble = 1) && (print STDERR "\nNo Org-ID given"));
39 :     ($org_id =~ m/^\d+\.\d+$/o) || (($trouble = 1) && (print STDERR "\nOrg-ID $org_id is malformed"));
40 : overbeek 1.1
41 :     for ($i=0; $i < @ARGV; ++$i)
42 :     {
43 :     if (!-e $ARGV[$i])
44 :     {
45 :     $trouble = 1;
46 : overbeek 1.4 print STDERR "\nERROR: File $ARGV[$i] does not exist";
47 : overbeek 1.1 }
48 :     }
49 : overbeek 1.4 die "\n\nAborting due to invalid args\n\n usage: $usage\n\n" if $trouble;
50 : overbeek 1.1
51 :     (($contigs_file = shift @ARGV) && (-s $contigs_file))
52 :     || die "Contigs file $contigs_file has zero size\n\n\tusage: $usage\n\n";
53 :    
54 :     ((@tbls = @ARGV) > 0) || die "\n\tusage: $usage\n\n";
55 :    
56 :     $len_of = &load_contig_lens($contigs_file);
57 :     $regions = &load_regions(@tbls);
58 :     # die Dumper($features);
59 :    
60 : overbeek 1.3 $gap_num = 0;
61 : overbeek 1.4 foreach $contig (sort keys %$len_of)
62 : overbeek 1.1 {
63 : overbeek 1.4 print STDERR "Scanning $contig ...\n" if $ENV{VERBOSE};
64 :    
65 : overbeek 1.1 $x = $regions->{$contig};
66 : overbeek 1.4 $x = defined($x) ? $x : [];
67 :    
68 :     # if ($x->[0] > 1) { unshift @$x, [0, 0]; }
69 :     # if ($x->[-1] > $len_of->{$contig}) { push @$x, [1+$len_of->{$contig}, 1+$len_of->{$contig}]; }
70 : overbeek 1.1
71 : overbeek 1.4 unshift @$x, [0, 0];
72 :     push @$x, [1+$len_of->{$contig}, 1+$len_of->{$contig}];
73 : overbeek 1.1
74 :     for ($i=1; $i < @$x; ++$i)
75 :     {
76 :     $gap_beg = $x->[$i-1]->[RIGHT] + 1;
77 :     $gap_end = $x->[$i]->[LEFT] - 1;
78 :     $gap_len = $gap_end - $gap_beg + 1;
79 :     $gap_loc = "$contig\_$gap_beg\_$gap_end";
80 :    
81 : overbeek 1.3 print "fig|$org_id.gap.".(++$gap_num) . "\t$gap_loc\n";
82 : overbeek 1.1 }
83 :     }
84 : overbeek 1.4 print STDERR "\n$this_tool done\n\n" if $ENV{VERBOSE};
85 : overbeek 1.1 exit(0);
86 :    
87 :    
88 :     sub load_contig_lens
89 :     {
90 :     my ($contigs_file) = @_;
91 :    
92 : overbeek 1.4 my $num_contigs = 0;
93 :     my $num_chars = 0;
94 :    
95 : overbeek 1.1 my $len_of = {};
96 :    
97 : overbeek 1.4 print STDERR "Loading $contigs_file ...\n" if $ENV{VERBOSE};
98 : overbeek 1.1 open (CONTIGS, "<$contigs_file") or die "could not open $contigs_file to read";
99 :     while (($id, $seqP) = &FIG::read_fasta_record(\*CONTIGS))
100 :     {
101 : overbeek 1.4 ++$num_contigs;
102 :     $num_chars += $len_of->{$id} = length($$seqP);
103 : overbeek 1.1 }
104 :     close(CONTIGS) or die "could not close $contigs_file";
105 : overbeek 1.4 print "Loaded $num_contigs contig lengths (total $num_chars chars) from $contigs_file\n\n" if $ENV{VERBOSE};
106 :    
107 : overbeek 1.1 return $len_of;
108 :     }
109 :    
110 :     sub load_regions
111 :     {
112 :     my (@tbl_files) = @_;
113 :     my ($fid, $org, $locus, $contig, $beg, $end, $len);
114 :    
115 :     my $regions = {};
116 :     foreach my $tbl_file (@tbl_files)
117 :     {
118 :     my $num_fids = 0;
119 :     open(TBL,"<$tbl_file") || die "could not open $tbl_file";
120 :     print STDERR "Loading $tbl_file ...\n" if $ENV{VERBOSE};
121 :    
122 :     while (defined($entry = <TBL>))
123 :     {
124 :     ++$num_fids;
125 :     chomp $entry;
126 :     if ($entry =~ /^(\S+)\s+(\S+)/)
127 :     {
128 :     ($fid, $locus) = ($1, $2);
129 :     if ($fid !~ m/^fig\|(\d+\.\d+)\.(rna|peg|orf)\.(\d+)$/)
130 :     {
131 :     die "Invalid FID $fid";
132 :     }
133 :    
134 :     ($contig, $beg, $end) = $fig->boundaries_of($locus);
135 :     $len = 1 + abs($end-$beg);
136 :    
137 :     unless (defined($regions->{$contig})) { $regions->{$contig} = []; }
138 :     push(@ { $regions->{$contig} }, [ &FIG::min($beg, $end), &FIG::max($beg, $end), $fid ] );
139 :     }
140 :     else
141 :     {
142 :     print STDERR "Skipping invalid entry $entry\n";
143 :     }
144 :     }
145 :     print STDERR "Loaded $num_fids features from $tbl_file\n\n" if $ENV{VERBOSE};
146 :     close(TBL);
147 :     }
148 :    
149 :     foreach my $contig (keys(%$regions))
150 :     {
151 :     my $x = [ sort { $a->[LEFT] <=> $b->[LEFT] } @ { $regions->{$contig} } ];
152 :     for (my $i=1; $i < @$x; ++$i)
153 :     {
154 :     while (($i < @$x) && ($x->[$i-1]->[RIGHT] >= $x->[$i]->[LEFT]))
155 :     {
156 :     print STDERR "Merging $i:\t$x->[$i-1]->[FID]\t$x->[$i]->[FID]\n"
157 :     if ($ENV{VERBOSE} && ($ENV{VERBOSE} > 1));
158 :     $x->[$i-1]->[RIGHT] = &FIG::max( $x->[$i-1]->[RIGHT], $x->[$i]->[RIGHT] );
159 :     splice @$x, $i, 1;
160 :     }
161 :     }
162 :    
163 :     $regions->{$contig} = $x;
164 :     print STDERR "After merger, $contig has ", (scalar @$x), " regions\n\n" if $ENV{VERBOSE};
165 :     }
166 :    
167 :     return $regions;
168 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3