| line |
stmt |
bran |
cond |
sub |
pod |
time |
code |
|
1
|
|
|
|
|
|
|
# |
|
2
|
|
|
|
|
|
|
# bioperl module for Bio::LiveSeq::Transcript |
|
3
|
|
|
|
|
|
|
# |
|
4
|
|
|
|
|
|
|
# Please direct questions and support issues to |
|
5
|
|
|
|
|
|
|
# |
|
6
|
|
|
|
|
|
|
# Cared for by Joseph Insana |
|
7
|
|
|
|
|
|
|
# |
|
8
|
|
|
|
|
|
|
# Copyright Joseph Insana |
|
9
|
|
|
|
|
|
|
# |
|
10
|
|
|
|
|
|
|
# You may distribute this module under the same terms as perl itself |
|
11
|
|
|
|
|
|
|
# |
|
12
|
|
|
|
|
|
|
# POD documentation - main docs before the code |
|
13
|
|
|
|
|
|
|
|
|
14
|
|
|
|
|
|
|
=head1 NAME |
|
15
|
|
|
|
|
|
|
|
|
16
|
|
|
|
|
|
|
Bio::LiveSeq::Transcript - Transcript class for LiveSeq |
|
17
|
|
|
|
|
|
|
|
|
18
|
|
|
|
|
|
|
=head1 SYNOPSIS |
|
19
|
|
|
|
|
|
|
|
|
20
|
|
|
|
|
|
|
# documentation needed |
|
21
|
|
|
|
|
|
|
|
|
22
|
|
|
|
|
|
|
=head1 DESCRIPTION |
|
23
|
|
|
|
|
|
|
|
|
24
|
|
|
|
|
|
|
This stores information about coding sequences (CDS). |
|
25
|
|
|
|
|
|
|
The implementation is that a Transcript object accesses a collection of |
|
26
|
|
|
|
|
|
|
Exon objects, inferring from them the nucleotide structure and sequence. |
|
27
|
|
|
|
|
|
|
|
|
28
|
|
|
|
|
|
|
=head1 AUTHOR - Joseph A.L. Insana |
|
29
|
|
|
|
|
|
|
|
|
30
|
|
|
|
|
|
|
Email: Insana@ebi.ac.uk, jinsana@gmx.net |
|
31
|
|
|
|
|
|
|
|
|
32
|
|
|
|
|
|
|
=head1 APPENDIX |
|
33
|
|
|
|
|
|
|
|
|
34
|
|
|
|
|
|
|
The rest of the documentation details each of the object |
|
35
|
|
|
|
|
|
|
methods. Internal methods are usually preceded with a _ |
|
36
|
|
|
|
|
|
|
|
|
37
|
|
|
|
|
|
|
=cut |
|
38
|
|
|
|
|
|
|
|
|
39
|
|
|
|
|
|
|
# Let the code begin... |
|
40
|
|
|
|
|
|
|
|
|
41
|
|
|
|
|
|
|
package Bio::LiveSeq::Transcript; |
|
42
|
|
|
|
|
|
|
|
|
43
|
2
|
|
|
2
|
|
10
|
use strict; |
|
|
2
|
|
|
|
|
3
|
|
|
|
2
|
|
|
|
|
47
|
|
|
44
|
|
|
|
|
|
|
# use Carp qw(carp cluck); |
|
45
|
2
|
|
|
2
|
|
10
|
use Bio::LiveSeq::Exon; # uses Exon to create new exon in case of deletion |
|
|
2
|
|
|
|
|
3
|
|
|
|
2
|
|
|
|
|
35
|
|
|
46
|
2
|
|
|
2
|
|
7
|
use base qw(Bio::LiveSeq::SeqI); |
|
|
2
|
|
|
|
|
3
|
|
|
|
2
|
|
|
|
|
5226
|
|
|
47
|
|
|
|
|
|
|
|
|
48
|
|
|
|
|
|
|
=head2 new |
|
49
|
|
|
|
|
|
|
|
|
50
|
|
|
|
|
|
|
Title : new |
|
51
|
|
|
|
|
|
|
Usage : $transcript = Bio::LiveSeq::Transcript->new(-exons => \@obj_refs); |
|
52
|
|
|
|
|
|
|
|
|
53
|
|
|
|
|
|
|
Function: generates a new Bio::LiveSeq::Transcript |
|
54
|
|
|
|
|
|
|
Returns : reference to a new object of class Transcript |
|
55
|
|
|
|
|
|
|
Errorcode -1 |
|
56
|
|
|
|
|
|
|
Args : reference to an array of Exon object references |
|
57
|
|
|
|
|
|
|
|
|
58
|
|
|
|
|
|
|
=cut |
|
59
|
|
|
|
|
|
|
|
|
60
|
|
|
|
|
|
|
sub new { |
|
61
|
6
|
|
|
6
|
1
|
20
|
my ($thing, %args) = @_; |
|
62
|
6
|
|
33
|
|
|
30
|
my $class = ref($thing) || $thing; |
|
63
|
6
|
|
|
|
|
14
|
my ($obj,%transcript); |
|
64
|
|
|
|
|
|
|
|
|
65
|
6
|
|
|
|
|
10
|
my @exons=@{$args{-exons}}; |
|
|
6
|
|
|
|
|
21
|
|
|
66
|
|
|
|
|
|
|
|
|
67
|
6
|
|
|
|
|
8
|
$obj = \%transcript; |
|
68
|
6
|
|
|
|
|
14
|
$obj = bless $obj, $class; |
|
69
|
|
|
|
|
|
|
|
|
70
|
6
|
50
|
|
|
|
16
|
unless (@exons) { |
|
71
|
0
|
|
|
|
|
0
|
$obj->warn("$class not initialised because exons array empty"); |
|
72
|
0
|
|
|
|
|
0
|
return(-1); |
|
73
|
|
|
|
|
|
|
} |
|
74
|
|
|
|
|
|
|
|
|
75
|
|
|
|
|
|
|
# now useless, after start and end methods have been overridden here |
|
76
|
6
|
|
|
|
|
14
|
my $firstexon = $exons[0]; |
|
77
|
|
|
|
|
|
|
#my $lastexon = $exons[-1]; |
|
78
|
|
|
|
|
|
|
#my $start = $firstexon->start; |
|
79
|
|
|
|
|
|
|
#my $end = $lastexon->end; |
|
80
|
6
|
|
|
|
|
33
|
my $strand = $firstexon->strand; |
|
81
|
6
|
|
|
|
|
14
|
my $seq = $firstexon->{'seq'}; |
|
82
|
6
|
|
|
|
|
29
|
$obj->alphabet('rna'); |
|
83
|
|
|
|
|
|
|
|
|
84
|
6
|
50
|
|
|
|
26
|
unless (_checkexons(\@exons)) { |
|
85
|
0
|
|
|
|
|
0
|
$obj->warn("$class not initialised because of problems in the exon structure"); |
|
86
|
0
|
|
|
|
|
0
|
return(-1); |
|
87
|
|
|
|
|
|
|
} |
|
88
|
6
|
|
|
|
|
19
|
$obj->{'strand'}=$strand; |
|
89
|
6
|
|
|
|
|
16
|
$obj->{'exons'}=\@exons; |
|
90
|
6
|
|
|
|
|
10
|
$obj->{'seq'}=$seq; |
|
91
|
|
|
|
|
|
|
|
|
92
|
|
|
|
|
|
|
# set Transcript into each Exon |
|
93
|
6
|
|
|
|
|
11
|
my $exon; |
|
94
|
6
|
|
|
|
|
12
|
foreach $exon (@exons) { |
|
95
|
34
|
|
|
|
|
56
|
$exon->{'transcript'}=$obj; |
|
96
|
|
|
|
|
|
|
} |
|
97
|
6
|
|
|
|
|
29
|
return $obj; |
|
98
|
|
|
|
|
|
|
} |
|
99
|
|
|
|
|
|
|
|
|
100
|
|
|
|
|
|
|
|
|
101
|
|
|
|
|
|
|
=head2 all_Exons |
|
102
|
|
|
|
|
|
|
|
|
103
|
|
|
|
|
|
|
Title : all_Exons |
|
104
|
|
|
|
|
|
|
Usage : $transcript_obj->all_Exons() |
|
105
|
|
|
|
|
|
|
Function: returns references to all Exon objects the Transcript is composed of |
|
106
|
|
|
|
|
|
|
Example : foreach $exon ($transcript->all_Exons()) { do_something } |
|
107
|
|
|
|
|
|
|
Returns : array of object references |
|
108
|
|
|
|
|
|
|
Args : none |
|
109
|
|
|
|
|
|
|
|
|
110
|
|
|
|
|
|
|
=cut |
|
111
|
|
|
|
|
|
|
|
|
112
|
|
|
|
|
|
|
sub all_Exons { |
|
113
|
306
|
|
|
306
|
1
|
457
|
my $self=shift; |
|
114
|
306
|
|
|
|
|
523
|
my $exonsref=$self->{'exons'}; |
|
115
|
306
|
|
|
|
|
514
|
my @exons=@{$exonsref}; |
|
|
306
|
|
|
|
|
649
|
|
|
116
|
306
|
|
|
|
|
508
|
my @newexons; |
|
117
|
|
|
|
|
|
|
my $exon; |
|
118
|
306
|
|
|
|
|
607
|
foreach $exon (@exons) { |
|
119
|
1642
|
50
|
|
|
|
2570
|
unless ($exon->obj_valid) { |
|
120
|
0
|
|
|
|
|
0
|
$self->warn("$exon no more valid, start or end label lost, skipping....",1); # ignorable |
|
121
|
|
|
|
|
|
|
} else { |
|
122
|
1642
|
|
|
|
|
2358
|
push(@newexons,$exon); |
|
123
|
|
|
|
|
|
|
} |
|
124
|
|
|
|
|
|
|
} |
|
125
|
306
|
50
|
|
|
|
763
|
if ($#exons != $#newexons) { |
|
126
|
|
|
|
|
|
|
# update exons field |
|
127
|
0
|
|
|
|
|
0
|
$self->{'exons'}=\@newexons; |
|
128
|
|
|
|
|
|
|
} |
|
129
|
306
|
|
|
|
|
729
|
return (@newexons); |
|
130
|
|
|
|
|
|
|
} |
|
131
|
|
|
|
|
|
|
|
|
132
|
|
|
|
|
|
|
=head2 downstream_seq |
|
133
|
|
|
|
|
|
|
|
|
134
|
|
|
|
|
|
|
Title : downstream_seq |
|
135
|
|
|
|
|
|
|
Usage : $transcript_obj->downstream_seq() |
|
136
|
|
|
|
|
|
|
: $transcript_obj->downstream_seq(64) |
|
137
|
|
|
|
|
|
|
Function: returns a string of nucleotides downstream of the end of the |
|
138
|
|
|
|
|
|
|
CDS. If there is some information of the real mRNA, from features in |
|
139
|
|
|
|
|
|
|
an attached Gene object, it will return up to those boundaries. |
|
140
|
|
|
|
|
|
|
Otherwise it will return 1000 nucleotides. |
|
141
|
|
|
|
|
|
|
If an argument is given it will override the default 1000 number |
|
142
|
|
|
|
|
|
|
and return instead /that/ requested number of nucleotides. |
|
143
|
|
|
|
|
|
|
But if a Gene object is attached, this argument will be ignored. |
|
144
|
|
|
|
|
|
|
Returns : string |
|
145
|
|
|
|
|
|
|
Args : an optional integer number of nucleotides to be returned instead of |
|
146
|
|
|
|
|
|
|
the default if no gene attached |
|
147
|
|
|
|
|
|
|
|
|
148
|
|
|
|
|
|
|
=cut |
|
149
|
|
|
|
|
|
|
|
|
150
|
|
|
|
|
|
|
sub downstream_seq { |
|
151
|
23
|
|
|
23
|
1
|
69
|
my ($self,$howmany)=@_; |
|
152
|
23
|
|
|
|
|
44
|
my $str; |
|
153
|
23
|
50
|
|
|
|
57
|
if (defined ($howmany)) { |
|
154
|
0
|
0
|
|
|
|
0
|
unless ($howmany > 0) { |
|
155
|
0
|
|
|
|
|
0
|
$self->throw("No sense in asking less than 1 downstream nucleotides!"); |
|
156
|
|
|
|
|
|
|
} |
|
157
|
|
|
|
|
|
|
} else { |
|
158
|
23
|
100
|
|
|
|
133
|
unless ($self->{'seq'}->alphabet eq 'rna') { # if rna retrieve until the end |
|
159
|
|
|
|
|
|
|
#$str=$DNAobj->labelsubseq($self->end,undef,undef,"unsecuremoderequested"); |
|
160
|
|
|
|
|
|
|
#return(substr($str,1)); # delete first nucleotide that is the last of Transcript |
|
161
|
19
|
50
|
|
|
|
96
|
if ($self->gene) { # if there is Gene object attached fetch relevant info |
|
162
|
19
|
|
|
|
|
85
|
$str=$self->{'seq'}->labelsubseq($self->end,undef,$self->gene->maxtranscript->end); # retrieve from end of this Transcript to end of the maxtranscript |
|
163
|
19
|
|
|
|
|
101
|
$str=substr($str,1); # delete first nucleotide that is the last of Transcript |
|
164
|
19
|
50
|
|
|
|
76
|
if (CORE::length($str) > 0) { |
|
165
|
19
|
|
|
|
|
81
|
return($str); |
|
166
|
|
|
|
|
|
|
} else { # if there was no downstream through the gene's maxtranscript, go the usual way |
|
167
|
0
|
|
|
|
|
0
|
$howmany = 1000; |
|
168
|
|
|
|
|
|
|
} |
|
169
|
|
|
|
|
|
|
} else { |
|
170
|
0
|
|
|
|
|
0
|
$howmany = 1000; |
|
171
|
|
|
|
|
|
|
} |
|
172
|
|
|
|
|
|
|
} |
|
173
|
|
|
|
|
|
|
} |
|
174
|
4
|
|
|
|
|
17
|
my @exons=$self->all_Exons; |
|
175
|
4
|
|
|
|
|
12
|
my $strand=$self->strand(); |
|
176
|
4
|
|
|
|
|
10
|
my $lastexon=$exons[-1]; |
|
177
|
4
|
|
|
|
|
27
|
my $lastexonlength=$lastexon->length; |
|
178
|
|
|
|
|
|
|
# $howmany nucs after end of last exon |
|
179
|
|
|
|
|
|
|
#my $downstream_seq=$lastexon->subseq($lastexonlength+1,undef,$howmany); |
|
180
|
4
|
|
|
|
|
15
|
my $downstream_seq; |
|
181
|
|
|
|
|
|
|
|
|
182
|
4
|
50
|
|
|
|
13
|
if ($howmany) { |
|
183
|
0
|
|
|
|
|
0
|
$downstream_seq=substr($lastexon->labelsubseq($self->end,$howmany,undef,"unsecuremoderequested"),1); |
|
184
|
|
|
|
|
|
|
} else { |
|
185
|
4
|
50
|
|
|
|
11
|
if ($strand == 1) { |
|
186
|
4
|
|
|
|
|
28
|
$downstream_seq=substr($lastexon->labelsubseq($self->end,undef,$self->{'seq'}->end,"unsecuremoderequested"),1); |
|
187
|
|
|
|
|
|
|
} else { |
|
188
|
0
|
|
|
|
|
0
|
$downstream_seq=substr($lastexon->labelsubseq($self->end,undef,$self->{'seq'}->start,"unsecuremoderequested"),1); |
|
189
|
|
|
|
|
|
|
} |
|
190
|
|
|
|
|
|
|
} |
|
191
|
4
|
|
|
|
|
25
|
return $downstream_seq; |
|
192
|
|
|
|
|
|
|
} |
|
193
|
|
|
|
|
|
|
|
|
194
|
|
|
|
|
|
|
=head2 upstream_seq |
|
195
|
|
|
|
|
|
|
|
|
196
|
|
|
|
|
|
|
Title : upstream_seq |
|
197
|
|
|
|
|
|
|
Usage : $transcript_obj->upstream_seq() |
|
198
|
|
|
|
|
|
|
: $transcript_obj->upstream_seq(64) |
|
199
|
|
|
|
|
|
|
Function: just like downstream_seq but returns nucleotides before the ATG |
|
200
|
|
|
|
|
|
|
Note : the default, if no Gene information present and no nucleotides |
|
201
|
|
|
|
|
|
|
number given, is to return up to 400 nucleotides. |
|
202
|
|
|
|
|
|
|
|
|
203
|
|
|
|
|
|
|
=cut |
|
204
|
|
|
|
|
|
|
|
|
205
|
|
|
|
|
|
|
sub upstream_seq { |
|
206
|
1
|
|
|
1
|
1
|
3
|
my ($self,$howmany)=@_; |
|
207
|
1
|
50
|
|
|
|
4
|
if (defined ($howmany)) { |
|
208
|
0
|
0
|
|
|
|
0
|
unless ($howmany > 0) { |
|
209
|
0
|
|
|
|
|
0
|
$self->throw("No sense in asking less than 1 upstream nucleotides!"); |
|
210
|
|
|
|
|
|
|
} |
|
211
|
|
|
|
|
|
|
} else { |
|
212
|
1
|
50
|
|
|
|
5
|
unless ($self->{'seq'}->alphabet eq 'rna') { # if rna retrieve from the start |
|
213
|
1
|
50
|
|
|
|
3
|
if ($self->gene) { # if there is Gene object attached fetch relevant info |
|
214
|
1
|
|
|
|
|
4
|
my $str=$self->{'seq'}->labelsubseq($self->gene->maxtranscript->start,undef,$self->start); # retrieve from start of maxtranscript to start of this Transcript |
|
215
|
1
|
|
|
|
|
5
|
chop $str; # delete last nucleotide that is the A of starting ATG |
|
216
|
1
|
50
|
|
|
|
4
|
if (length($str) > 0) { |
|
217
|
1
|
|
|
|
|
6
|
return($str); |
|
218
|
|
|
|
|
|
|
} else { # if there was no upstream through the gene's maxtranscript, go the usual way |
|
219
|
0
|
|
|
|
|
0
|
$howmany = 400; |
|
220
|
|
|
|
|
|
|
} |
|
221
|
|
|
|
|
|
|
} else { |
|
222
|
0
|
|
|
|
|
0
|
$howmany = 400; |
|
223
|
|
|
|
|
|
|
} |
|
224
|
|
|
|
|
|
|
} |
|
225
|
|
|
|
|
|
|
} |
|
226
|
0
|
|
|
|
|
0
|
my @exons=$self->all_Exons; |
|
227
|
0
|
|
|
|
|
0
|
my $firstexon=$exons[0]; |
|
228
|
|
|
|
|
|
|
|
|
229
|
0
|
|
|
|
|
0
|
my $upstream_seq; |
|
230
|
0
|
|
|
|
|
0
|
my $strand=$self->strand(); |
|
231
|
|
|
|
|
|
|
|
|
232
|
0
|
0
|
|
|
|
0
|
if ($howmany) {# $howmany nucs before begin of first exon |
|
233
|
0
|
|
|
|
|
0
|
my $labelbefore=$firstexon->label(-$howmany,$firstexon->start); |
|
234
|
0
|
0
|
|
|
|
0
|
if ($labelbefore < 1) { |
|
235
|
0
|
0
|
|
|
|
0
|
if ($strand == 1) { |
|
236
|
0
|
|
|
|
|
0
|
$labelbefore=$self->{'seq'}->start; |
|
237
|
|
|
|
|
|
|
} else { |
|
238
|
0
|
|
|
|
|
0
|
$labelbefore=$self->{'seq'}->end; |
|
239
|
|
|
|
|
|
|
} |
|
240
|
|
|
|
|
|
|
} |
|
241
|
0
|
|
|
|
|
0
|
$upstream_seq=$firstexon->labelsubseq($labelbefore,undef,$firstexon->start,"unsecuremoderequested"); |
|
242
|
0
|
|
|
|
|
0
|
chop $upstream_seq; |
|
243
|
|
|
|
|
|
|
} else { |
|
244
|
0
|
0
|
|
|
|
0
|
if ($strand == 1) { |
|
245
|
0
|
|
|
|
|
0
|
$upstream_seq=$firstexon->labelsubseq($self->{'seq'}->start,undef,$self->start,"unsecuremoderequested"); |
|
246
|
0
|
|
|
|
|
0
|
chop $upstream_seq; # delete last nucleotide that is the A of starting ATG |
|
247
|
|
|
|
|
|
|
} else { |
|
248
|
0
|
|
|
|
|
0
|
$upstream_seq=$firstexon->labelsubseq($self->{'seq'}->end,undef,$self->start,"unsecuremoderequested"); |
|
249
|
0
|
|
|
|
|
0
|
chop $upstream_seq; # delete last nucleotide that is the A of starting ATG |
|
250
|
|
|
|
|
|
|
} |
|
251
|
|
|
|
|
|
|
} |
|
252
|
0
|
|
|
|
|
0
|
return $upstream_seq; |
|
253
|
|
|
|
|
|
|
} |
|
254
|
|
|
|
|
|
|
|
|
255
|
|
|
|
|
|
|
# These get redefined here, overriding the SeqI one because they draw their |
|
256
|
|
|
|
|
|
|
# information from the Exons a Transcript is built of |
|
257
|
|
|
|
|
|
|
# optional argument: firstlabel. If not given, it checks coordinate_start |
|
258
|
|
|
|
|
|
|
# This is useful when called by Translation |
|
259
|
|
|
|
|
|
|
# also used by _delete |
|
260
|
|
|
|
|
|
|
sub label { |
|
261
|
23
|
|
|
23
|
1
|
66
|
my ($self,$position,$firstlabel)=@_; |
|
262
|
23
|
50
|
|
|
|
68
|
unless ($position) { # if position = 0 complain ? |
|
263
|
0
|
|
|
|
|
0
|
$self->warn("Position not given or position 0"); |
|
264
|
0
|
|
|
|
|
0
|
return (-1); |
|
265
|
|
|
|
|
|
|
} |
|
266
|
23
|
|
|
|
|
90
|
my ($start,$end,$strand)=($self->start(),$self->end(),$self->strand()); |
|
267
|
23
|
|
|
|
|
53
|
my ($label,@labels,$length,$arraypos); |
|
268
|
23
|
100
|
|
|
|
64
|
unless (defined ($firstlabel)) { |
|
269
|
3
|
|
|
|
|
21
|
$firstlabel=$self->coordinate_start; # this is inside Transcript obj |
|
270
|
|
|
|
|
|
|
} |
|
271
|
23
|
|
|
|
|
75
|
my $coord_pos=$self->_inside_position($firstlabel); |
|
272
|
23
|
|
|
|
|
144
|
$length=$self->length; |
|
273
|
|
|
|
|
|
|
#if ($strand == 1) { |
|
274
|
23
|
100
|
|
|
|
91
|
if ($position < 1) { |
|
275
|
10
|
|
|
|
|
25
|
$position++; # to account for missing of 0 position |
|
276
|
|
|
|
|
|
|
} |
|
277
|
23
|
|
|
|
|
70
|
$arraypos=$position+$coord_pos-2; |
|
278
|
|
|
|
|
|
|
#print "\n=-=-=-=-DEBUG: arraypos $arraypos, pos $position, coordpos: $coord_pos"; |
|
279
|
23
|
50
|
|
|
|
106
|
if ($arraypos < 0) { |
|
|
|
50
|
|
|
|
|
|
|
280
|
0
|
|
|
|
|
0
|
$label=$self->{'seq'}->label($arraypos,$start,$strand); #? |
|
281
|
|
|
|
|
|
|
} elsif ($arraypos >= $length) { |
|
282
|
0
|
|
|
|
|
0
|
$label=$self->{'seq'}->label($arraypos-$length+2,$end,$strand); #? |
|
283
|
|
|
|
|
|
|
} else { # inside the Transcript |
|
284
|
23
|
|
|
|
|
90
|
@labels=$self->all_labels; |
|
285
|
23
|
|
|
|
|
1162
|
$label=$labels[$arraypos]; |
|
286
|
|
|
|
|
|
|
} |
|
287
|
|
|
|
|
|
|
#} |
|
288
|
|
|
|
|
|
|
} |
|
289
|
|
|
|
|
|
|
|
|
290
|
|
|
|
|
|
|
# argument: label |
|
291
|
|
|
|
|
|
|
# returns: position of label according to coord_start |
|
292
|
|
|
|
|
|
|
# errorcode: 0 label not found |
|
293
|
|
|
|
|
|
|
# optional argument: firstlabel. If not given, it checks coordinate_start |
|
294
|
|
|
|
|
|
|
# This is useful when called by Translation |
|
295
|
|
|
|
|
|
|
sub position { |
|
296
|
16
|
|
|
16
|
1
|
57
|
my ($self,$label,$firstlabel)=@_; |
|
297
|
16
|
50
|
|
|
|
100
|
unless ($self->{'seq'}->valid($label)) { |
|
298
|
0
|
|
|
|
|
0
|
$self->warn("label is not valid"); |
|
299
|
0
|
|
|
|
|
0
|
return (0); |
|
300
|
|
|
|
|
|
|
} |
|
301
|
16
|
100
|
|
|
|
56
|
unless (defined ($firstlabel)) { |
|
302
|
1
|
|
|
|
|
4
|
$firstlabel=$self->coordinate_start; # this is inside Transcript obj |
|
303
|
|
|
|
|
|
|
} |
|
304
|
16
|
50
|
|
|
|
66
|
if ($label == $firstlabel) { |
|
305
|
0
|
|
|
|
|
0
|
return (1); |
|
306
|
|
|
|
|
|
|
} |
|
307
|
16
|
|
|
|
|
54
|
my ($start,$end,$strand)=($self->start(),$self->end(),$self->strand()); |
|
308
|
16
|
|
|
|
|
39
|
my ($position,$in_pos,$out_pos,$coord_pos); |
|
309
|
16
|
|
|
|
|
50
|
my $length=$self->length; |
|
310
|
16
|
|
|
|
|
88
|
$coord_pos=$self->_inside_position($firstlabel); |
|
311
|
16
|
50
|
|
|
|
104
|
if ($self->valid($label)) { # if label is inside the Transcript |
|
312
|
16
|
|
|
|
|
82
|
$in_pos=$self->_inside_position($label); |
|
313
|
16
|
|
|
|
|
63
|
$position=$in_pos-$coord_pos+1; |
|
314
|
16
|
50
|
|
|
|
62
|
if ($position <= 0) { |
|
315
|
0
|
|
|
|
|
0
|
return ($position-1); # accounts for the missing of the 0 position |
|
316
|
|
|
|
|
|
|
} |
|
317
|
|
|
|
|
|
|
} else { |
|
318
|
0
|
0
|
|
|
|
0
|
if ($self->follows($end,$label)) { # label after end of transcript |
|
|
|
0
|
|
|
|
|
|
|
319
|
0
|
|
|
|
|
0
|
$out_pos=$self->{'seq'}->position($label,$end,$strand); |
|
320
|
|
|
|
|
|
|
#print "\n+++++++++DEBUG label $label FOLLOWS end $end outpos $out_pos coordpos $coord_pos"; |
|
321
|
0
|
|
|
|
|
0
|
$position=$out_pos+$length-$coord_pos; |
|
322
|
|
|
|
|
|
|
} elsif ($self->follows($label,$start)) { # label before begin of transcript |
|
323
|
|
|
|
|
|
|
#print "\n+++++++++DEBUG label $label BEFORE start $start outpos $out_pos coordpos $coord_pos"; |
|
324
|
0
|
|
|
|
|
0
|
$out_pos=$self->{'seq'}->position($label,$start,$strand); |
|
325
|
0
|
|
|
|
|
0
|
$position=$out_pos-$coord_pos+1; |
|
326
|
|
|
|
|
|
|
} else { # label is in intron (not valid, not after, not before)! |
|
327
|
0
|
|
|
|
|
0
|
$self->warn("Cannot give position of label pointing to intron according to CDS numbering!",1); |
|
328
|
0
|
|
|
|
|
0
|
return (0); |
|
329
|
|
|
|
|
|
|
} |
|
330
|
|
|
|
|
|
|
} |
|
331
|
16
|
|
|
|
|
103
|
return ($position); |
|
332
|
|
|
|
|
|
|
} |
|
333
|
|
|
|
|
|
|
|
|
334
|
|
|
|
|
|
|
sub seq { |
|
335
|
99
|
|
|
99
|
1
|
278
|
my $self=shift; |
|
336
|
99
|
|
|
|
|
207
|
my ($exon,$str); |
|
337
|
99
|
|
|
|
|
375
|
my @exons=$self->all_Exons(); |
|
338
|
99
|
|
|
|
|
195
|
foreach $exon (@exons) { |
|
339
|
531
|
|
|
|
|
1540
|
$str .= $exon->seq(); |
|
340
|
|
|
|
|
|
|
} |
|
341
|
99
|
|
|
|
|
561
|
return $str; |
|
342
|
|
|
|
|
|
|
} |
|
343
|
|
|
|
|
|
|
|
|
344
|
|
|
|
|
|
|
sub length { |
|
345
|
39
|
|
|
39
|
1
|
82
|
my $self=shift; |
|
346
|
39
|
|
|
|
|
102
|
my ($exon,$length); |
|
347
|
39
|
|
|
|
|
124
|
my @exons=$self->all_Exons(); |
|
348
|
39
|
|
|
|
|
89
|
foreach $exon (@exons) { |
|
349
|
211
|
|
|
|
|
542
|
$length += $exon->length(); |
|
350
|
|
|
|
|
|
|
} |
|
351
|
39
|
|
|
|
|
137
|
return $length; |
|
352
|
|
|
|
|
|
|
} |
|
353
|
|
|
|
|
|
|
|
|
354
|
|
|
|
|
|
|
sub all_labels { |
|
355
|
130
|
|
|
130
|
1
|
188
|
my $self=shift; |
|
356
|
130
|
|
|
|
|
190
|
my ($exon,@labels); |
|
357
|
130
|
|
|
|
|
370
|
my @exons=$self->all_Exons(); |
|
358
|
130
|
|
|
|
|
217
|
foreach $exon (@exons) { |
|
359
|
710
|
|
|
|
|
1876
|
push (@labels,$exon->all_labels()); |
|
360
|
|
|
|
|
|
|
} |
|
361
|
130
|
|
|
|
|
23675
|
return @labels; |
|
362
|
|
|
|
|
|
|
} |
|
363
|
|
|
|
|
|
|
|
|
364
|
|
|
|
|
|
|
# redefined here so that it will retrieve effective subseq without introns |
|
365
|
|
|
|
|
|
|
# otherwise it would have retrieved an underlying DNA (possibly with introns) |
|
366
|
|
|
|
|
|
|
# subsequence |
|
367
|
|
|
|
|
|
|
# Drawback: this is really bulky, label->position and then a call to |
|
368
|
|
|
|
|
|
|
# subseq that will do the opposite position-> label |
|
369
|
|
|
|
|
|
|
# |
|
370
|
|
|
|
|
|
|
# one day this can be rewritten as the main one so that the normal subseq |
|
371
|
|
|
|
|
|
|
# will rely on this one and hence avoid this double (useless and lengthy) |
|
372
|
|
|
|
|
|
|
# conversion between labels and positions |
|
373
|
|
|
|
|
|
|
sub old_labelsubseq { |
|
374
|
0
|
|
|
0
|
0
|
0
|
my ($self,$start,$length,$end)=@_; |
|
375
|
0
|
|
|
|
|
0
|
my ($pos1,$pos2); |
|
376
|
0
|
0
|
|
|
|
0
|
if ($start) { |
|
377
|
0
|
0
|
|
|
|
0
|
unless ($self->valid($start)) { |
|
378
|
0
|
|
|
|
|
0
|
$self->warn("Start label not valid"); return (-1); |
|
|
0
|
|
|
|
|
0
|
|
|
379
|
|
|
|
|
|
|
} |
|
380
|
0
|
|
|
|
|
0
|
$pos1=$self->position($start); |
|
381
|
|
|
|
|
|
|
} |
|
382
|
0
|
0
|
|
|
|
0
|
if ($end) { |
|
383
|
0
|
0
|
|
|
|
0
|
if ($end == $start) { |
|
384
|
0
|
|
|
|
|
0
|
$length=1; |
|
385
|
|
|
|
|
|
|
} else { |
|
386
|
0
|
0
|
|
|
|
0
|
unless ($self->valid($end)) { |
|
387
|
0
|
|
|
|
|
0
|
$self->warn("End label not valid"); return (-1); |
|
|
0
|
|
|
|
|
0
|
|
|
388
|
|
|
|
|
|
|
} |
|
389
|
0
|
0
|
|
|
|
0
|
unless ($self->follows($start,$end) == 1) { |
|
390
|
0
|
|
|
|
|
0
|
$self->warn("End label does not follow Start label!"); return (-1); |
|
|
0
|
|
|
|
|
0
|
|
|
391
|
|
|
|
|
|
|
} |
|
392
|
0
|
|
|
|
|
0
|
$pos2=$self->position($end); |
|
393
|
0
|
|
|
|
|
0
|
undef $length; |
|
394
|
|
|
|
|
|
|
} |
|
395
|
|
|
|
|
|
|
} |
|
396
|
0
|
|
|
|
|
0
|
return ($self->subseq($pos1,$pos2,$length)); |
|
397
|
|
|
|
|
|
|
} |
|
398
|
|
|
|
|
|
|
|
|
399
|
|
|
|
|
|
|
# rewritten, eventually |
|
400
|
|
|
|
|
|
|
|
|
401
|
|
|
|
|
|
|
sub labelsubseq { |
|
402
|
12
|
|
|
12
|
1
|
39
|
my ($self,$start,$length,$end,$unsecuremode)=@_; |
|
403
|
12
|
50
|
33
|
|
|
48
|
unless (defined $unsecuremode && |
|
404
|
|
|
|
|
|
|
$unsecuremode eq "unsecuremoderequested") |
|
405
|
|
|
|
|
|
|
{ # to skip security checks (faster) |
|
406
|
12
|
50
|
|
|
|
29
|
if ($start) { |
|
407
|
12
|
50
|
|
|
|
49
|
unless ($self->valid($start)) { |
|
408
|
0
|
|
|
|
|
0
|
$self->warn("Start label not valid"); return (-1); |
|
|
0
|
|
|
|
|
0
|
|
|
409
|
|
|
|
|
|
|
} |
|
410
|
|
|
|
|
|
|
} else { |
|
411
|
0
|
|
|
|
|
0
|
$start=$self->start; |
|
412
|
|
|
|
|
|
|
} |
|
413
|
12
|
100
|
|
|
|
44
|
if ($end) { |
|
414
|
6
|
50
|
|
|
|
21
|
if ($end == $start) { |
|
415
|
0
|
|
|
|
|
0
|
$length=1; |
|
416
|
0
|
|
|
|
|
0
|
undef $end; |
|
417
|
|
|
|
|
|
|
} else { |
|
418
|
6
|
|
|
|
|
14
|
undef $length; # end argument overrides length argument |
|
419
|
6
|
50
|
|
|
|
31
|
unless ($self->valid($end)) { |
|
420
|
0
|
|
|
|
|
0
|
$self->warn("End label not valid"); return (-1); |
|
|
0
|
|
|
|
|
0
|
|
|
421
|
|
|
|
|
|
|
} |
|
422
|
6
|
50
|
|
|
|
61
|
unless ($self->follows($start,$end) == 1) { |
|
423
|
0
|
|
|
|
|
0
|
$self->warn("End label does not follow Start label!"); return (-1); |
|
|
0
|
|
|
|
|
0
|
|
|
424
|
|
|
|
|
|
|
} |
|
425
|
|
|
|
|
|
|
} |
|
426
|
|
|
|
|
|
|
} else { |
|
427
|
6
|
|
|
|
|
32
|
$end=$self->end; |
|
428
|
|
|
|
|
|
|
} |
|
429
|
|
|
|
|
|
|
} |
|
430
|
12
|
|
|
|
|
46
|
my ($seq,$exon,$startexon,$endexon); my @exonlabels; |
|
|
12
|
|
|
|
|
0
|
|
|
431
|
12
|
|
|
|
|
39
|
my @exons=$self->all_Exons; |
|
432
|
|
|
|
|
|
|
EXONCHECK: |
|
433
|
12
|
|
|
|
|
30
|
foreach $exon (@exons) { |
|
434
|
56
|
100
|
100
|
|
|
192
|
if ((!(defined($startexon)))&&($exon->valid($start))) { # checks only if not yet found |
|
435
|
12
|
|
|
|
|
22
|
$startexon=$exon; |
|
436
|
|
|
|
|
|
|
} |
|
437
|
56
|
100
|
|
|
|
143
|
if ($exon->valid($end)) { |
|
438
|
10
|
|
|
|
|
20
|
$endexon=$exon; |
|
439
|
|
|
|
|
|
|
} |
|
440
|
56
|
100
|
100
|
|
|
213
|
if ((!(defined($seq)) && (defined($startexon)))) { # initializes only once |
|
441
|
12
|
100
|
66
|
|
|
65
|
if ((defined($endexon)) && ($endexon eq $startexon)) { # then perfect, we are finished |
|
442
|
5
|
100
|
|
|
|
22
|
if ($length) { |
|
443
|
1
|
|
|
|
|
8
|
$seq = $startexon->labelsubseq($start,$length,undef,"unsecuremoderequested"); |
|
444
|
|
|
|
|
|
|
|
|
445
|
|
|
|
|
|
|
|
|
446
|
1
|
|
|
|
|
5
|
last EXONCHECK; |
|
447
|
|
|
|
|
|
|
} else { |
|
448
|
4
|
|
|
|
|
23
|
$seq = $startexon->labelsubseq($start,undef,$end,"unsecuremoderequested"); |
|
449
|
|
|
|
|
|
|
} |
|
450
|
4
|
|
|
|
|
16
|
last EXONCHECK; |
|
451
|
|
|
|
|
|
|
} else { # get up to the end of the exon |
|
452
|
7
|
|
|
|
|
44
|
$seq = $startexon->labelsubseq($start,undef,undef,"unsecuremoderequested"); |
|
453
|
|
|
|
|
|
|
} |
|
454
|
|
|
|
|
|
|
} |
|
455
|
51
|
100
|
100
|
|
|
177
|
if (($startexon)&&($exon ne $startexon)) { |
|
456
|
7
|
100
|
|
|
|
26
|
if (defined($endexon)) { # we arrived to the last exon |
|
|
|
50
|
|
|
|
|
|
|
457
|
5
|
|
|
|
|
25
|
$seq .= $endexon->labelsubseq(undef,undef,$end,"unsecuremoderequested"); # get from the start of the exon |
|
458
|
5
|
|
|
|
|
24
|
last EXONCHECK; |
|
459
|
|
|
|
|
|
|
|
|
460
|
|
|
|
|
|
|
} elsif (defined($startexon)) { # we are in a whole-exon-in-the-middle case |
|
461
|
2
|
|
|
|
|
18
|
$seq .= $exon->seq; # we add it completely to the seq |
|
462
|
|
|
|
|
|
|
} # else, we still have to reach the start point, exon useless, we move on |
|
463
|
2
|
50
|
|
|
|
10
|
if ($length) { # if length argument specified |
|
464
|
2
|
50
|
33
|
|
|
16
|
if (($seq && (CORE::length($seq) >= $length))) { |
|
465
|
2
|
|
|
|
|
6
|
last EXONCHECK; |
|
466
|
|
|
|
|
|
|
} |
|
467
|
|
|
|
|
|
|
} |
|
468
|
|
|
|
|
|
|
} |
|
469
|
|
|
|
|
|
|
} |
|
470
|
12
|
100
|
|
|
|
32
|
if ($length) { |
|
471
|
6
|
|
|
|
|
54
|
return (substr($seq,0,$length)); |
|
472
|
|
|
|
|
|
|
} else { |
|
473
|
6
|
|
|
|
|
43
|
return ($seq); |
|
474
|
|
|
|
|
|
|
} |
|
475
|
|
|
|
|
|
|
} |
|
476
|
|
|
|
|
|
|
|
|
477
|
|
|
|
|
|
|
|
|
478
|
|
|
|
|
|
|
# argument: label |
|
479
|
|
|
|
|
|
|
# returns: the objref and progressive number of the Exon containing that label |
|
480
|
|
|
|
|
|
|
# errorcode: -1 |
|
481
|
|
|
|
|
|
|
sub in_which_Exon { |
|
482
|
1
|
|
|
1
|
0
|
3
|
my ($self,$label)=@_; |
|
483
|
1
|
|
|
|
|
2
|
my ($count,$exon); |
|
484
|
1
|
|
|
|
|
113
|
my @exons=$self->all_Exons; |
|
485
|
1
|
|
|
|
|
2
|
foreach $exon (@exons) { |
|
486
|
7
|
|
|
|
|
9
|
$count++; # 1st exon is numbered "1" |
|
487
|
7
|
100
|
|
|
|
16
|
if ($exon->valid($label)) { |
|
488
|
1
|
|
|
|
|
5
|
return ($exon,$count) |
|
489
|
|
|
|
|
|
|
} |
|
490
|
|
|
|
|
|
|
} |
|
491
|
0
|
|
|
|
|
0
|
return (-1); # if nothing found |
|
492
|
|
|
|
|
|
|
} |
|
493
|
|
|
|
|
|
|
|
|
494
|
|
|
|
|
|
|
# recoded to exploit the new fast labelsubseq() |
|
495
|
|
|
|
|
|
|
# valid only inside Transcript |
|
496
|
|
|
|
|
|
|
sub subseq { |
|
497
|
0
|
|
|
0
|
1
|
0
|
my ($self,$pos1,$pos2,$length) = @_; |
|
498
|
0
|
|
|
|
|
0
|
my ($str,$startlabel,$endlabel); |
|
499
|
0
|
0
|
|
|
|
0
|
if (defined ($pos1)) { |
|
500
|
0
|
0
|
|
|
|
0
|
if ($pos1 == 0) { # if position = 0 complain |
|
501
|
0
|
|
|
|
|
0
|
$self->warn("Position cannot be 0!"); return (-1); |
|
|
0
|
|
|
|
|
0
|
|
|
502
|
|
|
|
|
|
|
} |
|
503
|
0
|
0
|
0
|
|
|
0
|
if ((defined ($pos2))&&($pos1>$pos2)) { |
|
504
|
0
|
|
|
|
|
0
|
$self->warn("1st position($pos1) cannot be > 2nd position($pos2)!"); return (-1); |
|
|
0
|
|
|
|
|
0
|
|
|
505
|
|
|
|
|
|
|
} |
|
506
|
0
|
|
|
|
|
0
|
$startlabel=$self->label($pos1); |
|
507
|
0
|
0
|
|
|
|
0
|
unless ($self->valid($startlabel)) { |
|
508
|
0
|
|
|
|
|
0
|
$self->warn("Start label not valid"); return (-1); |
|
|
0
|
|
|
|
|
0
|
|
|
509
|
|
|
|
|
|
|
} |
|
510
|
0
|
0
|
|
|
|
0
|
if ($startlabel < 1) { |
|
511
|
0
|
|
|
|
|
0
|
$self->warn("position $pos1 not valid as start of subseq!"); return (-1); |
|
|
0
|
|
|
|
|
0
|
|
|
512
|
|
|
|
|
|
|
} |
|
513
|
|
|
|
|
|
|
} else { |
|
514
|
0
|
|
|
|
|
0
|
$startlabel=$self->start; |
|
515
|
|
|
|
|
|
|
} |
|
516
|
0
|
0
|
|
|
|
0
|
if (defined ($pos2)) { |
|
517
|
0
|
0
|
|
|
|
0
|
if ($pos2 == 0) { # if position = 0 complain |
|
518
|
0
|
|
|
|
|
0
|
$self->warn("Position cannot be 0!"); return (-1); |
|
|
0
|
|
|
|
|
0
|
|
|
519
|
|
|
|
|
|
|
} |
|
520
|
0
|
|
|
|
|
0
|
undef $length; |
|
521
|
0
|
0
|
0
|
|
|
0
|
if ((defined ($pos1))&&($pos1>$pos2)) { |
|
522
|
0
|
|
|
|
|
0
|
$self->warn("1st position($pos1) cannot be > 2nd position($pos2)!"); return (-1); |
|
|
0
|
|
|
|
|
0
|
|
|
523
|
|
|
|
|
|
|
} |
|
524
|
0
|
|
|
|
|
0
|
$endlabel=$self->label($pos2); |
|
525
|
0
|
0
|
|
|
|
0
|
unless ($self->valid($endlabel)) { |
|
526
|
0
|
|
|
|
|
0
|
$self->warn("End label not valid"); return (-1); |
|
|
0
|
|
|
|
|
0
|
|
|
527
|
|
|
|
|
|
|
} |
|
528
|
0
|
0
|
|
|
|
0
|
if ($endlabel < 1) { |
|
529
|
0
|
|
|
|
|
0
|
$self->warn("position $pos2 not valid as end of subseq!"); return (-1); |
|
|
0
|
|
|
|
|
0
|
|
|
530
|
|
|
|
|
|
|
} |
|
531
|
|
|
|
|
|
|
} else { |
|
532
|
0
|
0
|
|
|
|
0
|
unless (defined ($length)) { |
|
533
|
0
|
|
|
|
|
0
|
$endlabel=$self->end; |
|
534
|
|
|
|
|
|
|
} |
|
535
|
|
|
|
|
|
|
} |
|
536
|
0
|
|
|
|
|
0
|
return ($self->labelsubseq($startlabel,$length,$endlabel,"unsecuremoderequested")); |
|
537
|
|
|
|
|
|
|
} |
|
538
|
|
|
|
|
|
|
|
|
539
|
|
|
|
|
|
|
# works only inside the transcript, complains if asked outside |
|
540
|
|
|
|
|
|
|
sub old_subseq { |
|
541
|
0
|
|
|
0
|
0
|
0
|
my ($self,$pos1,$pos2,$length) = @_; |
|
542
|
0
|
|
|
|
|
0
|
my ($str,$startcount,$endcount,$seq,$seqlength); |
|
543
|
0
|
0
|
|
|
|
0
|
if (defined ($length)) { |
|
544
|
0
|
0
|
|
|
|
0
|
if ($length < 1) { |
|
545
|
0
|
|
|
|
|
0
|
$self->warn("No sense asking for a subseq of length < 1"); |
|
546
|
0
|
|
|
|
|
0
|
return (-1); |
|
547
|
|
|
|
|
|
|
} |
|
548
|
|
|
|
|
|
|
} |
|
549
|
0
|
|
|
|
|
0
|
my $firstlabel=$self->coordinate_start; # this is inside Transcript obj |
|
550
|
0
|
|
|
|
|
0
|
my $coord_pos=$self->_inside_position($firstlabel); # TESTME old |
|
551
|
0
|
|
|
|
|
0
|
$seq=$self->seq; |
|
552
|
0
|
|
|
|
|
0
|
$seqlength=CORE::length($seq); |
|
553
|
0
|
0
|
|
|
|
0
|
unless (defined ($pos1)) { |
|
554
|
0
|
|
|
|
|
0
|
$startcount=1+$coord_pos-1; # i.e. coord_pos |
|
555
|
|
|
|
|
|
|
} else { |
|
556
|
0
|
0
|
|
|
|
0
|
if ($pos1 == 0) { # if position = 0 complain |
|
|
|
0
|
|
|
|
|
|
|
557
|
0
|
|
|
|
|
0
|
$self->warn("Position cannot be 0!"); return (-1); |
|
|
0
|
|
|
|
|
0
|
|
|
558
|
|
|
|
|
|
|
} elsif ($pos1 < 0) { |
|
559
|
0
|
|
|
|
|
0
|
$pos1++; |
|
560
|
|
|
|
|
|
|
} |
|
561
|
0
|
0
|
0
|
|
|
0
|
if ((defined ($pos2))&&($pos1>$pos2)) { |
|
562
|
0
|
|
|
|
|
0
|
$self->warn("1st position ($pos1) cannot be > 2nd position ($pos2)!"); |
|
563
|
0
|
|
|
|
|
0
|
return (-1); |
|
564
|
|
|
|
|
|
|
} |
|
565
|
0
|
|
|
|
|
0
|
$startcount=$pos1+$coord_pos-1; |
|
566
|
|
|
|
|
|
|
} |
|
567
|
0
|
0
|
|
|
|
0
|
unless (defined ($pos2)) { |
|
568
|
|
|
|
|
|
|
; |
|
569
|
|
|
|
|
|
|
} else { |
|
570
|
0
|
0
|
|
|
|
0
|
if ($pos2 == 0) { # if position = 0 complain |
|
|
|
0
|
|
|
|
|
|
|
571
|
0
|
|
|
|
|
0
|
$self->warn("Position cannot be 0!"); return (-1); |
|
|
0
|
|
|
|
|
0
|
|
|
572
|
|
|
|
|
|
|
} elsif ($pos2 < 0) { |
|
573
|
0
|
|
|
|
|
0
|
$pos2++; |
|
574
|
|
|
|
|
|
|
} |
|
575
|
0
|
0
|
0
|
|
|
0
|
if ((defined ($pos1))&&($pos1>$pos2)) { |
|
576
|
0
|
|
|
|
|
0
|
$self->warn("1st position ($pos1) cannot be > 2nd position ($pos2)!"); |
|
577
|
0
|
|
|
|
|
0
|
return (-1); |
|
578
|
|
|
|
|
|
|
} |
|
579
|
0
|
|
|
|
|
0
|
$endcount=$pos2+$coord_pos-1; |
|
580
|
0
|
0
|
|
|
|
0
|
if ($endcount > $seqlength) { |
|
581
|
|
|
|
|
|
|
#print "\n###DEBUG###: pos1 $pos1 pos2 $pos2 coordpos $coord_pos endcount $endcount seqln $seqlength\n"; |
|
582
|
0
|
|
|
|
|
0
|
$self->warn("Cannot access end position after the end of Transcript"); |
|
583
|
0
|
|
|
|
|
0
|
return (-1); |
|
584
|
|
|
|
|
|
|
} |
|
585
|
0
|
|
|
|
|
0
|
$length=$endcount-$startcount+1; |
|
586
|
|
|
|
|
|
|
} |
|
587
|
|
|
|
|
|
|
#print "\n###DEBUG pos1 $pos1 pos2 $pos2 start $startcount end $endcount length $length coordpos $coord_pos\n"; |
|
588
|
0
|
|
|
|
|
0
|
my $offset=$startcount-1; |
|
589
|
0
|
0
|
|
|
|
0
|
if ($offset < 0) { |
|
|
|
0
|
|
|
|
|
|
|
590
|
0
|
|
|
|
|
0
|
$self->warn("Cannot access startposition before the beginning of Transcript, returning from start",1); # ignorable |
|
591
|
0
|
|
|
|
|
0
|
return (substr($seq,0,$length)); |
|
592
|
|
|
|
|
|
|
} elsif ($offset >= $seqlength) { |
|
593
|
0
|
|
|
|
|
0
|
$self->warn("Cannot access startposition after the end of Transcript"); |
|
594
|
0
|
|
|
|
|
0
|
return (-1); |
|
595
|
|
|
|
|
|
|
} else { |
|
596
|
0
|
|
|
|
|
0
|
$str=substr($seq,$offset,$length); |
|
597
|
0
|
0
|
|
|
|
0
|
if (CORE::length($str) < $length) { |
|
598
|
0
|
0
|
|
|
|
0
|
$self->warn("Attention, cannot return the length requested ". |
|
599
|
|
|
|
|
|
|
"for subseq",1) if $self->verbose > 0; # ignorable |
|
600
|
|
|
|
|
|
|
} |
|
601
|
0
|
|
|
|
|
0
|
return $str; |
|
602
|
|
|
|
|
|
|
} |
|
603
|
|
|
|
|
|
|
} |
|
604
|
|
|
|
|
|
|
|
|
605
|
|
|
|
|
|
|
# redefined so that it doesn't require other methods (after deletions) to |
|
606
|
|
|
|
|
|
|
# reset it. |
|
607
|
|
|
|
|
|
|
sub start { |
|
608
|
131
|
|
|
131
|
1
|
208
|
my $self = shift; |
|
609
|
131
|
|
|
|
|
207
|
my $exonsref=$self->{'exons'}; |
|
610
|
131
|
|
|
|
|
152
|
my @exons=@{$exonsref}; |
|
|
131
|
|
|
|
|
268
|
|
|
611
|
131
|
|
|
|
|
364
|
return ($exons[0]->start); |
|
612
|
|
|
|
|
|
|
} |
|
613
|
|
|
|
|
|
|
|
|
614
|
|
|
|
|
|
|
sub end { |
|
615
|
140
|
|
|
140
|
1
|
177
|
my $self = shift; |
|
616
|
140
|
|
|
|
|
235
|
my $exonsref=$self->{'exons'}; |
|
617
|
140
|
|
|
|
|
162
|
my @exons=@{$exonsref}; |
|
|
140
|
|
|
|
|
276
|
|
|
618
|
140
|
|
|
|
|
321
|
return ($exons[-1]->end); |
|
619
|
|
|
|
|
|
|
} |
|
620
|
|
|
|
|
|
|
|
|
621
|
|
|
|
|
|
|
|
|
622
|
|
|
|
|
|
|
# internal methods begin here |
|
623
|
|
|
|
|
|
|
|
|
624
|
|
|
|
|
|
|
# returns: position of label in transcript's all_labels |
|
625
|
|
|
|
|
|
|
# with STARTlabel == 1 |
|
626
|
|
|
|
|
|
|
# errorcode 0 -> label not found |
|
627
|
|
|
|
|
|
|
# argument: label |
|
628
|
|
|
|
|
|
|
sub _inside_position { |
|
629
|
55
|
|
|
55
|
|
144
|
my ($self,$label)=@_; |
|
630
|
55
|
|
|
|
|
170
|
my ($start,$end,$strand)=($self->start(),$self->end(),$self->strand()); |
|
631
|
55
|
|
|
|
|
102
|
my ($position,$checkme); |
|
632
|
55
|
|
|
|
|
205
|
my @labels=$self->all_labels; |
|
633
|
55
|
|
|
|
|
465
|
foreach $checkme (@labels) { |
|
634
|
59926
|
|
|
|
|
44583
|
$position++; |
|
635
|
59926
|
100
|
|
|
|
66334
|
if ($label == $checkme) { |
|
636
|
55
|
|
|
|
|
2167
|
return ($position); |
|
637
|
|
|
|
|
|
|
} |
|
638
|
|
|
|
|
|
|
} |
|
639
|
0
|
|
|
|
|
0
|
return (0); |
|
640
|
|
|
|
|
|
|
} |
|
641
|
|
|
|
|
|
|
|
|
642
|
|
|
|
|
|
|
# returns 1 OK or 0 ERROR |
|
643
|
|
|
|
|
|
|
# arguments: reference to array of Exon object references |
|
644
|
|
|
|
|
|
|
sub _checkexons { |
|
645
|
6
|
|
|
6
|
|
11
|
my ($exon,$thisstart); |
|
646
|
6
|
|
|
|
|
11
|
my $self=$exon; |
|
647
|
6
|
|
|
|
|
10
|
my $exonsref=$_[0]; |
|
648
|
6
|
|
|
|
|
8
|
my @exons=@{$exonsref}; |
|
|
6
|
|
|
|
|
14
|
|
|
649
|
|
|
|
|
|
|
|
|
650
|
6
|
|
|
|
|
12
|
my $firstexon = $exons[0]; |
|
651
|
|
|
|
|
|
|
|
|
652
|
6
|
50
|
|
|
|
23
|
unless (ref($firstexon) eq "Bio::LiveSeq::Exon") { |
|
653
|
0
|
|
|
|
|
0
|
$self->warn("Object not of class Exon"); |
|
654
|
0
|
|
|
|
|
0
|
return (0); |
|
655
|
|
|
|
|
|
|
} |
|
656
|
6
|
|
|
|
|
15
|
my $strand = $firstexon->strand; |
|
657
|
|
|
|
|
|
|
|
|
658
|
6
|
|
|
|
|
25
|
my $prevend = $firstexon->end; |
|
659
|
6
|
|
|
|
|
14
|
shift @exons; # skip first one |
|
660
|
6
|
|
|
|
|
15
|
foreach $exon (@exons) { |
|
661
|
28
|
50
|
|
|
|
64
|
unless (ref($exon) eq "Bio::LiveSeq::Exon") { # object class check |
|
662
|
0
|
|
|
|
|
0
|
$self->warn("Object not of class Exon"); |
|
663
|
0
|
|
|
|
|
0
|
return (0); |
|
664
|
|
|
|
|
|
|
} |
|
665
|
28
|
50
|
|
|
|
53
|
if ($exon->strand != $strand) { # strand consistency check |
|
666
|
0
|
|
|
|
|
0
|
$self->warn("Exons' strands not consistent when trying to create Transcript"); |
|
667
|
0
|
|
|
|
|
0
|
return (0); |
|
668
|
|
|
|
|
|
|
} |
|
669
|
28
|
|
|
|
|
57
|
$thisstart = $exon->start; |
|
670
|
28
|
50
|
|
|
|
66
|
unless ($exon->{'seq'}->follows($prevend,$thisstart,$strand)) { |
|
671
|
0
|
|
|
|
|
0
|
$self->warn("Exons not in correct order when trying to create Transcript"); |
|
672
|
0
|
|
|
|
|
0
|
return (0); |
|
673
|
|
|
|
|
|
|
} |
|
674
|
28
|
|
|
|
|
82
|
$prevend = $exon->end; |
|
675
|
|
|
|
|
|
|
} |
|
676
|
6
|
|
|
|
|
20
|
return (1); |
|
677
|
|
|
|
|
|
|
} |
|
678
|
|
|
|
|
|
|
|
|
679
|
|
|
|
|
|
|
=head2 get_Translation |
|
680
|
|
|
|
|
|
|
|
|
681
|
|
|
|
|
|
|
Title : valid |
|
682
|
|
|
|
|
|
|
Usage : $translation = $obj->get_Translation() |
|
683
|
|
|
|
|
|
|
Function: retrieves the reference to the object of class Translation (if any) |
|
684
|
|
|
|
|
|
|
attached to a LiveSeq object |
|
685
|
|
|
|
|
|
|
Returns : object reference |
|
686
|
|
|
|
|
|
|
Args : none |
|
687
|
|
|
|
|
|
|
|
|
688
|
|
|
|
|
|
|
=cut |
|
689
|
|
|
|
|
|
|
|
|
690
|
|
|
|
|
|
|
sub get_Translation { |
|
691
|
21
|
|
|
21
|
1
|
51
|
my $self=shift; |
|
692
|
21
|
|
|
|
|
107
|
return ($self->{'translation'}); # this is set when Translation->new is called |
|
693
|
|
|
|
|
|
|
} |
|
694
|
|
|
|
|
|
|
|
|
695
|
|
|
|
|
|
|
# this checks so that deletion spanning multiple exons is |
|
696
|
|
|
|
|
|
|
# handled accordingly and correctly |
|
697
|
|
|
|
|
|
|
# arguments: begin and end label of a deletion |
|
698
|
|
|
|
|
|
|
# this is called BEFORE any deletion in the chain |
|
699
|
|
|
|
|
|
|
sub _deletecheck { |
|
700
|
0
|
|
|
0
|
|
0
|
my ($self,$startlabel,$endlabel)=@_; |
|
701
|
0
|
|
|
|
|
0
|
my $exonsref=$self->{'exons'}; |
|
702
|
0
|
|
|
|
|
0
|
my @exons=@{$exonsref}; |
|
|
0
|
|
|
|
|
0
|
|
|
703
|
0
|
|
|
|
|
0
|
my ($startexon,$endexon,$exon); |
|
704
|
0
|
|
|
|
|
0
|
$startexon=$endexon=0; |
|
705
|
0
|
|
|
|
|
0
|
foreach $exon (@exons) { |
|
706
|
0
|
0
|
0
|
|
|
0
|
if (($startexon == 0)&&($exon->valid($startlabel))) { |
|
707
|
0
|
|
|
|
|
0
|
$startexon=$exon; # exon containing start of deletion |
|
708
|
|
|
|
|
|
|
} |
|
709
|
0
|
0
|
0
|
|
|
0
|
if (($endexon == 0)&&($exon->valid($endlabel))) { |
|
710
|
0
|
|
|
|
|
0
|
$endexon=$exon; # exon containing end of deletion |
|
711
|
|
|
|
|
|
|
} |
|
712
|
0
|
0
|
0
|
|
|
0
|
if (($startexon)&&($endexon)) { |
|
713
|
0
|
|
|
|
|
0
|
last; # don't check further |
|
714
|
|
|
|
|
|
|
} |
|
715
|
|
|
|
|
|
|
} |
|
716
|
0
|
|
|
|
|
0
|
my $nextend=$self->label(2,$endlabel); # retrieve the next label |
|
717
|
0
|
|
|
|
|
0
|
my $prevstart=$self->label(-1,$startlabel); # retrieve the prev label |
|
718
|
|
|
|
|
|
|
|
|
719
|
0
|
0
|
|
|
|
0
|
if ($startexon eq $endexon) { # intra-exon deletion |
|
720
|
0
|
0
|
0
|
|
|
0
|
if (($startexon->start eq $startlabel) && ($startexon->end eq $endlabel)) { |
|
|
|
0
|
|
|
|
|
|
|
|
|
0
|
|
|
|
|
|
|
721
|
|
|
|
|
|
|
# let's delete the entire exon |
|
722
|
0
|
|
|
|
|
0
|
my @newexons; |
|
723
|
0
|
|
|
|
|
0
|
foreach $exon (@exons) { |
|
724
|
0
|
0
|
|
|
|
0
|
unless ($exon eq $startexon) { |
|
725
|
0
|
|
|
|
|
0
|
push(@newexons,$exon); |
|
726
|
|
|
|
|
|
|
} |
|
727
|
|
|
|
|
|
|
} |
|
728
|
0
|
|
|
|
|
0
|
$self->{'exons'}=\@newexons; |
|
729
|
|
|
|
|
|
|
} elsif ($startexon->start eq $startlabel) { # special cases |
|
730
|
0
|
|
|
|
|
0
|
$startexon->{'start'}=$nextend; # set a new start of exon |
|
731
|
|
|
|
|
|
|
} elsif ($startexon->end eq $endlabel) { |
|
732
|
0
|
|
|
|
|
0
|
$startexon->{'end'}=$prevstart; # set a new end of exon |
|
733
|
|
|
|
|
|
|
} else { |
|
734
|
0
|
|
|
|
|
0
|
return; # no problem |
|
735
|
|
|
|
|
|
|
} |
|
736
|
|
|
|
|
|
|
} else { # two new exons to be created, inter-exons deletion |
|
737
|
0
|
|
|
|
|
0
|
my @newexons; |
|
738
|
|
|
|
|
|
|
my $exonobj; |
|
739
|
0
|
|
|
|
|
0
|
my $dna=$self->{'seq'}; |
|
740
|
0
|
|
|
|
|
0
|
my $strand=$self->strand; |
|
741
|
0
|
|
|
|
|
0
|
my $notmiddle=1; # flag for skipping exons in the middle of deletion |
|
742
|
0
|
|
|
|
|
0
|
foreach $exon (@exons) { |
|
743
|
0
|
0
|
|
|
|
0
|
if ($exon eq $startexon) { |
|
|
|
0
|
|
|
|
|
|
|
744
|
0
|
|
|
|
|
0
|
$exonobj=Bio::LiveSeq::Exon->new('-seq'=>$dna,'-start'=>$exon->start,'-end'=>$prevstart,'-strand'=>$strand); # new partial exon |
|
745
|
0
|
|
|
|
|
0
|
push(@newexons,$exonobj); |
|
746
|
0
|
|
|
|
|
0
|
$notmiddle=0; # now we enter totally deleted exons |
|
747
|
|
|
|
|
|
|
} elsif ($exon eq $endexon) { |
|
748
|
0
|
|
|
|
|
0
|
$exonobj=Bio::LiveSeq::Exon->new('-seq'=>$dna,'-start'=>$nextend,'-end'=>$exon->end,'-strand'=>$strand); # new partial exon |
|
749
|
0
|
|
|
|
|
0
|
push(@newexons,$exonobj); |
|
750
|
0
|
|
|
|
|
0
|
$notmiddle=1; # exiting totally deleted exons |
|
751
|
|
|
|
|
|
|
} else { |
|
752
|
0
|
0
|
|
|
|
0
|
if ($notmiddle) { # if before or after exons with deletion |
|
753
|
0
|
|
|
|
|
0
|
push(@newexons,$exon); |
|
754
|
|
|
|
|
|
|
}# else skip them |
|
755
|
|
|
|
|
|
|
} |
|
756
|
|
|
|
|
|
|
} |
|
757
|
0
|
|
|
|
|
0
|
$self->{'exons'}=\@newexons; |
|
758
|
|
|
|
|
|
|
} |
|
759
|
|
|
|
|
|
|
} |
|
760
|
|
|
|
|
|
|
|
|
761
|
|
|
|
|
|
|
=head2 translation_table |
|
762
|
|
|
|
|
|
|
|
|
763
|
|
|
|
|
|
|
Title : translation_table |
|
764
|
|
|
|
|
|
|
Usage : $name = $obj->translation_table; |
|
765
|
|
|
|
|
|
|
: $name = $obj->translation_table(11); |
|
766
|
|
|
|
|
|
|
Function: Returns or sets the translation_table used for translating the |
|
767
|
|
|
|
|
|
|
transcript. |
|
768
|
|
|
|
|
|
|
If it has never been set, it will return undef. |
|
769
|
|
|
|
|
|
|
Returns : an integer |
|
770
|
|
|
|
|
|
|
|
|
771
|
|
|
|
|
|
|
=cut |
|
772
|
|
|
|
|
|
|
|
|
773
|
|
|
|
|
|
|
sub translation_table { |
|
774
|
46
|
|
|
46
|
1
|
107
|
my ($self,$value) = @_; |
|
775
|
46
|
100
|
|
|
|
110
|
if (defined $value) { |
|
776
|
1
|
|
|
|
|
4
|
$self->{'translation_table'} = $value; |
|
777
|
|
|
|
|
|
|
} |
|
778
|
46
|
100
|
|
|
|
143
|
unless (exists $self->{'translation_table'}) { |
|
779
|
42
|
|
|
|
|
299
|
return; |
|
780
|
|
|
|
|
|
|
} else { |
|
781
|
4
|
|
|
|
|
24
|
return $self->{'translation_table'}; |
|
782
|
|
|
|
|
|
|
} |
|
783
|
|
|
|
|
|
|
} |
|
784
|
|
|
|
|
|
|
|
|
785
|
|
|
|
|
|
|
=head2 frame |
|
786
|
|
|
|
|
|
|
|
|
787
|
|
|
|
|
|
|
Title : frame |
|
788
|
|
|
|
|
|
|
Usage : $frame = $transcript->frame($label); |
|
789
|
|
|
|
|
|
|
Function: Returns the frame of a particular nucleotide. |
|
790
|
|
|
|
|
|
|
Frame can be 0 1 or 2 and means the position in the codon triplet |
|
791
|
|
|
|
|
|
|
of the particulat nucleotide. 0 is the first codon_position. |
|
792
|
|
|
|
|
|
|
Codon_position (1 2 3) is simply frame+1. |
|
793
|
|
|
|
|
|
|
If the label asked for is not inside the Transcript, -1 will be |
|
794
|
|
|
|
|
|
|
returned. |
|
795
|
|
|
|
|
|
|
Args : a label |
|
796
|
|
|
|
|
|
|
Returns : 0 1 or 2 |
|
797
|
|
|
|
|
|
|
Errorcode -1 |
|
798
|
|
|
|
|
|
|
|
|
799
|
|
|
|
|
|
|
=cut |
|
800
|
|
|
|
|
|
|
|
|
801
|
|
|
|
|
|
|
# args: label |
|
802
|
|
|
|
|
|
|
# returns: frame of nucleotide (0 1 2) |
|
803
|
|
|
|
|
|
|
# errorcode: -1 |
|
804
|
|
|
|
|
|
|
sub frame { |
|
805
|
7
|
|
|
7
|
1
|
509
|
my ($self,$inputlabel)=@_; |
|
806
|
7
|
|
|
|
|
23
|
my @labels=$self->all_labels; |
|
807
|
7
|
|
|
|
|
57
|
my ($label,$frame,$count); |
|
808
|
7
|
|
|
|
|
19
|
foreach $label (@labels) { |
|
809
|
9447
|
100
|
|
|
|
9493
|
if ($inputlabel == $label) { |
|
810
|
7
|
|
|
|
|
310
|
return ($count % 3); |
|
811
|
|
|
|
|
|
|
} |
|
812
|
9440
|
|
|
|
|
7600
|
$count++; # 0 1 2 3 4.... |
|
813
|
|
|
|
|
|
|
} |
|
814
|
0
|
|
|
|
|
|
return (-1); # label not found amid Transcript labels |
|
815
|
|
|
|
|
|
|
} |
|
816
|
|
|
|
|
|
|
|
|
817
|
|
|
|
|
|
|
1; |