récupérer une sous-séquence

Version imprimable

Voir 40 message(s) de cette discussion en une page

26/11/2010, 12h25
Philou67430

Dans un tableau @a, l'élément $a[-5] est le 5e en partant de la fin...
Les éléments de -5 à -1 sont donc les 5 derniers éléments du tableau.
26/11/2010, 12h31
Jasmine80

Code:

@{[substr($sequence, 0, $X-1) =~ /([ATCG])/g]}[-5..-1]

Comment sait-on que les 5 derniers éléments sont bien des lettres et pas de points? Grâce à l'expression régulière? Je crois que j'ai compris mais j'aurais dû mal de l'écrire pas moi même.

Merci beaucoup pour ton aide
29/11/2010, 09h27
Philou67430
L'expression régulière capture "globalement" (/g => toutes les occurences) des caractères A, T, C ou G et les place dans un tableau anonyme (à l'aide de [ ... ].
Il suffit de déréférencer le résultat pour les éléments de -5 à -1 (du 5e depuis la fin au dernier).

Pour être plus clair, tu peux écrire :
Code:

1 2 @nuc_before = substr($sequence, 0, $X-1) =~ /([ATCG])/g; join "", @nuc_before[-5..-1];
29/11/2010, 09h31
Jasmine80

Merci beaucoup pour ton aide, j'ai bien compris maintenant :D
07/02/2011, 13h28
Jasmine80

Code:

my @subseq2 = @{[substr($hash_ali{$id2}, 0, $X) =~ /([A-Z])/gi]}[-$N..-1], @{[substr($hash_ali{$id2}, $X-1) =~ /([A-Z])/gi]}[0 .. $N];

j'obtiens l'erreur 'Useless use of array slice in void context at best_variable_region.pl line 76.' savez-vous pourquoi?

Merci
07/02/2011, 13h49
Philou67430

Quelles sont les valeurs de :
- $N
- [substr($hash_ali{$id2}, 0, $X) =~ /([A-Z])/gi]
- [substr($hash_ali{$id2}, $X-1) =~ /([A-Z])/gi]

merci de me répondre si rapidement

j'utilise une boucle qui analyse un fichier, pour la majorité des cas, cela fonctionne mais il doit y en avoir un qui coince.

voici le 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
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
#!/usr/local/bin/perl
 
 
 
use strict;
use warnings;
 
use Data::Dumper;
use Bio::SeqIO;
 
 
#-------------------------------- best_variable_region_per_seq.pl
 
 
 
# fichier d'entrées
my $infile = 'P:/Theorie/Leonid/pamgene/purA/test.fsa';
my $in  = Bio::SeqIO->new(-file => $infile , '-format' => 'fasta');
 
# taille de la région de part et d'autre du point central
my $N =  4;
my $l;
 
my %hash_ali;
 
my $start = 0;
 
my $gap_symbol = '';
 
while ( my $seq_IO = $in->next_seq() ) {
 
	$hash_ali{$seq_IO->primary_id} = $seq_IO->seq ;
	$l = length ($seq_IO->seq);
 
 
	# on recherche la position où toutes 
	# les séquences sont alignées
	my ($s) = $seq_IO->seq =~ m/^([^a-z]*)/i;
	my $l = length ($s) + 1;
 
	if ($start < $l){
		$start = $l
	}
 
	($gap_symbol) = $seq_IO->seq =~ m/([^a-z])/i;
}
 
 
 
my $k = $l - 2 * $N + 1;
 
