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

Annotation of /SeedViewer/MetagenomeAnalysis.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.21 - (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 : jared 1.21 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 :     sub key2taxa {
124 :     if(defined $_[1]) {
125 :     my $t = $_[0]->{key2taxa}->{$_[1]}->{str};
126 :     $t =~ s/\t+$//;
127 :     return $t;
128 :     }
129 :     return '';
130 :     }
131 :    
132 :    
133 :     sub split_taxstr {
134 :     my @r = split(':', $_[1]);
135 :     return \@r;
136 :     }
137 :    
138 :     sub join_taxstr {
139 :     # do I really want an ending colon?
140 :     return join(':', @{$_[1]}).':';
141 :     }
142 :    
143 :    
144 :     sub evalue2log {
145 :     return 10 * (log($_[1]) / log(10));
146 :     }
147 :    
148 :     sub log2evalue {
149 :     return 10**($_[1]/10);
150 :     }
151 :    
152 :    
153 : olson 1.13 #
154 :     # Determine the correct dbid for this job. Use sims.database_list
155 :     # to find the version that the analysis was run with.
156 :     #
157 : paarmann 1.1 sub get_dataset_id {
158 : olson 1.13 my($self, $dataset) = @_;
159 :    
160 :     my $id = $self->{dbid_cache}->{$dataset};
161 :     return $id if defined($id);
162 :    
163 :     my($dbname, $type) = split(/:/, $dataset);
164 :    
165 :     my $dbs = $self->job->metaxml->get_metadata('sims.database_list');
166 :    
167 :     my @this = grep { $_->{name} eq $dbname } @$dbs;
168 :     if (@this)
169 :     {
170 :     my $vers = $this[0]->{version};
171 :    
172 :     #
173 :     # Now we can find the dbid ni the database.
174 :     #
175 :     my $res = $self->dbh->selectcol_arrayref(qq(SELECT dbid
176 :     FROM seq_db
177 :     WHERE name = ? AND version = ?
178 :     AND tax_db_name = ?), undef,
179 :     $dbname, $vers, $type);
180 :     if (@$res)
181 :     {
182 : paczian 1.14 #print STDERR "Found @$res for $dbname $type $vers\n";
183 : olson 1.13 $id = $res->[0];
184 :     $self->{dbid_cache}->{$dataset} = $id;
185 :     return $id;
186 :     }
187 : dsouza 1.20 #print STDERR "Did not find anything for dataset='$dataset' '$dbname' '$type' '$vers'\n";
188 : olson 1.13 }
189 : dsouza 1.20 #print STDERR "did not find a vers for dataset='$dataset' $dbname $type\n" . Dumper($dbs);
190 : paarmann 1.1 }
191 :    
192 : paarmann 1.3 #******************************************************************************
193 :     #* MANAGING QUERY CRITERIA
194 :     #******************************************************************************
195 :    
196 :     =pod
197 :    
198 : parrello 1.6 =over 4
199 :    
200 : paarmann 1.3 =item * B<query_evalue> (I<evalue>)
201 :    
202 :     Set/get the expectation value which is currently used to query the database.
203 :     Parameter I<evalue> has to be a float or in '1e-5'-like format or undef.
204 :    
205 :     =cut
206 :    
207 :     sub query_evalue {
208 :     if(scalar(@_)>1) {
209 :     $_[0]->{query}->{evalue} = $_[1];
210 :     }
211 :     return $_[0]->{query}->{evalue};
212 :     }
213 :    
214 :    
215 :     =pod
216 :    
217 :     =item * B<query_bitscore> (I<score>)
218 :    
219 :     Set/get the bitscore which is currently used to query the database.
220 :     Parameter I<score> has to be a float or undef.
221 :    
222 :     =cut
223 :    
224 :     sub query_bitscore {
225 :     if(scalar(@_)>1) {
226 :     $_[0]->{query}->{bitscore} = $_[1];
227 :     }
228 :     return $_[0]->{query}->{bitscore};
229 :     }
230 :    
231 :    
232 :     =pod
233 :    
234 :     =item * B<query_align_len> (I<length>)
235 :    
236 :     Set/get the minimum alignment which is currently used to query the database.
237 :     Parameter I<length> has to be a positive integer or undef.
238 :    
239 :     =cut
240 :    
241 :     sub query_align_len {
242 :     if(scalar(@_)>1) {
243 :     if($_[1] and $_[1]<0) {
244 :     die "Alignment length has to be positive: ".$_[1];
245 :     }
246 :     $_[0]->{query}->{align_len} = $_[1];
247 :     }
248 :     return $_[0]->{query}->{align_len};
249 :     }
250 :    
251 :    
252 :     =pod
253 :    
254 :     =item * B<query_identity> (I<percent>)
255 :    
256 :     Set/get the minimum percent identity which is currently used to query the database.
257 :     Parameter I<percent> has to be a number in 0..100 or undef.
258 :    
259 :     =cut
260 :    
261 :     sub query_identity {
262 :     if(scalar(@_)>1) {
263 :     if($_[1] and ($_[1]<0 or $_[1]>100)) {
264 :     die "Identity has to be between 0 and 100: ".$_[1];
265 :     }
266 :     $_[0]->{query}->{identity} = $_[1];
267 :     }
268 :     return $_[0]->{query}->{identity};
269 :     }
270 :    
271 :    
272 :     =pod
273 :    
274 :     =item * B<query_load_from_cgi> (I<cgi>, [I<dataset>])
275 :    
276 :     Sets all query parameter to the values provided in the CGI query object I<cgi>.
277 :     This method recognises 'evalue', 'pvalue' (bitscore), 'alignment_length' and
278 :     'percent_identity' as query criteria. Any missing param will be set to undef.
279 :     If the optional parameter I<dataset> is set to one of the accepted datasets
280 :     (db types), the method will additionally load the defaults for this type into
281 :     the CGI object.
282 :    
283 :    
284 :     =cut
285 :    
286 :     sub query_load_from_cgi {
287 :     my ($self, $cgi, $dataset) = @_;
288 :    
289 :     unless(ref $cgi and $cgi->isa("CGI")) {
290 :     die "Query load from cgi requires a valid CGI object.";
291 :     }
292 :    
293 :     # load the defaults if necessary
294 :     if($dataset and $self->get_dataset_id($dataset)) {
295 :    
296 :     my $d = $self->get_dataset_id($dataset);
297 :    
298 :     my @v = qw( evalue bitscore align_len identity );
299 :     foreach my $v (@v) {
300 :     if(!defined($cgi->param($v)) and QUERY_DEFAULTS->{$d}->{$v}) {
301 :     $cgi->param($v, QUERY_DEFAULTS->{$d}->{$v});
302 :     }
303 :     }
304 :     }
305 :    
306 :     # set the query params
307 :     my $evalue = $cgi->param('evalue') || '';
308 :     $self->query_evalue($evalue);
309 :    
310 :     my $bitscore = $cgi->param('bitscore') || '';
311 :     $self->query_bitscore($bitscore);
312 :    
313 :     my $align_len = $cgi->param('align_len') || '';
314 :     $self->query_align_len($align_len);
315 :    
316 :     my $identity = $cgi->param('identity') || '';
317 :     $self->query_identity($identity);
318 :    
319 :     return $self;
320 :    
321 :     }
322 :    
323 :    
324 :     =pod
325 :    
326 :     =item * B<get_where_clause> ()
327 :    
328 :     Returns for the current query parameters the where clause as applicable to the
329 :     tax_sim_XYZ table SQL queries. The method will take care of all conversions to
330 :     eg the logscore evalues.
331 :    
332 :     =cut
333 :    
334 :     sub get_where_clause {
335 :     my ($self) = @_;
336 :    
337 :     my @params;
338 :    
339 :     if($self->{query}->{evalue}) {
340 :     push @params, "logpsc<=".$self->evalue2log($self->{query}->{evalue});
341 :     }
342 :    
343 :     if($self->{query}->{bitscore}) {
344 :     push @params, "bsc>=".$self->{query}->{bitscore};
345 :     }
346 :    
347 :     if($self->{query}->{align_len}) {
348 :     push @params, "ali_ln>=".$self->{query}->{align_len};
349 :     }
350 :    
351 :     if($self->{query}->{identity}) {
352 :     push @params, "iden>=".$self->{query}->{identity};
353 :     }
354 :    
355 :     return join(' and ', @params);
356 :    
357 :     }
358 :    
359 :    
360 :    
361 :     #******************************************************************************
362 :     #* OTHER
363 :     #******************************************************************************
364 :    
365 :    
366 :     =pod
367 :    
368 : paarmann 1.1 =item * B<get_sequence> (I<sequence_id>)
369 :    
370 :     Retrieve the sequence I<sequence_id> from the metagenome job directory.
371 :    
372 :     =cut
373 :    
374 :     sub get_sequence {
375 :     my ($self, $id) = @_;
376 :    
377 :     my $sequence_file = $self->job->org_dir.'/contigs';
378 :    
379 :     my $sequence = '';
380 :     open(FASTA, "<$sequence_file") or die "Unable to read metagenome sequences: $!";
381 :     while(<FASTA>) {
382 :     next unless /^\>$id/;
383 :    
384 :     while(<FASTA>) {
385 :     last if /^>/;
386 :     chomp;
387 :     $sequence .= $_;
388 :     }
389 :     }
390 :    
391 :     return $sequence;
392 :    
393 :     }
394 :    
395 : paczian 1.15 =pod
396 :    
397 :     =item * B<get_sequences_fasta> (I<sequence_ids>)
398 :    
399 :     Retrieve the sequences for a given list of I<sequence_ids> from the metagenome job directory in fasta format.
400 :    
401 :     =cut
402 :    
403 :     sub get_sequences_fasta {
404 :     my ($self, $ids) = @_;
405 :    
406 :     # get the path to the sequence file
407 :     my $sequence_file = $self->job->org_dir.'/contigs';
408 :    
409 :     # hash the ids
410 :     my %idh = map { $_ => 1 } @$ids;
411 :    
412 :     # store the result
413 :     my $sequence = '';
414 :    
415 :     # iterate over the sequence file
416 :     open(FASTA, "<$sequence_file") or die "Unable to read metagenome sequences: $!";
417 :     while(<FASTA>) {
418 :    
419 :     # look only at id lines
420 :     next unless /^\>/;
421 :    
422 :     # get the id
423 :     my ($id) = /^\>(.+)/;
424 :     chomp($id);
425 :    
426 :     # find the id in the requested id hash
427 :     next unless exists ($idh{$id});
428 :    
429 :     # id was found, delete it from the hash
430 :     delete $idh{$id};
431 :    
432 :     # print the header line
433 :     $sequence .= $_;
434 :    
435 :     # print all the following lines to the result
436 :     while(<FASTA>) {
437 :     last if /^>/;
438 :     $sequence .= $_;
439 :     }
440 :    
441 :     last unless (scalar(keys(%idh)));
442 :     }
443 :    
444 :     return $sequence;
445 :    
446 :     }
447 : paarmann 1.1
448 : paarmann 1.3 =pod
449 : paarmann 1.1
450 : paarmann 1.3 =item * B<get_hits_count> (I<dataset_name>)
451 :    
452 :     Given a dataset name (db_id), this method returns
453 : paarmann 1.1 the total number of sequences that contain a hit.
454 :    
455 :     =cut
456 :    
457 :     sub get_hits_count {
458 : paarmann 1.3 my ($self, $dataset) = @_;
459 : jared 1.10
460 :     my $table = $self->dbtable_best_psc;
461 : paarmann 1.1 my $dbid = $self->get_dataset_id($dataset);
462 : paarmann 1.3 my $where = $self->get_where_clause();
463 :     $where = ($where) ? "and $where" : '';
464 : paarmann 1.1
465 : jared 1.10 my $sth = $self->dbh->prepare("select count(distinct id1) from $table where dbid=$dbid $where");
466 : paarmann 1.1 $sth->execute;
467 :     my ($result) = $sth->fetchrow_array;
468 :    
469 :     return $result;
470 :    
471 :     }
472 :    
473 :    
474 : paarmann 1.3 =pod
475 :    
476 :     =item * B<get_group_counts> (I<dataset_name>, [I<group>, I<filter1>, I<filter2>])
477 :    
478 :     Given a dataset name (db_id), this method returns the total counts for all
479 :     taxonomy groups of a certain depth which are hit. If no group name I<group>
480 :     was given, the method returns counts for tax_group_1.
481 :     Optionally, I<group> may be 'tax_group_2' or 'tax_group_3' and in that case
482 :     any optional provided filters I<filter1> and I<filter2> will be applied to
483 :     the column 'tax_group_1' and 'tax_group_2' respectively.
484 :    
485 :     =cut
486 :    
487 :     sub get_group_counts {
488 :     my ($self, $dataset, $group, $filter1, $filter2) = @_;
489 :    
490 : jared 1.10 my $table = $self->dbtable_best_psc;
491 : paarmann 1.3 my $dbid = $self->get_dataset_id($dataset);
492 :     my $where = $self->get_where_clause();
493 :     $where = ($where) ? "and $where" : '';
494 :     $group = 'tax_group_1' unless($group);
495 :    
496 :     my @filters;
497 :     push @filters, "tax_group_1='$filter1'" if($filter1);
498 :     push @filters, "tax_group_2='$filter2'" if($filter2);
499 :     my $filter = (scalar(@filters)) ? 'and '.join(' and ', @filters) : '';
500 :    
501 : 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");
502 : paarmann 1.3 $sth->execute;
503 :     my $result = $sth->fetchall_arrayref();
504 :    
505 : dsouza 1.20 #print STDERR "get_group_counts: ds=$dataset group=$group filter1=$filter1 filter2=$filter2\n";
506 :     #print STDERR Dumper($result);
507 : paarmann 1.3 return $result;
508 :    
509 :     }
510 :    
511 :    
512 :     =pod
513 :    
514 :     =item * B<get_taxa_counts> (I<dataset_name>)
515 :    
516 :     Given a dataset name (db_id), this method returns the total counts for all
517 :     taxonomy strings which are hit.
518 :    
519 :     =cut
520 :    
521 :     sub get_taxa_counts {
522 :     my ($self, $dataset) = @_;
523 :    
524 : jared 1.9 my $table = $self->dbtable_best_psc;
525 : paarmann 1.3 my $dbid = $self->get_dataset_id($dataset);
526 :     my $where = $self->get_where_clause();
527 :     $where = ($where) ? "and $where" : '';
528 :    
529 : paczian 1.16 my $sth = $self->dbh->prepare("select tax_str as tax, count(*) from $table where dbid=$dbid $where group by tax");
530 : paarmann 1.3
531 :     $sth->execute;
532 :     my $result = $sth->fetchall_arrayref();
533 :    
534 :     return $result;
535 :    
536 :     }
537 :    
538 :    
539 : jared 1.21 sub get_tax2genomeid {
540 :     my ($self, $taxstring) = @_;
541 :    
542 :     my ($genomeid) = $self->dbh->selectrow_array(qq(SELECT seq_num
543 :     FROM rdp_to_tax
544 :     WHERE tax_str= ?), undef, $taxstring);
545 :    
546 :     return $genomeid;
547 :     }
548 :    
549 :    
550 : paarmann 1.3 =pod
551 :    
552 :     =item * B<get_subsystem_counts> (I<dataset_name>)
553 :    
554 :     Given a dataset name (db_id), this method returns the total counts for all
555 :     subsystems which are hit.
556 :    
557 :     =cut
558 :    
559 :     sub get_subsystem_counts {
560 :     my ($self, $dataset) = @_;
561 :    
562 : jared 1.9 my $table = $self->dbtable_best_psc;
563 : paarmann 1.3 my $dbid = $self->get_dataset_id($dataset);
564 :     my $where = $self->get_where_clause();
565 :     $where = ($where) ? "and $where" : '';
566 : paarmann 1.4
567 : 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");
568 : paarmann 1.3
569 :     $sth->execute;
570 :     my $result = $sth->fetchall_arrayref();
571 :    
572 :     return $result;
573 :    
574 :     }
575 :    
576 :    
577 :     =pod
578 :    
579 :     =item * B<get_sequence_subset> (I<dataset_name>, I<filter>)
580 :    
581 :     Given a dataset name (db_id), this method returns all sequence ids,
582 :     the alignment length, the match id and the taxonomy string for all
583 :     sequences which match the criteria and have their tax_str start with
584 :     the filter string I<filter>.
585 :    
586 :     =cut
587 :    
588 :     sub get_sequence_subset {
589 :     my ($self, $dataset, $filter) = @_;
590 :    
591 : jared 1.10 my $table = $self->dbtable_best_psc;
592 : paarmann 1.3 my $dbid = $self->get_dataset_id($dataset);
593 :     my $where = $self->get_where_clause();
594 :     $where = ($where) ? "and $where" : '';
595 :    
596 : paczian 1.15 $filter =~ s/'/''/g;
597 :    
598 : 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%'");
599 : paarmann 1.3 $sth->execute;
600 :     my $result = $sth->fetchall_arrayref();
601 :    
602 :     return $result;
603 :    
604 :     }
605 :    
606 : jared 1.19 =pod
607 :    
608 :     =item * B<get_sequence_subset_genome> (I<genome>)
609 :    
610 :     Given a dataset name (db_id), this method returns all sequence ids,
611 :     the alignment length, the match id and the taxonomy string for all
612 :     sequences which match the criteria and have their tax_str start with
613 :     the filter string I<filter>.
614 :    
615 :     =cut
616 :    
617 :     sub get_sequence_subset_genome {
618 :     my ($self, $genome) = @_;
619 :    
620 :     my $table = $self->dbtable_best_psc;
621 :     my $dbid = $self->get_dataset_id("SEED:seed_genome_tax");
622 :     my $where = $self->get_where_clause();
623 :     $where = ($where) ? "and $where" : '';
624 :    
625 :     my ($tax_id) = $self->dbh->selectrow_array(qq(SELECT tax_str
626 :     FROM rdp_to_tax
627 :     WHERE seq_num= ?), undef, $genome);
628 :    
629 :     if($tax_id =~ /(\S+)\s/){
630 :     $tax_id = $1;
631 :     }
632 :    
633 :     my $sth = $self->dbh->prepare(qq(SELECT id1, ali_ln, id2, tax_str
634 :     FROM $table
635 :     WHERE dbid=? $where and tax_str=?));
636 :    
637 :    
638 :     $sth->execute($dbid, $tax_id);
639 :     my $result = $sth->fetchall_arrayref();
640 :    
641 :     return $result;
642 :    
643 :     }
644 : paarmann 1.3
645 : paarmann 1.4 =pod
646 :    
647 : jared 1.5 =item * B<get_recruitment_plot_data> (I<genome>)
648 :    
649 :     Given a genome id (83333.1), this method returns all sequence ids,
650 :     the alignment length, the match id and the taxonomy string for all
651 :     sequences which match the criteria and have their tax_str start equal
652 :     the genome tax string I<filter>.
653 :    
654 :     =cut
655 :    
656 :     sub get_recruitment_plot_data {
657 :     my ($self, $genome) = @_;
658 :    
659 : jared 1.9 my $table = $self->dbtable_best_psc;
660 : olson 1.13 my $dbid = $self->get_dataset_id("SEED:seed_genome_tax");
661 : jared 1.5 my $where = $self->get_where_clause();
662 :     $where = ($where) ? "and $where" : '';
663 :    
664 : olson 1.13 my ($tax_id) = $self->dbh->selectrow_array(qq(SELECT tax_str
665 :     FROM rdp_to_tax
666 : jared 1.21 WHERE seq_num= ? and dbid= ?), undef, $genome, $dbid);
667 : jared 1.5
668 : jared 1.9 if($tax_id =~ /(\S+)\s/){
669 :     $tax_id = $1;
670 :     }
671 : jared 1.5
672 : olson 1.13 my $sth = $self->dbh->prepare(qq(SELECT id1, id2, b2, e2, logpsc
673 :     FROM $table
674 :     WHERE dbid=? $where and tax_str=?));
675 :    
676 :    
677 :     $sth->execute($dbid, $tax_id);
678 : jared 1.5 my $result = $sth->fetchall_arrayref();
679 :    
680 :     return $result;
681 :    
682 :     }
683 :    
684 :    
685 :    
686 :    
687 :     =pod
688 :    
689 : paarmann 1.4 =item * B<get_hits_for_sequence> (I<seq_id>, I<dataset>, I<limit>)
690 :    
691 :     Given a sequence id I<seq_id> (id1) and a dataset name (db_id), this method returns
692 :     the first I<limit> rows of hit data for this sequence. If no I<limit> is provided, it
693 :     will default to 10.
694 :     It returns (match id, taxonomy string, log evalue, bitscore, alignment length,
695 :     percent identity, start1, end1) per hit.
696 :    
697 :     =cut
698 :    
699 :     sub get_hits_for_sequence {
700 :     my ($self, $id, $dataset, $limit) = @_;
701 :    
702 :     my $table = $self->dbtable;
703 :     my $dbid = $self->get_dataset_id($dataset);
704 :     $limit = 10 unless ($limit);
705 : arodri7 1.18
706 : olson 1.13 my $sth = $self->dbh->prepare(qq(SELECT id2, tax_str, logpsc, bsc, ali_ln,
707 :     iden, b1, e1
708 :     FROM $table
709 :     WHERE id1=? AND dbid=? AND rank_psc < ?));
710 :     $sth->execute($id, $dbid, $limit);
711 : arodri7 1.18
712 :     my $result = $sth->fetchall_arrayref();
713 :    
714 :     return $result;
715 :    
716 :     }
717 :    
718 :     =item * B<get_best_hit_for_sequence> (I<seq_id>, I<dataset>)
719 :    
720 :     Given a sequence id I<seq_id> (id1) and a dataset name (db_id), this method returns
721 :     the first row of hit data for this sequence.
722 :     It returns (match id, taxonomy string, log evalue, bitscore, alignment length,
723 :     percent identity, start1, end1) per hit.
724 :    
725 :     =cut
726 :    
727 :     sub get_best_hit_for_sequence {
728 :     my ($self, $id, $dataset) = @_;
729 :    
730 :     my $table = $self->dbtable;
731 :     my $dbid = $self->get_dataset_id($dataset);
732 :    
733 :     my $sth = $self->dbh->prepare(qq(SELECT id2, tax_str, logpsc, bsc, ali_ln,
734 :     iden, b1, e1, b2, e2
735 :     FROM $table
736 :     WHERE id1=? AND dbid=?));
737 :     $sth->execute($id, $dbid);
738 :    
739 : paarmann 1.4 my $result = $sth->fetchall_arrayref();
740 :    
741 :     return $result;
742 :    
743 :     }
744 :    
745 :    
746 : paarmann 1.3 =pod
747 :    
748 : paarmann 1.1 =item * B<get_align_len_range> (I<dataset_name>)
749 :    
750 :     Given a dataset name (db_id), this method returns
751 :     the minimum and maximum alignment length.
752 :    
753 :     =cut
754 :    
755 :     sub get_align_len_range {
756 :     my ($self, $dataset) = @_;
757 :    
758 :     my $table = $self->dbtable;
759 :     my $dbid = $self->get_dataset_id($dataset);
760 :    
761 : paarmann 1.3 my $sth = $self->dbh->prepare("select min(ali_ln), max(ali_ln) from $table where dbid=$dbid");
762 : paarmann 1.1 $sth->execute;
763 :     my ($min, $max) = $sth->fetchrow_array;
764 :    
765 :     return ($min, $max);
766 :    
767 :     }
768 :    
769 : jared 1.11 =pod
770 :    
771 :     =item * B<get_genome_id> (I<tax_str>)
772 :    
773 :     =cut
774 :    
775 :     sub get_genome_id {
776 :     my ($self, $tax_str) = @_;
777 : paczian 1.14 $tax_str =~ s/'/''/g;
778 :     my $retval = $self->dbh->selectrow_array("select seq_num from rdp_to_tax where tax_str='". $tax_str . "'");
779 :     if (ref($retval) eq 'ARRAY') {
780 :     return $retval->[0];
781 :     } else {
782 :     return $retval;
783 :     }
784 : jared 1.11 }
785 :    
786 : parrello 1.6 =pod
787 :    
788 :     =back
789 : paarmann 1.1
790 : parrello 1.6 =cut

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3