1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
| while (my @row=$sth->fetchrow_array){
# $row[0] = id, $row[1] = organism, $row[2] = sequence
# séquence $row[2] recherche des sous-séquences $1 (motif $for_e3)
# clé1 $i valeur auto-incrémentée permettant une clé unique
# clé2.0 SOUS_SEQ : valeur de $1
# $hash_seq_for{1}{SOUS_SEQ} = TTGCGTAAGTCGTAGTCAT
# clé2.1 POSITTION position de départ de SOUS_SEQ dans $row[2] $-[1]
# $hash_seq_for{1}{POSITION} = 256
# clé2.2 SCORE valeur : le score observé > $score_threshold
# $hash_seq_for{1}{SCORE} = 75
my %hash_seq_for;
# idem pour motif $rev_e3
my %hash_seq_rev;
#____________________________________________________________________#
# 1 : RECUPERATION DES SOUS-SEQUENCES POUVANT S'APPARIER AUX AMORCES
# $hash_seq_for{ID}{SOUS_SEQ}
# AINSI QUE DE LEUR POSITION
# $hash_seq_for{ID}{POSITION}
#____________________________________________________________________#
# les $num_e3 nucléotides en 3' des 2 amorces s'apparient
# pour le reste des 2 amorces, on regarde si l'appariement
# est au moins de $score_threshold pourcents
# $i permet de créer une clé unique
my $i = 0;
# REGEXP
# $for_length_e5 nucléotides, suivi du motif $for_e3
# récupération des nucléotides précédant le motif $for_e3
while ($row[2] =~ /([A-Z]{0,$for_length_e5})$for_e3/g){
$hash_seq_for{$i}{POSITION} = $-[1];
$hash_seq_for{$i}{SOUS_SEQ} = $1;
$i++;
### for : $-[1]
}
$i = 0;
# motif $rev_e3 suivi de $rev_length_e5 nucléotides
# récupération des nucléotides suivant le motif $for_e3
while ($row[2] =~ /$rev_e3([A-Z]{0,$rev_length_e5})/g){
$hash_seq_rev{$i}{POSITION} = $-[1];
$hash_seq_rev{$i}{SOUS_SEQ} = $1;
$i++;
### rev : $-[1]
}
#____________________________________________________________________#
# 2 : CALCUL DES SCORES DES SOUS-SEQUENCES DE %hash_seq_for
# ON NE GARDE QUE LES SCORES >= $score_threshold
# $hash_seq_for{ID}{SCORE}
#____________________________________________________________________#
# calcul des scores pour %hash_seq_for
my $for_button = 0;
foreach my $i (keys %hash_seq_for){
my $for_score = score_calc($hash_seq_for{$i}{SOUS_SEQ}, \@a_for_e5);
### for score : $for_score
# si le score contient des lettres, c'est un message d'erreur
if ($for_score =~ /[a-z]/){
print "SENS id : $row[0] $for_score\n";
}
# si l'amorce sens s'apparie correctement à la
# séquence aspécifique de la table Fungi2
# on garde la valeur du score obtenu
elsif ($for_score >= $score_threshold){
$hash_seq_for{$i}{SCORE} = $for_score;
$for_button = 1;
}
}
#____________________________________________________________________#
# 3 : CALCUL DES SCORES DES SOUS-SEQUENCES DE %hash_seq_rev
# ON NE GARDE QUE LES SCORES >= $score_threshold
# $hash_seq_rev{ID}{SCORE}
#____________________________________________________________________#
# si au moins une des amorces sens s'apparie correctement à la
# séquence aspécifique de la table Fungi2 on analyse les amorces anti-sens
# si au moins une des sous-seq de %hash_seq_for a un score >= $score_threshold
# on va analyser les scores des sous-seq de %hash_seq_rev
my $rev_button = 0;
if ($for_button == 1){
foreach my $i (keys %hash_seq_rev){
my $rev_score = score_calc($hash_seq_rev{$i}{SOUS_SEQ}, \@a_rev_e5);
### rev score : $rev_score
# si le score contient des lettres, c'est un message d'erreur
if ($rev_score =~ /[a-z]/){
print "SENS id : $row[0] $rev_score\n";
}
elsif ($rev_score >= $score_threshold){
$hash_seq_rev{$i}{SCORE} = $rev_score;
$rev_button = 1;
}
}
}
# si $for_button est nul, cela signifie que %hash_seq_for
# n'a aucun score >= $score_threshold. Inutile d'aller plus loin
# on passe à l'entrée suivante de la DB
else{
next;
}
#____________________________________________________________________#
# 4 : RECHERCHE DES COUPLES
#____________________________________________________________________#
# recherche de tous les couples possibles de $db_seq_for/$db_seq_rev dont
# les 2 scores sont >= $score_threshold et $db_seq_for précède $db_seq_rev
# si $rev_button est nul, cela signifie que %hash_seq_rev
# n'a aucun score >= $score_threshold. Inutile d'aller plus loin
# on passe à l'entrée suivante de la DB
my @i_for = grep{ exists $hash_seq_for{$_}{SCORE}} keys %hash_seq_for;
my @i_rev = grep{ exists $hash_seq_rev{$_}{SCORE}} keys %hash_seq_rev;
}
my $rows = $sth->rows;
print "\n$rows entrées\n";
$sth->finish;
$dbh->disconnect();
sub score_calc {
my $db_seq_e5 = shift;
if ($db_seq_e5 !~ /^[ATCG]*$/){
return "ERREUR le motif $db_seq_e5 récupéré possède des nucléotides dégénérés\n";
}
my @a_db_seq_e5 = split(//, $db_seq_e5);
# array de l'extrémité 3' de
# l'amorce sens ou anti-sens
# même taille que @a_db_seq_e5
my $ref_a_e5 = shift;
if ($#{$ref_a_e5} != $#a_db_seq_e5){
return "ERREUR le motif récupéré a une taille $#a_db_seq_e5 incorrecte ($#{$ref_a_e5} attendus)\n";
}
# score d'identité entre l'extrémité 5' de l'amorce
# et la séquence de la DB
my $score = 0;
for my $i (0..$#a_db_seq_e5){
if (${$ref_a_e5}[$i] eq $a_db_seq_e5[$i]){
$score ++;
}
}
if( $score > 0 ){
$score = $score / ($#a_db_seq_e5 +1) * 100;
}
$score = sprintf ("%.3f", $score);
return $score;
} |
manipulation de hash de hash
Répondre avec citation
Partager