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

Annotation of /SeedViewer/MetagenomeAnalysis.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : paarmann 1.1 package SeedViewer::MetagenomeAnalysis;
2 :    
3 : parrello 1.6 # $Id: MetagenomeAnalysis.pm,v 1.5 2008/04/29 14:08:42 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 :    
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 : parrello 1.6 =over 4
153 :    
154 : paarmann 1.3 =item * B<query_evalue> (I<evalue>)
155 :    
156 :     Set/get the expectation value which is currently used to query the database.
157 :     Parameter I<evalue> has to be a float or in '1e-5'-like format or undef.
158 :    
159 :     =cut
160 :    
161 :     sub query_evalue {
162 :     if(scalar(@_)>1) {
163 :     $_[0]->{query}->{evalue} = $_[1];
164 :     }
165 :     return $_[0]->{query}->{evalue};
166 :     }
167 :    
168 :    
169 :     =pod
170 :    
171 :     =item * B<query_bitscore> (I<score>)
172 :    
173 :     Set/get the bitscore which is currently used to query the database.
174 :     Parameter I<score> has to be a float or undef.
175 :    
176 :     =cut
177 :    
178 :     sub query_bitscore {
179 :     if(scalar(@_)>1) {
180 :     $_[0]->{query}->{bitscore} = $_[1];
181 :     }
182 :     return $_[0]->{query}->{bitscore};
183 :     }
184 :    
185 :    
186 :     =pod
187 :    
188 :     =item * B<query_align_len> (I<length>)
189 :    
190 :     Set/get the minimum alignment which is currently used to query the database.
191 :     Parameter I<length> has to be a positive integer or undef.
192 :    
193 :     =cut
194 :    
195 :     sub query_align_len {
196 :     if(scalar(@_)>1) {
197 :     if($_[1] and $_[1]<0) {
198 :     die "Alignment length has to be positive: ".$_[1];
199 :     }
200 :     $_[0]->{query}->{align_len} = $_[1];
201 :     }
202 :     return $_[0]->{query}->{align_len};
203 :     }
204 :    
205 :    
206 :     =pod
207 :    
208 :     =item * B<query_identity> (I<percent>)
209 :    
210 :     Set/get the minimum percent identity which is currently used to query the database.
211 :     Parameter I<percent> has to be a number in 0..100 or undef.
212 :    
213 :     =cut
214 :    
215 :     sub query_identity {
216 :     if(scalar(@_)>1) {
217 :     if($_[1] and ($_[1]<0 or $_[1]>100)) {
218 :     die "Identity has to be between 0 and 100: ".$_[1];
219 :     }
220 :     $_[0]->{query}->{identity} = $_[1];
221 :     }
222 :     return $_[0]->{query}->{identity};
223 :     }
224 :    
225 :    
226 :     =pod
227 :    
228 :     =item * B<query_load_from_cgi> (I<cgi>, [I<dataset>])
229 :    
230 :     Sets all query parameter to the values provided in the CGI query object I<cgi>.
231 :     This method recognises 'evalue', 'pvalue' (bitscore), 'alignment_length' and
232 :     'percent_identity' as query criteria. Any missing param will be set to undef.
233 :     If the optional parameter I<dataset> is set to one of the accepted datasets
234 :     (db types), the method will additionally load the defaults for this type into
235 :     the CGI object.
236 :    
237 :    
238 :     =cut
239 :    
240 :     sub query_load_from_cgi {
241 :     my ($self, $cgi, $dataset) = @_;
242 :    
243 :     unless(ref $cgi and $cgi->isa("CGI")) {
244 :     die "Query load from cgi requires a valid CGI object.";
245 :     }
246 :    
247 :     # load the defaults if necessary
248 :     if($dataset and $self->get_dataset_id($dataset)) {
249 :    
250 :     my $d = $self->get_dataset_id($dataset);
251 :    
252 :     my @v = qw( evalue bitscore align_len identity );
253 :     foreach my $v (@v) {
254 :     if(!defined($cgi->param($v)) and QUERY_DEFAULTS->{$d}->{$v}) {
255 :     $cgi->param($v, QUERY_DEFAULTS->{$d}->{$v});
256 :     }
257 :     }
258 :     }
259 :    
260 :     # set the query params
261 :     my $evalue = $cgi->param('evalue') || '';
262 :     $self->query_evalue($evalue);
263 :    
264 :     my $bitscore = $cgi->param('bitscore') || '';
265 :     $self->query_bitscore($bitscore);
266 :    
267 :     my $align_len = $cgi->param('align_len') || '';
268 :     $self->query_align_len($align_len);
269 :    
270 :     my $identity = $cgi->param('identity') || '';
271 :     $self->query_identity($identity);
272 :    
273 :     return $self;
274 :    
275 :     }
276 :    
277 :    
278 :     =pod
279 :    
280 :     =item * B<get_where_clause> ()
281 :    
282 :     Returns for the current query parameters the where clause as applicable to the
283 :     tax_sim_XYZ table SQL queries. The method will take care of all conversions to
284 :     eg the logscore evalues.
285 :    
286 :     =cut
287 :    
288 :     sub get_where_clause {
289 :     my ($self) = @_;
290 :    
291 :     my @params;
292 :    
293 :     if($self->{query}->{evalue}) {
294 :     push @params, "logpsc<=".$self->evalue2log($self->{query}->{evalue});
295 :     }
296 :    
297 :     if($self->{query}->{bitscore}) {
298 :     push @params, "bsc>=".$self->{query}->{bitscore};
299 :     }
300 :    
301 :     if($self->{query}->{align_len}) {
302 :     push @params, "ali_ln>=".$self->{query}->{align_len};
303 :     }
304 :    
305 :     if($self->{query}->{identity}) {
306 :     push @params, "iden>=".$self->{query}->{identity};
307 :     }
308 :    
309 :     return join(' and ', @params);
310 :    
311 :     }
312 :    
313 :    
314 :    
315 :     #******************************************************************************
316 :     #* OTHER
317 :     #******************************************************************************
318 :    
319 :    
320 :    
321 :    
322 :    
323 :     =pod
324 :    
325 : paarmann 1.1 =item * B<get_sequence> (I<sequence_id>)
326 :    
327 :     Retrieve the sequence I<sequence_id> from the metagenome job directory.
328 :    
329 :     =cut
330 :    
331 :     sub get_sequence {
332 :     my ($self, $id) = @_;
333 :    
334 :     my $sequence_file = $self->job->org_dir.'/contigs';
335 :    
336 :     my $sequence = '';
337 :     open(FASTA, "<$sequence_file") or die "Unable to read metagenome sequences: $!";
338 :     while(<FASTA>) {
339 :     next unless /^\>$id/;
340 :    
341 :     while(<FASTA>) {
342 :     last if /^>/;
343 :     chomp;
344 :     $sequence .= $_;
345 :     }
346 :     }
347 :    
348 :     return $sequence;
349 :    
350 :     }
351 :    
352 :    
353 : paarmann 1.3 =pod
354 : paarmann 1.1
355 : paarmann 1.3 =item * B<get_hits_count> (I<dataset_name>)
356 :    
357 :     Given a dataset name (db_id), this method returns
358 : paarmann 1.1 the total number of sequences that contain a hit.
359 :    
360 :     =cut
361 :    
362 :     sub get_hits_count {
363 : paarmann 1.3 my ($self, $dataset) = @_;
364 : paarmann 1.1
365 :     my $table = $self->dbtable;
366 :     my $dbid = $self->get_dataset_id($dataset);
367 : paarmann 1.3 my $where = $self->get_where_clause();
368 :     $where = ($where) ? "and $where" : '';
369 : paarmann 1.1
370 : 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");
371 : paarmann 1.1 $sth->execute;
372 :     my ($result) = $sth->fetchrow_array;
373 :    
374 :     return $result;
375 :    
376 :     }
377 :    
378 :    
379 : paarmann 1.3 =pod
380 :    
381 :     =item * B<get_group_counts> (I<dataset_name>, [I<group>, I<filter1>, I<filter2>])
382 :    
383 :     Given a dataset name (db_id), this method returns the total counts for all
384 :     taxonomy groups of a certain depth which are hit. If no group name I<group>
385 :     was given, the method returns counts for tax_group_1.
386 :     Optionally, I<group> may be 'tax_group_2' or 'tax_group_3' and in that case
387 :     any optional provided filters I<filter1> and I<filter2> will be applied to
388 :     the column 'tax_group_1' and 'tax_group_2' respectively.
389 :    
390 :     =cut
391 :    
392 :     sub get_group_counts {
393 :     my ($self, $dataset, $group, $filter1, $filter2) = @_;
394 :    
395 :     my $table = $self->dbtable;
396 :     my $dbid = $self->get_dataset_id($dataset);
397 :     my $where = $self->get_where_clause();
398 :     $where = ($where) ? "and $where" : '';
399 :     $group = 'tax_group_1' unless($group);
400 :    
401 :     my @filters;
402 :     push @filters, "tax_group_1='$filter1'" if($filter1);
403 :     push @filters, "tax_group_2='$filter2'" if($filter2);
404 :     my $filter = (scalar(@filters)) ? 'and '.join(' and ', @filters) : '';
405 :    
406 :     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");
407 :     $sth->execute;
408 :     my $result = $sth->fetchall_arrayref();
409 :    
410 :     return $result;
411 :    
412 :     }
413 :    
414 :    
415 :     =pod
416 :    
417 :     =item * B<get_taxa_counts> (I<dataset_name>)
418 :    
419 :     Given a dataset name (db_id), this method returns the total counts for all
420 :     taxonomy strings which are hit.
421 :    
422 :     =cut
423 :    
424 :     sub get_taxa_counts {
425 :     my ($self, $dataset) = @_;
426 :    
427 :     my $table = $self->dbtable;
428 :     my $dbid = $self->get_dataset_id($dataset);
429 :     my $where = $self->get_where_clause();
430 :     $where = ($where) ? "and $where" : '';
431 :    
432 :     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");
433 :    
434 :     $sth->execute;
435 :     my $result = $sth->fetchall_arrayref();
436 :    
437 :     return $result;
438 :    
439 :     }
440 :    
441 :    
442 :     =pod
443 :    
444 :     =item * B<get_subsystem_counts> (I<dataset_name>)
445 :    
446 :     Given a dataset name (db_id), this method returns the total counts for all
447 :     subsystems which are hit.
448 :    
449 :     =cut
450 :    
451 :     sub get_subsystem_counts {
452 :     my ($self, $dataset) = @_;
453 :    
454 :     my $table = $self->dbtable;
455 :     my $dbid = $self->get_dataset_id($dataset);
456 :     my $where = $self->get_where_clause();
457 :     $where = ($where) ? "and $where" : '';
458 : paarmann 1.4
459 :     #print STDERR "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";
460 :     #die;
461 : paarmann 1.3 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");
462 :    
463 :     $sth->execute;
464 :     my $result = $sth->fetchall_arrayref();
465 :    
466 :     return $result;
467 :    
468 :     }
469 :    
470 :    
471 :     =pod
472 :    
473 :     =item * B<get_sequence_subset> (I<dataset_name>, I<filter>)
474 :    
475 :     Given a dataset name (db_id), this method returns all sequence ids,
476 :     the alignment length, the match id and the taxonomy string for all
477 :     sequences which match the criteria and have their tax_str start with
478 :     the filter string I<filter>.
479 :    
480 :     =cut
481 :    
482 :     sub get_sequence_subset {
483 :     my ($self, $dataset, $filter) = @_;
484 :    
485 :     my $table = $self->dbtable;
486 :     my $dbid = $self->get_dataset_id($dataset);
487 :     my $where = $self->get_where_clause();
488 :     $where = ($where) ? "and $where" : '';
489 :    
490 :     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%'");
491 :    
492 :     $sth->execute;
493 :     my $result = $sth->fetchall_arrayref();
494 :    
495 :     return $result;
496 :    
497 :     }
498 :    
499 :    
500 : paarmann 1.4 =pod
501 :    
502 : jared 1.5 =item * B<get_recruitment_plot_data> (I<genome>)
503 :    
504 :     Given a genome id (83333.1), this method returns all sequence ids,
505 :     the alignment length, the match id and the taxonomy string for all
506 :     sequences which match the criteria and have their tax_str start equal
507 :     the genome tax string I<filter>.
508 :    
509 :     =cut
510 :    
511 :     sub get_recruitment_plot_data {
512 :     my ($self, $genome) = @_;
513 :    
514 :     my $table = $self->dbtable;
515 :     my $dbid = $self->get_dataset_id("SEED");
516 :     my $where = $self->get_where_clause();
517 :     $where = ($where) ? "and $where" : '';
518 :    
519 :     my ($tax_id) = $self->dbh->selectrow_array("select tax_str from rdp_to_tax where seq_num='". $genome . "'");
520 :    
521 :     my $sth = $self->dbh->prepare("select s.id1, s.id2, s.b2, s.e2, s.logpsc 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='$tax_id'");
522 :    
523 :     $sth->execute;
524 :     my $result = $sth->fetchall_arrayref();
525 :    
526 :     return $result;
527 :    
528 :     }
529 :    
530 :    
531 :    
532 :    
533 :     =pod
534 :    
535 : paarmann 1.4 =item * B<get_hits_for_sequence> (I<seq_id>, I<dataset>, I<limit>)
536 :    
537 :     Given a sequence id I<seq_id> (id1) and a dataset name (db_id), this method returns
538 :     the first I<limit> rows of hit data for this sequence. If no I<limit> is provided, it
539 :     will default to 10.
540 :     It returns (match id, taxonomy string, log evalue, bitscore, alignment length,
541 :     percent identity, start1, end1) per hit.
542 :    
543 :     =cut
544 :    
545 :     sub get_hits_for_sequence {
546 :     my ($self, $id, $dataset, $limit) = @_;
547 :    
548 :     my $table = $self->dbtable;
549 :     my $dbid = $self->get_dataset_id($dataset);
550 :     $limit = 10 unless ($limit);
551 :    
552 :     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;");
553 :     $sth->execute;
554 :     my $result = $sth->fetchall_arrayref();
555 :    
556 :     return $result;
557 :    
558 :     }
559 :    
560 :    
561 : paarmann 1.3 =pod
562 :    
563 : paarmann 1.1 =item * B<get_align_len_range> (I<dataset_name>)
564 :    
565 :     Given a dataset name (db_id), this method returns
566 :     the minimum and maximum alignment length.
567 :    
568 :     =cut
569 :    
570 :     sub get_align_len_range {
571 :     my ($self, $dataset) = @_;
572 :    
573 :     my $table = $self->dbtable;
574 :     my $dbid = $self->get_dataset_id($dataset);
575 :    
576 : paarmann 1.3 my $sth = $self->dbh->prepare("select min(ali_ln), max(ali_ln) from $table where dbid=$dbid");
577 : paarmann 1.1 $sth->execute;
578 :     my ($min, $max) = $sth->fetchrow_array;
579 :    
580 :     return ($min, $max);
581 :    
582 :     }
583 :    
584 : parrello 1.6 =pod
585 :    
586 :     =back
587 : paarmann 1.1
588 : parrello 1.6 =cut

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3