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

Annotation of /FigKernelScripts/generate_intergenic_blocks.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : olson 1.3 #
2 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
3 :     # for Interpretations of Genomes. All Rights Reserved.
4 :     #
5 :     # This file is part of the SEED Toolkit.
6 :     #
7 :     # The SEED Toolkit is free software. You can redistribute
8 :     # it and/or modify it under the terms of the SEED Toolkit
9 :     # Public License.
10 :     #
11 :     # You should have received a copy of the SEED Toolkit Public License
12 :     # along with this program; if not write to the University of Chicago
13 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
14 :     # Genomes at veronika@thefig.info or download a copy from
15 :     # http://www.theseed.org/LICENSE.TXT.
16 :     #
17 :    
18 : overbeek 1.1
19 :     use FIG;
20 :     my $fig = new FIG;
21 :    
22 :     my $usage = "usage: generate_intergenic_blocks Dir";
23 :    
24 :     (
25 :     ($dir = shift @ARGV)
26 :     )
27 :     || die $usage;
28 :    
29 :     if (! -d $dir) { die "$dir must exist with the gene blocks in it" }
30 :    
31 :     open(REG,"<$dir/region.tbl") || die "could not open $dir/region.tbl";
32 :     while (defined($_ = <REG>))
33 :     {
34 :     ($peg,undef,$contig,$beg,$end,$blockid) = split(/\t/,$_);
35 :     $in{$peg} = [$blockid,$contig,$beg,$end];
36 :     push(@{$contains{$blockid}},$peg);
37 :     }
38 :     close(REG);
39 :    
40 :     open(BLOCK,"<$dir/block.tbl") || die "could not open $dir/block.tbl";
41 :     $max = 0;
42 :     while (defined($_ = <BLOCK>))
43 :     {
44 :     chop;
45 :     ($blockid,$gene) = split(/\t/,$_);
46 :     $gene_of{$blockid} = $gene;
47 :     $max = &FIG::max($max,$blockid);
48 :     }
49 :     close(BLOCK);
50 :     $blockid = $max+1;
51 :    
52 :     open(BLOCKS,">$dir/intergenic_block.tbl") || die "could not open intergenic_block.tbl";
53 :     open(REGIONS,">$dir/intergenic_region.tbl") || die "could not open intergenic_region.tbl";
54 :    
55 :    
56 : overbeek 1.2 open(NEIGH,">neigh") || die "aborted";
57 :     foreach $peg (keys(%in))
58 :     {
59 :     ($peg1,$peg2) = &conn($peg);
60 :    
61 :     $n = $in{$peg}->[0];
62 :     if ($peg1 && ($n1 = $in{$peg1})) { $precedes{$n}->{$n1->[0]}++ }
63 :     if ($peg2 && ($n1 = $in{$peg2})) { $follows{$n}->{$n1->[0]}++ }
64 :     $peg1 = $peg1 ? $peg1 : "";
65 :     $peg2 = $peg2 ? $peg2 : "";
66 :     print NEIGH "$peg\t$peg1\t$peg2\n";
67 :     }
68 :     #close(NEIGH);
69 :    
70 :     ####################
71 :     #open(NEIGH,"<neigh") || die "aborted";
72 :     #while (defined($_ = <NEIGH>))
73 : overbeek 1.1 #{
74 : overbeek 1.2 # if ($_ =~ /^(\S+)\t(\S+)\t(\S+)/)
75 :     # {
76 :     # ($peg,$peg1,$peg2) = ($1,$2,$3);
77 :     #
78 :     # $n = $in{$peg}->[0];
79 :     # if ($peg1 && ($n1 = $in{$peg1})) { $precedes{$n}->{$n1->[0]}++ }
80 :     # if ($peg2 && ($n1 = $in{$peg2})) { $follows{$n}->{$n1->[0]}++ }
81 :     # }
82 : overbeek 1.1 #}
83 :     #close(NEIGH);
84 :     ####################
85 :    
86 :     foreach $n (keys(%follows))
87 :     {
88 :     $x = $follows{$n};
89 :     @poss = sort { $x->{$b} <=> $x->{$a} } keys(%$x);
90 :     if ((@poss == 1) || ($x->{$poss[0]} > $x->{$poss[1]}) && ($x->{$poss[0]} > 2))
91 :     {
92 :     $key = join("\t",sort ($n,$poss[0]));
93 :     $pair{$key}++;
94 :     }
95 :     }
96 :    
97 :     foreach $n (keys(%precedes))
98 :     {
99 :     $x = $precedes{$n};
100 :     @poss = sort { $x->{$b} <=> $x->{$a} } keys(%$x);
101 :     if ((@poss == 1) || ($x->{$poss[0]} > $x->{$poss[1]}) && ($x->{$poss[0]} > 2))
102 :     {
103 :     $key = join("\t",sort ($n,$poss[0]));
104 :     $pair{$key}++;
105 :     }
106 :     }
107 :    
108 :     foreach $tuple (keys(%pair))
109 :     {
110 :     ($n1,$n2) = split(/\t/,$tuple);
111 :     $pegs1 = $contains{$n1};
112 :     $pegs2 = $contains{$n2};
113 :     @poss = ();
114 :    
115 :     foreach $peg1 (@$pegs1)
116 :     {
117 :     $genome1 = &FIG::genome_of($peg1);
118 :     (undef,$contig1,$beg1,$end1) = @{$in{$peg1}};
119 :     foreach $peg2 (@$pegs2)
120 :     {
121 :     $genome2 = &FIG::genome_of($peg2);
122 :     if ($genome1 eq $genome2)
123 :     {
124 :     (undef,$contig2,$beg2,$end2) = @{$in{$peg2}};
125 :     if ($contig1 eq $contig2)
126 :     {
127 :     # warn "$peg1 $peg2 $contig1 $beg1 $end1 $contig2 $beg2 $end2\n";
128 :     if (($beg1 < $beg2) && ($end1 < $beg2) && ($beg1 < $end2) && ($end1 < $end2))
129 :     {
130 :     $start = &FIG::max($beg1,$end1) + 1;
131 :     $end = &FIG::min($beg2,$end2) - 1;
132 :     if ($start < $end)
133 :     {
134 :     push(@poss,[$peg1,$peg2,$contig1,$start,$end,
135 :     $fig->dna_seq($genome1,join("_",($contig1,$start,$end)))]);
136 :     }
137 :     }
138 :     elsif (($beg1 > $beg2) && ($end1 > $beg2) && ($beg1 > $end2) && ($end1 > $end2))
139 :     {
140 :     $start = &FIG::min($beg1,$end1) - 1;
141 :     $end = &FIG::max($beg2,$end2) + 1;
142 :     if ($start > $end)
143 :     {
144 :     push(@poss,[$peg1,$peg2,$contig1,$start,$end,
145 :     $fig->dna_seq($genome1,join("_",($contig1,$start,$end)))]);
146 :     }
147 :     }
148 :     # print &Dumper(\@poss);
149 :     }
150 :     }
151 :     }
152 :     }
153 :    
154 :     @poss = sort { length($a->[5]) <=> length($b->[5]) } @poss;
155 :     if ((@poss > 1) &&
156 :     (((length($poss[$#poss]->[5]) - length($poss[$#poss]->[0])) / length($poss[$#poss]->[5])) < 0.1))
157 :     {
158 :     my @reordered = map { [&pegN($_->[0]) . ":" . &pegN($_->[1]),
159 :     &FIG::genome_of($_->[0]),$_->[5],join("_",($_->[2],$_->[3],$_->[4]))] } @poss;
160 :     my @aligned = &align_DNA_for_a_set_of_seqs($fig,\@reordered);
161 :     $name = $gene_of{$n1} . "-" . $gene_of{$n2};
162 :     print BLOCKS "$blockid\t$name\n";
163 :    
164 :     foreach $tuple (@aligned)
165 :     {
166 :     &write_tuple($blockid,\*REGIONS,$tuple);
167 :     }
168 :     $blockid++;
169 :     }
170 :     }
171 :    
172 :     sub pegN {
173 :     my($peg) = @_;
174 :    
175 :     if ($peg =~ /(\d+)$/) { return $1 }
176 :     return undef;
177 :     }
178 :    
179 :     sub write_tuple {
180 :     my($blockid,$fh,$tuple) = @_;
181 :    
182 :     my($pegs,$genome,$ali_seq,$loc) = @$tuple;
183 :     if ($loc =~ /^(\S+)_(\d+)_(\d+)$/)
184 :     {
185 :     print $fh join("\t",($pegs,$genome,$1,$2,$3,$blockid,$ali_seq)),"\n";
186 :     }
187 :     }
188 :    
189 :    
190 :     sub align_DNA_for_a_set_of_seqs {
191 :     my($fig,$tuples) = @_;
192 :     my($tuple,$id,$dseq,$hdr,$aseq,%aligned);
193 :    
194 :     my $tmpS = "/tmp/seq$$";
195 :    
196 :     open(SEQS,">$tmpS.fasta") || die "could not open $tmpS";
197 :     foreach $tuple (@$tuples)
198 :     {
199 :     (undef,$id,$dseq) = @$tuple;
200 :     print SEQS ">$id\n$dseq\n";
201 :     }
202 :     close(SEQS);
203 :     system "$FIG_Config::ext_bin/clustalw -infile=$tmpS.fasta -align -outorder=aligned > /dev/null";
204 :     my @ali = `$FIG_Config::bin/clustal_to_fasta < $tmpS.aln`;
205 :     while (($hdr = shift @ali) && ($aseq = shift @ali))
206 :     {
207 :     if ($hdr =~ /^>(\S+)/)
208 :     {
209 :     chop $aseq;
210 :     $aligned{$1} = $aseq;
211 :     }
212 :     }
213 :     return map { my($pegs,$genome,undef,@rest) = @$_; $aligned{$genome} ? [$pegs,$genome,$aligned{$genome},@rest] : () } @$tuples;
214 :     }
215 :    
216 :     sub conn {
217 :     my($peg) = @_;
218 :     my($beg,$end,$i,@with_pos,@close_genes);
219 :    
220 :     (undef,$beg,$end) = $fig->boundaries_of($fig->feature_location($peg));
221 :     @close_genes = $fig->close_genes($peg,4000);
222 :     @with_pos = sort { $a->[0] <=> $b->[0] }
223 :     map {
224 :     my $peg1 = $_;
225 :     my(undef,$beg,$end) = $fig->boundaries_of($fig->feature_location($peg1));
226 :     [int($beg+$end)/2,$peg1]
227 :     }
228 :     @close_genes;
229 :     for ($i=0; ($i < @with_pos) && ($with_pos[$i]->[1] ne $peg); $i++) {}
230 :    
231 :     my($after,$before);
232 :     if (($i < $#with_pos) && ($beg < $end)) { $after = $with_pos[$i+1]->[1] }
233 :     if (($i < $#with_pos) && ($beg > $end)) { $before = $with_pos[$i+1]->[1] }
234 :     if (($i > 0) && ($beg < $end)) { $before = $with_pos[$i-1]->[1] }
235 :     if (($i > 0) && ($beg > $end)) { $after = $with_pos[$i-1]->[1] }
236 :     return ($before,$after);
237 :     }
238 :    
239 :     sub before {
240 :     my($n,$strand) = @_;
241 :     my($n1,$strand1);
242 :    
243 :     if ($strand eq "+")
244 :     {
245 :     $n1 = $really_precedes{$n};
246 :     $strand1 = ($really_follows{$n1} == $n) ? "+" : "-";
247 :     }
248 :     else
249 :     {
250 :     $n1 = $really_follows{$n};
251 :     $strand1 = ($really_precedes{$n1} == $n) ? "-" : "+";
252 :     }
253 :     return ($n1,$strand1);
254 :     }
255 :    
256 :     sub after {
257 :     my($n,$strand) = @_;
258 :     my($n1,$strand1);
259 :    
260 :     if ($strand eq "+")
261 :     {
262 :     $n1 = $really_follows{$n};
263 :     $strand1 = ($n1 && ($really_precedes{$n1} == $n)) ? "+" : "-";
264 :     }
265 :     else
266 :     {
267 :     $n1 = $really_precedes{$n};
268 :     $strand1 = ($n1 && ($really_follows{$n1} == $n)) ? "-" : "+";
269 :     }
270 :     return ($n1,$strand1);
271 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3