couper une séquence fasta en petite séquence

Version imprimable

31/10/2014, 17h32
Isabella83

1 pièce(s) jointe(s)

couper une séquence fasta en petite séquence

Bonjour à tous,
Je suis face à un problème que je n'arrive pas à résoudre. Peut être puis je le faire en perl.
J'ai une séquence de taille N. Je cherche à la "couper" en petit fragment, mais la moyenne des tailles de ces petits fragments doit être égale à une certaine valeur (qui peut changer), mais pour l'exemple mettons 50 à +/- 2 nt près.
Les séquences doivent se suivre mais pas se chevaucher !
Je ne sais vraiment pas comment m'y prendre. Quelqu'un aurait une idée ?
j'ai fait une image pour illustrer mon problème.Pièce jointe 160348
Une idée pour me faire avancer ?
Merci d'avance !
31/10/2014, 18h58
Lolo78

Oui, ça paraît tout-à-fait possible en PerL

Pourrais-tu donner un échantillon de cette séquence? Je suppose qu'il s'agit d'une série de lettres bout à bout, du genre ATGGGTACCAAT. Et je suppose que la séquence figurera dans une variable (ou que l'on peut la lire dans une variable à partir d'un fichier).

Peut-on couper à n'importe quel endroit dans la séquence?

Enfin, si je comprends bien, tu veux varier la taille des fragments par rapport à la valeur pivot choisie pour que le dernier fragment ne soit pas beaucoup plus petit que les autres. C'est cela? Si oui, il faut trouver la taille totale de la séquence et, dans l'exemple d'une valeur pivot de 50, trouver la valeur comprise entre 48 et 52 qui divise exactement la taille de la séquence ou, à défaut, celle donnant le plus grand reste.

J'attends tes réponses à ces questions avant de te proposer une solution en Perl.
03/11/2014, 09h38
Isabella83
Citation:

Envoyé par Lolo78

Pourrais-tu donner un échantillon de cette séquence?

Par exemple :

Code:

CCGCAGTCGCCGCCTCCGCAGCAGCAGCAGCAGCAGCAACAACTACAACAACAACAACAACAGCCACCAGGGCAGCAACAACCACACCAGTTACAACAGCCCCAACACCAACAGCGATAACAACAACATATAGATCGCGGAATCTGCTGAAAGCCACCGCCAAAGTGAAGAAGAAGAAAAAGTGGAAGCAAATCTCCTTCGTCGACGACCAGCAACGCTTGAAGCAAGCAGCGACGTCGGCGGAGCGCAAACGCTCCACCAGCAGCGCCCCGCCCGCCCTCGAGAGCAACGTACCGCCGCGGTTGGGGCGCTTTCTGCGGCAAAAGAGCAGCGAAACTAACAGCGCTGCAACGTCAGCAGCAGCGGCGGCAGTGCAAACGGAAGCGGAAACAGGACGCCGGAAATCGTCGGCGGCATCGATAATCGCAGCGGCAGCGGGCGCTGCTGGTTTGCACAACAATAACAGCAGCAGCTGCGCAGCGACCGCCAGCGAGGAGCAAAACTCGAGCGGTACGCTCACCGGCAGCAGCGGCAATACTTCCACGTCGACGTGGTGCGGCAGCATCAACTTCCTGGTGTATATGGTGTGCTCGTTTTGCTGCTGTTGCTACAACTTTCGCAATTCGCCTTTCGATGTAGAAAATGGCCAGGGCGCTAGATCACCGCTCGAGGGCGGCTCACCCAGCGCCGGTTTGGTACTACAGAATTTGCCGCAGCGTCGCGAATCGTTTCTATATCGCTCCGATTCGGACTTTGAGATGTCACCCAAGTCCATGTCCCGAAACTCAAGCATCGCTAGCGAAAGGTTTAAAGAACAGGAAGCCTCTATACTCGTTGACCGATCCCATGGCGAGGATCTCATTGTGACGCCTTTCGCCCAAATTTTAGCCAGCTTACGTTCAGTGCGCAACAATTTATTAAGTCTGACAAATGTTCCAGCATCAAACAAATCCAGGCGGCCCAACCAATCGTCGTCGGCCTCGCGATCAGGAAATCCACCAGGAGCCCCGCTATCCCAAGGCGAGGAGGCATACACCCGCCTGGCCACCGACACAATCGAGGAGCTGGACTGGTGCCTCGACCAGCTGGAGACCATCCAGACCCATCGCAGCGTCTCCGACATGGCGTCGCTTAAGTTCAAACGCATGCTCAACAAGGAGCTGTCGCACTTTAGCGAGTCCAGCAGATCGGGAAATCAGATTTCCGAATATATATGTTCCACATTTTTGGACAAGCAACAGGAGTTCGACTTGCCATCGCTGCGCGTGGAGGATAATCCCGAGCTGGTCGCCGCCAATGCAGCCGCTGGTCAACAGTCCGCTGGACAGTATGCACGCTCCCGATCGCCGCGCGGTCCGCCCATGTCGCAGATCAGCGGCGTGAAGAGACCGCTGTCGCATACGAATAGCTTTACCGGCGAACGTTTGCCCACCTTCGGTGTGGAGACACCCAGGGAGAATGAGCTGGGCACGCTGCTCGGCGAACTGGACACCTGGGGTATTCAGATATTCAGCATCGGCGAGTTCAGTGTCAATCGACCGCTCACCTGTGTGGCATACACCATATTTCAGAGTAGAGAATTACTGACCAGTCTTATGATACCACCGAAAACTTTTCTTAACTTTATGTCTACTCTGGAGGACCACTACGTCAAAGACAATCCGTTTCACAATTCGCTGCATGCCGCCGATGTGACACAAAGCACTAATGTGCTACTCAATACACCGGCGCTGGAGGGCGTATTCACGCCGCTCGAAGTGGGCGGGGCGCTGTTCGCCGCTTGCATACACGATGTTGATCATCCCGGCTTAACCAATCAGTTCTTGGTTAACTCAAGTTCCGAACTAGCATTAATGTACAATGACGAATCTGTTTTGGAAAATCATCATTTAGCTGTTGCCTTTAAATTATTACAAAATCAAGGATGTGATATATTCTGTAATATGCAAAAGAAACAACGCCAAACATTGAGGAAAATGGTTATTGATATTGTGCTGTCCACGGACATGTCCAAGCACATGAGTCTGCTGGCCGACCTAAAGACAATGGTGGAAACCAAAAAGGTGGCCGGCTCCGGAGTACTGCTGCTGGACAACTACACCGATCGCATACAGGTGCTTGAGAATCTGGTGCACTGCGCCGATCTGAGCAATCCCACCAAGCCGCTGCCGCTTTACAAGCGCTGGGTAGCCCTGCTCATGGAGGAGTTCTTCCTGCAGGGCGATAAGGAACGCGAATCGGGCATGGACATTAGTCCCATGTGCGATCGCCATAATGCCACCATTGAGAAGTCGCAGGTGGGCTTCATCGACTACATCGTCCACCCGCTATGGGAGACCTGGGCGGACCTGGTGCATCCGGATGCCCAGGATATACTCGACACGCTTGAAGAGAACAGAGACTACTACCAGAGCATGATACCGCCTTCGCCGCCGCCATCGGGCGTCGATGAGAATCCGCAGGAGGACAGGATACGCTTTCAAGTAACCCTTGAGGAATCCGATCAGGAGAACCTCGCCGAACTGGAGGAGGGCGACGAGAGTGGTGGCGAGAGCACCACCACAGGCACAACCGGAACCACCGCTGCATCCGCGCTAAGCGGAGCTGGTGGCGGTGGCGGTGGAGGCGGGGGAATGGCACCCAGAACGGGTGGCTGCCAAAACCAACCGCAACACGGTGGAATGTGACGGAGAGTCGTGGGAATTTATCGCAAATTACAGGAAAAGGCTCACAACTTTTTCCTATTACGTTAGTGACTACTATTAACACAGAAAAAACCAAAACAAAAACCTAAAACGAAAAGCAAAAAAAAAAACAAAGAAAAACTCTTGGAAAAAATGC

