Discussion :

[Jasmine80]

Bonjour,


J'ai plusieurs séquences du même gène mais dont certaines sont dans le sens 5'->3' et d'autres dans le sens 3'->5'. Il faudrait un algorithme qui me permette de classer ces séquences en 2 groupes pour ensuite pouvoir effectuer la fonction revcom sur l'un d'eux. Quelle approche pourrais-je utiliser? A l'oeil, on voit clairement qu'il y a un mélange de 2 types de séquences mais existe-t'il un module capable de distinguer ces 2 groupes?

Une approche naive, serait de prendre la première séquence et de calculer le pourcentage d'identité qu'elle possède avec toutes les autres séquences. Au deçà d'un certain seuil, la séquence analysée serait considérée comme prise dans le mauvais sens et retournée.

Voyez-vous plus intelligent comme approche?


Merci,

[Jasmine80]

Que pensez-vous de ce programme ... je sais qu'il reste des choses à améliorer mais il fonctionne.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
#!/usr/local/bin/perl
 
 
use strict;
use warnings;
use Bio::SeqIO;
 
my $in  = Bio::SeqIO->new(-file => "carB.txt", '-format' => 'Fasta');
 
# séquences qui serviront de référence aux deux groupes
my $ref_seq1;
my $ref_seq2;
 
 
# liste des séquences du groupe1, du groupe2 et celles sans groupe
my (@group1, @group2, @undef);
 
# on récupère la première séquence qui servira de référence au groupe 1
while ( my $seq = $in->next_seq()){
 
	push (@group1, $seq->primary_id );
	$ref_seq1 = $seq->seq;
	last;
 
}
 
# clé : id    val : sequence
my %sequences_tab;
 
# décanucléotide contenu dans les séquences du groupe 1
my ($ref1_polynuc) = $ref_seq1 =~ m/\w{10}(\w{10})/;
 
# décanucléotide contenu dans les séquences du groupe 2
my $ref2_polynuc;
 
 
 
 
while ( my $seq = $in->next_seq()){
 
	$sequences_tab{$seq->primary_id} = $seq->seq;
	print "=> ".$seq->primary_id."\n";
	# la séuqnece contient le polynucléotide de référénce du groupe 1
	# mise dans le groupe 1
	if ($seq->seq =~ m/$ref1_polynuc/){
		push (@group1, $seq->primary_id );
	}
	# sinon est mise dans le groupe undef
	else{
		push (@undef, $seq->primary_id );
		# le polynuc de référence du second groupe sera contenu dans
		# la dernière séquence passant par ce else
		($ref2_polynuc) = $seq->seq =~ /\w{10}(\w{10})/; 
	}
}
 
# on vérifie que $ref2_polynuc n'appartient pas au groupe 1
if($ref_seq1 =~ /$ref2_polynuc/){
	# on doit changer de polynucléotide de référence du groupe 2
	die "erreur du polynuc de référence du groupe 2";
 
}
 
for  my $i (0..$#undef){
 
	if ($sequences_tab{$undef[$i]} =~ m/$ref2_polynuc/){
		push (@group2, $undef[$i]);
		delete $undef[$i];
	}
}
 
# si il reste des séquences non classées, on recommence avec un autre polynucléotide de référence
if (@undef > 0){
	my ($ref1_polynuc) = $ref_seq1 =~ m/\w{30}(\w{10})/;
	my ($ref2_polynuc) = $ref_seq2 =~ m/\w{30}(\w{10})/;
 
	for  my $i (0..$#undef){
 
		if ($sequences_tab{$undef[$i]} =~ m/$ref1_polynuc/){
			push (@group1, $undef[$i]);
			delete $undef[$i];
		}
		elsif ($sequences_tab{$undef[$i]} =~ m/$ref2_polynuc/){
			push (@group2, $undef[$i]);
			delete $undef[$i];
		}
	}
}
 
# si après 2 vérifications, toutes les séquences ne sont pas classées
# on arrête le programme
if (@undef){
	die "toutes les séquences n'ont pas été classées";
}
 
print "groupe 1 :\n";
map{print "\t$_\n";} @group1;
 
 
print "groupe 2 :\n";
map{print "\t$_\n";} @group2;

Merci beaucoup pour vos conseils,

[Zwiter]

Salut Jasmine80,
L'alignement que tu proposais dans ton 1er message n'etait pas une bonne idée ?
Est ce que tes séquences sont chevauchantes ?
est-ce plutot comme ca :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
----    ---------     --------  ------
  ---------  ------------   --------  -------
--------  ------  ---------  -------

ou bien :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
          -----------------
  -----------------
    ----------------------------------

Quelles sont leur tailles ?

Pourquoi ne pas constituer un grand fichier fasta contenant les séquences ET leur reverse.
Je suis sur qu'il existe des softs/scripts pour rassembler les contigs issus de séquencages.

Z.

[Jasmine80]

Le problème est que j'ai 163 séquences du même gène, d'un peu plus de 3000 nucléotides et séparées par espèces dans plusieurs fichiers. Je dois dire que les séquences sont très belles et ont toutes le même ordre de taille, il n'y a donc pas de problème de chevauchement.

Un autre problème est que quand j'ai mes 2 groupes, comment savoir lequel est du 5'->3' et lequel est du 3'->5'. Je vais considérer que le groupe majoritaire est le sens 5'->3' que je conserve. En cas d'égalité, ça sera au petit bonheur la chance.

Le but est de créer un consensus par espèce pour ensuite en faire un alignement.

Merci pour ton aide,

[Zwiter]

Le sens global est facile a trouver car tu peux reperer le sens de la traduction : Il suffit de compter le nombre de codons stop dans chaque brins, pour les 3 cadres de lectures, ou les 6 en considerant le reverse.

SI tout va bien, tu trouveras le 'bon sens', ou plutot le bon brin, et en plus, tu pourras separer tes 2 groupes.
Z.

[Jasmine80]

C'est encore plus simple que je ne pensais, toutes mes séquences possèdent le codon start ... je n'ai pas encore vérifié pour toutes mais si c'est le cas, je me suis vraiment cassée la tête inutilement .