[Perl] Mélanger une chaîne de caractères de grande longueur

[MaliciaR]

Hello, tout le monde

J'ai une pitite question pour laquelle je n'ai pas trop d'idées actuellement ce soir (j'espère juste ce soir, sinon bobo ).
Je cherche à faire la chose suivante : prendre un génome (ie, une chaîne de caractères de grande longueurs > 5*10^6 caractères) et la mélanger, càd créer un "génome aléatoire" contenant le même nombre de chaque base se trouvant dans le génome de départ mais enchaînées totalement au hasard. Je n'arrive pas à trouver comment faire ça en Perl, j'ai vu que PHP avait une fonction appelée str_shuffle. Y aurait-il un truc équivalent en Perl...? Peut-être le List::Util avec la sous-routine shuffle...?

Sinon, je pense que ça va prendre un pitit moment comme traitement vu la longueur de la chaîne. J'ai pensé à faire autrement, mais je ne sais pas si c'est mieux : il s'agit de compter le nombre de chaque base présente dans le génome que je veux "aléatoriser" (c'est tout c*n à faire), puis faire une boucle for qui inclut un random; pour différentes valeurs du random, je vais demander qu'il choisisse l'une des 4 bases (ATGC). Pour chaque base ajoutée, il décrémentera son compteur respectif, jusqu'à ce que toutes les bases soient utilisées. C'est pas super dur à faire non plus...

Qu'en pensez-vous? Laquelle des deux solutions serait meilleure et si c'est la 1re, auriez-vous une idée de comment faire...?


Merci d'avance
-----

[Philou67]

Avec cet uni-ligne, tu constitues une chaine aléatoire de tes 4 bases.
La construction de la chaine ne prend que quelques secondes sur ma machine (Core 2 duo).

time perl -e '$genome = join "", map { ("A", "T", "G", "C")[int(rand(4))] } 1 .. 46e5'

real    0m10.403s
user    0m10.249s
sys     0m0.093s

La commande time sert juste à donner une idée du temps d'exécution.

En revanche, comme la chaine est en mémoire, sa taille est limitée par cette dernière. Dans mon contexte particulier, je n'ai pu dépasser la taille de 47E5.
En revanche, en écrivant la chaine dans un fichier, il ne devrait y avoir aucune limite.

time perl -e 'open my $GENOME, ">", "genome.txt";foreach (1 .. 5e6) { print $GENOME ("A", "T", "G", "C")[int(rand(4))] }'

real    0m4.430s
user    0m3.609s
sys     0m0.280s

Tout ceci à ré-écrire proprement naturellement.

Explications :

• map exécute le code entre { } pour chaque élément de la liste fourni en paramètre (1 .. 5E6).
• ("A", "T", "G", "C") constitue un tableau anonyme dont je prend un élément au hasard par l'opérateur [int(rand(4))]
• 1 .. 5E6 est un tableau constitué d'éléments numériques de 1 à 5E6 (opérateur .. dans un contexte de liste).

[MaliciaR]

Coucou,

Merci
Le truc est que ce n'est pas ce que je souhaite faire Càd qu'ici tu crées un génome aléatoire from scratch, alors que je souhaite en faire un avec une composition donnée (ie, un nombre donné de chaque base A, C, G, T) mais que ces bases s'enchaînent aléatoirement. Spour ça que je cherche soit à "mélanger" un génome déjà existant, soit à le recréer mais en gardant sa constitution.

Je fais mon idée de boucle for, on verra bien ce que ça donne

Des idées et conseils supplémentaires sont bienvenus
Merci!

[Philou67]

Tu connais le nombre de chaque base ?

[A voir en vidéo sur Futura]

[MaliciaR]

Je peux l'avoir pour chaque génome donné, oui, c'est tout bête à coder

[Philou67]

En prenant des hypothèses fixes de nombre de base pour le genome, on peut écrire ceci qui semble fonctionner :

#!/usr/bin/perl

use strict;
use warnings;

my %genome = ( A => 1E6, T => 1.5E6, G => 1.5E6, C => 1E6 );
my @gene   = sort keys %genome;

print "Liste des gênes : ", @gene, "\n";

sub nb_gene {
  my $nb_gene = 0;

  $nb_gene += $genome{$_} foreach keys %genome;

  print "Nombre de gêne du genôme : $nb_gene\n";
  return $nb_gene;
}

sub extraire_gene {
  my ($numero_gene) = @_;

  my $cumul_numero = 0;

  foreach my $gene (@gene) {
    if ($numero_gene < $cumul_numero + $genome{$gene}) {
      $genome{$gene}--;

      return $gene;
    }
    $cumul_numero += $genome{$gene};
  }
}