Citation:

Envoyé par Lolo78

Peut-on couper à n'importe quel endroit dans la séquence?

Oui on peut couper à n'importe quel endroit dans la séquence.
Petite rectification, des fragments peuvent se chevaucher.
J'ai également oublié de préciser que j'ai un nombre attendu de fragments (ici on peut prendre par exemple 10 fragments)

Par exemple :
ACAACAATAACAGCAGCAGCTGCGCAGCGACCGCCAGCGAGGAGCAAAACTCGAGCGGTACGCTCACCGGCAGCAGCGGCAATACTTCCACGTCGACGTGGTGCGGCAGCATCAACTTCCTGGTGT
Peut donner :
Code:

1 2 3 ACAACAATAACAGCAGCAGCTGCGCAGCG ------------------GCTGCGCAGCGACCGCCAGCGAGGA ---------------GCAGCTGCGCAGCGACCGCCAGCGAGGAGCA
Etc ...
La seule condition : la moyenne des longueurs des fragments générés doit être comprise entre 48 et 52.

En fait, si tu ne mets pas de contraintes "biologique", ne pas couper après A ou pas entre G et C... il n'y a acune raison pour que tes morceaux aient une taille variable, tu choisis juste leur taille et il te découpe tes morceaux.

Code:

1
2
3
4
5
6
7
8
9
10
 
use strict;
 
