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

View of /SeedViewer/MetagenomeAnalysis.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (download) (as text) (annotate)
Thu Apr 10 14:46:51 2008 UTC (11 years, 8 months ago) by paarmann
Branch: MAIN
Changes since 1.1: +1 -2 lines
removed debug statement

package SeedViewer::MetagenomeAnalysis;

# $Id: MetagenomeAnalysis.pm,v 1.2 2008/04/10 14:46:51 paarmann Exp $

use strict;
use warnings;

use FIG_Config;
use DBI;

1;




sub new {
  my ($class, $job) = @_;
  
  # check job
  unless (ref $job and $job->isa('RAST::Job') and $job->metagenome) {
    return undef;
  }

  # connect to database
  my $dbh;
  eval {

    my $host     = $FIG_Config::mgrast_dbhost;
    my $database = $FIG_Config::mgrast_db;
    my $user     = $FIG_Config::mgrast_dbuser;
    my $password = '';
    
    $dbh = DBI->connect("DBI:mysql:database=$database;host=$host", $user, $password, 
			{ RaiseError => 1, AutoCommit => 0, PrintError => 0 }) ||
					  die "database connect error.";
  };
  if ($@) {
    warn "Unable to connect to metagenomics database: $@\n";
    return undef;
  }

  # create object
  my $self = { job => $job,
	       dbh => $dbh,
	       key2taxa => undef,
	     };
  bless $self, $class;

  # load key2taxa mapping
  $self->get_key2taxa_mapping();
  
  return $self;

}




sub job {
  return $_[0]->{job};
}


sub dbh {
  return $_[0]->{dbh};
}


sub dbtable {
  unless (defined $_[0]->{dbtable}) {
    $_[0]->{dbtable} = 'tax_sim_'.$_[0]->job->id; 
  }
  return $_[0]->{dbtable};
}


sub get_key2taxa_mapping {

  unless (defined($_[0]->{key2taxa})) {
    my $sth  = $_[0]->dbh->prepare("select dbkey, str from tax_item");
    $sth->execute;
    $_[0]->{key2taxa} = $sth->fetchall_hashref('dbkey');
  }

  return $_[0]->{key2taxa};

}


sub key2taxa {    
  if(defined $_[1]) {
    my $t = $_[0]->{key2taxa}->{$_[1]}->{str};
    $t =~ s/\t+$//;
    return $t;
  }
  return '';
}


sub split_taxstr {
  my @r = split(':', $_[1]);
  return \@r;
}

sub join_taxstr {
    # do I really want an ending colon?
  return join(':', @{$_[1]}).':';
}


sub evalue2log {
  return 10 * (log($_[1]) / log(10));
}

sub log2evalue {
  return 10**($_[1]/10);
}


sub get_dataset_id {
    # mapping contained in table db_type
    # unstable, changes every time 
  my $db_ids = { 'SEED' => 2, 
		 'Greengenes' => 3, 
		 'RDP' => 1, 
		 'LSU' => 4, 
		 'SSU' => 5,
		 'Subsystem' => 6,
	       };
  return $db_ids->{ $_[1] } || undef;
}


=item * B<get_sequence> (I<sequence_id>)

Retrieve the sequence I<sequence_id> from the metagenome job directory.

=cut 

sub get_sequence {
  my ($self, $id) = @_;
  
  my $sequence_file = $self->job->org_dir.'/contigs';

  my $sequence = '';
  open(FASTA, "<$sequence_file") or die "Unable to read metagenome sequences: $!";
  while(<FASTA>) {
    next unless /^\>$id/;

    while(<FASTA>) {
      last if /^>/;
      chomp;
      $sequence .= $_;
    }
  }
  
  return $sequence;

}


=item * B<get_hits_count> (I<dataset_name>, I<where_clause>)

Given a dataset name (db_id) and optional a where clause, this method returns
the total number of sequences that contain a hit. 

=cut 

sub get_hits_count {
  my ($self, $dataset, $where) = @_;

  my $table = $self->dbtable; 
  my $dbid  = $self->get_dataset_id($dataset);
  $where = '1' unless ($where);

  my $sth = $self->dbh->prepare("select count(*) from ( select id1, min(rank_psc) from $table where dbid=$dbid and $where group by id1) as b");
  $sth->execute;
  my ($result) = $sth->fetchrow_array;

  return $result;

}


=item * B<get_align_len_range> (I<dataset_name>)

Given a dataset name (db_id), this method returns
the minimum and maximum alignment length.

=cut 

sub get_align_len_range {
  my ($self, $dataset) = @_;

  my $table = $self->dbtable; 
  my $dbid  = $self->get_dataset_id($dataset);

  my $sth = $self->dbh->prepare("select min(ali_ln), max(ali_ln) from $table");
  $sth->execute;
  my ($min, $max) = $sth->fetchrow_array;

  return ($min, $max);

}



MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3