line |
stmt |
bran |
cond |
sub |
pod |
time |
code |
1
|
|
|
|
|
|
|
# module Bio::PopGen::TagHaplotype.pm |
2
|
|
|
|
|
|
|
# |
3
|
|
|
|
|
|
|
# Please direct questions and support issues to |
4
|
|
|
|
|
|
|
# |
5
|
|
|
|
|
|
|
# Cared for by Pedro M. Gomez-Fabre |
6
|
|
|
|
|
|
|
# |
7
|
|
|
|
|
|
|
# Copyright Pedro M. Gomez-Fabre |
8
|
|
|
|
|
|
|
# |
9
|
|
|
|
|
|
|
# You may distribute this module under the same term as perl itself |
10
|
|
|
|
|
|
|
# |
11
|
|
|
|
|
|
|
|
12
|
|
|
|
|
|
|
# POD documentation - main docs before the code |
13
|
|
|
|
|
|
|
|
14
|
|
|
|
|
|
|
=head1 NAME |
15
|
|
|
|
|
|
|
|
16
|
|
|
|
|
|
|
Bio::PopGen::TagHaplotype.pm - Haplotype tag object. |
17
|
|
|
|
|
|
|
|
18
|
|
|
|
|
|
|
=head1 SYNOPSIS |
19
|
|
|
|
|
|
|
|
20
|
|
|
|
|
|
|
use Bio::PopGen::TagHaplotype; |
21
|
|
|
|
|
|
|
|
22
|
|
|
|
|
|
|
my $obj = Bio::PopGen::TagHaplotype -> new($hap); |
23
|
|
|
|
|
|
|
|
24
|
|
|
|
|
|
|
=head1 DESCRIPTION |
25
|
|
|
|
|
|
|
|
26
|
|
|
|
|
|
|
This module take as input a haplotype and try toe get the minimal set |
27
|
|
|
|
|
|
|
of SNP that define the haplotype. This module can be use alone. But |
28
|
|
|
|
|
|
|
due to the tagging haplotype process is exponential one. My suggestion |
29
|
|
|
|
|
|
|
is that before to use this module you pass your data under Select.mp |
30
|
|
|
|
|
|
|
module also on this folder. In any case if, you provide an haplotype |
31
|
|
|
|
|
|
|
the module will try to find the answer to your question. |
32
|
|
|
|
|
|
|
|
33
|
|
|
|
|
|
|
=head1 CONSTRUCTORS |
34
|
|
|
|
|
|
|
|
35
|
|
|
|
|
|
|
my $obj = Bio::PopGen::TagHaplotype -> new($hap); |
36
|
|
|
|
|
|
|
|
37
|
|
|
|
|
|
|
were $hap is the reference to an array of array with the haplotype. |
38
|
|
|
|
|
|
|
|
39
|
|
|
|
|
|
|
$hap= [[0, 0, 0], |
40
|
|
|
|
|
|
|
[1, 0, 0], |
41
|
|
|
|
|
|
|
[0, 1, 1] |
42
|
|
|
|
|
|
|
]; |
43
|
|
|
|
|
|
|
|
44
|
|
|
|
|
|
|
=head1 FEEDBACK |
45
|
|
|
|
|
|
|
|
46
|
|
|
|
|
|
|
=head2 Mailing Lists |
47
|
|
|
|
|
|
|
|
48
|
|
|
|
|
|
|
User feedback is an integral part of the evolution of this and other |
49
|
|
|
|
|
|
|
Bioperl modules. Send your comments and suggestions preferably to |
50
|
|
|
|
|
|
|
the Bioperl mailing list. Your participation is much appreciated. |
51
|
|
|
|
|
|
|
|
52
|
|
|
|
|
|
|
bioperl-l@bioperl.org - General discussion |
53
|
|
|
|
|
|
|
http://bioperl.org/wiki/Mailing_lists - About the mailing lists |
54
|
|
|
|
|
|
|
|
55
|
|
|
|
|
|
|
=head2 Support |
56
|
|
|
|
|
|
|
|
57
|
|
|
|
|
|
|
Please direct usage questions or support issues to the mailing list: |
58
|
|
|
|
|
|
|
|
59
|
|
|
|
|
|
|
I |
60
|
|
|
|
|
|
|
|
61
|
|
|
|
|
|
|
rather than to the module maintainer directly. Many experienced and |
62
|
|
|
|
|
|
|
reponsive experts will be able look at the problem and quickly |
63
|
|
|
|
|
|
|
address it. Please include a thorough description of the problem |
64
|
|
|
|
|
|
|
with code and data examples if at all possible. |
65
|
|
|
|
|
|
|
|
66
|
|
|
|
|
|
|
=head2 Reporting Bugs |
67
|
|
|
|
|
|
|
|
68
|
|
|
|
|
|
|
Report bugs to the Bioperl bug tracking system to help us keep track |
69
|
|
|
|
|
|
|
of the bugs and their resolution. Bug reports can be submitted via |
70
|
|
|
|
|
|
|
the web: |
71
|
|
|
|
|
|
|
|
72
|
|
|
|
|
|
|
https://github.com/bioperl/bioperl-live/issues |
73
|
|
|
|
|
|
|
|
74
|
|
|
|
|
|
|
=head1 AUTHOR - Pedro M. Gomez-Fabre |
75
|
|
|
|
|
|
|
|
76
|
|
|
|
|
|
|
Email pgf18872-at-gsk-dot-com |
77
|
|
|
|
|
|
|
|
78
|
|
|
|
|
|
|
=cut |
79
|
|
|
|
|
|
|
|
80
|
|
|
|
|
|
|
|
81
|
|
|
|
|
|
|
# Let the code begin... |
82
|
|
|
|
|
|
|
|
83
|
|
|
|
|
|
|
package Bio::PopGen::TagHaplotype; |
84
|
1
|
|
|
1
|
|
454
|
use strict; |
|
1
|
|
|
|
|
1
|
|
|
1
|
|
|
|
|
25
|
|
85
|
|
|
|
|
|
|
|
86
|
1
|
|
|
1
|
|
4
|
use Data::Dumper; |
|
1
|
|
|
|
|
2
|
|
|
1
|
|
|
|
|
39
|
|
87
|
1
|
|
|
1
|
|
369
|
use Storable qw(dclone); |
|
1
|
|
|
|
|
2183
|
|
|
1
|
|
|
|
|
51
|
|
88
|
|
|
|
|
|
|
|
89
|
1
|
|
|
1
|
|
5
|
use base qw(Bio::Root::Root); |
|
1
|
|
|
|
|
2
|
|
|
1
|
|
|
|
|
300
|
|
90
|
|
|
|
|
|
|
|
91
|
|
|
|
|
|
|
my $USAGE = <
|
92
|
|
|
|
|
|
|
Usage: |
93
|
|
|
|
|
|
|
Bio::PopGen::TagHaplotype->new(-haplotype_block => \$hapblockref) |
94
|
|
|
|
|
|
|
|
95
|
|
|
|
|
|
|
EOF |
96
|
|
|
|
|
|
|
; |
97
|
|
|
|
|
|
|
|
98
|
|
|
|
|
|
|
=head2 new |
99
|
|
|
|
|
|
|
|
100
|
|
|
|
|
|
|
Title : new |
101
|
|
|
|
|
|
|
Function: constructor of the class. |
102
|
|
|
|
|
|
|
Returns : self hash |
103
|
|
|
|
|
|
|
Args : input haplotype (array of array) |
104
|
|
|
|
|
|
|
Status : public |
105
|
|
|
|
|
|
|
|
106
|
|
|
|
|
|
|
=cut |
107
|
|
|
|
|
|
|
|
108
|
|
|
|
|
|
|
#------------------------ |
109
|
|
|
|
|
|
|
sub new{ |
110
|
|
|
|
|
|
|
#------------------------ |
111
|
1
|
|
|
1
|
1
|
146
|
my ($class, @args) = @_; |
112
|
|
|
|
|
|
|
|
113
|
1
|
|
|
|
|
14
|
my $self = $class->SUPER::new(@args); |
114
|
|
|
|
|
|
|
|
115
|
1
|
|
|
|
|
13
|
my ($haplotype_block) = $self->_rearrange([qw(HAPLOTYPE_BLOCK)],@args); |
116
|
|
|
|
|
|
|
|
117
|
1
|
50
|
|
|
|
6
|
if ($haplotype_block) { |
118
|
1
|
|
|
|
|
5
|
$self->haplotype_block($haplotype_block); |
119
|
|
|
|
|
|
|
} |
120
|
|
|
|
|
|
|
else{ |
121
|
0
|
|
|
|
|
0
|
$self->throw("haplotype has not been supplied\n$USAGE"); |
122
|
|
|
|
|
|
|
} |
123
|
|
|
|
|
|
|
|
124
|
|
|
|
|
|
|
# check that the haplotype block is well formed. |
125
|
1
|
|
|
|
|
6
|
for (my $i=0; $i<$#$haplotype_block+1; $i++){ |
126
|
4
|
50
|
|
|
|
6
|
if ( $#{$haplotype_block->[0]} != |
|
4
|
|
|
|
|
6
|
|
127
|
4
|
|
|
|
|
14
|
$#{$haplotype_block->[$i]} ){ |
128
|
|
|
|
|
|
|
|
129
|
0
|
|
|
|
|
0
|
$self->throw("The haplotype matrix is not well formed (Not squared)"); |
130
|
|
|
|
|
|
|
} |
131
|
|
|
|
|
|
|
} |
132
|
|
|
|
|
|
|
|
133
|
|
|
|
|
|
|
# make the calculation |
134
|
1
|
|
|
|
|
4
|
my $tag_list = _scan_snp( $self ->haplotype_block ); |
135
|
|
|
|
|
|
|
|
136
|
1
|
50
|
|
|
|
4
|
if ($tag_list){ |
137
|
1
|
|
|
|
|
4
|
$self ->tag_list($tag_list); |
138
|
|
|
|
|
|
|
} |
139
|
|
|
|
|
|
|
else { |
140
|
0
|
|
|
|
|
0
|
$self ->tag_list(undef); |
141
|
|
|
|
|
|
|
} |
142
|
|
|
|
|
|
|
|
143
|
1
|
50
|
|
|
|
2
|
if ( defined $self->tag_list){ |
144
|
1
|
|
|
|
|
1
|
$self ->tag_length(scalar @{$self->tag_list}); |
|
1
|
|
|
|
|
2
|
|
145
|
|
|
|
|
|
|
} |
146
|
|
|
|
|
|
|
else { |
147
|
0
|
|
|
|
|
0
|
$self ->tag_length(0); #"NO TAGS FOUND!" |
148
|
|
|
|
|
|
|
} |
149
|
|
|
|
|
|
|
|
150
|
1
|
|
|
|
|
4
|
return $self; |
151
|
|
|
|
|
|
|
} |
152
|
|
|
|
|
|
|
|
153
|
|
|
|
|
|
|
=head2 haplotype_block |
154
|
|
|
|
|
|
|
|
155
|
|
|
|
|
|
|
Title : haplotype_block |
156
|
|
|
|
|
|
|
Usage : my $haplotype_block = $TagHaplotype->haplotype_block(); |
157
|
|
|
|
|
|
|
Function: Get the haplotype block for a haplotype tagging selection |
158
|
|
|
|
|
|
|
Returns : reference of array |
159
|
|
|
|
|
|
|
Args : reference of array with haplotype pattern |
160
|
|
|
|
|
|
|
|
161
|
|
|
|
|
|
|
|
162
|
|
|
|
|
|
|
=cut |
163
|
|
|
|
|
|
|
|
164
|
|
|
|
|
|
|
sub haplotype_block{ |
165
|
2
|
|
|
2
|
1
|
5
|
my ($self) =shift; |
166
|
2
|
100
|
|
|
|
8
|
return $self->{'_haplotype_block'} = shift if @_; |
167
|
1
|
|
|
|
|
5
|
return $self->{'_haplotype_block'}; |
168
|
|
|
|
|
|
|
} |
169
|
|
|
|
|
|
|
|
170
|
|
|
|
|
|
|
|
171
|
|
|
|
|
|
|
=head2 input_block |
172
|
|
|
|
|
|
|
|
173
|
|
|
|
|
|
|
Title : input_block |
174
|
|
|
|
|
|
|
Usage : $obj->input_block() |
175
|
|
|
|
|
|
|
Function: returns haplotype block. By now will produce the same output than |
176
|
|
|
|
|
|
|
$self->haplotype_block. but for compatibility, this method is kept. |
177
|
|
|
|
|
|
|
This method is deprecated. |
178
|
|
|
|
|
|
|
Returns : reference to array of array with the haplotype input value |
179
|
|
|
|
|
|
|
Args : none |
180
|
|
|
|
|
|
|
Status : public |
181
|
|
|
|
|
|
|
|
182
|
|
|
|
|
|
|
=cut |
183
|
|
|
|
|
|
|
|
184
|
|
|
|
|
|
|
#------------------------ |
185
|
|
|
|
|
|
|
sub input_block{ |
186
|
|
|
|
|
|
|
#------------------------ |
187
|
0
|
|
|
0
|
1
|
0
|
my $self = shift; |
188
|
|
|
|
|
|
|
|
189
|
0
|
|
|
|
|
0
|
$self->warn(ref($self). "::input_block - deprecated method. Use haplotype_block() instead."); |
190
|
0
|
|
|
|
|
0
|
return $self->haplotype_block; |
191
|
|
|
|
|
|
|
} |
192
|
|
|
|
|
|
|
|
193
|
|
|
|
|
|
|
=head2 tag_list |
194
|
|
|
|
|
|
|
|
195
|
|
|
|
|
|
|
Title : tag_list |
196
|
|
|
|
|
|
|
Usage : $obj->tag_list() |
197
|
|
|
|
|
|
|
Function: returns the list of SNPs combination that identify the |
198
|
|
|
|
|
|
|
haplotype. All combinations are displayed as arrays |
199
|
|
|
|
|
|
|
Returns : reference to array of array. |
200
|
|
|
|
|
|
|
Args : none |
201
|
|
|
|
|
|
|
Status : public |
202
|
|
|
|
|
|
|
|
203
|
|
|
|
|
|
|
=cut |
204
|
|
|
|
|
|
|
|
205
|
|
|
|
|
|
|
#------------------------ |
206
|
|
|
|
|
|
|
sub tag_list{ |
207
|
|
|
|
|
|
|
#------------------------ |
208
|
4
|
|
|
4
|
1
|
5
|
my ($self) = shift; |
209
|
4
|
100
|
|
|
|
8
|
return $self->{'_tag_list'}= shift if @_; |
210
|
3
|
|
|
|
|
12
|
return $self->{'_tag_list'}; |
211
|
|
|
|
|
|
|
} |
212
|
|
|
|
|
|
|
|
213
|
|
|
|
|
|
|
=head2 tag_length |
214
|
|
|
|
|
|
|
|
215
|
|
|
|
|
|
|
Title : tag_length |
216
|
|
|
|
|
|
|
Usage : $obj->tag_length() |
217
|
|
|
|
|
|
|
Function: returns the length of the tag. |
218
|
|
|
|
|
|
|
Returns : scalar |
219
|
|
|
|
|
|
|
Args : none |
220
|
|
|
|
|
|
|
Status : public |
221
|
|
|
|
|
|
|
|
222
|
|
|
|
|
|
|
=cut |
223
|
|
|
|
|
|
|
|
224
|
|
|
|
|
|
|
#------------------------ |
225
|
|
|
|
|
|
|
sub tag_length{ |
226
|
|
|
|
|
|
|
#------------------------ |
227
|
2
|
|
|
2
|
1
|
4
|
my ($self) =shift; |
228
|
2
|
100
|
|
|
|
4
|
return $self ->{'_tag_length'} = shift if @_; |
229
|
1
|
|
|
|
|
5
|
return $self ->{'_tag_length'}; |
230
|
|
|
|
|
|
|
} |
231
|
|
|
|
|
|
|
|
232
|
|
|
|
|
|
|
=head2 _scan_snp |
233
|
|
|
|
|
|
|
|
234
|
|
|
|
|
|
|
Title : _scan_snp |
235
|
|
|
|
|
|
|
Usage : internal |
236
|
|
|
|
|
|
|
Function: scan sets increasing the length until find a non degenerated |
237
|
|
|
|
|
|
|
pattern. |
238
|
|
|
|
|
|
|
Returns : scalar |
239
|
|
|
|
|
|
|
Args : none |
240
|
|
|
|
|
|
|
Status : private |
241
|
|
|
|
|
|
|
|
242
|
|
|
|
|
|
|
=cut |
243
|
|
|
|
|
|
|
|
244
|
|
|
|
|
|
|
#------------------------ |
245
|
|
|
|
|
|
|
sub _scan_snp{ |
246
|
|
|
|
|
|
|
#------------------------ |
247
|
1
|
|
|
1
|
|
4
|
my ($hap)=@_; |
248
|
|
|
|
|
|
|
|
249
|
1
|
|
|
|
|
2
|
my $hap_length = scalar @{$hap->[0]}; ## store the haplotype length |
|
1
|
|
|
|
|
2
|
|
250
|
|
|
|
|
|
|
|
251
|
1
|
|
|
|
|
5
|
for my $i(1..$hap_length){ |
252
|
|
|
|
|
|
|
|
253
|
2
|
|
|
|
|
10
|
my $list = _gen_comb($hap_length, $i); |
254
|
|
|
|
|
|
|
|
255
|
2
|
|
|
|
|
5
|
my $snp_collection = _scan_combinations($hap, $list); |
256
|
|
|
|
|
|
|
|
257
|
|
|
|
|
|
|
# if there is any element on the collection. |
258
|
|
|
|
|
|
|
# We have reached our goal and |
259
|
|
|
|
|
|
|
# we can stop the calculation. |
260
|
2
|
100
|
|
|
|
5
|
if($#$snp_collection>-1){ |
261
|
1
|
|
|
|
|
3
|
return $snp_collection; |
262
|
|
|
|
|
|
|
} |
263
|
|
|
|
|
|
|
} |
264
|
|
|
|
|
|
|
} |
265
|
|
|
|
|
|
|
|
266
|
|
|
|
|
|
|
=head2 _gen_comb |
267
|
|
|
|
|
|
|
|
268
|
|
|
|
|
|
|
Title : _gen_comb |
269
|
|
|
|
|
|
|
Usage : internal |
270
|
|
|
|
|
|
|
Function: we supply the length of the haplotype and the length of the |
271
|
|
|
|
|
|
|
word we need to find and the functions returns the possible |
272
|
|
|
|
|
|
|
list of combinations. |
273
|
|
|
|
|
|
|
Returns : scalar |
274
|
|
|
|
|
|
|
Args : none |
275
|
|
|
|
|
|
|
Status : private |
276
|
|
|
|
|
|
|
|
277
|
|
|
|
|
|
|
=cut |
278
|
|
|
|
|
|
|
|
279
|
|
|
|
|
|
|
#------------------------ |
280
|
|
|
|
|
|
|
sub _gen_comb{ |
281
|
|
|
|
|
|
|
#------------------------ |
282
|
|
|
|
|
|
|
|
283
|
2
|
|
|
2
|
|
3
|
my ($hap_length,$n) = @_; |
284
|
|
|
|
|
|
|
|
285
|
2
|
|
|
|
|
4
|
my @array = (); # list with all elements we have to combine |
286
|
|
|
|
|
|
|
|
287
|
|
|
|
|
|
|
|
288
|
2
|
|
|
|
|
4
|
for(0..$hap_length-1){ push @array, $_ }; |
|
6
|
|
|
|
|
7
|
|
289
|
|
|
|
|
|
|
|
290
|
|
|
|
|
|
|
# |
291
|
|
|
|
|
|
|
# we need some parameters to create the combination list. |
292
|
|
|
|
|
|
|
# This parameters can be changed if we can modify the list values |
293
|
|
|
|
|
|
|
# |
294
|
|
|
|
|
|
|
|
295
|
2
|
|
|
|
|
3
|
my $m = -1; # this parameter start the calculation at value |
296
|
|
|
|
|
|
|
# m+1 on the recursive cicle. |
297
|
|
|
|
|
|
|
|
298
|
2
|
|
|
|
|
3
|
my $value = []; ## seems to have not too much sense here, but is |
299
|
|
|
|
|
|
|
## needed on the recursion and need to be started |
300
|
|
|
|
|
|
|
## from here |
301
|
2
|
|
|
|
|
2
|
my $list = []; |
302
|
|
|
|
|
|
|
|
303
|
2
|
|
|
|
|
5
|
_generateCombinations ( \@array, \$m, \$n, $value, $list); |
304
|
|
|
|
|
|
|
|
305
|
2
|
|
|
|
|
4
|
return $list; |
306
|
|
|
|
|
|
|
|
307
|
|
|
|
|
|
|
} |
308
|
|
|
|
|
|
|
|
309
|
|
|
|
|
|
|
=head2 _generateCombinations |
310
|
|
|
|
|
|
|
|
311
|
|
|
|
|
|
|
Title : _generateCombinations |
312
|
|
|
|
|
|
|
Usage : internal |
313
|
|
|
|
|
|
|
Function: Recursive function that produce all combinations for a set |
314
|
|
|
|
|
|
|
|
315
|
|
|
|
|
|
|
i.e.: |
316
|
|
|
|
|
|
|
|
317
|
|
|
|
|
|
|
1, 2, 3, 4 |
318
|
|
|
|
|
|
|
|
319
|
|
|
|
|
|
|
and word of B<3> will produce: |
320
|
|
|
|
|
|
|
|
321
|
|
|
|
|
|
|
1, 2, 3 |
322
|
|
|
|
|
|
|
1, 2, 4 |
323
|
|
|
|
|
|
|
1, 3, 4 |
324
|
|
|
|
|
|
|
2, 3, 4 |
325
|
|
|
|
|
|
|
|
326
|
|
|
|
|
|
|
Returns : |
327
|
|
|
|
|
|
|
Args : none |
328
|
|
|
|
|
|
|
Status : private |
329
|
|
|
|
|
|
|
|
330
|
|
|
|
|
|
|
=cut |
331
|
|
|
|
|
|
|
|
332
|
|
|
|
|
|
|
#------------------------ |
333
|
|
|
|
|
|
|
sub _generateCombinations{ |
334
|
|
|
|
|
|
|
#------------------------ |
335
|
5
|
|
|
5
|
|
6
|
my ($rarr, $rm, $rn, $rvalue,$rlist)=@_; |
336
|
|
|
|
|
|
|
|
337
|
5
|
|
|
|
|
8
|
for (my $i = ($$rm+1); $i
|
338
|
9
|
|
|
|
|
13
|
push (my @value2,@$rvalue,$rarr->[$i]); |
339
|
9
|
100
|
|
|
|
14
|
if (scalar @value2<$$rn){ |
340
|
3
|
|
|
|
|
6
|
_generateCombinations($rarr,\$i, $rn, \@value2, $rlist); |
341
|
|
|
|
|
|
|
} |
342
|
9
|
100
|
|
|
|
11
|
if (scalar @value2==$$rn){ |
343
|
6
|
|
|
|
|
7
|
push @$rlist, [@value2]; |
344
|
|
|
|
|
|
|
} |
345
|
9
|
50
|
|
|
|
19
|
if(scalar @value2>$$rn){ |
346
|
0
|
|
|
|
|
0
|
last; |
347
|
|
|
|
|
|
|
} |
348
|
|
|
|
|
|
|
} |
349
|
|
|
|
|
|
|
} |
350
|
|
|
|
|
|
|
|
351
|
|
|
|
|
|
|
# take the list of combinations |
352
|
|
|
|
|
|
|
# i.e.: 1 2 3 |
353
|
|
|
|
|
|
|
# 1 2 4 |
354
|
|
|
|
|
|
|
# 1 3 4 |
355
|
|
|
|
|
|
|
# 2 3 4 |
356
|
|
|
|
|
|
|
# |
357
|
|
|
|
|
|
|
# generate a sub array from the haplotype with the snp tag for the combination |
358
|
|
|
|
|
|
|
# and check all haplotypes for these columns. |
359
|
|
|
|
|
|
|
# if two haplotypes have the same value. we can not define the haplotype |
360
|
|
|
|
|
|
|
# without ambiguity. |
361
|
|
|
|
|
|
|
# Will return a list of valid combinations (SNP Tags) |
362
|
|
|
|
|
|
|
# |
363
|
|
|
|
|
|
|
|
364
|
|
|
|
|
|
|
=head2 _scan_combinations |
365
|
|
|
|
|
|
|
|
366
|
|
|
|
|
|
|
Title : _scan_combinations |
367
|
|
|
|
|
|
|
Usage : internal |
368
|
|
|
|
|
|
|
Function: take the haplotype and a list of possible combination |
369
|
|
|
|
|
|
|
for that length. Generate a subset and scan it to find if |
370
|
|
|
|
|
|
|
the information is enough to define the haplotype set. |
371
|
|
|
|
|
|
|
Returns : |
372
|
|
|
|
|
|
|
Args : none |
373
|
|
|
|
|
|
|
Status : private |
374
|
|
|
|
|
|
|
|
375
|
|
|
|
|
|
|
=cut |
376
|
|
|
|
|
|
|
|
377
|
|
|
|
|
|
|
#------------------------ |
378
|
|
|
|
|
|
|
sub _scan_combinations { |
379
|
|
|
|
|
|
|
#------------------------ |
380
|
|
|
|
|
|
|
|
381
|
2
|
|
|
2
|
|
3
|
my($hap,$list) = @_; |
382
|
|
|
|
|
|
|
|
383
|
2
|
|
|
|
|
2
|
my $valid_combination = undef; |
384
|
|
|
|
|
|
|
|
385
|
|
|
|
|
|
|
# we have to check every snp combinations from the list |
386
|
2
|
|
|
|
|
3
|
for my $i (0..$#$list){ |
387
|
|
|
|
|
|
|
|
388
|
|
|
|
|
|
|
# extract from the big array the one we will use for tag calculations |
389
|
6
|
|
|
|
|
9
|
my $subArray = _get_subArray ($hap, $list->[$i]); |
390
|
|
|
|
|
|
|
|
391
|
6
|
|
|
|
|
10
|
my $degeneration = _deg_test($subArray); |
392
|
|
|
|
|
|
|
|
393
|
6
|
100
|
|
|
|
19
|
if(!$degeneration){ |
394
|
1
|
|
|
|
|
2
|
push @$valid_combination, [@{$list->[$i]}]; |
|
1
|
|
|
|
|
3
|
|
395
|
|
|
|
|
|
|
} |
396
|
|
|
|
|
|
|
} |
397
|
2
|
|
|
|
|
3
|
return $valid_combination; |
398
|
|
|
|
|
|
|
} |
399
|
|
|
|
|
|
|
|
400
|
|
|
|
|
|
|
# return 1 if two arrays are degenerated (same haplotype) |
401
|
|
|
|
|
|
|
#------------------------ |
402
|
|
|
|
|
|
|
sub _deg_test{ |
403
|
|
|
|
|
|
|
#------------------------ |
404
|
|
|
|
|
|
|
|
405
|
6
|
|
|
6
|
|
7
|
my ($hap)= @_; |
406
|
|
|
|
|
|
|
|
407
|
|
|
|
|
|
|
# for every sub array we compare each element with the rest |
408
|
6
|
|
|
|
|
7
|
for my $c1(0..$#$hap){ |
409
|
11
|
|
|
|
|
15
|
for my $c2($c1+1..$#$hap){ |
410
|
21
|
|
|
|
|
23
|
my $degeneration = compare_arrays($hap->[$c1], $hap->[$c2]); |
411
|
21
|
100
|
|
|
|
26
|
if ($degeneration){ |
412
|
|
|
|
|
|
|
# if the two arrays are the same |
413
|
5
|
|
|
|
|
7
|
return 1; |
414
|
|
|
|
|
|
|
} |
415
|
|
|
|
|
|
|
} |
416
|
|
|
|
|
|
|
} |
417
|
|
|
|
|
|
|
} |
418
|
|
|
|
|
|
|
|
419
|
|
|
|
|
|
|
#------------------------ |
420
|
|
|
|
|
|
|
sub _get_subArray { |
421
|
|
|
|
|
|
|
#------------------------ |
422
|
6
|
|
|
6
|
|
6
|
my($hap, $combination) =@_; |
423
|
|
|
|
|
|
|
|
424
|
6
|
|
|
|
|
7
|
my $out = []; # output array to be tested |
425
|
|
|
|
|
|
|
|
426
|
6
|
|
|
|
|
9
|
for my $i (0..$#$hap){ |
427
|
24
|
|
|
|
|
25
|
foreach(@$combination){ |
428
|
36
|
|
|
|
|
31
|
push @{$out->[$i]}, $hap->[$i][$_]; |
|
36
|
|
|
|
|
51
|
|
429
|
|
|
|
|
|
|
} |
430
|
|
|
|
|
|
|
} |
431
|
6
|
|
|
|
|
8
|
return $out; |
432
|
|
|
|
|
|
|
} |
433
|
|
|
|
|
|
|
|
434
|
|
|
|
|
|
|
# |
435
|
|
|
|
|
|
|
# take two arrays and compare their values |
436
|
|
|
|
|
|
|
# Returns : 1 if the two values are the same |
437
|
|
|
|
|
|
|
# 0 if the values are different |
438
|
|
|
|
|
|
|
# |
439
|
|
|
|
|
|
|
|
440
|
|
|
|
|
|
|
#------------------------ |
441
|
|
|
|
|
|
|
sub compare_arrays { |
442
|
|
|
|
|
|
|
#------------------------ |
443
|
21
|
|
|
21
|
0
|
23
|
my ($first, $second) = @_; |
444
|
21
|
50
|
|
|
|
26
|
return 0 unless @$first == @$second; |
445
|
21
|
|
|
|
|
26
|
for (my $i = 0; $i < @$first; $i++) { |
446
|
25
|
100
|
|
|
|
43
|
return 0 if $first->[$i] ne $second->[$i]; |
447
|
|
|
|
|
|
|
} |
448
|
5
|
|
|
|
|
6
|
return 1; |
449
|
|
|
|
|
|
|
} |
450
|
|
|
|
|
|
|
|
451
|
|
|
|
|
|
|
1; |