my $chaine="CCGCAGTCGCCGCCTCCGCAGCAGCAGCAGCAGCAGCAACAACTACAACAACAACAACAACAGCCACCAGGGCAGCAACAACCACACCAGTTACAACAGCCCCAACACCAACAGCGATAACAACAACATATAGATCGCGGAATCTGCTGAAAGCCACCGCCAAAGTGAAGAAGAAGAAAAAGTGGAAGCAAATCTCCTTCGTCGACGACCAGCAACGCTTGAAGCAAGCAGCGACGTCGGCGGAGCGCAAACGCTCCACCAGCAGCGCCCCGCCCGCCCTCGAGAGCAACGTACCGCCGCGGTTGGGGCGCTTTCTGCGGCAAAAGAGCAGCGAAACTAACAGCGCTGCAACGTCAGCAGCAGCGGCGGCAGTGCAAACGGAAGCGGAAACAGGACGCCGGAAATCGTCGGCGGCATCGATAATCGCAGCGGCAGCGGGCGCTGCTGGTTTGCACAACAATAACAGCAGCAGCTGCGCAGCGACCGCCAGCGAGGAGCAAAACTCGAGCGGTACGCTCACCGGCAGCAGCGGCAATACTTCCACGTCGACGTGGTGCGGCAGCATCAACTTCCTGGTGTATATGG"; 
my $pas = 50;
my $fin = length($chaine)-$pas-1; 
for (my $deb=0; $deb<=$fin; $deb++){
    my $morceau = substr($chaine, $deb, $pas);
    print "$morceau\n"
}

comme ça tu récupères tous les sous-chaines de longueur 50 dans ta chaine

03/11/2014, 17h01
Isabella83

Malheureusement c'est beaucoup plus compliqué que cela ...
Il faut couper aléatoirement, avec des tailles différentes (entre 16 nt et 50 nt) , et des séquences qui se chevauchent, comme dans mon exemple. et tout cela avec des fragments qui ont une moyenne de taille de 50 +-2nt.
j'ai commencé à faire ceci avec

Code:

