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

Annotation of /SeedViewer/MetagenomeAnalysis.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : paarmann 1.1 package SeedViewer::MetagenomeAnalysis;
2 :    
3 : jared 1.21 # $Id: MetagenomeAnalysis.pm,v 1.20 2008/08/15 22:46:48 dsouza Exp $
4 : paarmann 1.1
5 :     use strict;
6 :     use warnings;
7 :    
8 :     use FIG_Config;
9 :     use DBI;
10 : olson 1.13 use Data::Dumper;
11 : paarmann 1.1
12 :     1;
13 :    
14 : paarmann 1.3 use constant QUERY_DEFAULTS =>
15 :     { 1 => { evalue => '1e-05', align_len => 50 }, # RDP
16 :     2 => { evalue => '0.01' }, # SEED
17 :     3 => { evalue => '1e-05', align_len => 50 }, # Greengenes
18 :     4 => { evalue => '1e-05', align_len => 50 }, # LSU
19 :     5 => { evalue => '1e-05', align_len => 50 }, # SSU
20 :     6 => { evalue => '0.01' }, # Subsystem
21 :     };
22 :    
23 : paarmann 1.1 sub new {
24 :     my ($class, $job) = @_;
25 :    
26 :     # check job
27 :     unless (ref $job and $job->isa('RAST::Job') and $job->metagenome) {
28 :     return undef;
29 :     }
30 :    
31 :     # connect to database
32 :     my $dbh;
33 :     eval {
34 :    
35 : olson 1.7 my $dbms = $FIG_Config::mgrast_dbms;
36 : paarmann 1.1 my $host = $FIG_Config::mgrast_dbhost;
37 :     my $database = $FIG_Config::mgrast_db;
38 :     my $user = $FIG_Config::mgrast_dbuser;
39 :     my $password = '';
40 : olson 1.7
41 :     if ($dbms eq 'Pg')
42 :     {
43 : jared 1.8 $dbh = DBI->connect("DBI:Pg:dbname=$database;host=$host", $user, $password,
44 : olson 1.7 { RaiseError => 1, AutoCommit => 0, PrintError => 0 }) ||
45 :     die "database connect error.";
46 :     }
47 :     elsif ($dbms eq 'mysql' or $dbms eq '') # Default to mysql
48 :     {
49 :     $dbh = DBI->connect("DBI:mysql:database=$database;host=$host", $user, $password,
50 : paarmann 1.1 { RaiseError => 1, AutoCommit => 0, PrintError => 0 }) ||
51 : olson 1.7 die "database connect error.";
52 :     }
53 :     else
54 :     {
55 :     die "MetagenomeAnalysis: unknown dbms '$dbms'";
56 :     }
57 :    
58 : paarmann 1.1 };
59 :     if ($@) {
60 :     warn "Unable to connect to metagenomics database: $@\n";
61 :     return undef;
62 :     }
63 :    
64 :     # create object
65 :     my $self = { job => $job,
66 :     dbh => $dbh,
67 :     key2taxa => undef,
68 : paarmann 1.3 query => {},
69 : arodri7 1.22 dbid_cache => {},
70 : paarmann 1.1 };
71 :     bless $self, $class;
72 :    
73 :     # load key2taxa mapping
74 :     $self->get_key2taxa_mapping();
75 :    
76 :     return $self;
77 :    
78 :     }
79 :    
80 :     sub job {
81 :     return $_[0]->{job};
82 :     }
83 :    
84 :    
85 :     sub dbh {
86 :     return $_[0]->{dbh};
87 :     }
88 :    
89 :    
90 :     sub dbtable {
91 :     unless (defined $_[0]->{dbtable}) {
92 :     $_[0]->{dbtable} = 'tax_sim_'.$_[0]->job->id;
93 :     }
94 :     return $_[0]->{dbtable};
95 :     }
96 :    
97 : jared 1.9 sub dbtable_best_psc {
98 :     unless (defined $_[0]->{dbtable}) {
99 :     $_[0]->{dbtable} = 'tax_sim_best_by_psc_'.$_[0]->job->id;
100 :     }
101 :     return $_[0]->{dbtable};
102 :     }
103 :    
104 :     sub dbtable_best_iden {
105 :     unless (defined $_[0]->{dbtable}) {
106 :     $_[0]->{dbtable} = 'tax_sim_best_by_iden_'.$_[0]->job->id;
107 :     }
108 :     return $_[0]->{dbtable};
109 :     }
110 : paarmann 1.1
111 :     sub get_key2taxa_mapping {
112 :    
113 :     unless (defined($_[0]->{key2taxa})) {
114 :     my $sth = $_[0]->dbh->prepare("select dbkey, str from tax_item");
115 :     $sth->execute;
116 :     $_[0]->{key2taxa} = $sth->fetchall_hashref('dbkey');
117 :     }
118 :    
119 :     return $_[0]->{key2taxa};
120 :    
121 :     }
122 :    
123 : arodri7 1.22
124 : paarmann 1.1 sub key2taxa {
125 :     if(defined $_[1]) {
126 :     my $t = $_[0]->{key2taxa}->{$_[1]}->{str};
127 :     $t =~ s/\t+$//;
128 :     return $t;
129 :     }
130 :     return '';
131 :     }
132 :    
133 :    
134 :     sub split_taxstr {
135 :     my @r = split(':', $_[1]);
136 :     return \@r;
137 :     }
138 :    
139 :     sub join_taxstr {
140 :     # do I really want an ending colon?
141 :     return join(':', @{$_[1]}).':';
142 :     }
143 :    
144 :    
145 :     sub evalue2log {
146 :     return 10 * (log($_[1]) / log(10));
147 :     }
148 :    
149 :     sub log2evalue {
150 :     return 10**($_[1]/10);
151 :     }
152 :    
153 :    
154 : olson 1.13 #
155 :     # Determine the correct dbid for this job. Use sims.database_list
156 :     # to find the version that the analysis was run with.
157 :     #
158 : paarmann 1.1 sub get_dataset_id {
159 : olson 1.13 my($self, $dataset) = @_;
160 :    
161 :     my $id = $self->{dbid_cache}->{$dataset};
162 :     return $id if defined($id);
163 :    
164 :     my($dbname, $type) = split(/:/, $dataset);
165 :    
166 :     my $dbs = $self->job->metaxml->get_metadata('sims.database_list');
167 :    
168 :     my @this = grep { $_->{name} eq $dbname } @$dbs;
169 :     if (@this)
170 :     {
171 :     my $vers = $this[0]->{version};
172 :    
173 :     #
174 :     # Now we can find the dbid ni the database.
175 :     #
176 :     my $res = $self->dbh->selectcol_arrayref(qq(SELECT dbid
177 :     FROM seq_db
178 :     WHERE name = ? AND version = ?
179 :     AND tax_db_name = ?), undef,
180 :     $dbname, $vers, $type);
181 :     if (@$res)
182 :     {
183 : paczian 1.14 #print STDERR "Found @$res for $dbname $type $vers\n";
184 : olson 1.13 $id = $res->[0];
185 :     $self->{dbid_cache}->{$dataset} = $id;
186 :     return $id;
187 :     }
188 : dsouza 1.20 #print STDERR "Did not find anything for dataset='$dataset' '$dbname' '$type' '$vers'\n";
189 : olson 1.13 }
190 : dsouza 1.20 #print STDERR "did not find a vers for dataset='$dataset' $dbname $type\n" . Dumper($dbs);
191 : paarmann 1.1 }
192 :    
193 : paarmann 1.3 #******************************************************************************
194 :     #* MANAGING QUERY CRITERIA
195 :     #******************************************************************************
196 :    
197 :     =pod
198 :    
199 : parrello 1.6 =over 4
200 :    
201 : paarmann 1.3 =item * B<query_evalue> (I<evalue>)
202 :    
203 :     Set/get the expectation value which is currently used to query the database.
204 :     Parameter I<evalue> has to be a float or in '1e-5'-like format or undef.
205 :    
206 :     =cut
207 :    
208 :     sub query_evalue {
209 :     if(scalar(@_)>1) {
210 :     $_[0]->{query}->{evalue} = $_[1];
211 :     }
212 :     return $_[0]->{query}->{evalue};
213 :     }
214 :    
215 :    
216 :     =pod
217 :    
218 :     =item * B<query_bitscore> (I<score>)
219 :    
220 :     Set/get the bitscore which is currently used to query the database.
221 :     Parameter I<score> has to be a float or undef.
222 :    
223 :     =cut
224 :    
225 :     sub query_bitscore {
226 :     if(scalar(@_)>1) {
227 :     $_[0]->{query}->{bitscore} = $_[1];
228 :     }
229 :     return $_[0]->{query}->{bitscore};
230 :     }
231 :    
232 :    
233 :     =pod
234 :    
235 :     =item * B<query_align_len> (I<length>)
236 :    
237 :     Set/get the minimum alignment which is currently used to query the database.
238 :     Parameter I<length> has to be a positive integer or undef.
239 :    
240 :     =cut
241 :    
242 :     sub query_align_len {
243 :     if(scalar(@_)>1) {
244 :     if($_[1] and $_[1]<0) {
245 :     die "Alignment length has to be positive: ".$_[1];
246 :     }
247 :     $_[0]->{query}->{align_len} = $_[1];
248 :     }
249 :     return $_[0]->{query}->{align_len};
250 :     }
251 :    
252 :    
253 :     =pod
254 :    
255 :     =item * B<query_identity> (I<percent>)
256 :    
257 :     Set/get the minimum percent identity which is currently used to query the database.
258 :     Parameter I<percent> has to be a number in 0..100 or undef.
259 :    
260 :     =cut
261 :    
262 :     sub query_identity {
263 :     if(scalar(@_)>1) {
264 :     if($_[1] and ($_[1]<0 or $_[1]>100)) {
265 :     die "Identity has to be between 0 and 100: ".$_[1];
266 :     }
267 :     $_[0]->{query}->{identity} = $_[1];
268 :     }
269 :     return $_[0]->{query}->{identity};
270 :     }
271 :    
272 :    
273 :     =pod
274 :    
275 :     =item * B<query_load_from_cgi> (I<cgi>, [I<dataset>])
276 :    
277 :     Sets all query parameter to the values provided in the CGI query object I<cgi>.
278 :     This method recognises 'evalue', 'pvalue' (bitscore), 'alignment_length' and
279 :     'percent_identity' as query criteria. Any missing param will be set to undef.
280 :     If the optional parameter I<dataset> is set to one of the accepted datasets
281 :     (db types), the method will additionally load the defaults for this type into
282 :     the CGI object.
283 :    
284 :    
285 :     =cut
286 :    
287 :     sub query_load_from_cgi {
288 :     my ($self, $cgi, $dataset) = @_;
289 :    
290 :     unless(ref $cgi and $cgi->isa("CGI")) {
291 :     die "Query load from cgi requires a valid CGI object.";
292 :     }
293 :    
294 :     # load the defaults if necessary
295 :     if($dataset and $self->get_dataset_id($dataset)) {
296 :    
297 :     my $d = $self->get_dataset_id($dataset);
298 :    
299 :     my @v = qw( evalue bitscore align_len identity );
300 :     foreach my $v (@v) {
301 :     if(!defined($cgi->param($v)) and QUERY_DEFAULTS->{$d}->{$v}) {
302 :     $cgi->param($v, QUERY_DEFAULTS->{$d}->{$v});
303 :     }
304 :     }
305 :     }
306 :    
307 :     # set the query params
308 :     my $evalue = $cgi->param('evalue') || '';
309 :     $self->query_evalue($evalue);
310 :    
311 :     my $bitscore = $cgi->param('bitscore') || '';
312 :     $self->query_bitscore($bitscore);
313 :    
314 :     my $align_len = $cgi->param('align_len') || '';
315 :     $self->query_align_len($align_len);
316 :    
317 :     my $identity = $cgi->param('identity') || '';
318 :     $self->query_identity($identity);
319 :    
320 :     return $self;
321 :    
322 :     }
323 :    
324 :    
325 :     =pod
326 :    
327 :     =item * B<get_where_clause> ()
328 :    
329 :     Returns for the current query parameters the where clause as applicable to the
330 :     tax_sim_XYZ table SQL queries. The method will take care of all conversions to
331 :     eg the logscore evalues.
332 :    
333 :     =cut
334 :    
335 :     sub get_where_clause {
336 :     my ($self) = @_;
337 :    
338 :     my @params;
339 :    
340 :     if($self->{query}->{evalue}) {
341 :     push @params, "logpsc<=".$self->evalue2log($self->{query}->{evalue});
342 :     }
343 :    
344 :     if($self->{query}->{bitscore}) {
345 :     push @params, "bsc>=".$self->{query}->{bitscore};
346 :     }
347 :    
348 :     if($self->{query}->{align_len}) {
349 :     push @params, "ali_ln>=".$self->{query}->{align_len};
350 :     }
351 :    
352 :     if($self->{query}->{identity}) {
353 :     push @params, "iden>=".$self->{query}->{identity};
354 :     }
355 :    
356 :     return join(' and ', @params);
357 :    
358 :     }
359 :    
360 :    
361 :    
362 :     #******************************************************************************
363 :     #* OTHER
364 :     #******************************************************************************
365 :    
366 :    
367 :     =pod
368 :    
369 : paarmann 1.1 =item * B<get_sequence> (I<sequence_id>)
370 :    
371 :     Retrieve the sequence I<sequence_id> from the metagenome job directory.
372 :    
373 :     =cut
374 :    
375 :     sub get_sequence {
376 :     my ($self, $id) = @_;
377 :    
378 :     my $sequence_file = $self->job->org_dir.'/contigs';
379 :    
380 :     my $sequence = '';
381 :     open(FASTA, "<$sequence_file") or die "Unable to read metagenome sequences: $!";
382 :     while(<FASTA>) {
383 :     next unless /^\>$id/;
384 :    
385 :     while(<FASTA>) {
386 :     last if /^>/;
387 :     chomp;
388 :     $sequence .= $_;
389 :     }
390 :     }
391 :    
392 :     return $sequence;
393 :    
394 :     }
395 :    
396 : paczian 1.15 =pod
397 :    
398 :     =item * B<get_sequences_fasta> (I<sequence_ids>)
399 :    
400 :     Retrieve the sequences for a given list of I<sequence_ids> from the metagenome job directory in fasta format.
401 :    
402 :     =cut
403 :    
404 :     sub get_sequences_fasta {
405 :     my ($self, $ids) = @_;
406 :    
407 :     # get the path to the sequence file
408 :     my $sequence_file = $self->job->org_dir.'/contigs';
409 :    
410 :     # hash the ids
411 :     my %idh = map { $_ => 1 } @$ids;
412 :    
413 :     # store the result
414 :     my $sequence = '';
415 :    
416 :     # iterate over the sequence file
417 :     open(FASTA, "<$sequence_file") or die "Unable to read metagenome sequences: $!";
418 :     while(<FASTA>) {
419 :    
420 :     # look only at id lines
421 :     next unless /^\>/;
422 :    
423 :     # get the id
424 :     my ($id) = /^\>(.+)/;
425 :     chomp($id);
426 :    
427 :     # find the id in the requested id hash
428 :     next unless exists ($idh{$id});
429 :    
430 :     # id was found, delete it from the hash
431 :     delete $idh{$id};
432 :    
433 :     # print the header line
434 :     $sequence .= $_;
435 :    
436 :     # print all the following lines to the result
437 :     while(<FASTA>) {
438 :     last if /^>/;
439 :     $sequence .= $_;
440 :     }
441 :    
442 :     last unless (scalar(keys(%idh)));
443 :     }
444 :    
445 :     return $sequence;
446 :    
447 :     }
448 : paarmann 1.1
449 : paarmann 1.3 =pod
450 : paarmann 1.1
451 : paarmann 1.3 =item * B<get_hits_count> (I<dataset_name>)
452 :    
453 :     Given a dataset name (db_id), this method returns
454 : paarmann 1.1 the total number of sequences that contain a hit.
455 :    
456 :     =cut
457 :    
458 :     sub get_hits_count {
459 : paarmann 1.3 my ($self, $dataset) = @_;
460 : jared 1.10
461 :     my $table = $self->dbtable_best_psc;
462 : paarmann 1.1 my $dbid = $self->get_dataset_id($dataset);
463 : paarmann 1.3 my $where = $self->get_where_clause();
464 :     $where = ($where) ? "and $where" : '';
465 : paarmann 1.1
466 : jared 1.10 my $sth = $self->dbh->prepare("select count(distinct id1) from $table where dbid=$dbid $where");
467 : paarmann 1.1 $sth->execute;
468 :     my ($result) = $sth->fetchrow_array;
469 :    
470 :     return $result;
471 :    
472 :     }
473 :    
474 :    
475 : paarmann 1.3 =pod
476 :    
477 :     =item * B<get_group_counts> (I<dataset_name>, [I<group>, I<filter1>, I<filter2>])
478 :    
479 :     Given a dataset name (db_id), this method returns the total counts for all
480 :     taxonomy groups of a certain depth which are hit. If no group name I<group>
481 :     was given, the method returns counts for tax_group_1.
482 :     Optionally, I<group> may be 'tax_group_2' or 'tax_group_3' and in that case
483 :     any optional provided filters I<filter1> and I<filter2> will be applied to
484 :     the column 'tax_group_1' and 'tax_group_2' respectively.
485 :    
486 :     =cut
487 :    
488 :     sub get_group_counts {
489 :     my ($self, $dataset, $group, $filter1, $filter2) = @_;
490 :    
491 : jared 1.10 my $table = $self->dbtable_best_psc;
492 : paarmann 1.3 my $dbid = $self->get_dataset_id($dataset);
493 :     my $where = $self->get_where_clause();
494 :     $where = ($where) ? "and $where" : '';
495 :     $group = 'tax_group_1' unless($group);
496 :    
497 :     my @filters;
498 :     push @filters, "tax_group_1='$filter1'" if($filter1);
499 :     push @filters, "tax_group_2='$filter2'" if($filter2);
500 :     my $filter = (scalar(@filters)) ? 'and '.join(' and ', @filters) : '';
501 :    
502 : arodri7 1.22 print STDERR "select $group as tax, count(*) as num from $table where dbid=$dbid $where $filter group by tax";
503 : 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");
504 : paarmann 1.3 $sth->execute;
505 :     my $result = $sth->fetchall_arrayref();
506 :    
507 : dsouza 1.20 #print STDERR "get_group_counts: ds=$dataset group=$group filter1=$filter1 filter2=$filter2\n";
508 :     #print STDERR Dumper($result);
509 : paarmann 1.3 return $result;
510 :    
511 :     }
512 :    
513 :    
514 :     =pod
515 :    
516 :     =item * B<get_taxa_counts> (I<dataset_name>)
517 :    
518 :     Given a dataset name (db_id), this method returns the total counts for all
519 :     taxonomy strings which are hit.
520 :    
521 :     =cut
522 :    
523 :     sub get_taxa_counts {
524 :     my ($self, $dataset) = @_;
525 :    
526 : jared 1.9 my $table = $self->dbtable_best_psc;
527 : paarmann 1.3 my $dbid = $self->get_dataset_id($dataset);
528 :     my $where = $self->get_where_clause();
529 :     $where = ($where) ? "and $where" : '';
530 :    
531 : paczian 1.16 my $sth = $self->dbh->prepare("select tax_str as tax, count(*) from $table where dbid=$dbid $where group by tax");
532 : paarmann 1.3
533 :     $sth->execute;
534 :     my $result = $sth->fetchall_arrayref();
535 :    
536 :     return $result;
537 :    
538 :     }
539 :    
540 : jared 1.21 sub get_tax2genomeid {
541 :     my ($self, $taxstring) = @_;
542 :    
543 :     my ($genomeid) = $self->dbh->selectrow_array(qq(SELECT seq_num
544 : arodri7 1.22 FROM rdp_to_tax
545 :     WHERE tax_str= ?), undef, $taxstring);
546 :    
547 : jared 1.21 return $genomeid;
548 :     }
549 :    
550 :    
551 : paarmann 1.3 =pod
552 :    
553 :     =item * B<get_subsystem_counts> (I<dataset_name>)
554 :    
555 :     Given a dataset name (db_id), this method returns the total counts for all
556 :     subsystems which are hit.
557 :    
558 :     =cut
559 :    
560 :     sub get_subsystem_counts {
561 :     my ($self, $dataset) = @_;
562 :    
563 : jared 1.9 my $table = $self->dbtable_best_psc;
564 : paarmann 1.3 my $dbid = $self->get_dataset_id($dataset);
565 :     my $where = $self->get_where_clause();
566 :     $where = ($where) ? "and $where" : '';
567 : paarmann 1.4
568 : paczian 1.17 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 $where group by tax_group_1, tax_group_2, tax_group_3, tax_str");
569 : paarmann 1.3
570 :     $sth->execute;
571 :     my $result = $sth->fetchall_arrayref();
572 :    
573 :     return $result;
574 :    
575 :     }
576 :    
577 :    
578 :     =pod
579 :    
580 :     =item * B<get_sequence_subset> (I<dataset_name>, I<filter>)
581 :    
582 :     Given a dataset name (db_id), this method returns all sequence ids,
583 :     the alignment length, the match id and the taxonomy string for all
584 :     sequences which match the criteria and have their tax_str start with
585 :     the filter string I<filter>.
586 :    
587 :     =cut
588 :    
589 :     sub get_sequence_subset {
590 :     my ($self, $dataset, $filter) = @_;
591 :    
592 : jared 1.10 my $table = $self->dbtable_best_psc;
593 : paarmann 1.3 my $dbid = $self->get_dataset_id($dataset);
594 :     my $where = $self->get_where_clause();
595 :     $where = ($where) ? "and $where" : '';
596 :    
597 : paczian 1.15 $filter =~ s/'/''/g;
598 :    
599 : arodri7 1.22 my $sth = $self->dbh->prepare("select id1, ali_ln, id2, tax_str, logpsc, bsc, iden, b1, e1, b2, e2 from $table where dbid=$dbid $where and tax_str like '$filter%'");
600 : paarmann 1.3 $sth->execute;
601 :     my $result = $sth->fetchall_arrayref();
602 :    
603 :     return $result;
604 :    
605 :     }
606 :    
607 : jared 1.19 =pod
608 :    
609 :     =item * B<get_sequence_subset_genome> (I<genome>)
610 :    
611 :     Given a dataset name (db_id), this method returns all sequence ids,
612 :     the alignment length, the match id and the taxonomy string for all
613 :     sequences which match the criteria and have their tax_str start with
614 :     the filter string I<filter>.
615 :    
616 :     =cut
617 :    
618 :     sub get_sequence_subset_genome {
619 :     my ($self, $genome) = @_;
620 :    
621 :     my $table = $self->dbtable_best_psc;
622 :     my $dbid = $self->get_dataset_id("SEED:seed_genome_tax");
623 :     my $where = $self->get_where_clause();
624 :     $where = ($where) ? "and $where" : '';
625 :    
626 :     my ($tax_id) = $self->dbh->selectrow_array(qq(SELECT tax_str
627 :     FROM rdp_to_tax
628 :     WHERE seq_num= ?), undef, $genome);
629 :    
630 :     if($tax_id =~ /(\S+)\s/){
631 :     $tax_id = $1;
632 :     }
633 :    
634 : arodri7 1.22 my $sth = $self->dbh->prepare(qq(SELECT id1, ali_ln, id2, tax_str, logpsc, bsc, iden, b1, e1, b2, e2
635 : jared 1.19 FROM $table
636 :     WHERE dbid=? $where and tax_str=?));
637 :    
638 :    
639 :     $sth->execute($dbid, $tax_id);
640 :     my $result = $sth->fetchall_arrayref();
641 :    
642 :     return $result;
643 :    
644 :     }
645 : paarmann 1.3
646 : paarmann 1.4 =pod
647 :    
648 : jared 1.5 =item * B<get_recruitment_plot_data> (I<genome>)
649 :    
650 :     Given a genome id (83333.1), this method returns all sequence ids,
651 :     the alignment length, the match id and the taxonomy string for all
652 :     sequences which match the criteria and have their tax_str start equal
653 :     the genome tax string I<filter>.
654 :    
655 :     =cut
656 :    
657 :     sub get_recruitment_plot_data {
658 :     my ($self, $genome) = @_;
659 :    
660 : jared 1.9 my $table = $self->dbtable_best_psc;
661 : olson 1.13 my $dbid = $self->get_dataset_id("SEED:seed_genome_tax");
662 : jared 1.5 my $where = $self->get_where_clause();
663 :     $where = ($where) ? "and $where" : '';
664 :    
665 : olson 1.13 my ($tax_id) = $self->dbh->selectrow_array(qq(SELECT tax_str
666 :     FROM rdp_to_tax
667 : arodri7 1.22 WHERE seq_num= ?), undef, $genome);
668 : jared 1.5
669 : jared 1.9 if($tax_id =~ /(\S+)\s/){
670 :     $tax_id = $1;
671 :     }
672 : jared 1.5
673 : olson 1.13 my $sth = $self->dbh->prepare(qq(SELECT id1, id2, b2, e2, logpsc
674 :     FROM $table
675 :     WHERE dbid=? $where and tax_str=?));
676 :    
677 :    
678 :     $sth->execute($dbid, $tax_id);
679 : jared 1.5 my $result = $sth->fetchall_arrayref();
680 :    
681 :     return $result;
682 :    
683 :     }
684 :    
685 :    
686 :    
687 :    
688 :     =pod
689 :    
690 : paarmann 1.4 =item * B<get_hits_for_sequence> (I<seq_id>, I<dataset>, I<limit>)
691 :    
692 :     Given a sequence id I<seq_id> (id1) and a dataset name (db_id), this method returns
693 :     the first I<limit> rows of hit data for this sequence. If no I<limit> is provided, it
694 :     will default to 10.
695 :     It returns (match id, taxonomy string, log evalue, bitscore, alignment length,
696 :     percent identity, start1, end1) per hit.
697 :    
698 :     =cut
699 :    
700 :     sub get_hits_for_sequence {
701 :     my ($self, $id, $dataset, $limit) = @_;
702 :    
703 :     my $table = $self->dbtable;
704 :     my $dbid = $self->get_dataset_id($dataset);
705 :     $limit = 10 unless ($limit);
706 : arodri7 1.18
707 : olson 1.13 my $sth = $self->dbh->prepare(qq(SELECT id2, tax_str, logpsc, bsc, ali_ln,
708 :     iden, b1, e1
709 :     FROM $table
710 :     WHERE id1=? AND dbid=? AND rank_psc < ?));
711 :     $sth->execute($id, $dbid, $limit);
712 : arodri7 1.18
713 :     my $result = $sth->fetchall_arrayref();
714 :    
715 :     return $result;
716 :    
717 :     }
718 :    
719 :     =item * B<get_best_hit_for_sequence> (I<seq_id>, I<dataset>)
720 :    
721 :     Given a sequence id I<seq_id> (id1) and a dataset name (db_id), this method returns
722 :     the first row of hit data for this sequence.
723 :     It returns (match id, taxonomy string, log evalue, bitscore, alignment length,
724 :     percent identity, start1, end1) per hit.
725 :    
726 :     =cut
727 :    
728 :     sub get_best_hit_for_sequence {
729 :     my ($self, $id, $dataset) = @_;
730 :    
731 :     my $table = $self->dbtable;
732 :     my $dbid = $self->get_dataset_id($dataset);
733 :    
734 :     my $sth = $self->dbh->prepare(qq(SELECT id2, tax_str, logpsc, bsc, ali_ln,
735 :     iden, b1, e1, b2, e2
736 :     FROM $table
737 :     WHERE id1=? AND dbid=?));
738 :     $sth->execute($id, $dbid);
739 :    
740 : paarmann 1.4 my $result = $sth->fetchall_arrayref();
741 :    
742 :     return $result;
743 :    
744 :     }
745 :    
746 :    
747 : paarmann 1.3 =pod
748 :    
749 : paarmann 1.1 =item * B<get_align_len_range> (I<dataset_name>)
750 :    
751 :     Given a dataset name (db_id), this method returns
752 :     the minimum and maximum alignment length.
753 :    
754 :     =cut
755 :    
756 :     sub get_align_len_range {
757 :     my ($self, $dataset) = @_;
758 :    
759 :     my $table = $self->dbtable;
760 :     my $dbid = $self->get_dataset_id($dataset);
761 :    
762 : paarmann 1.3 my $sth = $self->dbh->prepare("select min(ali_ln), max(ali_ln) from $table where dbid=$dbid");
763 : paarmann 1.1 $sth->execute;
764 :     my ($min, $max) = $sth->fetchrow_array;
765 :    
766 :     return ($min, $max);
767 :    
768 :     }
769 :    
770 : jared 1.11 =pod
771 :    
772 :     =item * B<get_genome_id> (I<tax_str>)
773 :    
774 :     =cut
775 :    
776 :     sub get_genome_id {
777 :     my ($self, $tax_str) = @_;
778 : paczian 1.14 $tax_str =~ s/'/''/g;
779 :     my $retval = $self->dbh->selectrow_array("select seq_num from rdp_to_tax where tax_str='". $tax_str . "'");
780 :     if (ref($retval) eq 'ARRAY') {
781 :     return $retval->[0];
782 :     } else {
783 :     return $retval;
784 :     }
785 : jared 1.11 }
786 :    
787 : parrello 1.6 =pod
788 :    
789 : arodri7 1.22 =item * B<get_mg_comparison_table_metadata>
790 :    
791 :     returns the metadata (column names, headers) for the metagenome tables used in comparing metagenomes
792 :    
793 :     =cut
794 :     sub get_mg_comparison_table_metadata{
795 :     my ($self) = @_;
796 :     my $column_metadata = {};
797 :     my $desc = $self->data('dataset_desc');
798 :     my $metagenome = $self->application->cgi->param('metagenome') || '';
799 :     my $next_col;
800 :    
801 :     if (dataset_is_phylo($desc)){
802 :     $column_metadata->{Domain} = {'value'=>'Domain',
803 :     'header' => { name => 'Domain', filter => 1, operator => 'combobox',
804 :     visible => 0, show_control => 1 },
805 :     'order' => 1,
806 :     'visible' => 1,
807 :     'group' => 'permanent'
808 :     };
809 :     $column_metadata->{level1} = {'value'=>'Taxa Level 1',
810 :     'header' => { name => '', filter => 1, operator => 'combobox',
811 :     sortable => 1, width => 150 },
812 :     'order' => 2,
813 :     'visible' => 1,
814 :     'group' => 'permanent'
815 :     };
816 :     $column_metadata->{level2} = {'value'=>'Taxa Level 2',
817 :     'header' => { name => '', filter => 1, operator => 'combobox',
818 :     width => 150 },
819 :     'order' => 3,
820 :     'visible' => 1,
821 :     'group' => 'permanent'
822 :     };
823 :     $column_metadata->{level3} = {'value'=>'Taxa Level 3',
824 :     'header' => { name => '', filter => 1, operator => 'combobox',
825 :     width => 150 },
826 :     'order' => 4,
827 :     'visible' => 1,
828 :     'group' => 'permanent'
829 :     };
830 :     $column_metadata->{organism} = {'value'=>'Organism',
831 :     'header' => { name => 'Organism Name', filter => 1 },
832 :     'order' => 5,
833 :     'visible' => 1,
834 :     'group' => 'permanent'};
835 :     $next_col = 6;
836 :     }
837 :     elsif (dataset_is_metabolic($desc)){
838 :     $column_metadata->{hierarchy1} = {'value'=>'Subsystem Hierarchy 1',
839 :     'header' => { name => 'Subsystem Hierarchy 1', filter => 1,
840 :     operator => 'combobox', width => 150, sortable => 1 },
841 :     'order' => 1,
842 :     'visible' => 1,
843 :     'group' => 'permanent'
844 :     };
845 :     $column_metadata->{hierarchy2} = {'value'=>'Subsystem Hierarchy 1',
846 :     'header' => { name => 'Subsystem Hierarchy 2', filter => 1,
847 :     width => 150 },
848 :     'order' => 2,
849 :     'visible' => 1,
850 :     'group' => 'permanent'
851 :     };
852 :     $column_metadata->{hierarchy3} = {'value'=>'Subsystem Hierarchy 1',
853 :     'header' => { name => 'Subsystem Name', filter => 1,
854 :     sortable => 1, width => 150 },
855 :     'order' => 3,
856 :     'visible' => 1,
857 :     'group' => 'permanent'
858 :     };
859 :     $next_col = 4;
860 :     }
861 :    
862 :     # add your metagenome to permanent and add the other possible metagenomes to the select listbox
863 :     # check for available metagenomes
864 :     my $rast = $self->application->data_handle('RAST');
865 :     my $available = {};
866 :     if (ref($rast)) {
867 :     my $public_metagenomes = &get_public_metagenomes($self->app->dbmaster, $rast);
868 :     foreach my $pmg (@$public_metagenomes) {
869 :     $column_metadata->{$pmg->[0]} = {'value' => 'Public - ' . $pmg->[1],
870 :     'header' => { name => $pmg->[0],
871 :     filter => 1,
872 :     operators => ['equal', 'unequal', 'less', 'more'],
873 :     sortable => 1,
874 :     width => 150,
875 :     tooltip => $pmg->[1] . '(' . $pmg->[0] . ')'
876 :     },
877 :     };
878 :     if ($pmg->[0] eq $metagenome){
879 :     $column_metadata->{$pmg->[0]}->{order} = $next_col;
880 :     $column_metadata->{$pmg->[0]}->{visible} = 1;
881 :     $column_metadata->{$pmg->[0]}->{group} = 'permanent';
882 :     }
883 :     else{
884 :     $column_metadata->{$pmg->[0]}->{visible} = 0;
885 :     $column_metadata->{$pmg->[0]}->{group} = 'metagenomes';
886 :     }
887 :     }
888 :    
889 :     if ($self->application->session->user) {
890 :    
891 :     my $mgs = $rast->Job->get_jobs_for_user($self->application->session->user, 'view', 1);
892 :    
893 :     # build hash from all accessible metagenomes
894 :     foreach my $mg_job (@$mgs) {
895 :     $column_metadata->{$mg_job->genome_id} = {'value' => 'Private - ' . $mg_job->genome_name,
896 :     'header' => { name => $mg_job->genome_id,
897 :     filter => 1,
898 :     operators => ['equal', 'unequal', 'less', 'more'],
899 :     sortable => 1,
900 :     width => 150,
901 :     tooltip => $mg_job->genome_name . '(' . $mg_job->genome_id . ')'
902 :     },
903 :     };
904 :     if ( ($mg_job->metagenome) && ($mg_job->genome_id eq $metagenome) ) {
905 :     $column_metadata->{$mg_job->genome_id}->{order} = $next_col;
906 :     $column_metadata->{$mg_job->genome_id}->{visible} = 1;
907 :     $column_metadata->{$mg_job->genome_id}->{group} = 'permanent';
908 :     }
909 :     else{
910 :     $column_metadata->{$mg_job->genome_id}->{visible} = 0;
911 :     $column_metadata->{$mg_job->genome_id}->{group} = 'metagenomes';
912 :     }
913 :     }
914 :     }
915 :     }
916 :     else {
917 :     # no rast/user, no access to metagenomes
918 :     }
919 :    
920 :     return $column_metadata;
921 :     }
922 :    
923 :    
924 :     =pod
925 :    
926 : parrello 1.6 =back
927 : paarmann 1.1
928 : parrello 1.6 =cut

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3