foreach my $id (sort keys %hash_ali){
 
	my %score_id;
 
	my @array_score;
 
	for my $X ( $N + 1 + $start ..$k){
 
		# my $subseq = join "", @{[substr($hash_ali{$id}, 0, $X-1) =~ /([A-Z])/gi]}[-$N..-1], @{[substr($hash_ali{$id}, $X-1) =~ /([A-Z])/gi]}[0 .. $N];
 
		my @subseq = @{[substr($hash_ali{$id}, 0, $X) =~ /([A-Z])/gi]}[-$N..-1], @{[substr($hash_ali{$id}, $X-1) =~ /([A-Z])/gi]}[0 .. $N];
 
=h
		my $rx = join ('\\'.$gap_symbol.'*', @subseq );
		my ($subseq_ali) =  $hash_ali{$id} =~m/($rx)/;	
 
		my $pos = $X - 2 * $N;
 
		my $score = 0;
 
		foreach my $id2 (keys %hash_ali){
 
			if ($id2 ne $id){
 
 
				my @subseq2 = @{[substr($hash_ali{$id2}, 0, $X) =~ /([A-Z])/gi]}[-$N..-1], @{[substr($hash_ali{$id2}, $X-1) =~ /([A-Z])/gi]}[0 .. $N];
 
				# print join "$id\t".('', @subseq)."\t$id2".('', @subseq2)."\n";				
				for my $j (0..$#subseq){
 
					if ($subseq[$j] ne $subseq2[$j]){
						$score ++;
					}				
				}
 
			}
 
		}
 
		print "$id\t$pos\t$score\t$subseq_ali\n";
 
		push @{$array_score[$score]}, $score;
 
=cut	
 
	}
 
 
	 # $score_id{$id} = \(sort {$b<=>$a} @array_score)[0..10];
 
	# print Dumper %score_id ;
 
}

et voici le fichier analysé

Citation:

>SHAE
----------ACTGAGAAGATTG--AAATACGTACATTA-CATAATTCAG
>SSAP
----GGGCCAACTGAGAAGATTG--AAATCTTCACGTCA-CATAATTCAG
>Staphylococcus_epidermidis
AGCAGAGCAAGCAGACGTAATTGCTAGATTTTCTGGTGGTAACAATGCGG
>Staphylococcus_oralis
ATCTGGCCCAACTGAGAAAGTTG--AGATACGAACACCA-ACCAACTCGC
>Staphylococcus_carnosus
---GGGGCCAACTGAGAATGTTG--AAATGCGAACACCA-ACGAGTTCAC
>Staphylococcus_haemolyticus
----------ACTGAGAAGATTG--AAATACGTACATTA-CATAATTCAG
>Staphylococcus_hominis
GCCTGGACCTACTGAGAAGATTG--AAATATGAACGCCA-CATAATTCAG
>Staphylococcus_aureus
ATCTGGACCAACTGAGAAGATAG--AAATTTGTACATTA-CATAATTCTG

le but étant de trouver la séquence de taille (2 * $N) la plus variable pour chaque séquence (par rapport à l'ensemble des autres séquences)

merci :)

07/02/2011, 14h02
Jasmine80

... join

j'ai oublié le join "", :?:? ... désolée

:ccool: ça fonctionne maintenant
07/02/2011, 14h08
Jasmine80

il y a toujours un problème

a mais non, c'était bon
07/02/2011, 15h09
Philou67430

Je pense avoir trouvé d'après ta réponse :

Code:

@table = $element1, $element2;

ne fait pas ce que l'on voudrait que cela fasse.
En effet, dans ce contexte, la priorité des opérateurs fait que cette expression s'exécute en fait comme ceci :

Code:

(@table = $element1), $element2;

Résultat, on n'a qu'un élément dans @table, et $element2 est inutilisé (d'où l'erreur).

Pour faire ce que l'on attends (mettre $element1 et $element2 dans @table) il faut écrire :

Code:

@table = ($element1, $element2);

merci, cela fonctionne très bien mais j'ai encore un petit problème
Code:

1 2 my @subseq = (@{[substr($hash_ali{$id}, 0, $X) =~ /([A-Z])/gi]}[-$N..-1], @{[substr($hash_ali{$id}, $X-1) =~ /([A-Z])/gi]}[1 .. $N]); print "@{[substr($hash_ali{$id}, 0, $X) =~ /([A-Z])/gi]}[-$N..-1]\n@{[substr($hash_ali{$id}, $X-1) =~ /([A-Z])/gi]}[1 .. $N])\n@subseq\n";
on obtient bien

Citation:

T G A G
A A G A
T G A G A A G A

et dans l'alignement on a TGAGAAGA => donc, c'est bon

par contre pour

Citation:

A T T G
A A T A
A T T G A A T A

