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

Annotation of /SeedViewer/MetagenomeAnalysis.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3