[Bio] / FigKernelPackages / TBLstuff.pm Repository:
ViewVC logotype

Annotation of /FigKernelPackages/TBLstuff.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 package TBLstuff;
2 :    
3 :     use Data::Dumper;
4 :     use Carp;
5 :    
6 :     sub current_features_of {
7 :     my($fig,$genome,@types) = @_;
8 :    
9 :     if (! @types) { @types = ('peg','rna') }
10 :     my @curr_features = ();
11 :     foreach my $type (@types)
12 :     {
13 :     foreach my $fid ($fig->all_features($genome,$type))
14 :     {
15 :     my $loc = $fig->feature_location($fid);
16 :     my @aliases = $fig->feature_aliases($fid);
17 :     push(@curr_features,[$type,$loc,{ id => $fid, aliases => \@aliases }]);
18 :     }
19 :     }
20 :     return wantarray ? @curr_features : \@curr_features;
21 :     }
22 :    
23 :     sub flat_tbl_to_features {
24 :     my(@files) = @_;
25 :    
26 :     my @features = ();
27 :     foreach my $file (@files)
28 :     {
29 :     if (open(TMP,"<$file"))
30 :     {
31 :     while (defined($_ = <TMP>))
32 :     {
33 :     chomp;
34 :     my($fid,$loc,@aliases) = split(/\t/,$_);
35 :     if ($fid && ($fid =~ /^fig\|\d+\.\d+\.([a-zA-Z]+)/))
36 :     {
37 :     push(@features,[$1,$loc,{ id => $fid, aliases => \@aliases }]);
38 :     }
39 :     }
40 :     close(TMP);
41 :     }
42 :     }
43 :     return wantarray ? @features : \@features;
44 :     }
45 :    
46 :     sub compare_sets_of_features {
47 :     my($fig,$set1,$set2) = @_;
48 :    
49 :     my @set1 = map { my($contig,$beg,$end) = $fig->boundaries_of($_->[1]);
50 :     my($left,$right) = sort { $a <=> $b } ($beg,$end);
51 :     [1,$_,$contig,$beg,$end,$left,$right,&frame($_->[0],$beg,$end)] } @$set1;
52 :    
53 :     my @set2 = map { my($contig,$beg,$end) = $fig->boundaries_of($_->[1]);
54 :     my($left,$right) = sort { $a <=> $b } ($beg,$end);
55 :     [2,$_,$contig,$beg,$end,$left,$right,&frame($_->[0],$beg,$end)] } @$set2;
56 :    
57 :     my @merged = sort { ($a->[1]->[0] cmp $b->[1]->[0]) or # type
58 :     ($a->[2] cmp $b->[2]) or # contig
59 :     ($a->[7] <=> $b->[7]) or # frame
60 :     ($a->[5] <=> $b->[5]) # left coord
61 :     } (@set1,@set2);
62 :    
63 :     my @output;
64 :     my $x = shift @merged;
65 :     while ($x)
66 :     {
67 :     print STDERR "Looking at ",&Dumper($x);
68 :    
69 :     if (! (@merged && &comparable($x,$merged[0])))
70 :     {
71 :     print STDERR "not comparable\n";
72 :     push(@output,['',
73 :     $x->[2], # contig
74 :     $x->[3] + $x->[4], # mid*2
75 :     ($x->[0] == 1) ? ($x,undef) : (undef,$x)
76 :     ]
77 :     );
78 :     $x = shift @merged;
79 :     }
80 :     else
81 :     {
82 :     my $y = shift @merged;
83 :     print STDERR "comparable",&Dumper($y);
84 :    
85 :     if ($x->[1]->[0] eq 'peg')
86 :     {
87 :     if ($x->[4] == $y->[4]) # if ends match
88 :     {
89 :     my $relation = ($x->[3] == $y->[3]) ? 'same' : 'diff-start';
90 :     push(@output,[$relation,
91 :     $x->[2], # contig
92 :     $x->[3] + $x->[4], # mid*2
93 :     ($x->[0] < $y->[0]) ? ($x,$y) : ($y,$x) # entries
94 :     ]);
95 :     $x = shift @merged;
96 :     }
97 :     else
98 :     {
99 :     ($x,$y) = ($y,$x) if ($x->[6] > $y->[6]); # flip if $y is longer
100 :    
101 :     push(@output,['',
102 :     $x->[2], # contig
103 :     $x->[3] + $x->[4], # mid*2
104 :     ($x->[0] == 1) ? ($x,undef) : (undef,$x)
105 :     ]
106 :     );
107 :     $x = $y;
108 :     }
109 :     }
110 :     elsif (&correspond($x,$y))
111 :     {
112 :     my $relation = (($x->[3] == $y->[3]) && ($x->[4] == $y->[4])) ? 'same' : 'diff-ends';
113 :     push(@output,[$relation,
114 :     $x->[2], # contig
115 :     $x->[3] + $x->[4], # mid*2
116 :     ($x->[0] < $y->[0]) ? ($x,$y) : ($y,$x) # entries
117 :     ]);
118 :     $x = shift @merged;
119 :     }
120 :     else
121 :     {
122 :     ($x,$y) = ($y,$x) if ($x->[6] > $y->[6]); # flip if $y is longer
123 :    
124 :     push(@output,['',
125 :     $x->[2], # contig
126 :     $x->[3] + $x->[4], # mid*2
127 :     ($x->[0] == 1) ? ($x,undef) : (undef,$x)
128 :     ]
129 :     );
130 :     $x = $y;
131 :     }
132 :     }
133 :     }
134 :     return wantarray ? @output : \@output;
135 :     }
136 :    
137 :     sub comparable {
138 :     my($x,$y) = @_;
139 :    
140 :     return (($x->[0] != $y->[0]) && # diff sets
141 :     ($x->[1]->[0] eq $y->[1]->[0]) && # same type
142 :     ($x->[2] eq $y->[2]) && # same contig
143 :     ($x->[7] == $y->[7])); # same frame
144 :     }
145 :    
146 :     sub correspond {
147 :     my($x,$y) = @_; # testing if two non-pegs overlap enough
148 :    
149 :     my $overlap = 1 + &FIG::min($x->[6],$y->[6]) -
150 :     &FIG::max($x->[5],$y->[5]);
151 :     my $maxln = &FIG::max(($x->[6] - $x->[5])+1, (($y->[6] - $y->[5])+1));
152 :     return ($overlap >= (0.8 * $maxln));
153 :     }
154 :    
155 :     sub frame {
156 :     my($type,$beg,$end) = @_;
157 :    
158 :     if ($type eq 'peg')
159 :     {
160 :     return (($end % 3)+1) * ($end <=> $beg);
161 :     }
162 :     else
163 :     {
164 :     return ($end <=> $beg);
165 :     }
166 :     }
167 :    
168 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3