et dans l'alignement on a ATTG--AAATA => donc, c'est mauvais

il manque donc un nucléotide

Si j'utilise
Code:

1 2 my @subseq = (@{[substr($hash_ali{$id}, 0, $X) =~ /([A-Z])/gi]}[-$N..-1], @{[substr($hash_ali{$id}, $X-1) =~ /([A-Z])/gi]}[0 .. $N]); print "@{[substr($hash_ali{$id}, 0, $X) =~ /([A-Z])/gi]}[-$N..-1]\n@{[substr($hash_ali{$id}, $X-1) =~ /([A-Z])/gi]}[0 .. $N])\n@subseq\n";
c'est l'inverse, on trouve le second mais pas le premier

est-ce clair? C'est ce que j'obtiens avec l'algorithme utilisé plus haut

merci pour ton aide

voici le code adapté afin de ne pas devoir lire de fichier

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;
use Bio::SeqIO;
use List::MoreUtils;
 
 
#-------------------------------- best_variable_region_per_seq.pl
 
 
 
# fichier d'entrées
my $infile = 'P:/Theorie/Leonid/pamgene/purA/test.fsa';
my $in  = Bio::SeqIO->new(-file => $infile , '-format' => 'fasta');
 
# taille de la région de part et d'autre du point central
my $N =  4;
 
# nombre des meilleurs hits à trouver
my $nb_hit = 4;
 
my $l;
 
my $start = 0;
 
my $gap_symbol = '';
 
my %hash_ali = (
 
	SHAE	=>	'----------ACTGAGAAGATTG--AAATACGTACATTA-CATAATTCAG',
	SSAP	=>	'----GGGCCAACTGAGAAGATTG--AAATCTTCACGTCA-CATAATTCAG',
	Staphylococcus_epidermidis	=>	'AGCAGAGCAAGCAGACGTAATTGCTAGATTTTCTGGTGGTAACAATGCGG',
	Staphylococcus_oralis	=>	'ATCTGGCCCAACTGAGAAAGTTG--AGATACGAACACCA-ACCAACTCGC',
	Staphylococcus_carnosus	=>	'---GGGGCCAACTGAGAATGTTG--AAATGCGAACACCA-ACGAGTTCAC',
	Staphylococcus_haemolyticus	=>	'----------ACTGAGAAGATTG--AAATACGTACATTA-CATAATTCAG',
	Staphylococcus_hominis	=>	'GCCTGGACCTACTGAGAAGATTG--AAATATGAACGCCA-CATAATTCAG',
	Staphylococcus_aureus	=>	'ATCTGGACCAACTGAGAAGATAG--AAATTTGTACATTA-CATAATTCTG',
 
 
);
 
 
foreach my $id  ( keys %hash_ali ) {
 
	$l = length ($hash_ali{$id});
 
 
	# on recherche la position où toutes 
	# les séquences sont alignées
	my ($s) = $hash_ali{$id}=~ m/^([^a-z]*)/i;
	my $l = length ($s) + 1;
 
	if ($start < $l){
		$start = $l
	}
 
	($gap_symbol) = $hash_ali{$id} =~ m/([^a-z])/i;
}
 
