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
180
181
182
183
184
185
186
187
| #!/usr/local/bin/perl
use strict;
use warnings;
use FileHandle;
use Bio::SeqIO;
# clé : id valeurs : array contenant la séquence
my %h_seq;
# position dans l'alignement du 1ier nucléotide de la séquence
my %h_start;
# position dans l'alignement du dernier nucléotide de la séquence
my %h_end;
# liste des organismes à analyser
my %h_orga = (
'Escherichia_coli' => 1,
);
#------------------------------------------------#
# récupération des sequences alignées #
#------------------------------------------------#
my $file_in = 'P:/Theorie/PCR_Bact_Hybridation/sequences/Test.fsa';
my $in = Bio::SeqIO->new(-file => => $file_in , -format => 'fasta');
while ( my $seq = $in->next_seq() ) {
my $s = $seq->seq;
$s =~ s/[-~]/\./g;
if($s =~ /^\.*([a-z])[.a-z]+?([a-z])\.*$/i){
# position dans l'alignement du 1ier nucléotide de la séquence
$h_start{$seq->primary_id} = $-[1];
# position dans l'alignement du dernier nucléotide de la séquence
$h_end{$seq->primary_id} = $-[2];
}
# clé : id valeurs : array contenant la séquence
@{$h_seq{$seq->primary_id}} = split(//, $s);
}
my $window = 5;
#------------------------------------------------#
# récupération des motifs par séquences #
#------------------------------------------------#
# clé : organisme valeurs array contenant les index de fin et de début des motifs de $window nucléotides de cette séquence
my %h_pattern;
foreach my $k (keys %h_seq){
if($h_orga{$k}){
# $h_start{$k} = position du premier nucléotide dans la séquence alignée
# $h_end{$k} = position du dernier nucléotide dans la séquence alignée
# position dans la séquence
my $index = $h_start{$k};
# motif de $window nucléotide
my $pattern= '';
# index de départ de ce motif dans la séquence alignée
my $pattern_start;
# index de fin de ce motif dans la séquence alignée
my $pattern_end;
# découpe en tronçons de $window nucléotides
while ( $index < $h_end{$k}+1) {
if(!$pattern){
$pattern_start = $index;
}
if((length($pattern) < $window) && ($h_seq{$k}[$index] !~ /\./)){
$pattern .= $h_seq{$k}[$index];
}
if(length($pattern)==$window){
$pattern_end = $index;
my $score_rel = &pattern_analyze (\%h_seq, $k, $pattern, $pattern_start, $pattern_end);
print "$pattern, $pattern_start, $pattern_end, $score_rel\n";
$pattern = '';
# on recule l'index
$index = $pattern_start;
}
$index++;
}
}
}
#------------------------------------------------#
# analyse des motifs par séquences #
#------------------------------------------------#
# sub pattern_analyze
#----------------------------
# Entrée : - une référence vers un hash contenant les séquences
# ------- - l'organisme de référence de ce hash (à comparer aux autres organismes de ce hash)
# - le motif à comparer
# - son index de départ et de fin dans l'alignement
# Retour : - le pourcentage de nucléotide différent de la séquence de référence par rapport aux autres
sub pattern_analyze{
my ($ref_seq, $orga_ref, $pattern, $pattern_start, $pattern_end) = @_;
my $seq_ref;
my @a_seq_compare;
my $ref_seq_compare = \@a_seq_compare;
foreach my $id (keys %{$ref_seq}){
if($id =~ $orga_ref){
$seq_ref = '';
for(my $l = $pattern_start; $l<$pattern_end; $l++){
$seq_ref .= $h_seq{$orga_ref}[$l];
}
}
else{
my $seq = '';
for(my $l = $pattern_start; $l<=$pattern_end; $l++){
$seq .= $h_seq{$id}[$l];
}
push(@a_seq_compare, $seq);
}
}
my $score_rel = &CompareSequencesScore($seq_ref , $ref_seq_compare);
}
# sub CompareSequencesScore
#----------------------------
# Entrée : - une séquence de référence
# ------- - une référence vers une liste de séquence à comparer
# => toutes ces séquences doivent provenir d'un même alignement
# et donc avoir la même longueur
# Retour : - le pourcentage de nucléotide différent de la séquence de référence par rapport aux autres
sub CompareSequencesScore
{
my $Seq_Ref = $_[0];
my $Ref_A_SeqCompare = $_[1];
# Test que toutes les séquences de l'alignement ait la même longueur
$Seq_Ref = uc($Seq_Ref);
foreach (@{$Ref_A_SeqCompare}){$_=uc($_);}
# Toutes les séquences mises en majuscules
my @A_Long;
my %H_Long;
map{push(@A_Long, length($_))}@{$Ref_A_SeqCompare};
map{$H_Long{$_}=1;}@A_Long;
my $NBrDiffTaille = keys(%H_Long);
if($NBrDiffTaille!=1){print "ERREUR LES SEQUENCES DE L'ALIGNEMENT N'ONT PAS TOUTES LA MEME TAILLE!!!!!\n\n";}
# création du tableau des séquences à comparer
my @A_Tableau;
for (my $a=0; $a<@{$Ref_A_SeqCompare}; $a++)
{
my @Array = split('', ${$Ref_A_SeqCompare}[$a]);
for ($b=0; $b<@Array; $b++)
{
$A_Tableau[$b][$a]=$Array[$b];
}
}
my @A_Seq_Ref = split('', $Seq_Ref);
my @PositionsCles;
my $score = 0;
my $NbrNucTot = 0;
map {$NbrNucTot+= length($_);} @{$Ref_A_SeqCompare};
# print "Nombre de nucléotides comparés : ".$NbrNucTot."\n";
for (my $b =0; $b<@A_Seq_Ref; $b++)
{
foreach my $a (@{$A_Tableau[$b]}){if($A_Seq_Ref[$b] ne $a){$score++;}}
}
# print "Score absolu : $score\n";
my $ScoreRel = sprintf( "%.3f" ,$score/$NbrNucTot);
return($ScoreRel);
} |
découper une séquence en tronçons
Répondre avec citation
Envoyé par Jasmine80
Partager