Analyse d'un alignement de séquences

Bonjour,

Voici le problème. J'ai un alignement de 53 séquences de 1498 nucléotides. Je voudrais retrouver pour chaque séquence, les tronçons de 5 nucléotides dont au minimum 3 nucléotides différents pour les 52 autres séquences alignées.
Quelle est la façon la plus simple de traiter les données?
Je pourrais placer l'ensemble des séquences dans un tableau à 2 dimensions, une matrice, et m'y déplacer nucléotide par nucléotide via une fenêtre de 5 nucléotides. Cela me permettrait de faire mes comparaisons par 5 colonnes ... mais cela sera un travail énorme vu le nombre de comparaisons à faire. N'existerait-il pas des modules pouvant alléger cette tâche?
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 0 | 1 | 2 | 3 ... 1497 0 A . T C ( . = gap = absence de nucléotide = nucléotide différent) 1 A T T C 2 A C T C 3 G T T A ... 51
Avec vous une idée à suggérer?

Pour ce qui est de la comparaison des pentanucléotides, existe-il une fonction renvoyant le nombre de lettre différentes entre 2 mots? Ou devrait-je comparer nucléotide par nucléotide? Il y a plusieurs mois, j'avais regardé dans les module destinés à la lexicographie mais je n'avais rien trouvé.

Merci pour votre aide.

Voici le début de mon programme :

Code:

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
#!/usr/local/bin/perl
 
use strict;
use warnings;
use Data::Dumper;
 
 
my %h_sequences;
$h_sequences{"A"} = "ATCG.TTGC";
$h_sequences{"B"} = "ATGG.TCCT";
$h_sequences{"C"} = "ATGGATTGG";
 
my $length_ali = length($h_sequences{"A"});
my $nbr_seq = (keys %h_sequences);
 
# création de la matrice
my @a_matrix;
my $i = 0;
 
