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

Annotation of /FigKernelPackages/Contigs.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 package Contigs;
2 :    
3 :     #
4 :     # Build and manage information about a set of contigs. They can be
5 :     # supplied as file names (in which case the sequence locations are
6 :     # indexed), or as id/sequence pairs ( [ id, seq ], ... ) or sequence
7 :     # entries ( [ id, def, seq ] ), in which case an index of references
8 :     # is built. When a sequence on disk is used, it is cached in memory.
9 :     # (This should be replaced by indexing subsequences on the disk, but
10 :     # for now ....)
11 :     #
12 :     # New contigs object:
13 :     #
14 :     # $contigs = Contigs->new( $file, ... ) # files
15 :     # $contigs = Contigs->new( [ $file, ... ] ) # ref to file list
16 :     # $contigs = Contigs->new( [ $id, $seq ], ... ) # id-seq pairs
17 :     # $contigs = Contigs->new( [ [ $id, $seq ], ... ] ) # ref to pair list
18 :     # $contigs = Contigs->new( [ $id, $def, $seq ], ... ) # id-def-seq triples
19 :     # $contigs = Contigs->new( [ [ $id, $def, $seq ], ... ] ) # ref to triple list
20 :     #
21 :     # Contigs access functions:
22 :     #
23 :     # @ids = $contigs->ids( );
24 :     # \@ids = $contigs->ids( );
25 :     # $length = $contigs->length( $id );
26 :     # $sequence = $contigs->sequence( $id );
27 :     # \$sequence = $contigs->sequenceP( $id );
28 :     # $subseq = $contigs->subseq( $id, $beg, $end );
29 :     #
30 :     # Test scripts:
31 :     #
32 :     # perl -e 'use Contigs; $contigs=Contigs->new("Vibrio/contigs"); print join "\n", $contigs->ids(), ""'
33 :     #
34 :     # perl -e 'use Contigs; $contigs=Contigs->new("Vibrio/contigs"); $c="Contig83"; print $contigs->sequence($c), "\n", ${$contigs->sequenceP($c)}, "\n"'
35 :     #
36 :     # perl -e 'use Contigs; $contigs=Contigs->new("Vibrio/contigs"); $c="Contig83"; print $contigs->length($c), "\n", $contigs->subseq($c,1,50), "\n", $contigs->subseq($c,50,1), "\n"'
37 :     # 1097
38 :     # ggtTTACTCAAAATATTAAATCAGCAGGTCAATTTGAAAGCTACgaagaT
39 :     # AtcttcGTAGCTTTCAAATTGACCTGCTGATTTAATATTTTGAGTAAacc
40 :     #
41 :     # perl -e 'use gjoseqlib; use Contigs; $contigs=Contigs->new(read_fasta("Vibrio/contigs")); $c="Contig83"; print $contigs->length($c), "\n", $contigs->subseq($c,1,50), "\n", $contigs->subseq($c,50,1), "\n"'
42 :     # 1097
43 :     # ggtTTACTCAAAATATTAAATCAGCAGGTCAATTTGAAAGCTACgaagaT
44 :     # AtcttcGTAGCTTTCAAATTGACCTGCTGATTTAATATTTTGAGTAAacc
45 :     #
46 :    
47 :     use strict;
48 :     use gjoseqlib; # read_fasta(), DNA_subseq()
49 :     use Data::Dumper;
50 :    
51 :     #
52 :     # $contigs = Contigs->new( $file, ... )
53 :     # $contigs = Contigs->new( [ $file, ... ] )
54 :     # $contigs = Contigs->new( [ $id, $seq ], ... )
55 :     # $contigs = Contigs->new( [ [ $id, $seq ], ... ] )
56 :     #
57 :    
58 :     sub new
59 :     {
60 :     my $class = shift;
61 :     my $self = {};
62 :    
63 :     my $type;
64 :     my @files;
65 :     my @seqs;
66 :    
67 :     # Scalars are file names
68 :     if ( ! ref $_[0] )
69 :     {
70 :     @files = grep { -f } @_;
71 :     if ( @files != @_ )
72 :     {
73 :     print STDERR "The following file(s) were not found:\n";
74 :     foreach ( grep { ! -f } @_ ) { print STDERR " '$_'\n" }
75 :     }
76 :     $type = 'file';
77 :     }
78 :    
79 :     # Scalars in an array ref can be file names
80 :     elsif ( @_ == 1
81 :     && ! ref $_[0]->[0]
82 :     && ( @files = grep { -f } @{$_[0]} )
83 :     && ( @files == @{$_[0]} )
84 :     ) { $type = 'file' }
85 :    
86 :     # Array refs of scalars can be sequence entries
87 :     elsif ( ! ref $_[0]->[0]
88 :     && ( @seqs = grep { @$_ == 2 || @$_ == 3 } @_ )
89 :     && ( @seqs == @_ )
90 :     ) { $type = 'seqs' }
91 :    
92 :     # Arrays of array refs can be sequence entries
93 :     elsif ( ref $_[0]->[0]
94 :     && ( @seqs = grep { @$_ == 2 || @$_ == 3 } @{$_[0]} )
95 :     && ( @seqs = @{$_[0]} )
96 :     ) { $type = 'seqs' }
97 :    
98 :     else
99 :     {
100 :     print STDERR "Contigs->new called wtih arguments that could no be interpretted as files or sequences.\n";
101 :     print STDERR Dumper( @_ );
102 :     return undef;
103 :     }
104 :    
105 :     $self->{ type } = $type;
106 :    
107 :     if ( $type eq 'file' )
108 :     {
109 :     my ( $c_index, $c_length ) = index_contig_file( @files );
110 :     $self->{ contig_index } = $c_index;
111 :     $self->{ contig_len } = $c_length;
112 :     $self->{ contigs } = {};
113 :     }
114 :     elsif ( $type eq 'seqs' )
115 :     {
116 :     my %contigs = map { $_->[0] => $_ } @seqs;
117 :     $self->{ contigs } = \%contigs;
118 :    
119 :     my %contig_len = map { $_->[0] => length $_->[-1] } @seqs;
120 :     $self->{ contig_len } = \%contig_len;
121 :     }
122 :    
123 :     bless $self, $class;
124 :     }
125 :    
126 :    
127 :     sub ids
128 :     {
129 :     my ( $self ) = @_;
130 :     my @ids = sort { lc $a cmp lc $b } keys %{ $self->{ contig_len } };
131 :     wantarray ? @ids : \@ids;
132 :     }
133 :    
134 :    
135 :     sub length
136 :     {
137 :     my ( $self, $contig ) = @_;
138 :     my $len = $self->{ contig_len }->{ $contig };
139 :     print STDERR "Unable to evalute length( '$contig' ).\n" if ! defined $len;
140 :     return $len;
141 :     }
142 :    
143 :    
144 :     sub sequence
145 :     {
146 :     my ( $self, $contig ) = @_;
147 :    
148 :     my $c_data = $self->{ contigs }->{ $contig };
149 :    
150 :     if ( ! $c_data && $self->{ type } eq 'file' )
151 :     {
152 :     my $c_index = $self->{ contig_index }->{ $contig };
153 :     $c_data = read_contig_from_file( @$c_index ) if $c_index;
154 :     $self->{ contigs }->{ $contig } = $c_data if $c_data;
155 :     }
156 :    
157 :     if ( ! $c_data )
158 :     {
159 :     print STDERR "Unable to evalute sequence( '$contig' ).\n";
160 :     return undef;
161 :     }
162 :    
163 :     $c_data->[-1]
164 :     }
165 :    
166 :    
167 :     sub sequenceP
168 :     {
169 :     my ( $self, $contig ) = @_;
170 :    
171 :     my $c_data = $self->{ contigs }->{ $contig };
172 :    
173 :     if ( ! $c_data && $self->{ type } eq 'file' )
174 :     {
175 :     my $c_index = $self->{ contig_index }->{ $contig };
176 :     $c_data = read_contig_from_file( @$c_index ) if $c_index;
177 :     $self->{ contigs }->{ $contig } = $c_data if $c_data;
178 :     }
179 :    
180 :     if ( ! $c_data )
181 :     {
182 :     print STDERR "Unable to evalute sequenceP( '$contig' ).\n";
183 :     return undef;
184 :     }
185 :    
186 :     \$c_data->[-1]
187 :     }
188 :    
189 :    
190 :     sub subseq
191 :     {
192 :     my ( $self, $contig, $beg, $end ) = @_;
193 :    
194 :     my $c_data = $self->{ contigs }->{ $contig };
195 :    
196 :     if ( ! $c_data && $self->{ type } eq 'file' )
197 :     {
198 :     my $c_index = $self->{ contig_index }->{ $contig };
199 :     $c_data = read_contig_from_file( @$c_index ) if $c_index;
200 :     $self->{ contigs }->{ $contig } = $c_data if $c_data;
201 :     }
202 :    
203 :     if ( ! $c_data )
204 :     {
205 :     print STDERR "Unable to evalute subseq( '$contig', $beg, $end ).\n";
206 :     return undef;
207 :     }
208 :    
209 :     DNA_subseq( \$c_data->[-1], $beg, $end )
210 :     }
211 :    
212 :    
213 :     #===============================================================================
214 :     # Internal functions:
215 :     #
216 :     #
217 :     # ( \%records, \%lengths ) = index_contig_file( $file, ... )
218 :     #
219 :     # $record = [ $id, $def, $file, $offset, $seqbytes ];
220 :     #
221 :     #
222 :     # \@seq_entry = read_contig_from_file( $id, $def, $file, $offset, $seqbytes )
223 :     #
224 :     # $seq_entry = [ $id, $def, $seq ]
225 :     #
226 :     #===============================================================================
227 :    
228 :     sub index_contig_file
229 :     {
230 :     my %c_index;
231 :     my %c_length;
232 :     foreach my $file ( @_ )
233 :     {
234 :     my $prev_id = "";
235 :     open GREP, "grep -b '^>' < '$file' |";
236 :     while ( $_ = <GREP> )
237 :     {
238 :     s/^(\d+)\://;
239 :     my $off = $1 + 0;
240 :     my $prev_ind = $c_index{ $prev_id };
241 :     $prev_ind->[4] = $off - $prev_ind->[3] if $prev_ind;
242 :    
243 :     $off += CORE::length($_);
244 :     chomp;
245 :     my ( $id, $def ) = /^>(\S+)(\s+(\S.*)?)?$/;
246 :     $c_index{ $id } = [ $id, $def, $file, $off, undef ];
247 :     $prev_id = $id;
248 :     }
249 :     close GREP;
250 :    
251 :     my $prev_ind = $c_index{ $prev_id };
252 :     $prev_ind->[4] = (-s $file) - $prev_ind->[3] if $prev_ind;
253 :    
254 :     open SIZE, "fasta_size < '$file' |";
255 :     while ( $_ = <SIZE> )
256 :     {
257 :     my ( $size, $id ) = /^\s*(\d+)\s+(\S+)/;
258 :     $c_length{ $id } = $size if $id && $c_index{ $id };
259 :     }
260 :     close SIZE;
261 :     }
262 :    
263 :     ( \%c_index, \%c_length );
264 :     }
265 :    
266 :    
267 :     sub read_contig_from_file
268 :     {
269 :     my ( $id, $def, $file, $offset, $seqbytes ) = @_;
270 :    
271 :     my $seq;
272 :     open FILE, "<$file" or return undef;
273 :     seek FILE, $offset, 0;
274 :     read FILE, $seq, $seqbytes;
275 :     close FILE;
276 :     $seq =~ tr/ \n\r\t//d;
277 :    
278 :     [ $id, $def, $seq ]
279 :     }
280 :    
281 :    
282 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3