# Ouvrir le fichier du genome en ecriture
open my $GENOME, ">", "genome.txt" or die "Impossible d'ouvrir genome.txt ($!)\n";

# Ecrire le tableau en ordre aleatoire
my $nb_gene = nb_gene();
do {
  print $GENOME extraire_gene(int(rand(--$nb_gene)));
} until $nb_gene == 0;

En gros, on choisi aléatoirement le numéro d'un gêne parmi le nombre total de gêne restant à écrire ($nb_gene qui est décrémenté) en prenant soin à chaque fois que l'on a utilisé un gêne, de le supprimer du %genome ($genome{$gene}-- dans extraire_gene())

Pour vérifier le nombre de chaque base dans le génôme, voici un script un peu bourrin pour le connaitre :

#!/usr/bin/perl

use strict;
use warnings;

my @bases = ( 'A', 'T', 'G', 'C' );

foreach my $gene (@bases) {
  open my $GENOME, "<", "genome.txt" or die "Can't open genome.txt ($!)\n";
  print "$gene : ", scalar(@{[<$GENOME> =~ /$gene/g]}), "\n";
}

Je ne suis pas sûr d'employer les bons termes métiers au bon endroit, mais le code me semble répondre à ton besoin, et il est rapide (17s sur ma machine).

[Towl]

Un peu en retard, mais ca me parait extrement lent.
De mon point de vue, j'aurais plutot fait un code qui ressemble à cela (en python et codé à la porc, je connais pas le perl )

base = [A, T, G, C]
genome = {A: 1e6, T:1.5e6, G:1.5e6, C:1e6 )
genome_size = genome['A'] + genome['T'] + genome['G'] + genome['C'] 

file = open('output')
while genome_size != 0:
    g = base[rand(4)]
    if genome[g] == 0:
         continue
    file.print g
    genome[g] -= 1
    genome_size -= 1
file.close()

En gros l'idée c'est de prendre une base, de l'afficher et de décrementer son nombre d'occurence. Sur les boucles random comme cela, il est généralement moins couteux de faire plus de tour (ici on peut tomber sur une base qui n'est plus présente dans la liste à afficher) que de faire un moyen de contournement. Cela est vrai surtout quand il n'y a pas trop d'écart dans la distribution, comme cela est le cas dans ton exemple.

[MaliciaR]

Hello,

Merci pour le coup de main, mais ce n'est pas encore tout à fait ce que je voulais Cela dit, je trouve intréressant ce que tu as proposé et je vais l'utiliser pour certains autres tests

Voilà ce que j'ai fait moi (avec mes maigres connaissances, on ne se moque pas! ) :

#!/usr/bin/perl

#nom du programme : test-genome.pl 
#but : récupérer un génome, le découper en fragments de taille entre 800 et 1200 pb et rechercher des profiles de Tases dedans. Puis, le mélanger => #génome aléatoire et recommencer la même procédure. Chaque fileout est dans un nouveau fichier texte.


use strict;
use warnings;

###################################################
#Part 1: traitement d'un génome de référence (découpage en morceaux et test).


my $fich_genome = "/home/rayna/Desktop/metanosarcina.gen";
my $fragm_genome = "> /home/rayna/Desktop/metanoFragm.gen";

chomp $fich_genome;

open (FICH_GENOME, $fich_genome) or die ("Error 1: $!\n");
open (FRAGM_GENOME, $fragm_genome) or die ("Error 2: $!\n");


#Découpage en fragments de 1kb :
my $genome;
my $fragment = '';
my $last_fragment = '';
my $i = 0;

#Je le mets en une seule ligne : 
while ( my $ligne = <FICH_GENOME> ) {
	$ligne =~ s/\s//g;
	$genome .= $ligne;
}
for (my $position = 0; $position < length $genome; $position += 1000) {
	my $reste = (length $genome) - $position;
	if ( $reste > 1000 ) {
		$fragment = substr ($genome, $position, 1000);
		print FRAGM_GENOME ">fragm".$i."\n".$fragment."\n"; 
		$i++;
	} else {
		$last_fragment = substr ($genome, $position, $reste);
		print FRAGM_GENOME ">last\n".$last_fragment."\n";
	}
}

close FRAGM_GENOME;
close FICH_GENOME;


####################################################
#Part 2: dénombrement des types de bases dans le génome de référence, création d'un génome aléatoire, découpage et traitement.

#Compter le nombre de chaque type de bases dans le génome :
my ($nb_A, $nb_T, $nb_G, $nb_C, $nb_erreurs) = 0; 

$genome = uc($genome);

for ( my $each = 0; $each < length $genome; $each++ ) {
	my $base = substr($genome, $each, 1);
	if ($base eq "A") {
		$nb_A++;
	} elsif ($base eq "T") {
		$nb_T++;
	} elsif ($base eq "G") {
		$nb_G++;
	} elsif ($base eq "C") {
		$nb_C++;
	} else {
		print "Erreur : base inconnue\n\n";
		$nb_erreurs++;
	}
}
print "A = $nb_A\n";
print "T = $nb_T\n";
print "G = $nb_G\n";
print "C = $nb_C\n";



#Création du génome aléatoire :
my $random_genome = "/home/rayna/Desktop/randomMetano.gen";
open (RANDOM, "> $random_genome") or die ("Error 3: $!\n");

#D'abord, je stocke les bases dans une liste :
my @base_array = ("A", "T", "G", "C");
my $array_lentgh = @base_array; #je calcule la longueur de la liste

while (@base_array != 0 ) { #tant que la liste n'est pas vide
	my $random_number = int (rand(@base_array));
	my $base = $base_array[$random_number];
	print RANDOM $base_array[$random_number];
	if ($base eq "A") {
		$nb_A--;
		if ($nb_A == 0) {
			splice(@base_array, $random_number, 1); #supprimer l'élément A de la liste
		}
	} elsif ($base eq "T") {
		$nb_T--;
		if ($nb_T == 0) {
			splice(@base_array, $random_number, 1);
		}
	} elsif ($base eq "G") {
		$nb_G--;
		if ($nb_G == 0) {
			splice(@base_array, $random_number, 1);
		}
	} elsif ($base eq "C") {
		$nb_C--;
		if ($nb_C == 0) {
			splice(@base_array, $random_number, 1);
		}
	}
}


close RANDOM;

exit;

Le résultat de la Part 1 est un génome découpé en bouts de 1000 caractères + le dernier bout de < de 1000; celui de la Part 2 est un génome totalement aléatoire. Il faut juste que je m'amuse à pondérer les probas de tomber sur les bases qui sont fréquemment représentées (si j'ai un type de base représenté 2 fois plus qu'un autre, on est d'accord que ça va pas être super aléatoire à un moment donné ).

