[Bio] / FigKernelPackages / RAST_submission.pm Repository:
ViewVC logotype

Diff of /FigKernelPackages/RAST_submission.pm

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.3, Wed Oct 21 17:19:31 2009 UTC revision 1.4, Wed Oct 21 20:05:36 2009 UTC
# Line 1  Line 1 
1    
2  package RAST_submission;  package RAST_submission;
3    
4    
5  use strict;  use strict;
6  use Job48;  use Job48;
7  use JobUpload;  use JobUpload;
8  use Data::Dumper;  use Data::Dumper;
9    use FIG;
10    use FIG_Config;
11    use gjoseqlib;
12    
13    use LWP::UserAgent;
14    use Bio::DB::RefSeq;
15    use Bio::SeqIO;
16    
17  use base 'Class::Accessor';  use base 'Class::Accessor';
18    
19  __PACKAGE__->mk_accessors(qw(rast_dbmaster user_dbmaster user_obj));  __PACKAGE__->mk_accessors(qw(rast_dbmaster user_dbmaster user_obj project_cache_dir
20                                 contig_cache_dir max_cache_age ua));
21    
22  sub new  sub new
23  {  {
# Line 18  Line 27 
27          rast_dbmaster => $rast_dbmaster,          rast_dbmaster => $rast_dbmaster,
28          user_dbmaster => $user_dbmaster,          user_dbmaster => $user_dbmaster,
29          user_obj => $user_obj,          user_obj => $user_obj,
30            project_cache_dir => "$FIG_Config::var/ncbi_project_cache",
31            contig_cache_dir => "$FIG_Config::var/ncbi_contig_cache",
32            max_cache_age => 86400,
33            ua => LWP::UserAgent->new(),
34      };      };
35    
36        &FIG::verify_dir($self->{project_cache_dir});
37        &FIG::verify_dir($self->{contig_cache_dir});
38    
39    
40      return bless $self, $class;      return bless $self, $class;
41  }  }
42    
43    sub get_contig_ids_in_project_from_entrez
44    {
45        my($self, $params) = @_;
46    
47        #
48        # Determine the project ID to use. Which one we take depends on if
49        # we were passed a project id, a tax id, or a contig id.
50        #
51    
52        my $proj;
53        if ($params->{-tax_id})
54        {
55        }
56        elsif ($params->{-contig_id})
57        {
58            $proj = $self->determine_project_of_contig($params->{-contig_id});
59        }
60        elsif ($params->{-project_id})
61        {
62            $proj = $params->{-project_id};
63        }
64    
65        print STDERR "project is $proj\n";
66        my $project_data = $self->retrieve_project_data($proj);
67    
68        return $self->check_project_for_redundancy($project_data);
69    }
70    
71    sub determine_project_of_contig
72    {
73        my($self, $contig_id) = @_;
74    
75        my $file = $self->retrieve_contig_data($contig_id);
76        open(F, "<", $file) or die "cannot open contig data $file: $!";
77    
78        my $proj;
79        while (<F>)
80        {
81            if (/DBLINK\s+Project:(\d+)/)
82            {
83                $proj = $1;
84                last;
85            }
86        }
87        close(F);
88        return $proj;
89    
90    }
91    
92    sub check_project_for_redundancy
93    {
94        my($self, $file) = @_;
95    
96        my $seqio_object = Bio::SeqIO->new(
97                                           -file => $file ,
98                                           -format => "genbank",
99                                          );
100    
101        my @seqs;
102        my @ids;
103        while ( my $seq = $seqio_object->next_seq ) {
104            push(@seqs, [$seq->accession_number, $seq->seq]);
105            push(@ids, $seq->accession_number);
106        }
107    
108        my @redundancy = $self->test_for_redundancy(\@seqs);
109        return { ids => \@ids, redundancy_report => \@redundancy };
110    }
111    
112    sub test_for_redundancy {
113        my($self, $seqs) = @_;
114    
115        if (@$seqs < 2)
116        {
117            return ();
118        }
119    
120        my %lens = map { $_->[0] => length($_->[1]) } @$seqs;
121        my $tmp = "tmp.$$.fasta";
122        &gjoseqlib::print_alignment_as_fasta($tmp,$seqs);
123        system "formatdb -i $tmp -pF";
124        my @blastout = `blastall -m8 -i $tmp -d $tmp -p blastn -FF -e 1.0e-100`;
125        system "rm $tmp $tmp\.*";
126        my @tuples = ();
127        my %seen;
128        foreach my $hit (map { chomp; [split(/\t/,$_)] } @blastout)
129        {
130            my($id1,$id2,$iden,undef,undef,undef,$b1,$e1,$b2,$e2) = @$hit;
131            if ((! $seen{"$id1/$id2"}) && ($id1 ne $id2))
132            {
133                $seen{"$id1/$id2"} = 1;
134                if (($iden >= 98) &&
135                    (abs($e1 - $b1) > (0.9 * $lens{$id1})))
136                {
137                    push(@tuples,[$id1,$lens{$id1},$id2,$lens{$id2}]);
138                }
139            }
140        }
141    
142        return @tuples;
143    }
144    
145    sub retrieve_project_data
146    {
147        my($self, $project) = @_;
148    
149        my $cached_file = $self->project_cache_dir() . "/$project.gbff";
150        if (my(@stat) = stat($cached_file))
151        {
152            my $last_mod = $stat[9];
153            if (time - $last_mod < $self->max_cache_age)
154            {
155                return $cached_file;
156            }
157        }
158    
159        my $url = "http://www.ncbi.nlm.nih.gov/sites/entrez?Db=genomeprj&Cmd=Retrieve&list_uids=";
160        my $res = $self->ua->get($url.$project);
161        if (!$res->is_success)
162        {
163            die "error retrieving project data: " . $res->status_line;
164        }
165        my $search_result = $res->content;
166    
167        my @lines = split ( "\n" , $search_result);
168    
169        my $nr_seq = 0;
170        my $nr_proj = 0;
171        my $url_seq = "";
172        my $url_proj = "";
173        my $genome_name = "";
174    
175        my $found_genome_information_table = 0;
176    
177        my $next = "";
178        my $id_list  = "";
179        foreach my $line ( @lines )
180        {
181    
182            if ($line =~/Genome information:/){
183                $found_genome_information_table = 1;
184                next; # skip table line
185            };
186    
187            $found_genome_information_table = 0 if ( $found_genome_information_table and $line =~ /<\/table>/);
188            $id_list .= $line if ( $found_genome_information_table ); # collect id entries
189    
190            # print $line , "\n"  if (  $found_genome_information_table );
191        }
192    
193        my @ids;
194        my @blocks = split "<\/tr>" , $id_list ;
195        foreach my $block (@blocks)
196        {
197            my @local_ids = $block =~/([^>]+)<\/a><\/td>/gc ;
198            # print join "\t" , @local_ids  , "\n";
199            push @ids , $local_ids[0] if ($local_ids[0]);
200        }
201    
202        my $id_list = join(",", @ids);
203        my $query = "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nucleotide&id=" . $id_list . "&rettype=gb" ;
204        print STDERR $query , "\n";
205        my $resp = $self->ua->get($query);
206        if ($resp->is_success())
207        {
208            open(F, ">", $cached_file) or die "Cannot open $cached_file for writing: $!";
209            print F $resp->content;
210            close(F);
211            return $cached_file;
212        }
213        else
214        {
215            die "Error retrieving data: " . $resp->status_line;
216        }
217    }
218    
219    sub retrieve_contig_data
220    {
221        my($self, $contig) = @_;
222    
223        my $cached_file = $self->contig_cache_dir() . "/$contig.gbff";
224        if (my(@stat) = stat($cached_file))
225        {
226            my $last_mod = $stat[9];
227            if (time - $last_mod < $self->max_cache_age)
228            {
229                return $cached_file;
230            }
231        }
232    
233        my $query = "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nucleotide&id=" . $contig . "&rettype=gb" ;
234        print STDERR $query , "\n";
235        my $resp = $self->ua->get($query);
236        if ($resp->is_success())
237        {
238            open(F, ">", $cached_file) or die "Cannot open $cached_file for writing: $!";
239            print F $resp->content;
240            close(F);
241            return $cached_file;
242        }
243        else
244        {
245            die "Error retrieving data: " . $resp->status_line;
246        }
247    }
248    
249    
250  sub submit_RAST_job  sub submit_RAST_job
251  {  {
252      my($self, $params) = @_;      my($self, $params) = @_;

Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3