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

Annotation of /FigKernelScripts/make_probes_to_genes.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 use strict;
2 :     use SeedEnv;
3 :     use gjoseqlib;
4 :     use Data::Dumper;
5 : olson 1.7 # Fig was only used for testing deleted pegs; we now require
6 :     # the calling code to elide removed deleted genomes from any tbl files provided.
7 :     #use FIG;
8 :     #my $fig = new FIG;
9 : overbeek 1.1
10 : olson 1.6 my $usage = "usage: make_probes_to_genes Probes Contigs Tbl peg.probe.table probe.occ.table [more tbl files]";
11 : overbeek 1.1 my($probesF,$contigsF,$tblF,$peg_probeF,$probe_occF);
12 :    
13 :     (
14 :     ($probesF = shift @ARGV) &&
15 :     ($contigsF = shift @ARGV) &&
16 :     ($tblF = shift @ARGV) &&
17 :     ($peg_probeF = shift @ARGV) &&
18 :     ($probe_occF = shift @ARGV)
19 :     )
20 :     || die $usage;
21 : olson 1.6 my @other_tblF = @ARGV;
22 : overbeek 1.1
23 :     my @contigs = map { $_->[2] = lc $_->[2]; $_ } &gjoseqlib::read_fasta($contigsF);
24 :    
25 :     my @probes = ();
26 : olson 1.5 my %probelens;
27 :    
28 : olson 1.6 my $fh;
29 :     open($fh, "<", $probesF) or die "Cannot read $probesF: $!";
30 :    
31 :     while (<$fh>)
32 : overbeek 1.1 {
33 :     if ($_ =~ /^(\S+)\t(\S+)/)
34 :     {
35 :     my $probe_id = $1;
36 :     my $seq = lc $2;
37 :     my $seqR = lc &SeedUtils::rev_comp($seq);
38 :     push(@probes,[$probe_id,$seq,$seqR]);
39 : olson 1.5 $probelens{length($seq)}++;
40 :     }
41 :     }
42 : olson 1.6 close($fh);
43 : olson 1.5 my %contig_index;
44 :     my $probe_size;
45 :    
46 :     #
47 :     # If we have a single probe length, create an index of the contigs
48 :     # and use for later lookups.
49 :     #
50 :     if (keys %probelens == 1)
51 :     {
52 :     $probe_size = (keys %probelens)[0];
53 :     for my $contig (@contigs)
54 :     {
55 :     my($id, undef, $seq) = @$contig;
56 :     for (my $i = 0; $i < length($seq) - $probe_size + 1; $i++)
57 :     {
58 :     my $s = substr($seq, $i, $probe_size);
59 :     push(@{$contig_index{$s}}, [$id, $i]);
60 :     }
61 : overbeek 1.1 }
62 :     }
63 : overbeek 1.4 #print STDERR "got probes\n";
64 : olson 1.5
65 : overbeek 1.1 #my @pegs = map { ($_ =~ /^(\S+)\t(\S+)_(\d+)_(\d+)\s/) ? [$1,[$2,$3,$4]] : () } `cat $tblF`;
66 :     my @pegs;
67 : olson 1.6 for my $tbl_file ($tblF, @other_tblF)
68 : overbeek 1.1 {
69 : olson 1.6 open($fh, "<", $tbl_file) or die "Cannot open $tbl_file: $!";
70 :     while (<$fh>)
71 : overbeek 1.1 {
72 : olson 1.6 if ($_ =~ /^(fig\|\d+\.\d+\.(peg|rna)\.\d+)\t(\S+)/)
73 : overbeek 1.4 {
74 : olson 1.6 my $peg = $1;
75 :     my($c1,$b1,$e1) = &boundaries($3);
76 : olson 1.7
77 :     # It is responsibility of calling code to provide tbl files with
78 :     # deleted genomes removed.
79 :     # if (! $fig->is_deleted_fid($peg))
80 :     # {
81 :     # push(@pegs,[$peg,[$c1,$b1,$e1]]);
82 :     # }
83 :     push(@pegs,[$peg,[$c1,$b1,$e1]]);
84 : overbeek 1.4 }
85 : overbeek 1.1 }
86 : olson 1.6 close($fh);
87 : overbeek 1.1 }
88 :    
89 :     my %multiple;
90 :     my @all_hits = ();
91 :     foreach my $probe (@probes)
92 :     {
93 :     my($probe_id,$seq,$seqR) = @$probe;
94 : olson 1.5 my @matches = &matches_in_contigs($seq,\@contigs,'+', \%contig_index);
95 :     push(@matches,&matches_in_contigs($seqR,\@contigs,'-', \%contig_index));
96 : overbeek 1.1 if (@matches == 0)
97 :     {
98 :     print STDERR "$probe_id does not occur\n";
99 :     }
100 :     elsif (@matches > 1)
101 :     {
102 :     $multiple{$probe_id} = scalar @matches;
103 :     print STDERR "$probe_id occurs ",scalar @matches," times: ",join(",",map {join(":",@$_) } @matches),"\n";
104 :     }
105 :     push(@all_hits,map { [$probe_id,$_] } @matches);
106 :     }
107 :     @all_hits = sort { ($a->[1]->[0] cmp $b->[1]->[0]) or ($a->[1]->[1] <=> $b->[1]->[1]) } @all_hits;
108 : overbeek 1.4 @pegs = sort { ($a->[1]->[0] cmp $b->[1]->[0]) or
109 :     (&min($a->[1]->[1],$a->[1]->[2]) <=> &min($b->[1]->[1],$b->[1]->[2])) }
110 :     @pegs;
111 : overbeek 1.1 # print STDERR &Dumper(\@all_hits,\@pegs);
112 : overbeek 1.4
113 : overbeek 1.1 my @corr = ();
114 : overbeek 1.4 foreach my $hit_tuple (@all_hits)
115 : overbeek 1.1 {
116 : overbeek 1.4 my($probe_id,$hit) = @$hit_tuple;
117 : overbeek 1.1 my($contigH,$bH,$eH,$strandH) = @$hit;
118 : overbeek 1.4 foreach my $peg_tuple (@pegs)
119 : overbeek 1.1 {
120 : overbeek 1.4 my($peg,$loc) = @$peg_tuple;
121 :     my($contigP,$bP,$eP) = @$loc;
122 : dejongh 1.8 if ($contigH eq $contigP && &SeedUtils::between($bP,$bH,$eP))
123 : overbeek 1.1 {
124 : overbeek 1.4 my $n;
125 :     if ((($bP < $eP) && (($bH - $bP) > 3)) ||
126 :     (($bP > $eP) && (($bH - $eP) > 3)))
127 :     {
128 :     $n = $multiple{$probe_id} ? $multiple{$probe_id} : 0;
129 :     }
130 :     my $strandP;
131 :     if ($bP > $eP)
132 :     {
133 :     ($bP,$eP) = ($eP,$bP);
134 :     $strandP = '-';
135 :     }
136 :     else
137 :     {
138 :     $strandP = '+';
139 :     }
140 :     push(@corr,[$peg,$probe_id,$n,$strandP,$strandH,$contigH,$bH,$eH]);
141 : overbeek 1.1 }
142 :     }
143 :     }
144 :    
145 :     open(P2P,">$peg_probeF") || die "could not open $peg_probeF";
146 :     open(OCC,">$probe_occF") || die "could not open $probe_occF";
147 :     my %seen;
148 : overbeek 1.4 my %probes_in_peg;
149 : overbeek 1.1 foreach my $x (sort { &SeedUtils::by_fig_id($a->[0],$b->[0]) or ($a->[1] cmp $b->[1]) } @corr)
150 :     {
151 :     my($peg,$probe,$n,$strandH,$strandP,$contigH,$bH,$eH) = @$x;
152 :     if ($strandH eq $strandP)
153 :     {
154 :     print P2P "$peg\t$probe\n";
155 :     if (! $seen{$probe})
156 :     {
157 :     $seen{$probe} = 1;
158 :     print OCC "$probe\t$n\t$contigH,$bH,$eH\n";
159 :     }
160 :     $probes_in_peg{$x->[0]} = 1;
161 :     }
162 :     }
163 :    
164 :     foreach my $peg (map { $_->[0] } @pegs)
165 :     {
166 :     if (! $probes_in_peg{$peg})
167 :     {
168 :     print STDERR "$peg is not matched by probes\n";
169 :     }
170 :     }
171 :    
172 :     sub matches_in_contigs {
173 : olson 1.5 my($seq,$contigs,$strand, $idx) = @_;
174 : overbeek 1.1 my @hits;
175 : olson 1.5
176 :     if (%$idx)
177 :     {
178 :     my $hits = $idx->{$seq};
179 :     if ($hits)
180 :     {
181 :     for my $hit (@$hits)
182 :     {
183 :     my($ctg, $loc) = @$hit;
184 :     push(@hits, [$ctg, $loc + 1, $loc + length($seq), $strand]);
185 :     }
186 :     }
187 :     return @hits;
188 :     }
189 : overbeek 1.1 foreach my $contig (@$contigs)
190 :     {
191 :     my($idC,undef,$seqC) = @$contig;
192 :     my $off = 0;
193 :     while (($_ = index($seqC,$seq,$off)) >= 0)
194 :     {
195 :     push(@hits,[$idC,$_+1,$_+length($seq),$strand]);
196 :     $off = $_ + 1;
197 :     }
198 :     }
199 :     return @hits;
200 :     }
201 :    
202 :     sub boundaries {
203 :     if (!UNIVERSAL::isa($_[0],__PACKAGE__)) {
204 :     my ($package, $filename, $line) = caller;
205 :     # warn "Deprecated boundaries_of called from $filename line $line.";
206 :     } else {
207 :     shift;
208 :     }
209 :     my($location) = (@_ == 1) ? $_[0] : $_[1];
210 :     my($contigQ);
211 :    
212 :     if (defined($location))
213 :     {
214 :     my @exons = split(/\s*,\s*/,$location);
215 :     my($contig, $beg, $end);
216 :    
217 :     if (($exons[0] =~ /^(\S+)_(\d+)_\d+$/))
218 :     {
219 :     ($contig, $beg) = ($1,$2);
220 :     $contigQ = quotemeta $contig;
221 :    
222 :     if ($exons[$#exons] =~ /^$contigQ\_\d+_(\d+)$/)
223 :     {
224 :     $end = $1;
225 :     if ($beg > 0 && $end > 0)
226 :     {
227 :     my $strand = (($beg < $end) ? qq(+) : qq(-));
228 :     return ($contig, $beg, $end, $strand);
229 :     }
230 :     }
231 :     }
232 :     Cluck("Could not parse loc=$location.") if T(0);
233 :     }
234 :     return ();
235 :     }
236 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3