Conseils et suggestions d'améliorations bienvenus Merci!

[MaliciaR]

Je n'avais pas vu le message de Towl
Tu connais pas le perl? Chanceux!

Je regarderai comment ça fait et comment l'intégrer à quelque chose Merci!

[Philou67]

Towl, j'ai deux remarques à faire :
- la lenteur est toute relative : s'il faut 17s pour calculer un génome par jour, c'est tout à fait acceptable, en revanche, c'est pénalisant s'il faut calculer 100000 génomes par jour.
- le problème de l'algo que tu proposes est que la répartition n'est pas aléatoire (j'avais également codé en perl cet algo). En effet, lorsqu'une base est épuisée, le reste du génome ne contient plus cette base. Si la répartition des bases présente un grande disparité, les bases en faible quantité seront regroupées au début du génome. En procédant selon mon algorithme, qui est effectivement nettement plus long, on s'assure que toutes les bases sont réparties aléatoirement sur tout le génome.

[Philou67]

Malicia, pourrais-tu préciser le format de ton fichier génome ?
Il est constitué de lignes contenant des "\n" ? Avec combien de gêne pas ligne ? Si tu maitrises ce format, le plus simple est d'avoir une liste de caractère sans espaces ni retour à la ligne... que des gènes. La lecture de 1000 gène est alors extrêmement aisée :

$/ = \1000; # défini la taille d'un enregistrement en lecture
while (my $fragment = <FILE>) {
  # length $fragment <= 1000
}

Sinon, j'ai des doutes sur la possibilité de mettre dans une seule chaine, un génome très grand. Dans mon cas particulier, je n'ai pas pu créer de chaine de plus de 5E6 caractères (dépendant de la mémoire vive disponible au moment de l'exécution). Cela peut marcher un jour, et échouer un autre.

Pour la fragmentation, je te propose la simplification suivante :

...
my $fich_genome = "/home/rayna/Desktop/metanosarcina.gen";
my $fragm_genome = "> /home/rayna/Desktop/metanoFragm.gen";

chomp $fich_genome;

open (FICH_GENOME, $fich_genome) or die ("Error 1: $!\n");
open (FRAGM_GENOME, $fragm_genome) or die ("Error 2: $!\n");

#Découpage en fragments de 1kb :
my $genome;
my $i = 0;

#Je le mets en une seule ligne : 
while ( my $ligne = <FICH_GENOME> ) {
	$ligne =~ s/\s//g;
	$genome .= $ligne;
}
# s'écrit aussi :
# $genome = join "", map { s/\s//g; $_ } <FICH_GENOME>;

my $i = 0;
map { print FRAGM_GENOME ">fragm",$i++,"\n",$_, "\n" } $genome =~ /.{,1000}/g;
# s'écrit aussi :
# foreach my $fragment ($genome =~ /.{0,1000}/g) {
#   print FRAGM_GENOME ">fragm", $i++, "\n", $fragment, "\n"
# }

close FRAGM_GENOME;
close FICH_GENOME;

J'ai grossièrement testé, et cela semble correct. La seule chose que ne fait pas mon code, c'est de mettre "last" à la place de "fragm" pour le dernier fragment. Si c'est vraiment important, il convient alors de tester $i par rapport à (length $genome / 1000) pour afficher le texte adéquat, mais je doute que cela ait un intérêt réel.

Je regarde le reste de ton programme.

[Philou67]

Pour le dénombrement, il est souvent utile d'utiliser les tables de hachages, qui permettent un accès direct à une valeur nommée, plutôt que une suite de test (parfois interminable).
Je te propose la chose suivante :

#Part 2: dénombrement des types de bases dans le génome de référence, création d'un génome aléatoire, découpage et traitement.

#Compter le nombre de chaque type de bases dans le génome :
my %nb = { "A" => 0, "T" => 0, "G" => 0, "C" => 0 };
my $nb_erreur = 0;

foreach my $gene (split //, uc $genome) {
  if (exists $nb{$gene}) {
    $nb{$gene}++;
  }
  else {
    print "Erreur : [$gene] base inconnue\n\n";
    $nb_erreur++;
  }

print "$_ = $nb{$_}" foreach keys %nb;

[Philou67]

Pour la partie de génération aléatoire du génome, j'ai la même remarque que celle pour l'algorithme de Towl, la répartition n'est pas homogène, notamment si la disparité des bases est importante.

Il me semble préférable pour cette partie de l'algorithme, que tu utilises celui que j'ai fourni, qui s'adaptera très facilement à ton cas, surtout si tu utilises la table de hachage pour le nombre de base que j'ai présenté dans mon post précédent.

[MaliciaR]

Merci pour tes remarques Je regarde ça aujourd'hui.

Pour la partie de génération aléatoire du génome, j'ai la même remarque que celle pour l'algorithme de Towl, la répartition n'est pas homogène, notamment si la disparité des bases est importante.

Il me semble préférable pour cette partie de l'algorithme, que tu utilises celui que j'ai fourni, qui s'adaptera très facilement à ton cas, surtout si tu utilises la table de hachage pour le nombre de base que j'ai présenté dans mon post précédent.

Oui, je me suis rendue compte que question homogénéité c'était pas encore ça quand j'ai testé mon script hier Et quand je t'avais dit que ce n'était pas tout à fait ça, je me suis fiée aux noms de tes variables : du coup, pour moi, tu prenais les gènes (des sous-strings du génome, quoi) pour ton algo. C'est après que je me suis rendue compte que c'était les bases dont tu parlais. Faut pas oublier, je suis intelligente de temps à autre et hier j'avais déjà donné
Donc, j'ai revu les choses et j'ai un algo très similaire au tien, mais sans hash Z'aime pas les hash... On verra ce que ça donnera quand je l'aurai fait

[Philou67]

Malicia, j'avais bien précisé :

Je ne suis pas sûr d'employer les bons termes métiers au bon endroit, mais le code me semble répondre à ton besoin, et il est rapide (17s sur ma machine).

Pour les hash... ce rebut vient de quelle appréhension ?
Selon moi, c'est une des bases de la programmation en perl.

[MaliciaR]

Malicia, j'avais bien précisé

Oui, oui, je sais bien, c'est qu'en ce moment je ne sais plus trop où donner la tête alors parfois, je rate des trucs clairs...

Pour les hash... ce rebut vient de quelle appréhension ?
Selon moi, c'est une des bases de la programmation en perl.

Beh ça vient surtout d'un manque total de connaissances J'ai beaucoup de mal du coup pour les gérer, même si je comprends le principe. Mais j'essaie un max de les éviter parce que ça me prend des plombes pour dérouler l'algo et voir tout ce que ça fait. Et comme il faut absolument que j'avance vite, beh voilà... Tu as ta réponse