|
0
|
1 # $Id: BioPerl.pm,v 1.15 2001/12/14 16:40:15 heikki Exp $
|
|
|
2 #
|
|
|
3 # bioperl module for Bio::LiveSeq::IO::BioPerl
|
|
|
4 #
|
|
|
5 # Cared for by Joseph Insana <insana@ebi.ac.uk> <jinsana@gmx.net>
|
|
|
6 #
|
|
|
7 # Copyright Joseph Insana
|
|
|
8 #
|
|
|
9 # You may distribute this module under the same terms as perl itself
|
|
|
10 #
|
|
|
11 # POD documentation - main docs before the code
|
|
|
12
|
|
|
13 =head1 NAME
|
|
|
14
|
|
|
15 Bio::LiveSeq::IO::BioPerl - Loader for LiveSeq from EMBL entries with BioPerl
|
|
|
16
|
|
|
17 =head1 SYNOPSIS
|
|
|
18
|
|
|
19 my $db="EMBL";
|
|
|
20 my $file="../data/M20132";
|
|
|
21 my $id="HSANDREC";
|
|
|
22
|
|
|
23 my $loader=Bio::LiveSeq::IO::BioPerl->load(-db=>"$db", -file=>"$file");
|
|
|
24 or
|
|
|
25 my $loader=Bio::LiveSeq::IO::BioPerl->load(-db=>"$db", -id=>"$id");
|
|
|
26
|
|
|
27 my @translationobjects=$loader->entry2liveseq();
|
|
|
28
|
|
|
29 my $genename="AR";
|
|
|
30 my $gene=$loader->gene2liveseq(-gene_name => "$genename",
|
|
|
31 -getswissprotinfo => 0);
|
|
|
32
|
|
|
33 NOTE1: The only -db now supported is EMBL. Hence it defaults to EMBL.
|
|
|
34 NOTE2: -file requires a filename (and path if necessary) containing an
|
|
|
35 EMBL entry
|
|
|
36 -id will use Bio::DB::EMBL.pm to fetch the sequence from the web,
|
|
|
37 (bioperl wraparound to [w]getz from SRS)
|
|
|
38 NOTE3: To retrieve the swissprot (if possible) attached to the embl entry
|
|
|
39 (to get protein domains at dna level), only Bio::DB::EMBL.pm
|
|
|
40 is supported under BioPerl. Refer to Bio::LiveSeq::IO::SRS
|
|
|
41 otherwise.
|
|
|
42 NOTE4: NOTE3 is not implemented yet for bioperl, working on it
|
|
|
43
|
|
|
44
|
|
|
45 =head1 DESCRIPTION
|
|
|
46
|
|
|
47 This package uses BioPerl (SeqIO) to fetch a sequence database entry,
|
|
|
48 analyse it and create LiveSeq objects out of it.
|
|
|
49
|
|
|
50 A filename (or an ID that will fetch entry through the web) has to be passed
|
|
|
51 to this package which will return references to all translation objects
|
|
|
52 created from the EMBL entry. References to Transcription, DNA and Exon
|
|
|
53 objects can all be retrieved departing from these.
|
|
|
54
|
|
|
55 Alternatively, a specific "gene" name can be specified, together with
|
|
|
56 the embl-acc ID. This will create a LiveSeq::Gene object with all
|
|
|
57 relevant gene features attached/created.
|
|
|
58
|
|
|
59 ATTENTION: if web fetching is requested, the package HTTP::Request needs
|
|
|
60 to be installed.
|
|
|
61
|
|
|
62
|
|
|
63 =head1 AUTHOR - Joseph A.L. Insana
|
|
|
64
|
|
|
65 Email: Insana@ebi.ac.uk, jinsana@gmx.net
|
|
|
66
|
|
|
67 Address:
|
|
|
68
|
|
|
69 EMBL Outstation, European Bioinformatics Institute
|
|
|
70 Wellcome Trust Genome Campus, Hinxton
|
|
|
71 Cambs. CB10 1SD, United Kingdom
|
|
|
72
|
|
|
73 =head1 APPENDIX
|
|
|
74
|
|
|
75 The rest of the documentation details each of the object
|
|
|
76 methods. Internal methods are usually preceded with a _
|
|
|
77
|
|
|
78 =cut
|
|
|
79
|
|
|
80 # Let the code begin...
|
|
|
81
|
|
|
82 package Bio::LiveSeq::IO::BioPerl;
|
|
|
83 $VERSION=2.42;
|
|
|
84
|
|
|
85 # Version history:
|
|
|
86 # Thu Apr 6 00:25:46 BST 2000 v 1.0 begun
|
|
|
87 # Thu Apr 6 03:40:04 BST 2000 v 1.25 added Division
|
|
|
88 # Thu Apr 6 03:40:36 BST 2000 v 2.0 working
|
|
|
89 # Thu Apr 20 02:17:28 BST 2000 v 2.1 mRNA added to valid_feature_names
|
|
|
90 # Tue Jul 4 14:07:52 BST 2000 v 2.11 note&number added in val_qual_names
|
|
|
91 # Fri Sep 15 15:41:02 BST 2000 v 2.22 novelaasequence2gene now works without SRS
|
|
|
92 # Mon Jan 29 17:40:06 EST 2001 v 2.3 made it work with the new split_location of BioPerl 0.7
|
|
|
93 # Tue Apr 10 17:00:18 BST 2001 v 2.41 started work on support of DB::EMBL.pm
|
|
|
94 # Tue Apr 10 17:22:26 BST 2001 v 2.42 -id should work now
|
|
|
95
|
|
|
96 # TODO->TOCHECK
|
|
|
97 # each_secondary_access not working
|
|
|
98 # why array from each_tag_value($qual) ? When will there be more than one
|
|
|
99 # element in such array?
|
|
|
100 # what is the annotation object? ($seqobj->annotation)
|
|
|
101 # unsatisfied by both BioPerl binomial and SRS "org" to retrieve Organism info
|
|
|
102
|
|
|
103 use strict;
|
|
|
104 use Carp qw(cluck croak carp);
|
|
|
105 use vars qw($VERSION @ISA);
|
|
|
106 use Bio::SeqIO; # for -file entry loading
|
|
|
107
|
|
|
108 # Note, the following requires HTTP::Request. If the modules are not installed
|
|
|
109 # uncomment the following and use only -filename and don't request swissprotinfo
|
|
|
110 use Bio::DB::EMBL; # for -id entry loading
|
|
|
111
|
|
|
112 use Bio::LiveSeq::IO::Loader 2.0;
|
|
|
113
|
|
|
114 @ISA=qw(Bio::LiveSeq::IO::Loader);
|
|
|
115
|
|
|
116 # This package can in the future host other databases loading subroutines.
|
|
|
117 # e.g. ensembl2hash
|
|
|
118
|
|
|
119 =head2 load
|
|
|
120
|
|
|
121 Title : load
|
|
|
122 Usage : my $filename="../data/M20132";
|
|
|
123 $loader=Bio::LiveSeq::IO::BioPerl->load(-db=>"EMBL", -file=>"$filename");
|
|
|
124 or
|
|
|
125 $loader=Bio::LiveSeq::IO::BioPerl->load(-db=>"EMBL", -id=>"HSANDREC");
|
|
|
126
|
|
|
127 Function: loads an entry with BioPerl from a database into a hash
|
|
|
128 Returns : reference to a new object of class IO::BioPerl holding an entry
|
|
|
129 Errorcode 0
|
|
|
130 Args : an filename containing an EMBL entry OR an ID or ACCESSION code
|
|
|
131
|
|
|
132 =cut
|
|
|
133
|
|
|
134 sub load {
|
|
|
135 my ($thing, %args) = @_;
|
|
|
136 my $class = ref($thing) || $thing;
|
|
|
137 my ($obj,%loader);
|
|
|
138
|
|
|
139 my ($db,$filename,$id)=($args{-db},$args{-file},$args{-id});
|
|
|
140
|
|
|
141 if (defined($db)) {
|
|
|
142 unless ($db eq "EMBL") {
|
|
|
143 carp "Note: only EMBL now supported!";
|
|
|
144 return(0);
|
|
|
145 }
|
|
|
146 } else {
|
|
|
147 $db="EMBL";
|
|
|
148 }
|
|
|
149
|
|
|
150 if (defined($id) && defined($filename)) {
|
|
|
151 carp "You can either specify a -id or a -filename!";
|
|
|
152 return(0);
|
|
|
153 }
|
|
|
154
|
|
|
155 unless (defined($id) || defined($filename)) {
|
|
|
156 carp "You must specify either a -id or a -filename!";
|
|
|
157 return(0);
|
|
|
158 }
|
|
|
159
|
|
|
160 my $hashref;
|
|
|
161 if ($db eq "EMBL") {
|
|
|
162 my $test_transl=0; # change to 0 to avoid comparison of translation
|
|
|
163
|
|
|
164 # these can be changed for future needs
|
|
|
165 my @embl_valid_feature_names=qw(CDS CDS_span exon prim_transcript intron repeat_unit repeat_region mRNA);
|
|
|
166 my @embl_valid_qual_names=qw(gene codon_start db_xref product note number rpt_family transl_table);
|
|
|
167
|
|
|
168 # dunno yet how to implement test_transl again....
|
|
|
169 # probably on a one-on-one basis with each translation?
|
|
|
170 if ($test_transl) {
|
|
|
171 push (@embl_valid_qual_names,"translation"); # needed for test_transl
|
|
|
172 }
|
|
|
173
|
|
|
174 my $seqobj; # bioperl sequence object, to be passed to embl2hash
|
|
|
175
|
|
|
176 if (defined($filename)) {
|
|
|
177 my $stream = Bio::SeqIO->new('-file' => $filename, '-format' => 'EMBL');
|
|
|
178 $seqobj = $stream->next_seq();
|
|
|
179 } else { # i.e. if -id
|
|
|
180 my $embl = new Bio::DB::EMBL;
|
|
|
181 $seqobj = $embl->get_Seq_by_id($id); # EMBL ID or ACC
|
|
|
182 }
|
|
|
183
|
|
|
184 $hashref=&embl2hash($seqobj,\@embl_valid_feature_names,\@embl_valid_qual_names);
|
|
|
185 }
|
|
|
186 unless ($hashref) { return (0); }
|
|
|
187
|
|
|
188 %loader = (db => $db, filename => $filename, id => $id, hash => $hashref);
|
|
|
189 $obj = \%loader;
|
|
|
190 $obj = bless $obj, $class;
|
|
|
191 return $obj;
|
|
|
192 }
|
|
|
193
|
|
|
194 =head2 embl2hash
|
|
|
195
|
|
|
196 Title : embl2hash
|
|
|
197 Function: retrieves with BioPerl an EMBL entry, parses it and creates
|
|
|
198 a hash that contains all the information.
|
|
|
199 Returns : a reference to a hash
|
|
|
200 Errorcode: 0
|
|
|
201 Args : a BioPerl Sequence Object (from file or web fetching)
|
|
|
202 two array references to skip features and qualifiers (for
|
|
|
203 performance)
|
|
|
204 Example: @valid_features=qw(CDS exon prim_transcript mRNA);
|
|
|
205 @valid_qualifiers=qw(gene codon_start db_xref product rpt_family);
|
|
|
206 $hashref=&embl2hash($seqobj,\@valid_features,\@valid_qualifiers);
|
|
|
207
|
|
|
208 =cut
|
|
|
209
|
|
|
210 # arguments: Bioperl $seqobj
|
|
|
211 # to skip features and qualifiers (for performance), two array
|
|
|
212 # references must be passed (this can change into string arguments to
|
|
|
213 # be passed....)
|
|
|
214 # returns: a reference to a hash containing the important features requested
|
|
|
215 sub embl2hash {
|
|
|
216 my $seqobj=$_[0];
|
|
|
217 my %valid_features; my %valid_names;
|
|
|
218 if ($_[1]) {
|
|
|
219 %valid_features = map {$_, 1} @{$_[1]}; # to skip features
|
|
|
220 }
|
|
|
221 if ($_[2]) {
|
|
|
222 %valid_names = map {$_, 1} @{$_[2]}; # to skip qualifiers
|
|
|
223 }
|
|
|
224
|
|
|
225 my $annobj = $seqobj->annotation(); # what's this?
|
|
|
226
|
|
|
227 my $entry_Sequence = lc($seqobj->seq()); # SRS returns lowercase
|
|
|
228
|
|
|
229 my $entry_ID = $seqobj->display_id;
|
|
|
230 my $entry_AccNumber = $seqobj->accession; # or maybe accession_number ?
|
|
|
231 my $secondary_acc; # to fetch the other acc numbers
|
|
|
232 foreach $secondary_acc ($seqobj->get_secondary_accessions) { # not working!
|
|
|
233 $entry_AccNumber .= " $secondary_acc";
|
|
|
234 }
|
|
|
235 my $entry_Molecule = $seqobj->molecule; # this alone returns molec+division
|
|
|
236 my $entry_Division = $seqobj->division;
|
|
|
237 # fixed: now Molecule works in BioPerl, no need for next lines
|
|
|
238 #my @Molecule=split(" ",$entry_Molecule);
|
|
|
239 #my $entry_Division = pop(@Molecule); # only division
|
|
|
240 #$entry_Molecule = join(" ",@Molecule); # only molecule
|
|
|
241 my $entry_Description = $seqobj->desc;
|
|
|
242
|
|
|
243 my $speciesobj = $seqobj->species;
|
|
|
244 my $entry_Organism = $speciesobj->binomial;
|
|
|
245
|
|
|
246 my $entry_SeqLength = $seqobj->length;
|
|
|
247
|
|
|
248 # put into the hash
|
|
|
249 my %entryhash;
|
|
|
250 $entryhash{ID}=$entry_ID;
|
|
|
251 $entryhash{AccNumber}=$entry_AccNumber;
|
|
|
252 $entryhash{Molecule}=$entry_Molecule;
|
|
|
253 $entryhash{Division}=$entry_Division;
|
|
|
254 $entryhash{Description}=$entry_Description;
|
|
|
255 $entryhash{Organism}=$entry_Organism;
|
|
|
256 $entryhash{Sequence}=$entry_Sequence;
|
|
|
257 $entryhash{SeqLength}=$entry_SeqLength;
|
|
|
258
|
|
|
259 my @topfeatures=$seqobj->top_SeqFeatures();
|
|
|
260 # create features array
|
|
|
261 my $featuresnumber= scalar(@topfeatures);
|
|
|
262 $entryhash{FeaturesNumber}=$featuresnumber;
|
|
|
263 my $feature_name;
|
|
|
264 my @feature_qual_names; my @feature_qual_value;
|
|
|
265 my ($feature_qual_name,$feature_qual_number);
|
|
|
266 my @features;
|
|
|
267
|
|
|
268 my ($feat,$qual,$subfeat);
|
|
|
269 my @subfeat;
|
|
|
270 my $i=0;
|
|
|
271 foreach $feat (@topfeatures) {
|
|
|
272 my %feature;
|
|
|
273 $feature_name = $feat->primary_tag;
|
|
|
274 unless ($valid_features{$feature_name}) {
|
|
|
275 #print "skipping $feature_name\n";
|
|
|
276 next;
|
|
|
277 }
|
|
|
278 # works ok with 0.6.2
|
|
|
279 # if ($feature_name eq "CDS_span") { # case of CDS with various exons 0.6.2
|
|
|
280 # $feature_name="CDS"; # 0.6.2
|
|
|
281 my $featlocation=$feat->location; # 0.7
|
|
|
282 if (($feature_name eq "CDS")&&($featlocation->isa('Bio::Location::SplitLocationI'))) { # case of CDS with various exons BioPerl 0.7
|
|
|
283 # @subfeat=$feat->sub_SeqFeature; # 0.6.2
|
|
|
284 @subfeat=$featlocation->sub_Location(); # 0.7
|
|
|
285 my @transcript;
|
|
|
286 foreach $subfeat (@subfeat) {
|
|
|
287 my @range;
|
|
|
288 if ($subfeat->strand == -1) {
|
|
|
289 @range=($subfeat->end,$subfeat->start,$subfeat->strand);
|
|
|
290 } else {
|
|
|
291 @range=($subfeat->start,$subfeat->end,$subfeat->strand);
|
|
|
292 }
|
|
|
293 push (@transcript,\@range);
|
|
|
294 }
|
|
|
295 $feature{range}=\@transcript;
|
|
|
296 } else {
|
|
|
297 my @range;
|
|
|
298 ($feat->strand == -1) ? (@range = ($feat->end, $feat->start, $feat->strand) ) :
|
|
|
299 (@range = ( $feat->start,$feat->end,$feat->strand) );
|
|
|
300 # works ok with 0.6.2
|
|
|
301 if ($feature_name eq "CDS") { # case of single exon CDS (CDS name but not split location)
|
|
|
302 my @transcript=(\@range);
|
|
|
303 $feature{range}=\@transcript;
|
|
|
304 } else { # all other range features
|
|
|
305 $feature{range}=\@range;
|
|
|
306 }
|
|
|
307 }
|
|
|
308 $feature{location}="deprecated";
|
|
|
309
|
|
|
310 $feature{position}=$i;
|
|
|
311 $feature{name}=$feature_name;
|
|
|
312
|
|
|
313 @feature_qual_names= $feat->all_tags();
|
|
|
314 $feature_qual_number= scalar(@feature_qual_names);
|
|
|
315
|
|
|
316 $feature{qual_number}=$feature_qual_number;
|
|
|
317
|
|
|
318 my %feature_qualifiers;
|
|
|
319 for $qual (@feature_qual_names) {
|
|
|
320 $feature_qual_name=$qual;
|
|
|
321 unless ($valid_names{$feature_qual_name}) {
|
|
|
322 next;
|
|
|
323 }
|
|
|
324 @feature_qual_value=$feat->each_tag_value($qual);
|
|
|
325 #print "$qual => @feature_qual_value \n";
|
|
|
326 $feature_qualifiers{$feature_qual_name}=$feature_qual_value[0]; # ?
|
|
|
327 # maybe the whole array should be entered, not just the 1st element?
|
|
|
328 # what could be the other elements? TOCHECK!
|
|
|
329 }
|
|
|
330 $feature{qualifiers}=\%feature_qualifiers;
|
|
|
331 push (@features,\%feature); # array of features
|
|
|
332 $i++;
|
|
|
333 }
|
|
|
334 $entryhash{Features}=\@features; # put this also into the hash
|
|
|
335
|
|
|
336 my @cds; # array just of CDSs
|
|
|
337 for $i (0..$#features) {
|
|
|
338 if ($features[$i]->{'name'} eq "CDS") {
|
|
|
339 push(@cds,$features[$i]);
|
|
|
340 }
|
|
|
341 }
|
|
|
342 $entryhash{CDS}=\@cds; # put this also into the hash
|
|
|
343 return (\%entryhash);
|
|
|
344 }
|
|
|
345
|
|
|
346 =head2 novelaasequence2gene
|
|
|
347
|
|
|
348 Title : novelaasequence2gene
|
|
|
349 Usage : $gene=Bio::LiveSeq::IO::BioPerl->novelaasequence2gene(-aasequence => "MGLAAPTRS*");
|
|
|
350 : $gene=Bio::LiveSeq::IO::BioPerl->novelaasequence2gene(-aasequence => "MGLAAPTRS*",
|
|
|
351 -cusg_data => "58 44 7 29 3 3 480 267 105 143 122 39 144 162 14 59 53 25 233 292 19 113 88 246 28 68 161 231 27 102 128 151 67 60 138 131 48 61 153 19 233 73 150 31 129 38 147 71 138 43 181 81 44 15 255 118 312 392 236 82 20 10 14 141");
|
|
|
352 : $gene=Bio::LiveSeq::IO::BioPerl->novelaasequence2gene(-aasequence => "MGLAAPTRS*",
|
|
|
353 -cusg_data => "58 44 7 29 3 3 480 267 105 143 122 39 144 162 14 59 53 25 233 292 19 113 88 246 28 68 161 231 27 102 128 151 67 60 138 131 48 61 153 19 233 73 150 31 129 38 147 71 138 43 181 81 44 15 255 118 312 392 236 82 20 10 14 141",
|
|
|
354 -translation_table => "2",
|
|
|
355 -gene_name => "tyr-kinase");
|
|
|
356
|
|
|
357 Function: creates LiveSeq objects from a novel amino acid sequence,
|
|
|
358 using codon usage information (loaded from a file) to choose
|
|
|
359 codons according to relative frequencies.
|
|
|
360 If a codon_usage information is not specified,
|
|
|
361 the default is to use Homo sapiens data (taxonomy ID 9606).
|
|
|
362 If a translation_table ID is not specified, it will default to 1
|
|
|
363 (standard code).
|
|
|
364 Returns : reference to a Gene object containing references to LiveSeq objects
|
|
|
365 Errorcode 0
|
|
|
366 Args : string containing an amino acid sequence
|
|
|
367 string (optional) with codon usage data (64 integer numbers)
|
|
|
368 string (optional) specifying a gene_name
|
|
|
369 integer (optional) specifying a translation_table ID
|
|
|
370
|
|
|
371 =cut
|
|
|
372
|
|
|
373 sub novelaasequence2gene {
|
|
|
374 my ($self, %args) = @_;
|
|
|
375 my ($gene_name,$cusg_data,$aasequence,$ttabid)=($args{-gene_name},$args{-cusg_data},$args{-aasequence},$args{-translation_table});
|
|
|
376
|
|
|
377 my @species_codon_usage;
|
|
|
378 unless ($aasequence) {
|
|
|
379 carp "aasequence not given";
|
|
|
380 return (0);
|
|
|
381 }
|
|
|
382 unless ($gene_name) {
|
|
|
383 $gene_name="Novel Unknown";
|
|
|
384 }
|
|
|
385 unless ($ttabid) {
|
|
|
386 $ttabid=1;
|
|
|
387 }
|
|
|
388 unless ($cusg_data) {
|
|
|
389 @species_codon_usage=
|
|
|
390 qw(68664 118404 126679 51100 125600 123646 75667 210903 435317
|
|
|
391 139009 79303 135218 128429 192616 49456 161556 211962 131222
|
|
|
392 162837 213626 69346 140780 182506 219428 76684 189374 173010
|
|
|
393 310626 82647 202329 180955 250410 180001 118798 76398 160764
|
|
|
394 317359 119013 262630 359627 218376 186915 130857 377006 162826
|
|
|
395 113684 317703 441298 287040 245435 174805 133427 134523 108740
|
|
|
396 225633 185619 78463 240138 174021 244236 142435 8187 5913
|
|
|
397 14381); # updated 21Jul2000
|
|
|
398 } else {
|
|
|
399 @species_codon_usage=split(/ /,$cusg_data);
|
|
|
400 }
|
|
|
401
|
|
|
402 my $gene=Bio::LiveSeq::IO::Loader::_common_novelaasequence2gene(\@species_codon_usage,$ttabid,$aasequence,$gene_name);
|
|
|
403 return ($gene);
|
|
|
404 }
|
|
|
405
|
|
|
406 1;
|