my $k = $l - 2 * $N + 1;
 
 
foreach my $id (sort keys %hash_ali){
 
 
	my %score_id;
 
	my @array_score;
 
	for my $X ( $N + 1 + $start ..$k){
 
		# my $subseq = join "", @{[substr($hash_ali{$id}, 0, $X-1) =~ /([A-Z])/gi]}[-$N..-1], @{[substr($hash_ali{$id}, $X-1) =~ /([A-Z])/gi]}[0 .. $N];
 
		my @subseq = (@{[substr($hash_ali{$id}, 0, $X) =~ /([A-Z])/gi]}[-$N..-1], @{[substr($hash_ali{$id}, $X-1) =~ /([A-Z])/gi]}[1 .. $N]);
 
		my $rx = join ('\\'.$gap_symbol.'*', @subseq );
		my ($subseq_ali) =  $hash_ali{$id} =~m/($rx)/;	
 
		my $pos = $-[1];
 
		my $score = 0;
# print "@{[substr($hash_ali{$id}, 0, $X) =~ /([A-Z])/gi]}[-$N..-1]\n@{[substr($hash_ali{$id}, $X-1) =~ /([A-Z])/gi]}[0 .. $N])\n@subseq\n";	
# print "$pos\n$rx\n";
 
 
		foreach my $id2 (keys %hash_ali){
 
			if ($id2 ne $id){
 
 
				my @subseq2 = (@{[substr($hash_ali{$id2}, 0, $X) =~ /([A-Z])/gi]}[-$N..-1], @{[substr($hash_ali{$id2}, $X-1) =~ /([A-Z])/gi]}[1 .. $N]);
 
				for my $j (0..$#subseq){
 
					if ($subseq[$j] ne $subseq2[$j]){
						$score ++;
					}				
				}
 
			}
 
		}
 
		# print "$id\t$pos\t$score\t$subseq_ali\n";
 
		$array_score[$pos] = $score;
	}
 
 
	my @index = ();
 
	foreach my $s (0 .. $#array_score) {
	  @index = (sort { $array_score[$b] <=> $array_score[$a] } grep defined, @index, $s)[0.. $nb_hit];
	}
 
	print "$id\n";
	print map "$array_score[$_] (pos $_)\n", @index;	
	print "\n\n";
 
}

Citation:

SHAE
24 (pos 30)
23 (pos 29)
22 (pos 32)
21 (pos 31)
21 (pos 34)

SSAP
33 (pos 29)
30 (pos 28)
30 (pos 30)
27 (pos 31)
26 (pos 26)

on obtient donc bien les 4 meilleurs scores par séquence mais il y a toujours l'erreur citée plus haut, le fait que certaines sous-séquences erronées ne soient pas retrouvées dans l'alignement (cf : dans l'alignement on a ATTG--AAATA => donc, c'est mauvais)

Citation:

Use of uninitialized value in numeric comparison (<=>) at test.pl line 118.

but

je vais rappeler le but du code

dans l'alignement :

Code:

1
2
3
4
5
6
7
8
	SHAE				=>	'----------ACTGAGAAGATTG--AAATACGTACATTA-CATAATTCAG',
	SSAP				=>	'----GGGCCAACTGAGAAGATTG--AAATCTTCACGTCA-CATAATTCAG',
	Staphylococcus_epidermidis	=>	'AGCAGAGCAAGCAGACGTAATTGCTAGATTTTCTGGTGGTAACAATGCGG',
	Staphylococcus_oralis		=>	'ATCTGGCCCAACTGAGAAAGTTG--AGATACGAACACCA-ACCAACTCGC',
	Staphylococcus_carnosus		=>	'---GGGGCCAACTGAGAATGTTG--AAATGCGAACACCA-ACGAGTTCAC',
	Staphylococcus_haemolyticus	=>	'----------ACTGAGAAGATTG--AAATACGTACATTA-CATAATTCAG',
	Staphylococcus_hominis		=>	'GCCTGGACCTACTGAGAAGATTG--AAATATGAACGCCA-CATAATTCAG',
	Staphylococcus_aureus		=>	'ATCTGGACCAACTGAGAAGATAG--AAATTTGTACATTA-CATAATTCTG',

il faut , pour chaque séquence, trouver la région de 2*$N+1 nucléotides la plus variable par rapport au reste des séquences (avec $score comptant le nombre de différences totales)

est-ce plus clair?

voici le code corrigé

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
126
127
128
#!/usr/local/bin/perl
 
 
 
use strict;
use warnings;
 
use Data::Dumper;
use Bio::SeqIO;
use List::MoreUtils;
 
 
#-------------------------------- best_variable_region_per_seq.pl
 
 
 
# fichier d'entrées
my $infile = 'P:/Theorie/Leonid/pamgene/purA/test.fsa';
my $in  = Bio::SeqIO->new(-file => $infile , '-format' => 'fasta');
 
# taille de la région de part et d'autre du point central
my $N =  4;
 
# nombre des meilleurs hits à trouver
my $nb_hit = 4;
 
my $l;
 
my $start = 0;
 
my $gap_symbol = '';
 
