récupérer une sous-séquence

**Philou67430** · 26/11/2010, 13h25

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.

**Jasmine80** · 26/11/2010, 13h31

Code :

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

@{[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

**Philou67430** · 29/11/2010, 10h27

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 :

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

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

**Jasmine80** · 29/11/2010, 10h31

Merci beaucoup pour ton aide, j'ai bien compris maintenant

**Jasmine80** · 07/02/2011, 14h28

Code :

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

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

**Philou67430** · 07/02/2011, 14h49

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

**Jasmine80** · 07/02/2011, 14h53

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 :

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
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é

>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

**Jasmine80** · 07/02/2011, 15h02

j'ai oublié le join "",

... désolée

ça fonctionne maintenant

**Jasmine80** · 07/02/2011, 15h08

a mais non, c'était bon

**Philou67430** · 07/02/2011, 16h09

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

Code :

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

@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 :

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

(@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 :

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

@table = ($element1, $element2);

**Jasmine80** · 08/02/2011, 08h48

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

Code :

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

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

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

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 :

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

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

**Jasmine80** · 08/02/2011, 09h08

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

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
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";
 
}

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)

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

**Jasmine80** · 08/02/2011, 09h11

je vais rappeler le but du code

dans l'alignement :

Code :

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

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?

**Jasmine80** · 08/02/2011, 09h28

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
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";
 
}

super, ça fonctionne. Merci beaucoup pour ton aide

**Philou67430** · 08/02/2011, 11h25

Euh... qu'ai-je fais ?

récupérer une sous-séquence

Langage Perl

Discussions similaires

Partager

Partager