nucleotide_cut = int(rand(length(ma_sequence))

Mais je n'ai pas encore fini. Je pense calculer la taille des séquences générées et l'ajouter à une variable, et choisir un nucléotide aléatoire dans ma séquence jusqu’à ce que la moyenne soit comprise entre 48 et 52 nt. Il faut aussi que je prenne en compte le nombre de fragments attendu mais je ne sais pas trop comment ...

Bonjour,

tu peux essayer ceci:

Code:

1
2
3
4
5
6
7
8
9
10
11
use strict;
use warnings;
use Data::Dumper;
 
my $len = shift;
 
my $seq = "ACAACAATAACAGCAGCAGCTGCGCAGCGACCGCCAGCGAGGAGCAAAACTCGAGCGGTACGCTCACCGGCAGCAGCGGCAATACTTCCACGTCGACGTGGTGCGGCAGCATCAACTTCCTGGTGT";
 
my @fragments = unpack ("(a$len)*" , $seq); # faire des fragments de taille $len
$fragments[-1] = $1 if ($seq =~ /(\w{$len})$/); # mettre le dernier fragment en chevauchement pour qu'il ait la taille $len
print Dumper \@fragments;

Ce qui donne:

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
$ perl fasta.pl 9
$VAR1 = [
          'ACAACAATA',
          'ACAGCAGCA',
          'GCTGCGCAG',
          'CGACCGCCA',
          'GCGAGGAGC',
          'AAAACTCGA',
          'GCGGTACGC',
          'TCACCGGCA',
          'GCAGCGGCA',
          'ATACTTCCA',
          'CGTCGACGT',
          'GGTGCGGCA',
          'GCATCAACT',
          'TCCTGGTGT'
        ];
 
$ perl fasta.pl 20
$VAR1 = [
          'ACAACAATAACAGCAGCAGC',
          'TGCGCAGCGACCGCCAGCGA',
          'GGAGCAAAACTCGAGCGGTA',
          'CGCTCACCGGCAGCAGCGGC',
          'AATACTTCCACGTCGACGTG',
          'GTGCGGCAGCATCAACTTCC',
          'CAGCATCAACTTCCTGGTGT'
        ];

EDIT: je n'avais pas vu les deux derniers posts quand j'ai posté ce qui précède. La solution ci-dessus ne fonctionne évidemment pas avec ces nouvelles contraintes. Plus trop le temps maintenant, je verrai plus tard si je peux trouver quelque chose correspondant à tes nouvelles contraintes.

03/11/2014, 17h28
Isabella83

Je vais voir si je peux adapter ceci à mon problème de départ.
03/11/2014, 19h32
Lolo78

Citation:

Envoyé par Isabella83

Il faut couper aléatoirement, avec des tailles différentes (entre 16 nt et 50 nt) , et des séquences qui se chevauchent, comme dans mon exemple. et tout cela avec des fragments qui ont une moyenne de taille de 50 +-2nt.

Cela me paraît contradictoire. si tu choisis des tailles aléatoires comprises entre 16 et 50 nt, la probabilité d'obtenir une moyenne de 51 ou 52 est nulle, et celle d'obtenir 48, 49 ou 50, sans être nulle, est proche de 0, surtout si ta séquence d'origine est un tant soit peu longue et mène à un nombre important de fragments.
04/11/2014, 10h15
Isabella83

Pensez vous que cela soit possible en enlevant la contrainte de taille ? des fragments de 1 nt à 50 nt ?
04/11/2014, 11h25
6ril23

il est surement possible de faire un $taille=rand(50) et un $deb=rand(length(tailleChaine))

puis de faire un substr($chaine, $deb, $taille); inclus dans une boucle afin de limiter le nombre de morceaux créés

Je regarde ça de plus près cette après midi

PS: juste par pure curiosité, pourrais-tu nous dire le but de ce découpage car j'ai du mal à en voir l'intérêt? (à part peut être pour tester de l'assemblage)
04/11/2014, 11h28
Isabella83

Je cherche à reproduire une RNAse in silico !
Les biologistes de mon équipe ont fait une manip, et la taille moyennes des fragments après séquençage fait 50 nt. Mais effectivement, il ne faut surtout pas mettre de contrainte seuil sur la taille des fragments, la seule contrainte, c'est la moyenne des tailles des fragments = 50 nt +/-2 nt !
C'est la raison pour laquelle je cherche à couper aléatoirement une séquence ( des séquences de transcripts) en respectant ces conditions ! Ensuite je sélectionnerai aléatoirement un certain nombre de ces fragments pour faire différents test et supprimer le bruit de fond de nos données.

Si j'ai bien compris tu cherches à faire quelque chose de taille 50 (moyenne dans l'expérience) mais avec une découpe aléatoire, ce qui n'est pas le cas de la RNAse afin d'enlever le bruit de l'aléatoire

Dans ce cas, je pense que si tu fais tous les "morceaux" de longueurs 50 mais avec peu de chevauchement

Code:

1
2
3
4
5
6
7
8
9
10
use strict;
 
my $chaine="CCGCAGTCGCCGCCTCCGCAGCAGCAGCAGCAGCAGCAACAACTACAACAACAACAACAACAGCCACCAGGGCAGCAACAACCACACCAGTTACAACAGCCCCAACACCAACAGCGATAACAACAACATATAGATCGCGGAATCTGCTGAAAGCCACCGCCAAAGTGAAGAAGAAGAAAAAGTGGAAGCAAATCTCCTTCGTCGACGACCAGCAACGCTTGAAGCAAGCAGCGACGTCGGCGGAGCGCAAACGCTCCACCAGCAGCGCCCCGCCCGCCCTCGAGAGCAACGTACCGCCGCGGTTGGGGCGCTTTCTGCGGCAAAAGAGCAGCGAAACTAACAGCGCTGCAACGTCAGCAGCAGCGGCGGCAGTGCAAACGGAAGCGGAAACAGGACGCCGGAAATCGTCGGCGGCATCGATAATCGCAGCGGCAGCGGGCGCTGCTGGTTTGCACAACAATAACAGCAGCAGCTGCGCAGCGACCGCCAGCGAGGAGCAAAACTCGAGCGGTACGCTCACCGGCAGCAGCGGCAATACTTCCACGTCGACGTGGTGCGGCAGCATCAACTTCCTGGTGTATATGG"; 
my $pas = 50;
my $chevauchement = 3;
my $fin = length($chaine)-$pas-1; 
for ( my $deb=0; $deb<=$fin; $deb+=($pas-$chevauchement) ){
    my $morceau = substr($chaine, $deb, $pas);
	print ">$morceau\n"
}

04/11/2014, 15h27
Isabella83

Non je ne peux pas prendre que des morceaux de taille 50 nt, ni utiliser un "pas" pour parcourir la séquence.
Il faut vraiment que la coupure soit aléatoire.
Je pense que je vais choisir des positions aléatoires, calculer la taille de fragments, puis la moyenne a chaque fois, et prendre des positions aléatoire tant que ma moyenne n'est pas de 50 +/- 2 nt.
Je suis en train d'essayer ....
04/11/2014, 15h55
__dardanos__
Salut à tous.
J’ai bien lu ce qui a été écrit. Moi, je procéderais de la façon suivante :

Tirer un fragment :
- Tirage aléatoire gaussien (moyenne = 50, sigma = 2) pour la largeur du fragment,
- Tirage aléatoire uniforme sur toute la largeur de la séquence pour la position du centre du fragment,
Continuer à ajouter un nouveau fragment tant qu’un critère de recouvrement n’est pas satisfait (nt délaissés < 5%, par exemple).

Les largeurs seront bien aléatoires. Mais le nombre de fragments sera élevé. Et il y aura des nt non utilisés.

ok.
Alors voilà une solution pour avoir des morceaux de taille aléatoire et de démarrage aléatoire:

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
use strict;
 
my $chaine="CCGCAGTCGCCGCCTCCGCAGCAGCAGCAGCAGCAGCAACAACTACAACAACAACAACAACAGCCACCAGGGCAGCAACAACCACACCAGTTACAACAGCCCCAACACCAACAGCGATAACAACAACATATAGATCGCGGAATCTGCTGAAAGCCACCGCCAAAGTGAAGAAGAAGAAAAAGTGGAAGCAAATCTCCTTCGTCGACGACCAGCAACGCTTGAAGCAAGCAGCGACGTCGGCGGAGCGCAAACGCTCCACCAGCAGCGCCCCGCCCGCCCTCGAGAGCAACGTACCGCCGCGGTTGGGGCGCTTTCTGCGGCAAAAGAGCAGCGAAACTAACAGCGCTGCAACGTCAGCAGCAGCGGCGGCAGTGCAAACGGAAGCGGAAACAGGACGCCGGAAATCGTCGGCGGCATCGATAATCGCAGCGGCAGCGGGCGCTGCTGGTTTGCACAACAATAACAGCAGCAGCTGCGCAGCGACCGCCAGCGAGGAGCAAAACTCGAGCGGTACGCTCACCGGCAGCAGCGGCAATACTTCCACGTCGACGTGGTGCGGCAGCATCAACTTCCTGGTGTATATGG"; 
my $pas = 50;
 
for ( my $i=0; $i<=30; $i++ ){           # a toi de choisir le nombre de séquences que tu souhaites générer
	my $taille=rand($pas);
	my $deb=rand(length($chaine))-1;
	while (  $taille+$deb >= length($chaine) ){
		$taille=rand(50);
		$deb=rand(length($chaine))-1;
	}
        my $morceau = substr($chaine, $deb, $taille);
	print ">$morceau\n"
}

Résultat

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
>perl test.pl
>GGACGCCGGAAATCGTCGGCGGCATCGATAAT
>CACGTC
>TCGATAATCGCAGCGGCAGC
>ACCGCCAAAGTGAAGAAGAAGAA
>AGTG
>CACACC
>CACCGGCAGCAGCG
>GCGCCCCGCCCGCCCTCG
>CAGGGCAGCAACAACCACACCAGT
>CAACAACAACAGCCACCAGGGCAG
>ACGTCGGCGGAGCGCAAACGCTCCACCAGCAGCGCCCCGCCCGCCCT
>CGACGTCGGCGGAGCGCAAA
>TCTCCTTCGTCGACGACCAGCAACGCTTGAAGCAAG
>CAATAACAGCAGCAGCTGCGCAGCGACCGCCAGCGAGGAGCAAAAC
>AAGCAGCGACGTCGGCGGAGCGCAAACGCTCCACCAGCAGCGCCCCGCC
>GCCAAAGTGAAGAAGAAGAAAAAGTGGAAGCAAATCTCCTTCGTCG
>CAACCACACCAGTTACAA
>TGGTGCGGCAGCATCAACTTCCTGGTGTATA
>TGGGGCGCTTT
>AGTGCAAACGGAAGCGGAAACAGGACGCCGGAAATCGTCGGCGGCATCG
>CACCAACAGCGATAACAACAACATATAGAT
>GGTTTGCACAACAATAACAGCAGCAGCTGCGCAGCGACCGC
>CGTCGACGTG
>ACAACAACAACAACAACAGC
>CAACAACTACAAC
>AACAGGACGCC
>GCGGAATCTGCTGAAAGCCACCGCCAAAGTG
>GCCGGAAATCGTCGGC
>GCAGCAGCA
>CGACCGCCAGCGAG
>GCAGCA

04/11/2014, 18h59
Lolo78

Isabella, tu ne semble pas comprendre pas le problème que j'ai expliqué précédemment (et Cyril non plus, semble-t-il), c'est un problème de mathématiques (de probabilité), pas de programmation informatique pour l'instant.

Si tu veux avoir seulement une chance que tes fragments aient une taille moyenne de l'ordre de 50, il ne faut pas choisir des tailles aléatoires comprises entre 1 et 50, ni même entre 16 et 50. Le moyenne ne pourra jamais être 50. Le problème est manifeste dans le cas du programme de Cyril: il obtient des fragments d'une longueur compprise entre 4 et 48 ou 49, il n'y a aucune chance que la moyenne soit proche de 50. En fait, cette moyenne, dans ce cas, doit être aux alentours de 25.

Pour avoir une chance d'obtenir une moyenne de l'ordre de 50, il faut que 50 soit (à peu près) au milieu de l'intervalle dans lequel tu choisis des longueurs aléatoires. Si su choisis des longueurs aléatoires entre 48 et 52, tu es certaine de satisfaire le critère de la moyenne. Si tu veut "étaler" plus l'intervalle de longueurs aléatoires, il faut par exemple choisir aléatoirement tes valeurs entre 40 et 60, ou 30 et 70, ou 1 et 100, etc. Plus tu choisis un intervalle étroit, plus tu auras de bonnes chances d'avoir la bonne moyenne, et plus tu choisis un intervalle grand, moins il sera probable d'obtenir la bonne moyenne, et plus tu devras multiplier les essais avant d'obtenir ce que tu cherches.

La correction minimale au programme de Cyril serait d'écrire:

Code:

my $taille=rand($pas) + int ($pas/2);

Il faut donc que tu résolves cette incohérence dans les contraintes que tu donnes avant même de songer à coder quelque chose.
04/11/2014, 19h21
Isabella83

Citation:

Envoyé par Lolo78

Si tu veux avoir seulement une chance que tes fragments aient une taille moyenne de l'ordre de 50, il ne faut pas choisir des tailles aléatoires comprises entre 1 et 50, ni même entre 16 et 50.

Oui j'ai rectifié ce point là dans un précédent post (#11).
04/11/2014, 19h53
Lolo78

Citation:

Envoyé par Isabella83

Je cherche à reproduire une RNAse in silico !
Les biologistes de mon équipe ont fait une manip, et la taille moyennes des fragments après séquençage fait 50 nt. Mais effectivement, il ne faut surtout pas mettre de contrainte seuil sur la taille des fragments, la seule contrainte, c'est la moyenne des tailles des fragments = 50 nt +/-2 nt !
C'est la raison pour laquelle je cherche à couper aléatoirement une séquence ( des séquences de transcripts) en respectant ces conditions ! Ensuite je sélectionnerai aléatoirement un certain nombre de ces fragments pour faire différents test et supprimer le bruit de fond de nos données.

Si je comprends bien tu te moques même que tes fragments couvrent absolument toute la séquence, il pourrait y avoir un (ou deux) nucléotide(s) manquant à une extrémité ou l'autre de la chaîne (voir au milieu, mais c'est moins probable) sans que ça gêne.