foreach my $id (keys %h_sequences){
	@{$a_matrix[$i]} = split (//, $h_sequences{$id});
	$i++;	
}
 
# parcours position par position
for (my $a = 0; $a<$length_ali; $a++){
	# parcours espèce par espèce
	for (my $b = 0; $b<$nbr_seq-4; $b++){
		$a_matrix[$b][$a];
		$a_matrix[$b][$a+1];
		$a_matrix[$b][$a+2];
		$a_matrix[$b][$a+3];
		$a_matrix[$b][$a+4];
	}
 
}

Cela vous semble t'il bien?

Merci,

Voila, le programme est terminé et fonctionne bien même si il est loin d'être optimal.

Code:

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
#!/usr/local/bin/perl
 
use strict;
use warnings;
use Data::Dumper;
 
 
my %h_sequences;		 
$h_sequences{"A"} = "AGCG.TTGC";   # $i = 0
$h_sequences{"B"} = "CTGGTTCCT";   # $i = 1
$h_sequences{"C"} = "CGGGATTGG";   # $i = 2
 
=h
 
 
 
$h_sequences{"A"} = "AGCG.TTGC";   # $i = 0
$h_sequences{"B"} = "CTGGTTCCT";   # $i = 1
$h_sequences{"C"} = "CGGGATTGG";   # $i = 2
 
 
=> seq B pos 3 GTTCC     B correspond à la séquence $s1 comparée ($i) $s1 = 1				G T T  C C
# $scores_pentanuc[$p][$s1][$s2] 
# $scores_pentanuc[3][1][0]   score doit être > ou = à 3        # pentanucléotide de la séq A    $s2 = 0		G  . T  T G    scores_pentanuc =  3  (= nbr nucléotides différents entre GTTCC et G.TTG)
# $scores_pentanuc[3][1][2]  score doit être > ou = à 3         # pentanucléotide de la séq C   $s2 = 2		G A T T  G    scores_pentanuc = 3   (= nbr nucléotides différents entre GTTCC et GATTG)
# NE PAS COMPARER A SON PROPRE SCORE $scores_pentanuc[0][1][1] qui est d'office nul
 
=> tous lles pentanucléotides comparés à  GTTCC on donc bien au minimum un score de 3
$test = 1;
on garde ce pentanucléotide
 
 
=cut
 
# 3 séquences de 9 nucléotides => 27 nucléotides
# nombres de comparaisons  = nbr nuc * nbr seq * nbr seq = 3 * 3 * 9 = 81  
 
my $length_ali = length($h_sequences{"A"});
my $nbr_seq = (keys %h_sequences);
 
# création de la matrice
my @a_matrix;
my $i = 0;
 
my @a_id;
foreach my $id (sort keys %h_sequences){
	@{$a_matrix[$i]} = split (//, $h_sequences{$id});
	$a_id[$i]=$id;
	$i++;	
}
 
 
# Création de la matrice des scores mononucléotides vs mononucléotides
#--------------------------------------------------------------------------------
# array des scores
# $score[$n][$a]  séquence $n vs $m
# score = nb de nucléotides différents de toutes les autres séquences
my @scores;
 
my $nbr= 0;
# parcours seq par seq 
for (my $s1 = 0; $s1<$nbr_seq; $s1++){
	# parcours seq par seq  
	for (my $s2 = 0; $s2<$nbr_seq; $s2++){
		# parcours position par position
		for (my $p = 0; $p<$length_ali; $p++){	
			if ($a_matrix[$s1][$p] ne $a_matrix[$s2][$p]){
				$scores[$p][$s1][$s2]= 1;					
			}
			else{
				$scores[$p][$s1][$s2] = 0;
			}
		}
	}
}
 
# Création de la matrice des scores des pentanucléotides vs pentanucléotides
#----------------------------------------------------------------------------------------
# $scores_pentanuc[$p][$s1][$s2]
# score du pentanucléotide commençant à la position $p
# pour des séquences $s1 et $s2
my @scores_pentanuc;
 
# parcours seq par seq 
for (my $s1 = 0; $s1<$nbr_seq; $s1++){
	# parcours seq par seq  
	for (my $s2 = 0; $s2<$nbr_seq; $s2++){
		# parcours position par position
		for (my $p = 0; $p<$length_ali-4; $p++){	
			$scores_pentanuc[$p][$s1][$s2] = $scores[$p][$s1][$s2] + $scores[$p+1][$s1][$s2] + $scores[$p+2][$s1][$s2] + $scores[$p+3][$s1][$s2] + $scores[$p+4][$s1][$s2];
		}
	}
}
 
 
# Pour un pentanucléotide, on regarde si tous les autres
# pentanucléotides de l'alignement sont au minimum
# différents de 3 nucléotides (score > ou = 3)
 
# exemple, analyse du pentanucléotide de la séq $s1 (1)
# commençant à la position $p 0 (0)
# $scores_pentanuc[$p][$s1][$s2] 
# $scores_pentanuc[0][1][0]   score doit être > ou = à 3
# $scores_pentanuc[0][1][2]  score doit être > ou = à 3
# NE PAS COMPARER A SON PROPRE SCORE $scores_pentanuc[0][1][1] qui est d'office nul
 
print "Pentanucléotides possédant au minimum 3 nucléotides différents avec toutes les autres séquences\n";
# parcours position par position, pentanucléotide par pentanucléotide (renseigné par la position de son 1ier nuc)
for (my $p = 0; $p<$length_ali-4; $p++){	
	# parcours seq par seq 
	for (my $s1 = 0; $s1<$nbr_seq; $s1++){
		my $test = 1;
		# parcours seq par seq  
		for (my $s2 = 0; $s2<$nbr_seq; $s2++){
			if(($s1 != $s2) && ($scores_pentanuc[$p][$s1][$s2] < 3)){
				$test = 0;
			}
		}
		# test réussi si $test reste à 1
		if($test){
			my $pentanuc_succeed = $a_matrix[$s1][$p].$a_matrix[$s1][$p+1].$a_matrix[$s1][$p+2].$a_matrix[$s1][$p+3].$a_matrix[$s1][$p+4]; 
			print "=> seq $a_id[$s1] pos $p $pentanuc_succeed\n";
		}
	}
}