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

Annotation of /SeedViewer/MetagenomeAnalysis.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : paarmann 1.1 package SeedViewer::MetagenomeAnalysis;
2 :    
3 : paarmann 1.2 # $Id: MetagenomeAnalysis.pm,v 1.1 2008/04/07 21:20:45 paarmann 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 :    
16 :     sub new {
17 :     my ($class, $job) = @_;
18 :    
19 :     # check job
20 :     unless (ref $job and $job->isa('RAST::Job') and $job->metagenome) {
21 :     return undef;
22 :     }
23 :    
24 :     # connect to database
25 :     my $dbh;
26 :     eval {
27 :    
28 :     my $host = $FIG_Config::mgrast_dbhost;
29 :     my $database = $FIG_Config::mgrast_db;
30 :     my $user = $FIG_Config::mgrast_dbuser;
31 :     my $password = '';
32 :    
33 :     $dbh = DBI->connect("DBI:mysql:database=$database;host=$host", $user, $password,
34 :     { RaiseError => 1, AutoCommit => 0, PrintError => 0 }) ||
35 :     die "database connect error.";
36 :     };
37 :     if ($@) {
38 :     warn "Unable to connect to metagenomics database: $@\n";
39 :     return undef;
40 :     }
41 :    
42 :     # create object
43 :     my $self = { job => $job,
44 :     dbh => $dbh,
45 :     key2taxa => undef,
46 :     };
47 :     bless $self, $class;
48 :    
49 :     # load key2taxa mapping
50 :     $self->get_key2taxa_mapping();
51 :    
52 :     return $self;
53 :    
54 :     }
55 :    
56 :    
57 :    
58 :    
59 :     sub job {
60 :     return $_[0]->{job};
61 :     }
62 :    
63 :    
64 :     sub dbh {
65 :     return $_[0]->{dbh};
66 :     }
67 :    
68 :    
69 :     sub dbtable {
70 :     unless (defined $_[0]->{dbtable}) {
71 :     $_[0]->{dbtable} = 'tax_sim_'.$_[0]->job->id;
72 :     }
73 :     return $_[0]->{dbtable};
74 :     }
75 :    
76 :    
77 :     sub get_key2taxa_mapping {
78 :    
79 :     unless (defined($_[0]->{key2taxa})) {
80 :     my $sth = $_[0]->dbh->prepare("select dbkey, str from tax_item");
81 :     $sth->execute;
82 :     $_[0]->{key2taxa} = $sth->fetchall_hashref('dbkey');
83 :     }
84 :    
85 :     return $_[0]->{key2taxa};
86 :    
87 :     }
88 :    
89 :    
90 :     sub key2taxa {
91 :     if(defined $_[1]) {
92 :     my $t = $_[0]->{key2taxa}->{$_[1]}->{str};
93 :     $t =~ s/\t+$//;
94 :     return $t;
95 :     }
96 :     return '';
97 :     }
98 :    
99 :    
100 :     sub split_taxstr {
101 :     my @r = split(':', $_[1]);
102 :     return \@r;
103 :     }
104 :    
105 :     sub join_taxstr {
106 :     # do I really want an ending colon?
107 :     return join(':', @{$_[1]}).':';
108 :     }
109 :    
110 :    
111 :     sub evalue2log {
112 :     return 10 * (log($_[1]) / log(10));
113 :     }
114 :    
115 :     sub log2evalue {
116 :     return 10**($_[1]/10);
117 :     }
118 :    
119 :    
120 :     sub get_dataset_id {
121 :     # mapping contained in table db_type
122 :     # unstable, changes every time
123 :     my $db_ids = { 'SEED' => 2,
124 :     'Greengenes' => 3,
125 :     'RDP' => 1,
126 :     'LSU' => 4,
127 :     'SSU' => 5,
128 :     'Subsystem' => 6,
129 :     };
130 :     return $db_ids->{ $_[1] } || undef;
131 :     }
132 :    
133 :    
134 :     =item * B<get_sequence> (I<sequence_id>)
135 :    
136 :     Retrieve the sequence I<sequence_id> from the metagenome job directory.
137 :    
138 :     =cut
139 :    
140 :     sub get_sequence {
141 :     my ($self, $id) = @_;
142 :    
143 :     my $sequence_file = $self->job->org_dir.'/contigs';
144 :    
145 :     my $sequence = '';
146 :     open(FASTA, "<$sequence_file") or die "Unable to read metagenome sequences: $!";
147 :     while(<FASTA>) {
148 :     next unless /^\>$id/;
149 :    
150 :     while(<FASTA>) {
151 :     last if /^>/;
152 :     chomp;
153 :     $sequence .= $_;
154 :     }
155 :     }
156 :    
157 :     return $sequence;
158 :    
159 :     }
160 :    
161 :    
162 :     =item * B<get_hits_count> (I<dataset_name>, I<where_clause>)
163 :    
164 :     Given a dataset name (db_id) and optional a where clause, this method returns
165 :     the total number of sequences that contain a hit.
166 :    
167 :     =cut
168 :    
169 :     sub get_hits_count {
170 :     my ($self, $dataset, $where) = @_;
171 :    
172 :     my $table = $self->dbtable;
173 :     my $dbid = $self->get_dataset_id($dataset);
174 :     $where = '1' unless ($where);
175 :    
176 :     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");
177 :     $sth->execute;
178 :     my ($result) = $sth->fetchrow_array;
179 :    
180 :     return $result;
181 :    
182 :     }
183 :    
184 :    
185 :     =item * B<get_align_len_range> (I<dataset_name>)
186 :    
187 :     Given a dataset name (db_id), this method returns
188 :     the minimum and maximum alignment length.
189 :    
190 :     =cut
191 :    
192 :     sub get_align_len_range {
193 :     my ($self, $dataset) = @_;
194 :    
195 :     my $table = $self->dbtable;
196 :     my $dbid = $self->get_dataset_id($dataset);
197 :    
198 :     my $sth = $self->dbh->prepare("select min(ali_ln), max(ali_ln) from $table");
199 :     $sth->execute;
200 :     my ($min, $max) = $sth->fetchrow_array;
201 :    
202 :     return ($min, $max);
203 :    
204 :     }
205 :    
206 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3