my %hash_ali = (
 
	SHAE				=>	'----------ACTGAGAAGATTG--AAATACGTACATTA-CATAATTCAG',
	SSAP				=>	'----GGGCCAACTGAGAAGATTG--AAATCTTCACGTCA-CATAATTCAG',
	Staphylococcus_epidermidis	=>	'AGCAGAGCAAGCAGACGTAATTGCTAGATTTTCTGGTGGTAACAATGCGG',
	Staphylococcus_oralis		=>	'ATCTGGCCCAACTGAGAAAGTTG--AGATACGAACACCA-ACCAACTCGC',
	Staphylococcus_carnosus		=>	'---GGGGCCAACTGAGAATGTTG--AAATGCGAACACCA-ACGAGTTCAC',
	Staphylococcus_haemolyticus	=>	'----------ACTGAGAAGATTG--AAATACGTACATTA-CATAATTCAG',
	Staphylococcus_hominis		=>	'GCCTGGACCTACTGAGAAGATTG--AAATATGAACGCCA-CATAATTCAG',
	Staphylococcus_aureus		=>	'ATCTGGACCAACTGAGAAGATAG--AAATTTGTACATTA-CATAATTCTG',
 
 
);
 
 
foreach my $id  ( keys %hash_ali ) {
 
	$l = length ($hash_ali{$id});
 
 
	# on recherche la position où toutes 
	# les séquences sont alignées
	my ($s) = $hash_ali{$id}=~ m/^([^a-z]*)/i;
	my $l = length ($s) + 1;
 
	if ($start < $l){
		$start = $l
	}
 
	($gap_symbol) = $hash_ali{$id} =~ m/([^a-z])/i;
}
 
my $k = $l - 2 * $N + 1;
 
 
foreach my $id (sort keys %hash_ali){
 
 
	my %score_id;
 
	my @array_score;
 
	for my $X ( $N + 1 + $start ..$k){
 
		my @subseq = (@{[substr($hash_ali{$id}, 0, $X) =~ /([A-Z])/gi]}[-$N..-1], @{[substr($hash_ali{$id}, $X) =~ /([A-Z])/gi]}[0 .. $N]);
 
		my $rx = join ('\\'.$gap_symbol.'*', @subseq );
		my ($subseq_ali) =  $hash_ali{$id} =~m/($rx)/;	
 
		my $pos = $-[1];
 
		my $score = 0;
# print "@{[substr($hash_ali{$id}, 0, $X) =~ /([A-Z])/gi]}[-$N..-1]\n@{[substr($hash_ali{$id}, $X) =~ /([A-Z])/gi]}[0 .. $N])\n@subseq\n";	
# print "$pos\n$rx\n";
 
 
		foreach my $id2 (keys %hash_ali){
 
			if ($id2 ne $id){
 
 
				my @subseq2 = (@{[substr($hash_ali{$id2}, 0, $X) =~ /([A-Z])/gi]}[-$N..-1], @{[substr($hash_ali{$id2}, $X) =~ /([A-Z])/gi]}[0 .. $N]);
 
				for my $j (0..$#subseq){
 
					if ($subseq[$j] ne $subseq2[$j]){
						$score ++;
					}				
				}
 
			}
 
		}
 
		# print "$start\t$id\t$pos\t$score\t$subseq_ali\n";
 
		$array_score[$pos] = $score;
	}
 
 
	my @index = ();
	$start += 1;
 
	# attention, toutes les indexs de @array-score n'ont pas une valeur
	# pas de valeur si un gap se trouve à cette position
	foreach my $s ($start .. $#array_score) {
		if (defined $array_score[$s]){	
		  @index = (sort { $array_score[$b] <=> $array_score[$a] } grep defined, @index, $s)[0.. $nb_hit];
		}
	}
 
	print "$id\n";
	print map "$array_score[$_] (pos $_)\n", @index;	
	print "\n\n";
 
}

:ccool: super, ça fonctionne. Merci beaucoup pour ton aide

08/02/2011, 10h25
Philou67430

Euh... qu'ai-je fais ? 8O

Voir 40 message(s) de cette discussion en une page