récupérer une sous-séquence

**Jasmine80** · 25/11/2010, 14h16

Code :

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

my $sequence = 'ATC.GATGAT....ACAGTAGATA..CAGATATG..TTCGCG...ATATAGGCAAG...ATGGATA.CCAGTA.T';

Problème :

A partir du Xième caractère de la séquence X = 30 (2ème A de CAGATATG)
Récupérer N nucléotides de part et d'autre N = 5

Résultat attendu TA..CAGATATG..T

Quelle est la manière la plus simple de procéder?

Merci pour votre aide,

**djibril** · 25/11/2010, 14h36

Avec un substr, non ?

**Jasmine80** · 25/11/2010, 14h40

Envoyé par djibril

Avec un substr, non ?

Je ne sais pas j'ai commencé quelque chose de compliqué :

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
 
my @sequence_list = split (//, $sequence);
 
my $x = 30;
my $n = 5;
 
# $x - 1 pour avoir l'indice de l'array
$x --;
 
 
my $subseq .= $sequence_list[$x];
 
 
my $ok = 0;
 
 
# récupération de n nucléotides à droite à partir de $xième
while ($ok < $n){
 
	# on progresse 1 à 1
	$x ++;
 
	if ($sequence_list[$x] =~ m/[ATCG]/i ){		
 
		$ok ++;
		$subseq .= $sequence_list[$x];
	}
 
}

**Jasmine80** · 25/11/2010, 14h45

Peut-on dans une expression régulière donner la position? en utilisant $* ?

Code :

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

=~ [\w\.*\w\.*\w\.*\w\.*\w\.*[\w\.]\.\w\.*\w\.*\w\.*\w\.*\w]

Récupérer 5 lettres à gauche et 5 à droite en disant que la lettre rouge (ou le point) est à la 30ième position

... le problème, c'est que le nombre de lettres à récupérer est une variable $n

**djibril** · 25/11/2010, 14h55

Code :

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

1
2
3
4
5
6
7
my $sequence = 'ATC.GATGAT....ACAGTAGATA..CAGATATG..TTCGCG...ATATAGGCAAG...ATGGATA.CCAGTA.T';
my $motif = 'CAGATATG';
my $N = 5;
if ( my ($avant, $apres) = $sequence =~ m/(.{$N})$motif(.{$N})/ ) {
  print "Avant : $avant\n";
  print "Apres : $apres\n";
}

Avant : ATA..
Apres : ..TTC

**djibril** · 25/11/2010, 15h01

Je ne sais pas si j'ai bien saisi car je ne comprends pas ton résultat obtenu par rapport à ton énoncé. Si tu veux juste 5 caractères avant le X=30 et 5 après, pourquoi ne pas utiliser un subtr comme ceci :

Code :

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

1
2
my $X = 30;
print substr $sequence,$X-5, 2*$N;

.CAGATATG.

**Philou67430** · 25/11/2010, 16h28

Je suis d'accord avec djibril, substr est la meilleur solution.
L'équivalent regexp pourrait être :

Code :

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

/.{15}(.{10})/

mais il sera forcément moins performant.

**Jasmine80** · 26/11/2010, 07h09

Envoyé par djibril

Je ne sais pas si j'ai bien saisi car je ne comprends pas ton résultat obtenu par rapport à ton énoncé. Si tu veux juste 5 caractères avant le X=30 et 5 après

Je ne veux pas 5 caractères, je veux 5 lettres (A, T, C ou G) quelques soit le nombre de points intermédiaires.

**Jasmine80** · 26/11/2010, 07h40

Voici un programme faisant ce que je veux, mais il devrait y avoir moyen de faire un meilleur code :

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
 
#!/usr/local/bin/perl
 
use strict;
use warnings;
 
 
my $sequence = 'ATC.GATGAT....ACAGTAGATA..CAGATATG..TTCGCG...ATATAGGCAAG...ATGGATA.CCAGTA.T';
 
=h
A partier du Xième caractère (2ème A de CAGATATG)		X = 30
Récupérer N nucléotides de part et d'autre			N = 5
 
Résultat TA..CAGATATG..T
 
 
=cut
 
 
my @sequence_list = split (//, $sequence);
 
my $x = 30;
# $x - 1 pour avoir l'indice de l'array
$x --;
 
my $n = 5;
my $ok = 0;
 
my $x2 = $x;
 
 
my $subseq .= $sequence_list[$x];
 
 
 
# récupération de n nucléotides à droite à partir de $xième
while ($ok < $n){
 
	# on progresse 1 à 1
	$x ++;
 
	if ($sequence_list[$x] =~ m/[ATCG]/i ){		
 
		$ok ++;
		$subseq .= $sequence_list[$x];
	}
 
}
 
$ok = 0;
 
 
# récupération de n nucléotides à gauche à partir de $xième
while ($ok < $n){
 
	# on progresse 1 à 1
	$x2 --;
 
	if ($sequence_list[$x2] =~ m/[ATCG]/i ){		
 
		$ok ++;
		$subseq = $sequence_list[$x2].$subseq ;
 
		print "$x2 => $sequence_list[$x2]\n";
	}
 
}
 
 
 
print "$subseq\n";

**Philou67430** · 26/11/2010, 10h09

Je ne comprends pas Jasmine. Tu dis cela

Envoyé par Jasmine80

Je ne veux pas 5 caractères, je veux 5 lettres (A, T, C ou G) quelques soit le nombre de points intermédiaires.

et tu attends ce résultat :

Envoyé par Jasmine80

Code :

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

1
2
3
4
5
6
...
=h
A partier du Xième caractère (2ème A de CAGATATG)		X = 30
Récupérer N nucléotides de part et d'autre			N = 5
 
Résultat TA..CAGATATG..T

Pourrais-tu préciser, car j'ai une autre solution, mais il faudrait savoir si :
- les X caractères représentent X lettres ou X caractères quelconques (incluant le .)
- les N nucléotides peuvent-elle contenir des . comme dans ton exemple (donc ce serait des caractères et non des lettres) ?

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

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, 13h49

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, 13h53

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