|
0
|
1 package Bio::EnsEMBL::Compara::Homology;
|
|
|
2
|
|
|
3 use strict;
|
|
|
4 use Bio::EnsEMBL::Utils::Exception qw( deprecate throw warning );
|
|
|
5
|
|
|
6 use base ('Bio::EnsEMBL::Compara::AlignedMemberSet');
|
|
|
7
|
|
|
8 =head1 NAME
|
|
|
9
|
|
|
10 Bio::EnsEMBL::Compara::Homology - Homology between two proteins
|
|
|
11
|
|
|
12 =head1 SYNOPSIS
|
|
|
13
|
|
|
14 use Bio::EnsEMBL::Registry;
|
|
|
15
|
|
|
16 my $homology_adaptor = $reg->get_adaptor("Multi", "compara", "Homology");
|
|
|
17
|
|
|
18 ## Please, refer to Bio::EnsEMBL::Compara::DBSQL::MemberAdaptor
|
|
|
19 ## to see how to get a Member from the database. Also, you can
|
|
|
20 ## find alternative ways to fetch homologies in the POD for the
|
|
|
21 ## Bio::EnsEMBL::Compara::DBSQL::HomologyAdaptor module.
|
|
|
22
|
|
|
23 my $homologies = $homology_adaptor->fetch_all_by_Member($member);
|
|
|
24
|
|
|
25 foreach my $this_homology (@$homologies) {
|
|
|
26 my $homologue_genes = $this_homology->gene_list();
|
|
|
27 print join(" and ", @$homologue_genes), " are ",
|
|
|
28 $this_homology->description, "\n";
|
|
|
29 }
|
|
|
30
|
|
|
31 =head1 AUTHOR
|
|
|
32
|
|
|
33 Ensembl Team
|
|
|
34
|
|
|
35 =head1 COPYRIGHT
|
|
|
36
|
|
|
37 Copyright (c) 1999-2012. Ensembl Team
|
|
|
38
|
|
|
39 You may distribute this module under the same terms as perl itself
|
|
|
40
|
|
|
41 =head1 CONTACT
|
|
|
42
|
|
|
43 This modules is part of the EnsEMBL project (http://www.ensembl.org)
|
|
|
44
|
|
|
45 Questions can be posted to the ensembl-dev mailing list:
|
|
|
46 dev@ensembl.org
|
|
|
47
|
|
|
48 =head1 INHERITANCE
|
|
|
49
|
|
|
50 This class inherits all the methods and attributes from Bio::EnsEMBL::DBSQL::BaseAdaptor
|
|
|
51
|
|
|
52 =head1 APPENDIX
|
|
|
53
|
|
|
54 The rest of the documentation details each of the object methods. Internal methods are usually preceded with a _
|
|
|
55
|
|
|
56 =cut
|
|
|
57
|
|
|
58
|
|
|
59 =head2 subtype
|
|
|
60
|
|
|
61 Arg [1] : string $subtype (optional)
|
|
|
62 Example : $subtype = $homology->subtype();
|
|
|
63 $homology->subtype($subtype);
|
|
|
64 Description: getter/setter of string description of homology subtype.
|
|
|
65 Examples: 'Chordata', 'Euteleostomi', 'Homo sapiens'
|
|
|
66 Returntype : string
|
|
|
67 Exceptions : none
|
|
|
68 Caller : general
|
|
|
69
|
|
|
70 =cut
|
|
|
71
|
|
|
72 sub subtype {
|
|
|
73 my $self = shift;
|
|
|
74 # deprecate("Use taxonomy_level() instead.");
|
|
|
75 return $self->taxonomy_level(@_);
|
|
|
76 }
|
|
|
77
|
|
|
78
|
|
|
79 =head2 taxonomy_level
|
|
|
80
|
|
|
81 Arg [1] : string $taxonomy_level (optional)
|
|
|
82 Example : $taxonomy_level = $homology->taxonomy_level();
|
|
|
83 $homology->taxonomy_level($taxonomy_level);
|
|
|
84 Description: getter/setter of string description of homology taxonomy_level.
|
|
|
85 Examples: 'Chordata', 'Euteleostomi', 'Homo sapiens'
|
|
|
86 Returntype : string
|
|
|
87 Exceptions : none
|
|
|
88 Caller : general
|
|
|
89
|
|
|
90 =cut
|
|
|
91
|
|
|
92 sub taxonomy_level {
|
|
|
93 my $self = shift;
|
|
|
94 $self->{'_subtype'} = shift if(@_);
|
|
|
95 $self->{'_subtype'} = '' unless($self->{'_subtype'});
|
|
|
96 return $self->{'_subtype'};
|
|
|
97 }
|
|
|
98
|
|
|
99
|
|
|
100 =head2 taxonomy_alias
|
|
|
101
|
|
|
102 Arg [1] : string $taxonomy_alias (optional)
|
|
|
103 Example : $taxonomy_alias = $homology->taxonomy_alias();
|
|
|
104 $homology->taxonomy_alias($taxonomy_alias);
|
|
|
105 Description: get string description of homology taxonomy_alias.
|
|
|
106 Examples: 'Chordates', 'Bony vertebrates', 'Homo sapiens'
|
|
|
107 Defaults to taxonomy_level if alias is not in the db
|
|
|
108 Returntype : string
|
|
|
109 Exceptions : none
|
|
|
110 Caller : general
|
|
|
111
|
|
|
112 =cut
|
|
|
113
|
|
|
114 sub taxonomy_alias {
|
|
|
115
|
|
|
116 my $self = shift;
|
|
|
117
|
|
|
118 my $ancestor_node_id = $self->ancestor_node_id;
|
|
|
119 return unless $ancestor_node_id;
|
|
|
120
|
|
|
121 my $ancestor_node = $self->adaptor->db->get_GeneTreeNodeAdaptor->fetch_node_by_node_id($ancestor_node_id);
|
|
|
122 return unless $ancestor_node;
|
|
|
123
|
|
|
124 my $taxon_id = $ancestor_node->get_tagvalue('taxon_id');
|
|
|
125 return unless $taxon_id;
|
|
|
126
|
|
|
127 my $taxon = $self->adaptor->db->get_NCBITaxonAdaptor->fetch_node_by_taxon_id($taxon_id);
|
|
|
128 return unless $taxon;
|
|
|
129
|
|
|
130 return $taxon->ensembl_alias();
|
|
|
131 }
|
|
|
132
|
|
|
133
|
|
|
134 =head2 n
|
|
|
135
|
|
|
136 Arg [1] : float $n (optional)
|
|
|
137 Example : $n = $homology->n();
|
|
|
138 $homology->n(3);
|
|
|
139 Description: getter/setter of number of nonsynonymous positions for the homology.
|
|
|
140 Returntype : float
|
|
|
141 Exceptions : none
|
|
|
142 Caller : general
|
|
|
143
|
|
|
144 =cut
|
|
|
145
|
|
|
146 sub n {
|
|
|
147 my $self = shift;
|
|
|
148 $self->{'_n'} = shift if(@_);
|
|
|
149 return $self->{'_n'};
|
|
|
150 }
|
|
|
151
|
|
|
152
|
|
|
153 =head2 s
|
|
|
154
|
|
|
155 Arg [1] : float $s (optional)
|
|
|
156 Example : $s = $homology->s();
|
|
|
157 $homology->s(4);
|
|
|
158 Description: getter/setter of number of synonymous positions for the homology.
|
|
|
159 Returntype : float
|
|
|
160 Exceptions : none
|
|
|
161 Caller : general
|
|
|
162
|
|
|
163 =cut
|
|
|
164
|
|
|
165 sub s {
|
|
|
166 my $self = shift;
|
|
|
167 $self->{'_s'} = shift if(@_);
|
|
|
168 return $self->{'_s'};
|
|
|
169 }
|
|
|
170
|
|
|
171
|
|
|
172 =head2 lnl
|
|
|
173
|
|
|
174 Arg [1] : float $lnl (optional)
|
|
|
175 Example : $lnl = $homology->lnl();
|
|
|
176 $homology->lnl(-1234.567);
|
|
|
177 Description: getter/setter of number of the negative log likelihood for the dnds homology calculation.
|
|
|
178 Returntype : float
|
|
|
179 Exceptions : none
|
|
|
180 Caller : general
|
|
|
181
|
|
|
182 =cut
|
|
|
183
|
|
|
184 sub lnl {
|
|
|
185 my $self = shift;
|
|
|
186 $self->{'_lnl'} = shift if(@_);
|
|
|
187 return $self->{'_lnl'};
|
|
|
188 }
|
|
|
189
|
|
|
190 =head2 threshold_on_ds
|
|
|
191
|
|
|
192 Arg [1] : float $threshold_ond_ds (optional)
|
|
|
193 Example : $lnl = $homology->threshold_on_ds();
|
|
|
194 $homology->threshold_on_ds(1.01340);
|
|
|
195 Description: getter/setter of the threshold on ds for which the dnds ratio still makes sense.
|
|
|
196 Returntype : float
|
|
|
197 Exceptions : none
|
|
|
198 Caller : general
|
|
|
199
|
|
|
200 =cut
|
|
|
201
|
|
|
202 sub threshold_on_ds {
|
|
|
203 my $self = shift;
|
|
|
204 $self->{'_threshold_on_ds'} = shift if(@_);
|
|
|
205 return $self->{'_threshold_on_ds'};
|
|
|
206 }
|
|
|
207
|
|
|
208 =head2 dn
|
|
|
209
|
|
|
210 Arg [1] : floating $dn (can be undef)
|
|
|
211 Arg [2] : boolean $apply_threshold_on_ds (optional, default = 1)
|
|
|
212 Can be 0 or 1.
|
|
|
213 Example : $homology->dn or $homology->dn(0.1209)
|
|
|
214 if you want to retrieve dn without applying threshold_on_ds, the right call
|
|
|
215 is $homology->dn(undef,0).
|
|
|
216 Description: set/get the non synonymous subtitution rate
|
|
|
217 Returntype : floating
|
|
|
218 Exceptions :
|
|
|
219 Caller :
|
|
|
220
|
|
|
221 =cut
|
|
|
222
|
|
|
223
|
|
|
224 sub dn {
|
|
|
225 my ($self, $dn, $apply_threshold_on_ds) = @_;
|
|
|
226
|
|
|
227 $self->{'_dn'} = $dn if (defined $dn);
|
|
|
228 $apply_threshold_on_ds = 1 unless (defined $apply_threshold_on_ds);
|
|
|
229
|
|
|
230 unless (defined $self->ds(undef, $apply_threshold_on_ds)) {
|
|
|
231 return undef;
|
|
|
232 }
|
|
|
233
|
|
|
234 return $self->{'_dn'};
|
|
|
235 }
|
|
|
236
|
|
|
237 =head2 ds
|
|
|
238
|
|
|
239 Arg [1] : floating $ds (can be undef)
|
|
|
240 Arg [2] : boolean $apply_threshold_on_ds (optional, default = 1)
|
|
|
241 Can be 0 or 1.
|
|
|
242 Example : $homology->ds or $homology->ds(0.9846)
|
|
|
243 if you want to retrieve ds without applying threshold_on_ds, the right call
|
|
|
244 is $homology->dn(undef,0).
|
|
|
245 Description: set/get the synonymous subtitution rate
|
|
|
246 Returntype : floating
|
|
|
247 Exceptions :
|
|
|
248 Caller :
|
|
|
249
|
|
|
250 =cut
|
|
|
251
|
|
|
252
|
|
|
253 sub ds {
|
|
|
254 my ($self, $ds, $apply_threshold_on_ds) = @_;
|
|
|
255
|
|
|
256 $self->{'_ds'} = $ds if (defined $ds);
|
|
|
257 $apply_threshold_on_ds = 1 unless (defined $apply_threshold_on_ds);
|
|
|
258
|
|
|
259 if (defined $self->{'_ds'} &&
|
|
|
260 defined $self->{'_threshold_on_ds'} &&
|
|
|
261 $self->{'_ds'} > $self->{'_threshold_on_ds'}) {
|
|
|
262
|
|
|
263 if ($apply_threshold_on_ds) {
|
|
|
264 return undef;
|
|
|
265 } else {
|
|
|
266 warning("Threshold on ds values is switched off. Be aware that you may obtain saturated ds values that are not to be trusted, neither the dn/ds ratio\n");
|
|
|
267 }
|
|
|
268 }
|
|
|
269
|
|
|
270 return $self->{'_ds'};
|
|
|
271 }
|
|
|
272
|
|
|
273 =head2 dnds_ratio
|
|
|
274
|
|
|
275 Arg [1] : boolean $apply_threshold_on_ds (optional, default = 1)
|
|
|
276 Can be 0 or 1.
|
|
|
277 Example : $homology->dnds_ratio or
|
|
|
278 $homology->dnds_ratio(0) if you want to obtain a result
|
|
|
279 even when the dS is above the threshold on dS.
|
|
|
280 Description: return the ratio of dN/dS
|
|
|
281 Returntype : floating
|
|
|
282 Exceptions :
|
|
|
283 Caller :
|
|
|
284
|
|
|
285 =cut
|
|
|
286
|
|
|
287
|
|
|
288 sub dnds_ratio {
|
|
|
289 my $self = shift;
|
|
|
290 my $apply_threshold_on_ds = shift;
|
|
|
291
|
|
|
292 $apply_threshold_on_ds = 1 unless (defined $apply_threshold_on_ds);
|
|
|
293
|
|
|
294 my $ds = $self->ds(undef, $apply_threshold_on_ds);
|
|
|
295 my $dn = $self->dn(undef, $apply_threshold_on_ds);
|
|
|
296
|
|
|
297 unless (defined $dn &&
|
|
|
298 defined $ds &&
|
|
|
299 $ds !=0) {
|
|
|
300 return undef;
|
|
|
301 }
|
|
|
302
|
|
|
303 unless (defined $self->{'_dnds_ratio'}) {
|
|
|
304 $self->{'_dnds_ratio'} = sprintf("%.5f",$dn/$ds);
|
|
|
305 }
|
|
|
306
|
|
|
307 return $self->{'_dnds_ratio'};
|
|
|
308 }
|
|
|
309
|
|
|
310
|
|
|
311
|
|
|
312 =head2 print_homology
|
|
|
313
|
|
|
314 Example : $homology->print_homology
|
|
|
315 Description: This method prints a short descriptor of the homology
|
|
|
316 USE ONLY FOR DEBUGGING not for data output since the
|
|
|
317 format of this output may change as need dictates.
|
|
|
318
|
|
|
319 =cut
|
|
|
320
|
|
|
321 sub print_homology {
|
|
|
322 my $self = shift;
|
|
|
323
|
|
|
324 printf("Homology %d,%s,%s : ", $self->dbID, $self->description, $self->subtype);
|
|
|
325 foreach my $member (@{$self->gene_list}) {
|
|
|
326 printf("%s(%d)\t", $member->stable_id, $member->dbID);
|
|
|
327 }
|
|
|
328 print("\n");
|
|
|
329 }
|
|
|
330
|
|
|
331
|
|
|
332 =head2 get_all_PeptideAlignFeature
|
|
|
333
|
|
|
334 Description: returns a reference to an empty array as we don't have any
|
|
|
335 PeptideAlignFeatures associated to the homologies
|
|
|
336 Returntype : array ref
|
|
|
337 Exceptions :
|
|
|
338 Caller :
|
|
|
339
|
|
|
340 =cut
|
|
|
341
|
|
|
342 sub get_all_PeptideAlignFeature {
|
|
|
343
|
|
|
344 deprecated("Homologies don't have PeptideAlignFeatures any more. Use DBSQL::PeptideAlignFeatureAdaptor::fetch_all_by_qmember_id_hmember_id() instead");
|
|
|
345 return [];
|
|
|
346 }
|
|
|
347
|
|
|
348
|
|
|
349
|
|
|
350
|
|
|
351 =head2 homology_key
|
|
|
352
|
|
|
353 Example : my $key = $homology->homology_key;
|
|
|
354 Description: returns a string uniquely identifying this homology in world space.
|
|
|
355 uses the gene_stable_ids of the members and orders them by taxon_id
|
|
|
356 and concatonates them together.
|
|
|
357 Returntype : string
|
|
|
358 Exceptions :
|
|
|
359 Caller :
|
|
|
360
|
|
|
361 =cut
|
|
|
362
|
|
|
363 sub homology_key {
|
|
|
364 my $self = shift;
|
|
|
365 return $self->{'_homology_key'} if(defined($self->{'_homology_key'}));
|
|
|
366
|
|
|
367 my @genes = sort {$a->taxon_id <=> $b->taxon_id || $a->stable_id cmp $b->stable_id} @{$self->gene_list};
|
|
|
368 @genes = map ($_->stable_id, @genes);
|
|
|
369
|
|
|
370 my $homology_key = join('_', @genes);
|
|
|
371 return $homology_key;
|
|
|
372 }
|
|
|
373
|
|
|
374
|
|
|
375 =head2 node_id
|
|
|
376
|
|
|
377 Arg [1] : int $node_id (optional)
|
|
|
378 Example : $node_id = $homology->node_id();
|
|
|
379 $homology->subtype($node_id);
|
|
|
380 Description: getter/setter of integer that refer to a node_id in the protein_tree data.
|
|
|
381 Returntype : int
|
|
|
382 Exceptions : none
|
|
|
383 Caller : general
|
|
|
384
|
|
|
385 =cut
|
|
|
386
|
|
|
387 sub node_id {
|
|
|
388 my $self = shift;
|
|
|
389
|
|
|
390 $self->{'_ancestor_node_id'} = shift if(@_);
|
|
|
391 $self->{'_ancestor_node_id'} = '' unless($self->{'_ancestor_node_id'});
|
|
|
392 return $self->{'_ancestor_node_id'};
|
|
|
393
|
|
|
394 }
|
|
|
395
|
|
|
396 =head2 ancestor_node_id
|
|
|
397
|
|
|
398 Arg [1] : int $ancestor_node_id (optional)
|
|
|
399 Example : $ancestor_node_id = $homology->ancestor_node_id();
|
|
|
400 $homology->subtype($ancestor_node_id);
|
|
|
401 Description: getter/setter of integer that refer to the ancestor_node_id in the protein_tree data.
|
|
|
402 Returntype : int
|
|
|
403 Exceptions : none
|
|
|
404 Caller : general
|
|
|
405
|
|
|
406 =cut
|
|
|
407
|
|
|
408 sub ancestor_node_id {
|
|
|
409 my $self = shift;
|
|
|
410
|
|
|
411 $self->{'_ancestor_node_id'} = shift if(@_);
|
|
|
412 $self->{'_ancestor_node_id'} = '' unless($self->{'_ancestor_node_id'});
|
|
|
413 return $self->{'_ancestor_node_id'};
|
|
|
414
|
|
|
415 }
|
|
|
416
|
|
|
417
|
|
|
418 =head2 tree_node_id
|
|
|
419
|
|
|
420 Arg [1] : int $tree_node_id (optional)
|
|
|
421 Example : $tree_node_id = $homology->tree_node_id();
|
|
|
422 $homology->subtype($tree_node_id);
|
|
|
423 Description: getter/setter of integer that refer to the tree node_id in the protein_tree data.
|
|
|
424 Returntype : int
|
|
|
425 Exceptions : none
|
|
|
426 Caller : general
|
|
|
427
|
|
|
428 =cut
|
|
|
429
|
|
|
430 sub tree_node_id {
|
|
|
431 my $self = shift;
|
|
|
432
|
|
|
433 $self->{'_tree_node_id'} = shift if(@_);
|
|
|
434 $self->{'_tree_node_id'} = '' unless($self->{'_tree_node_id'});
|
|
|
435 return $self->{'_tree_node_id'};
|
|
|
436
|
|
|
437 }
|
|
|
438
|
|
|
439
|
|
|
440 1;
|
|
|
441
|