Du coup, je te suggères de choisir à chaque fois aléatoirement un pivot (rand sur la longueur de la séquence) et une longueur de chaîne (rand sur un intervalle centré sur 50). Tu divises la longueur N obtenue par 2 et prends les N/2 caractères précédant le pivot et les N/2 caractères suivants (en t'arrêtant naturellement aux extrémités de la séquence si le pivot est éloigné de cette extrémité de moins d'une demi-longueur.
05/11/2014, 09h55
Isabella83

Très bien,je vais essayer comme ceci et je reviendrai vers vous !
Merci beaucoup !
05/11/2014, 10h22
6ril23

Citation:

Envoyé par Lolo78

Isabella, tu ne semble pas comprendre pas le problème que j'ai expliqué précédemment (et Cyril non plus, semble-t-il), c'est un problème de mathématiques (de probabilité), pas de programmation informatique pour l'instant.

non j'avais juste lu tous les posts :mouarf::mouarf::mouarf::mouarf:

Citation:

Envoyé par Lolo78

Du coup, je te suggères de choisir à chaque fois aléatoirement un pivot (rand sur la longueur de la séquence) et une longueur de chaîne (rand sur un intervalle centré sur 50). Tu divises la longueur N obtenue par 2 et prends les N/2 caractères précédant le pivot et les N/2 caractères suivants (en t'arrêtant naturellement aux extrémités de la séquence si le pivot est éloigné de cette extrémité de moins d'une demi-longueur.

heu ça me parait pas foncièrement différent de ce que je fais: le désaventage de ta méthode, tu vas devoir tester que la sous séquence ne dépasse pas de chaque côté de la séquence.
05/11/2014, 11h51
vidici

je ne comprends pas vraiment l'énoncé du problème, mais s'il sagit d'essayer d'imiter le coté aléatoire d'un phénomène naturel, alors il faudrait couper selon une plage d'aléatoire et recommencer l'opération autant de fois que nécessaire, jusqu'a obtenir la moyenne voulue.
je pense
05/11/2014, 14h25
6ril23

Citation:

Envoyé par vidici

je ne comprends pas vraiment l'énoncé du problème, mais s'il sagit d'essayer d'imiter le coté aléatoire d'un phénomène naturel, alors il faudrait couper selon une plage d'aléatoire et recommencer l'opération autant de fois que nécessaire, jusqu'a obtenir la moyenne voulue.
je pense

en fait, si j'ai bien compris, Isabella cherche à créer des séquences aléatoires afin de pouvoir faire les même traitements que sur son jeu de données réelles.
Ainsi tu peux trouver les résultats sur/sous-exprimés par rapports à quelque chose d'aléatoire et donc potentiellement intéressant à étudier
05/11/2014, 19h07
Lolo78

Citation:

Envoyé par 6ril23

heu ça me parait pas foncièrement différent de ce que je fais: le désaventage de ta méthode, tu vas devoir tester que la sous séquence ne dépasse pas de chaque côté de la séquence.

Oui, j'avais mal lu ton code, ça ressemble sans doute. Sauf que pour la longueur du fragment tu fais un rand 50, ce qui donne un nombre aléatoire entre 0 et 50, et empêche que la moyenne soit de 50 (ça va générer une moyenne de l'ordre de 25). Il faudrait donc ajouter 25 au nombres aléatoires générés pour avoir les longueurs de fragments susceptibles d'avoir une moyenne proche de 25.

Citation:

Envoyé par 6ril23

en fait, si j'ai bien compris, Isabella cherche à créer des séquences aléatoires afin de pouvoir faire les même traitements que sur son jeu de données réelles.
Ainsi tu peux trouver les résultats sur/sous-exprimés par rapports à quelque chose d'aléatoire et donc potentiellement intéressant à étudier

mouvais..
en tout cas je me suis pris au jeu et brouillonner mon petit truc
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 #!/usr/bin/perl5.16 use strict; use warnings; my $chaine="CCGCAGTCGCCGCCTCCGCAGCAGCAGCAGCAGCAGCAACAACTACAACAACAACAACAACAGCCACCAGGGCAGCAACAACCACACCAGTTACAACAGCCCCAACACCAACAGCGATAACAACAACATATAGATCGCGGAATCTGCTGAAAGCCACCGCCAAAGTGAAGAAGAAGAAAAAGTGGAAGCAAATCTCCTTCGTCGACGACCAGCAACGCTTGAAGCAAGCAGCGACGTCGGCGGAGCGCAAACGCTCCACCAGCAGCGCCCCGCCCGCCCTCGAGAGCAACGTACCGCCGCGGTTGGGGCGCTTTCTGCGGCAAAAGAGCAGCGAAACTAACAGCGCTGCAACGTCAGCAGCAGCGGCGGCAGTGCAAACGGAAGCGGAAACAGGACGCCGGAAATCGTCGGCGGCATCGATAATCGCAGCGGCAGCGGGCGCTGCTGGTTTGCACAACAATAACAGCAGCAGCTGCGCAGCGACCGCCAGCGAGGAGCAAAACTCGAGCGGTACGCTCACCGGCAGCAGCGGCAATACTTCCACGTCGACGTGGTGCGGCAGCATCAACTTCCTGGTGTATATGG"; #To generate a random number between, for example, 100 and 150 my $range = 40; my $minimum = 30; my @hasard; my $dep=0; my $long_chaine = length($chaine); my $cpt = 0; my $moyen = 0; while ( 1 ){ my $random_number = int(rand($range)) + $minimum; # my $taille=$random_number; my $morceau = substr($chaine, $dep, $taille); $dep = $taille+$dep; print ">$morceau\n"; print ">$taille\n"; $cpt = $taille +$cpt; my $compte =$cpt; push(@hasard, $compte); if ($dep>$long_chaine) { #correction d une erreur #$moyen=$hasard[-1]/@hasard; #my $lmorceau=length($morceau); my $totalhasard = $hasard[-2]+length($morceau); $moyen=$totalhasard/@hasard; print "moyen >$moyen\n"; print "total chaine >$hasard[-1]\n"; print "long >@hasard\n"; my $longdernschaine = length($morceau); #print "morceau >$morceau\n"; if (($moyen > 48 && $moyen < 52) && ($longdernschaine > 30 && $longdernschaine < 60 )) { last; }else{ @hasard=(); $dep=0; $cpt = 0; $moyen = 0; } } }
ce qui donne par exemple
>CCGCAGTCGCCGCCTCCGCAGCAGCAGCAGCAGCA
>35
>GCAACAACTACAACAACAACAACAACAGCCACCAGGGCAGCAACAACCACACCAGTTAC
>59
>AACAGCCCCAACACCAACAGCGATAACAACAACATATAGATCGCGGAATCTG
>52
>CTGAAAGCCACCGCCAAAGTGAAGAAGAAGAAAAAGTGGAAGCAA
>45
>ATCTCCTTCGTCGACGACCAGCAACGCTTGAAGCAAGCAGCGACGTCGGCGGAGCG
>56
>CAAACGCTCCACCAGCAGCGCCCCGCCCGCCCTCGAGAGCAACGTACCGCCGCGG
>55
>TTGGGGCGCTTTCTGCGGCAAAAGAGCAGCGAAACTAACAGCGCTGCAACGTC
>53
>AGCAGCAGCGGCGGCAGTGCAAACGGAAGCGGAAACA
>37
>GGACGCCGGAAATCGTCGGCGGCATCGATAATCGCAGCGGCAGCGGGCGC
>50
>TGCTGGTTTGCACAACAATAACAGCAGCAGCTGCGCAGCGAC
>42
>CGCCAGCGAGGAGCAAAACTCGAGCGGTACGCTCACCGGCAG
>42
>CAGCGGCAATACTTCCACGTCGACGTGGTGCGGCAGCATCAACTTCCTGGTGTATATGG
>69
moyen >48.75

#au bout de la troisième boucle.
au bout d'un nombre un peu plus important de boucle

Salut.
Les solutions proposées fonctionnent, mais elles tiennent du bidouillage. Pour faire suite à mon premier post, je propose une solution :

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
#!/usr/bin/perl
# -*- coding:utf-8 -*-
 
use warnings;
use strict;
use locale;
use 5.010;
use Math::Random qw(random_uniform random_normal);
use List::Util qw(sum);
 
 
sub split_seq
{
    return -1 if (@_ != 1);
    (my $seq) = @_;
    my $N = length($seq);
 
    my $mean = 50;
    my $sigma = 2;
    my $full = 0.95;
 
    my $iter = 0;
    my $total = 0;
    my $occuped = 0;
    my @occupation ;
    for (my $i = 0; $i < $N; $i++) { push @occupation, 0; }
 
    while ($occuped / $N < $full)
    {
        my $start = int(random_uniform(1, 0, $N) + 0.5);
        my $width = int(random_normal(1, $mean, $sigma) + 0.5);
        my $finish = $start + $width;
        next if ($finish > $N);
        for (my $i = $start; $i < $finish; $i++) { $occupation[$i] = 1; }
        $occuped = sum(@occupation);
        my $snippet = substr($seq, $start, $width);
        $total += $width;
        say join(" ", ++($iter), $width, $snippet);
    }
    say "Moyenne : ", $total / $iter;
    return 0;
}
 
 
my $sequence =
"CCGCAGTCGCCGCCTCCGCAGCAGCAGCAGCAGCAGCAACAACTACAACAACAACAACAACAGCCACCAGGGCA" .
"GCAACAACCACACCAGTTACAACAGCCCCAACACCAACAGCGATAACAACAACATATAGATCGCGGAATCTGCT" .
"GAAAGCCACCGCCAAAGTGAAGAAGAAGAAAAAGTGGAAGCAAATCTCCTTCGTCGACGACCAGCAACGCTTGA" .
"AGCAAGCAGCGACGTCGGCGGAGCGCAAACGCTCCACCAGCAGCGCCCCGCCCGCCCTCGAGAGCAACGTACCG" .
"CCGCGGTTGGGGCGCTTTCTGCGGCAAAAGAGCAGCGAAACTAACAGCGCTGCAACGTCAGCAGCAGCGGCGGC" .
"AGTGCAAACGGAAGCGGAAACAGGACGCCGGAAATCGTCGGCGGCATCGATAATCGCAGCGGCAGCGGGCGCTG" .
"CTGGTTTGCACAACAATAACAGCAGCAGCTGCGCAGCGACCGCCAGCGAGGAGCAAAACTCGAGCGGTACGCTC" .
"ACCGGCAGCAGCGGCAATACTTCCACGTCGACGTGGTGCGGCAGCATCAACTTCCTGGTGTATATGG";
 
split_seq($sequence)

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
1 50 CAGCGAAACTAACAGCGCTGCAACGTCAGCAGCAGCGGCGGCAGTGCAAA
2 51 GAAGCAAGCAGCGACGTCGGCGGAGCGCAAACGCTCCACCAGCAGCGCCCC
3 50 GCAGCAGCAGCAGCAACAACTACAACAACAACAACAACAGCCACCAGGGC
4 48 AACAGCGATAACAACAACATATAGATCGCGGAATCTGCTGAAAGCCAC
5 51 CAGTTACAACAGCCCCAACACCAACAGCGATAACAACAACATATAGATCGC
6 51 TCGTCGGCGGCATCGATAATCGCAGCGGCAGCGGGCGCTGCTGGTTTGCAC
7 49 GCCACCGCCAAAGTGAAGAAGAAGAAAAAGTGGAAGCAAATCTCCTTCG
8 48 CCCTCGAGAGCAACGTACCGCCGCGGTTGGGGCGCTTTCTGCGGCAAA
9 49 AACGCTTGAAGCAAGCAGCGACGTCGGCGGAGCGCAAACGCTCCACCAG
10 51 CAGCAGCTGCGCAGCGACCGCCAGCGAGGAGCAAAACTCGAGCGGTACGCT
11 52 CGACGACCAGCAACGCTTGAAGCAAGCAGCGACGTCGGCGGAGCGCAAACGC
12 50 CAGCAGCGGCGGCAGTGCAAACGGAAGCGGAAACAGGACGCCGGAAATCG
13 53 CCGCAGCAGCAGCAGCAGCAGCAACAACTACAACAACAACAACAACAGCCACC
14 48 GGAGCAAAACTCGAGCGGTACGCTCACCGGCAGCAGCGGCAATACTTC
15 46 TGCGCAGCGACCGCCAGCGAGGAGCAAAACTCGAGCGGTACGCTCA
16 48 TTCGTCGACGACCAGCAACGCTTGAAGCAAGCAGCGACGTCGGCGGAG
17 54 CGCCCTCGAGAGCAACGTACCGCCGCGGTTGGGGCGCTTTCTGCGGCAAAAGAG
18 53 AAGAGCAGCGAAACTAACAGCGCTGCAACGTCAGCAGCAGCGGCGGCAGTGCA
19 47 AATAACAGCAGCAGCTGCGCAGCGACCGCCAGCGAGGAGCAAAACTC
20 49 TGGTTTGCACAACAATAACAGCAGCAGCTGCGCAGCGACCGCCAGCGAG
21 54 CAGCGGCAATACTTCCACGTCGACGTGGTGCGGCAGCATCAACTTCCTGGTGTA
22 50 AACTACAACAACAACAACAACAGCCACCAGGGCAGCAACAACCACACCAG
Moyenne : 50.0909090909091

08/11/2014, 16h28
Lolo78

La tentative est intéressante, mais je reste un peu perplexe, car rien n'indique dans ce qu'a dit Isabella que la taille des fragments doive suivre une distribution gaussienne (normale). Certes, cela facilite l'obtention d'une moyenne située dans les bornes voulues, surtout en prenant un delta de 2, mais ce qu'a dit Isabella faisait plutôt penser à une distribution aléatoire uniforme des tailles de fragments.

Même en supposant qu'il faille une distribution normale, un écart-type de 2 ne semble pas du tout correspondre à la description d'Isabella, si ce n'est qu'il garantit presque à coup sûr de tomber dans l'intervalle voulu.
09/11/2014, 15h16
vidici

Citation:

Envoyé par Lolo78

si ce n'est qu'il garantit presque à coup sûr de tomber dans l'intervalle voulu.

c'est sur. en partant avec une plage de 10 a 60. Il me faut 12000, 17509 etc .. boucles pour obtenir la fourchette voulue.
de 0 à 50 j'ai fait Controle C. :D