Mercurial > repos > jjohnson > crest

comparison Gene.pm @ 0:acc8d8bfeb9a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author jjohnson
date Wed, 08 Feb 2012 16:59:24 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:acc8d8bfeb9a
1 package Gene;
2 use strict;
3 use Transcript;
4 use Carp;
5 use Data::Dumper;
6
7 # a light weight gene structure is used here
8
9 my @Gene_slots;
10 BEGIN {
11 @Gene_slots = qw(NAME CHR START END STRAND EXONS TRANSCRIPTS TYPE);
12 }
13 use enum @Gene_slots;
14
15 my %attribute = (
16 name => NAME,
17 chr => CHR,
18 start => START,
19 end => END,
20 strand => STRAND,
21 transcripts => TRANSCRIPTS,
22 type => TYPE,
23 );
24
25 #using an array instead of a hash for the node
26 sub _accessor {
27 my $index = shift;
28 return sub {
29 my $self = shift;
30 return undef unless $self;
31 if (@_) {
32 $self->[$index] = shift;
33 }
34 return $self->[$index];
35 };
36 }
37
38 while(my($at, $idx) = each %attribute) {
39 no strict 'refs';
40 *$at = _accessor($idx);
41 }
42
43 sub new {
44 my $class = shift;
45 my $obj = [];
46 $obj->[TRANSCRIPTS] = [];
47 if (@_) {
48 my %arg = @_;
49 $obj->[NAME] = $arg{-NAME} if($arg{-NAME});
50 $obj->[CHR] = $arg{-CHR} if($arg{-CHR});
51 $obj->[START] = $arg{-START} if($arg{-START});
52 $obj->[END] = $arg{-END} if($arg{-END});
53 $obj->[STRAND] = $arg{-STRAND} if($arg{-STRAND});
54 $obj->[TRANSCRIPTS] = $arg{-TRANSCRIPTS} if($arg{-TRANSCRIPTS});
55 }
56 return bless $obj, $class;
57 }
58
59 sub add_transcript {
60 my ($self, $fea) = @_;
61 croak "You must add a Transcript type into a gene"
62 unless ($fea->isa('Transcript'));
63 if($self->[STRAND] && $self->[STRAND] ne $fea->strand) {
64 croak "The transcript has different orientation with the gene";
65 }
66 if($self->[CHR] && $self->[CHR] ne $fea->chr) {
67 croak "The transcript is on different chr with the gene";
68 }
69 # if($self->[TYPE] && $fea->type ne $fea->type) {
70 # croak "The type of the transcript are different from the gene";
71 # }
72 $self->[STRAND] = $fea->strand;
73 $self->[CHR] = $fea->chr;
74 # $self->[TYPE] = $fea->type;
75
76 push @{$self->[TRANSCRIPTS]}, $fea;
77
78 $self->[NAME] = $self->[NAME] ? $self->[NAME] . "," . $fea->name : $fea->name;
79
80
81 #update the start and end of the gene
82 $self->[START] = $fea->start if(!$self->[START] || $self->[START] > $fea->start);
83 $self->[END] = $fea->end if(!$self->[END] || $self->[END] < $fea->end);
84 }
85
86 sub get_start {
87 my ($self, $pos, $ext) = @_;
88 my $rtn = $pos;
89 foreach my $t (@{$self->[TRANSCRIPTS]}) {
90 my $tmp = $t->get_start($pos, $ext);
91 $rtn = $tmp if($tmp < $rtn);
92 }
93 return $rtn;
94 }
95
96 sub get_end {
97 my ($self, $pos, $ext) = @_;
98 my $rtn = $pos;
99 foreach my $t (@{$self->[TRANSCRIPTS]}) {
100 my $tmp = $t->get_end($pos, $ext);
101 $rtn = $tmp if($tmp > $rtn);
102 }
103 return $rtn;
104 }
105
106 sub overlap {
107 my ($self, $fea) = @_;
108
109 if(ref($fea) eq 'ARRAY') {
110 foreach my $t ( @{$self->[TRANSCRIPTS]} ) {
111 return 1 if($t->overlap($fea));
112 }
113 }
114 elsif($fea->isa('Transcript')) {
115 return if($fea->strand && $self->[STRAND] ne $fea->strand );
116 return if($fea->chr && $self->[CHR] ne $fea->chr) ;
117 #return if($fea->type && $self->[TYPE] ne $fea->type);
118 foreach my $e ( @{$fea->exons}) {
119 foreach my $t ( @{$self->[TRANSCRIPTS]} ) {
120 return 1 if($t->overlap($e));
121 }
122 }
123 }
124 else {
125 croak "Not implemented overlap";
126 }
127 return 0;
128 }
129
130 1;