Discussion :

[Isabella83]

Bonjour à tous,
J'ai des millions de séquences à analyser (taille variant entre 25 et 30 nt).
Je souhaite en fait identifier, pour chaque séquence, si il y a des sequences complémentaires à ma séquence, sachant que la complémentarité peut être sur toute la longueur de ma séquence, mais aussi bien sur les 10 premiers nucléotides, comme sur les 3 premiers, les 5 premier etc ...

Par exemple :
>seqA
UACGCAGAGGCCUAAGUAAAUAGUC
>seqB
GAAUUUCAUUUACCGGAUGCGUCUCC

Voici un exemple de complémentarité que je peux avoir :

Code :

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

1
2
----------------UACGCAGAGGCCUAAGUAAAUAGUC
GAAUUUCAUUUACCGGAUGCGUCUCC

Je ne vois pas du tout comment m'y prendre, quelqu'un aurait une idée à me proposer ?

Merci pour votre aide!

[Gardyen]

hum cela revient à aligner les séquences et les séquences complémentaires non ?
en recherchant des alignements locaux en priorité

[Isabella83]

oui c'est ca, j'ai recherché un module bioperl, mais je n'ai pas trouvé ...
Je pense peut etre utiliser la fonction index(), en cherchant toutes les sous chaines dans seq2 qui se trouvent dans seq1 ...

[Gardyen]

Transformer une séquence en sa séquence complémentaire est simple:

Code :

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

1
2
3
4
5
6
7
sub reverse_complement($){
	my $string = shift;
	my $res = $string;
	$res =~ tr/TAGCtagc/ATCGatcg/;
	$res = scalar reverse $res;
	return $res;
}

à adapter pour l'arn.

Ensuite il suffit de savoir quel logiciel d'alignement est suffisamment sensible pour obtenir de bons résultats... je n'ai pas testé plus en détail, mais en terme de sensibilité BLAST est pas mal il me semble, après il faut jouer avec les paramètres. Le prix à payer sera certainement en temps de calcul.

[Isabella83]

J'ai fait quelque chose comme ca :

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
#!/usr/bin/perl
use warnings;
use strict;
 
sub complement {
    $_[0] =~ y/CGATcgat/GCTAgcta/;
    return $_[0];
}
 
sub match {
	my ($s1, $s2) = @_; # TACGCAGAGGCCTAAGTAAATAGTC, GAATTTCATTTACCGGATGCGTCTCC
	#$s2 = reverse_complement $s2 ;
	complement $s2;
	print "$s1\n";
	my $s2l = length $s2;
	for (my $length = $s2l; $length; $length--) { # start from the longest possible substring
		for my $start (0 .. $s2l - $length) {     # starting position of the matching substring
			my $substr = substr $s2, $start, $length;
			#print "$substr\n";
			my $pos = index $s1, $substr; # searches for one string within another
			if ($pos + 1) {
				my $length_sub = length $substr;
				print "$pos\t$length_sub\n";
			}
		}
	}
}
 
 
print match('TACGCAGAGGCCTAAGTAAATAGTC',
'GAATTTCATTTACCGGATGCGTCTCC')
    ,"\n";

J'ai oublié de préciser que ce que je souhaite obtenir c'est la distance entre le début de seq1 et le premier nucléotide complémentaire à la seq1 de la seq2;
Mais je voudrais m'assurer que la complémentarité se produit bien entre la fin de la seq2 et le début de la seq1 pour toute séquence seq1 et seq2, comment puis je m'y prendre ?

Par exemple :
>seq1
UACGCAGAGGCCUAAGUAAAUAGUC
>seq2
GAAUUUCAUUUACCGGAUGCGUCUCC

Voici un exemple de complémentarité que je veux avoir :

Code :

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

1
2
3
 
----------------UACGCAGAGGCCUAAGUAAAUAGUC
GAAUUUCAUUUACCGGAUGCGUCUCC

Dans cet exemple, j'obtiens distance = 10

Et un exemple de complémentarité que je ne veux pas avoir :

>seq3
ATGCAGACAGACAGTAGACAGTAGAACCAGCAGAC
>seq4
ATCTTGGTCGTCTGGCGATGACAG

Code :

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

1
2
ATGCAGACAGACAGTAGACAGTAGAACCAGCAGAC
---------------------ATCTTGGTCGTCT

[trex7g2]

sur des millions de lignes, ton script risque d'être assez long... Tu as déjà deux boucles imbriquées dans ta sub, qui risque elle meme d'etre dans 1 voir 2 boucles pour croiser toutes tes séquences. Des logiciels comme BLAST utilisent des algorithmes très efficaces afin de trouver les alignements de séquence. Je te conseille de les utiliser, et éventuellement de trier les séquences qui t'interessent ou pas.

[Isabella83]

Ok comment faire un blast de mon fichier séquence contre mon autre fichier séquence avec bioperl ?