[Bio] / SeedViewer / MetagenomeAnalysis.pm Repository:
ViewVC logotype

Annotation of /SeedViewer/MetagenomeAnalysis.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : paarmann 1.1 package SeedViewer::MetagenomeAnalysis;
2 :    
3 : paarmann 1.3 # $Id: MetagenomeAnalysis.pm,v 1.2 2008/04/10 14:46:51 paarmann Exp $
4 : paarmann 1.1
5 :     use strict;
6 :     use warnings;
7 :    
8 :     use FIG_Config;
9 :     use DBI;
10 :    
11 :     1;
12 :    
13 :    
14 :    
15 : paarmann 1.3 use constant QUERY_DEFAULTS =>
16 :     { 1 => { evalue => '1e-05', align_len => 50 }, # RDP
17 :     2 => { evalue => '0.01' }, # SEED
18 :     3 => { evalue => '1e-05', align_len => 50 }, # Greengenes
19 :     4 => { evalue => '1e-05', align_len => 50 }, # LSU
20 :     5 => { evalue => '1e-05', align_len => 50 }, # SSU
21 :     6 => { evalue => '0.01' }, # Subsystem
22 :     };
23 :    
24 :    
25 :    
26 : paarmann 1.1
27 :     sub new {
28 :     my ($class, $job) = @_;
29 :    
30 :     # check job
31 :     unless (ref $job and $job->isa('RAST::Job') and $job->metagenome) {
32 :     return undef;
33 :     }
34 :    
35 :     # connect to database
36 :     my $dbh;
37 :     eval {
38 :    
39 :     my $host = $FIG_Config::mgrast_dbhost;
40 :     my $database = $FIG_Config::mgrast_db;
41 :     my $user = $FIG_Config::mgrast_dbuser;
42 :     my $password = '';
43 :    
44 :     $dbh = DBI->connect("DBI:mysql:database=$database;host=$host", $user, $password,
45 :     { RaiseError => 1, AutoCommit => 0, PrintError => 0 }) ||
46 :     die "database connect error.";
47 :     };
48 :     if ($@) {
49 :     warn "Unable to connect to metagenomics database: $@\n";
50 :     return undef;
51 :     }
52 :    
53 :     # create object
54 :     my $self = { job => $job,
55 :     dbh => $dbh,
56 :     key2taxa => undef,
57 : paarmann 1.3 query => {},
58 : paarmann 1.1 };
59 :     bless $self, $class;
60 :    
61 :     # load key2taxa mapping
62 :     $self->get_key2taxa_mapping();
63 :    
64 :     return $self;
65 :    
66 :     }
67 :    
68 :    
69 :    
70 :    
71 :     sub job {
72 :     return $_[0]->{job};
73 :     }
74 :    
75 :    
76 :     sub dbh {
77 :     return $_[0]->{dbh};
78 :     }
79 :    
80 :    
81 :     sub dbtable {
82 :     unless (defined $_[0]->{dbtable}) {
83 :     $_[0]->{dbtable} = 'tax_sim_'.$_[0]->job->id;
84 :     }
85 :     return $_[0]->{dbtable};
86 :     }
87 :    
88 :    
89 :     sub get_key2taxa_mapping {
90 :    
91 :     unless (defined($_[0]->{key2taxa})) {
92 :     my $sth = $_[0]->dbh->prepare("select dbkey, str from tax_item");
93 :     $sth->execute;
94 :     $_[0]->{key2taxa} = $sth->fetchall_hashref('dbkey');
95 :     }
96 :    
97 :     return $_[0]->{key2taxa};
98 :    
99 :     }
100 :    
101 :    
102 :     sub key2taxa {
103 :     if(defined $_[1]) {
104 :     my $t = $_[0]->{key2taxa}->{$_[1]}->{str};
105 :     $t =~ s/\t+$//;
106 :     return $t;
107 :     }
108 :     return '';
109 :     }
110 :    
111 :    
112 :     sub split_taxstr {
113 :     my @r = split(':', $_[1]);
114 :     return \@r;
115 :     }
116 :    
117 :     sub join_taxstr {
118 :     # do I really want an ending colon?
119 :     return join(':', @{$_[1]}).':';
120 :     }
121 :    
122 :    
123 :     sub evalue2log {
124 :     return 10 * (log($_[1]) / log(10));
125 :     }
126 :    
127 :     sub log2evalue {
128 :     return 10**($_[1]/10);
129 :     }
130 :    
131 :    
132 :     sub get_dataset_id {
133 : paarmann 1.3 # mapping contained in table db_type
134 :     # unstable, changes every time
135 : paarmann 1.1 my $db_ids = { 'SEED' => 2,
136 :     'Greengenes' => 3,
137 :     'RDP' => 1,
138 :     'LSU' => 4,
139 :     'SSU' => 5,
140 :     'Subsystem' => 6,
141 :     };
142 :     return $db_ids->{ $_[1] } || undef;
143 :     }
144 :    
145 :    
146 : paarmann 1.3 #******************************************************************************
147 :     #* MANAGING QUERY CRITERIA
148 :     #******************************************************************************
149 :    
150 :     =pod
151 :    
152 :     =item * B<query_evalue> (I<evalue>)
153 :    
154 :     Set/get the expectation value which is currently used to query the database.
155 :     Parameter I<evalue> has to be a float or in '1e-5'-like format or undef.
156 :    
157 :     =cut
158 :    
159 :     sub query_evalue {
160 :     if(scalar(@_)>1) {
161 :     $_[0]->{query}->{evalue} = $_[1];
162 :     }
163 :     return $_[0]->{query}->{evalue};
164 :     }
165 :    
166 :    
167 :     =pod
168 :    
169 :     =item * B<query_bitscore> (I<score>)
170 :    
171 :     Set/get the bitscore which is currently used to query the database.
172 :     Parameter I<score> has to be a float or undef.
173 :    
174 :     =cut
175 :    
176 :     sub query_bitscore {
177 :     if(scalar(@_)>1) {
178 :     $_[0]->{query}->{bitscore} = $_[1];
179 :     }
180 :     return $_[0]->{query}->{bitscore};
181 :     }
182 :    
183 :    
184 :     =pod
185 :    
186 :     =item * B<query_align_len> (I<length>)
187 :    
188 :     Set/get the minimum alignment which is currently used to query the database.
189 :     Parameter I<length> has to be a positive integer or undef.
190 :    
191 :     =cut
192 :    
193 :     sub query_align_len {
194 :     if(scalar(@_)>1) {
195 :     if($_[1] and $_[1]<0) {
196 :     die "Alignment length has to be positive: ".$_[1];
197 :     }
198 :     $_[0]->{query}->{align_len} = $_[1];
199 :     }
200 :     return $_[0]->{query}->{align_len};
201 :     }
202 :    
203 :    
204 :     =pod
205 :    
206 :     =item * B<query_identity> (I<percent>)
207 :    
208 :     Set/get the minimum percent identity which is currently used to query the database.
209 :     Parameter I<percent> has to be a number in 0..100 or undef.
210 :    
211 :     =cut
212 :    
213 :     sub query_identity {
214 :     if(scalar(@_)>1) {
215 :     if($_[1] and ($_[1]<0 or $_[1]>100)) {
216 :     die "Identity has to be between 0 and 100: ".$_[1];
217 :     }
218 :     $_[0]->{query}->{identity} = $_[1];
219 :     }
220 :     return $_[0]->{query}->{identity};
221 :     }
222 :    
223 :    
224 :     =pod
225 :    
226 :     =item * B<query_load_from_cgi> (I<cgi>, [I<dataset>])
227 :    
228 :     Sets all query parameter to the values provided in the CGI query object I<cgi>.
229 :     This method recognises 'evalue', 'pvalue' (bitscore), 'alignment_length' and
230 :     'percent_identity' as query criteria. Any missing param will be set to undef.
231 :     If the optional parameter I<dataset> is set to one of the accepted datasets
232 :     (db types), the method will additionally load the defaults for this type into
233 :     the CGI object.
234 :    
235 :    
236 :     =cut
237 :    
238 :     sub query_load_from_cgi {
239 :     my ($self, $cgi, $dataset) = @_;
240 :    
241 :     unless(ref $cgi and $cgi->isa("CGI")) {
242 :     die "Query load from cgi requires a valid CGI object.";
243 :     }
244 :    
245 :     # load the defaults if necessary
246 :     if($dataset and $self->get_dataset_id($dataset)) {
247 :    
248 :     my $d = $self->get_dataset_id($dataset);
249 :    
250 :     my @v = qw( evalue bitscore align_len identity );
251 :     foreach my $v (@v) {
252 :     if(!defined($cgi->param($v)) and QUERY_DEFAULTS->{$d}->{$v}) {
253 :     $cgi->param($v, QUERY_DEFAULTS->{$d}->{$v});
254 :     }
255 :     }
256 :     }
257 :    
258 :     # set the query params
259 :     my $evalue = $cgi->param('evalue') || '';
260 :     $self->query_evalue($evalue);
261 :    
262 :     my $bitscore = $cgi->param('bitscore') || '';
263 :     $self->query_bitscore($bitscore);
264 :    
265 :     my $align_len = $cgi->param('align_len') || '';
266 :     $self->query_align_len($align_len);
267 :    
268 :     my $identity = $cgi->param('identity') || '';
269 :     $self->query_identity($identity);
270 :    
271 :     return $self;
272 :    
273 :     }
274 :    
275 :    
276 :     =pod
277 :    
278 :     =item * B<get_where_clause> ()
279 :    
280 :     Returns for the current query parameters the where clause as applicable to the
281 :     tax_sim_XYZ table SQL queries. The method will take care of all conversions to
282 :     eg the logscore evalues.
283 :    
284 :     =cut
285 :    
286 :     sub get_where_clause {
287 :     my ($self) = @_;
288 :    
289 :     my @params;
290 :    
291 :     if($self->{query}->{evalue}) {
292 :     push @params, "logpsc<=".$self->evalue2log($self->{query}->{evalue});
293 :     }
294 :    
295 :     if($self->{query}->{bitscore}) {
296 :     push @params, "bsc>=".$self->{query}->{bitscore};
297 :     }
298 :    
299 :     if($self->{query}->{align_len}) {
300 :     push @params, "ali_ln>=".$self->{query}->{align_len};
301 :     }
302 :    
303 :     if($self->{query}->{identity}) {
304 :     push @params, "iden>=".$self->{query}->{identity};
305 :     }
306 :    
307 :     return join(' and ', @params);
308 :    
309 :     }
310 :    
311 :    
312 :    
313 :     #******************************************************************************
314 :     #* OTHER
315 :     #******************************************************************************
316 :    
317 :    
318 :    
319 :    
320 :    
321 :     =pod
322 :    
323 : paarmann 1.1 =item * B<get_sequence> (I<sequence_id>)
324 :    
325 :     Retrieve the sequence I<sequence_id> from the metagenome job directory.
326 :    
327 :     =cut
328 :    
329 :     sub get_sequence {
330 :     my ($self, $id) = @_;
331 :    
332 :     my $sequence_file = $self->job->org_dir.'/contigs';
333 :    
334 :     my $sequence = '';
335 :     open(FASTA, "<$sequence_file") or die "Unable to read metagenome sequences: $!";
336 :     while(<FASTA>) {
337 :     next unless /^\>$id/;
338 :    
339 :     while(<FASTA>) {
340 :     last if /^>/;
341 :     chomp;
342 :     $sequence .= $_;
343 :     }
344 :     }
345 :    
346 :     return $sequence;
347 :    
348 :     }
349 :    
350 :    
351 : paarmann 1.3 =pod
352 : paarmann 1.1
353 : paarmann 1.3 =item * B<get_hits_count> (I<dataset_name>)
354 :    
355 :     Given a dataset name (db_id), this method returns
356 : paarmann 1.1 the total number of sequences that contain a hit.
357 :    
358 :     =cut
359 :    
360 :     sub get_hits_count {
361 : paarmann 1.3 my ($self, $dataset) = @_;
362 : paarmann 1.1
363 :     my $table = $self->dbtable;
364 :     my $dbid = $self->get_dataset_id($dataset);
365 : paarmann 1.3 my $where = $self->get_where_clause();
366 :     $where = ($where) ? "and $where" : '';
367 : paarmann 1.1
368 : paarmann 1.3 my $sth = $self->dbh->prepare("select count(*) from ( select id1, min(rank_psc) from $table where dbid=$dbid $where group by id1) as b");
369 : paarmann 1.1 $sth->execute;
370 :     my ($result) = $sth->fetchrow_array;
371 :    
372 :     return $result;
373 :    
374 :     }
375 :    
376 :    
377 : paarmann 1.3 =pod
378 :    
379 :     =item * B<get_group_counts> (I<dataset_name>, [I<group>, I<filter1>, I<filter2>])
380 :    
381 :     Given a dataset name (db_id), this method returns the total counts for all
382 :     taxonomy groups of a certain depth which are hit. If no group name I<group>
383 :     was given, the method returns counts for tax_group_1.
384 :     Optionally, I<group> may be 'tax_group_2' or 'tax_group_3' and in that case
385 :     any optional provided filters I<filter1> and I<filter2> will be applied to
386 :     the column 'tax_group_1' and 'tax_group_2' respectively.
387 :    
388 :     =cut
389 :    
390 :     sub get_group_counts {
391 :     my ($self, $dataset, $group, $filter1, $filter2) = @_;
392 :    
393 :     my $table = $self->dbtable;
394 :     my $dbid = $self->get_dataset_id($dataset);
395 :     my $where = $self->get_where_clause();
396 :     $where = ($where) ? "and $where" : '';
397 :     $group = 'tax_group_1' unless($group);
398 :    
399 :     my @filters;
400 :     push @filters, "tax_group_1='$filter1'" if($filter1);
401 :     push @filters, "tax_group_2='$filter2'" if($filter2);
402 :     my $filter = (scalar(@filters)) ? 'and '.join(' and ', @filters) : '';
403 :    
404 :     my $sth = $self->dbh->prepare("select s.$group as tax, count(*) as num from ( select id1, min(rank_psc) as rank from $table where dbid=$dbid $where group by id1) as b inner join $table as s on b.id1=s.id1 and b.rank=s.rank_psc where dbid=$dbid $filter group by s.$group");
405 :     $sth->execute;
406 :     my $result = $sth->fetchall_arrayref();
407 :    
408 :     return $result;
409 :    
410 :     }
411 :    
412 :    
413 :     =pod
414 :    
415 :     =item * B<get_taxa_counts> (I<dataset_name>)
416 :    
417 :     Given a dataset name (db_id), this method returns the total counts for all
418 :     taxonomy strings which are hit.
419 :    
420 :     =cut
421 :    
422 :     sub get_taxa_counts {
423 :     my ($self, $dataset) = @_;
424 :    
425 :     my $table = $self->dbtable;
426 :     my $dbid = $self->get_dataset_id($dataset);
427 :     my $where = $self->get_where_clause();
428 :     $where = ($where) ? "and $where" : '';
429 :    
430 :     my $sth = $self->dbh->prepare("select s.tax_str as tax, count(*) as num from ( select id1, min(rank_psc) as rank from $table where dbid=$dbid $where group by id1) as b inner join $table as s on b.id1=s.id1 and b.rank=s.rank_psc where dbid=$dbid group by s.tax_str");
431 :    
432 :     $sth->execute;
433 :     my $result = $sth->fetchall_arrayref();
434 :    
435 :     return $result;
436 :    
437 :     }
438 :    
439 :    
440 :     =pod
441 :    
442 :     =item * B<get_subsystem_counts> (I<dataset_name>)
443 :    
444 :     Given a dataset name (db_id), this method returns the total counts for all
445 :     subsystems which are hit.
446 :    
447 :     =cut
448 :    
449 :     sub get_subsystem_counts {
450 :     my ($self, $dataset) = @_;
451 :    
452 :     my $table = $self->dbtable;
453 :     my $dbid = $self->get_dataset_id($dataset);
454 :     my $where = $self->get_where_clause();
455 :     $where = ($where) ? "and $where" : '';
456 :    
457 :     my $sth = $self->dbh->prepare("select s.tax_group_1, s.tax_group_2, s.tax_group_3, s.tax_str, count(*) as num from ( select id1, min(rank_psc) as rank from $table where dbid=$dbid $where group by id1) as b inner join $table as s on b.id1=s.id1 and b.rank=s.rank_psc where dbid=$dbid group by s.tax_group_1, s.tax_group_2, s.tax_group_3");
458 :    
459 :     $sth->execute;
460 :     my $result = $sth->fetchall_arrayref();
461 :    
462 :     return $result;
463 :    
464 :     }
465 :    
466 :    
467 :     =pod
468 :    
469 :     =item * B<get_sequence_subset> (I<dataset_name>, I<filter>)
470 :    
471 :     Given a dataset name (db_id), this method returns all sequence ids,
472 :     the alignment length, the match id and the taxonomy string for all
473 :     sequences which match the criteria and have their tax_str start with
474 :     the filter string I<filter>.
475 :    
476 :     =cut
477 :    
478 :     sub get_sequence_subset {
479 :     my ($self, $dataset, $filter) = @_;
480 :    
481 :     my $table = $self->dbtable;
482 :     my $dbid = $self->get_dataset_id($dataset);
483 :     my $where = $self->get_where_clause();
484 :     $where = ($where) ? "and $where" : '';
485 :    
486 :     my $sth = $self->dbh->prepare("select s.id1, s.ali_ln, s.id2, s.tax_str from ( select id1, min(rank_psc) as rank from $table where dbid=$dbid $where group by id1) as b inner join $table as s on b.id1=s.id1 and b.rank=s.rank_psc where dbid=$dbid and s.tax_str like '$filter%'");
487 :    
488 :     $sth->execute;
489 :     my $result = $sth->fetchall_arrayref();
490 :    
491 :     return $result;
492 :    
493 :     }
494 :    
495 :    
496 :     =pod
497 :    
498 : paarmann 1.1 =item * B<get_align_len_range> (I<dataset_name>)
499 :    
500 :     Given a dataset name (db_id), this method returns
501 :     the minimum and maximum alignment length.
502 :    
503 :     =cut
504 :    
505 :     sub get_align_len_range {
506 :     my ($self, $dataset) = @_;
507 :    
508 :     my $table = $self->dbtable;
509 :     my $dbid = $self->get_dataset_id($dataset);
510 :    
511 : paarmann 1.3 my $sth = $self->dbh->prepare("select min(ali_ln), max(ali_ln) from $table where dbid=$dbid");
512 : paarmann 1.1 $sth->execute;
513 :     my ($min, $max) = $sth->fetchrow_array;
514 :    
515 :     return ($min, $max);
516 :    
517 :     }
518 :    
519 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3