maximum d'une entité dans un fichier [Résolu]

[Isabella83]

Bonjour à tous,

J'ai un fichier de ce type :

>rho-RA
ATGCGTAGC
>rho-RB
ATGCGTAGCATGCGTAGC
>rho-RC
ATGCGTAGCATGCGTAGCATGCGTAGCATGCGTAGC
>pax6-RA
GCTGCATGATAG
>pax6-RB
GCTGCATGATAGGCTGCATGATAG

donc, je veux en premier temps determiner la taille de chacune des séquences, j'ai donc écrit cela :

Code :

1
2
3
4
5
6
7
8
9
10
11
12
13
my $fichier = 'test.fa';
open my $fh, '<', $fichier or die "Impossible de lire le fichier $fichier\n";
while ( my $ligne = <$fh> ) {
	chomp $ligne;
	if ($ligne =~ m{^>} ) {
    		print "$ligne\t";
	}
	else{
		my $length = length($ligne);
		print "$length\n";
	}
}
close $fichier;

Mais j'aimerai qu'il ne m'affiche pour chaque entité (ie dans l'exemple rho et pax6) seulement le transcript (RA,RB...) qui a la plus grande sequence, dans l'exemple rho-RC et pax6-RB

je ne sais pas comment faire, quelqu'un pourrait m'aider ?
merci d'avance,
et bonne fête de fin d'année

[Beniou]

Bonjour,

Le plus simple est d'utiliser des tables de hachages avec comme clé le nome de l'entité (rho et pax6 dans ton exemple) pour sauvegarder le plus long transcrit et sa longueur.

Quand on arrive sur une ligne commençant par '>', on sauvegarde l'entité et le transcrit puis on calcule la longueur de sa séquence sur la ligne suivante.

Voici un exemple qui fonctionne mais qui ne fait aucun contrôle sur le format du fichier : s'il y a une ligne vide ou bien si la séquence est sur plusieurs lignes cela ne fonctionnera pas...

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
use strict;
use warnings;
 
my $fasta_file = "sequences.fa";
 
my $current_entity = q{};
my $current_transcript = q{};
 
my %max_length_entity;
my %max_transcript;
 
open(my $FASTA, '<', $fasta_file) or die "$fasta_file : $!\n\n";
while(<$FASTA>){
    chomp;
    if(/^>(.+)[-](.+)/){
        $current_entity = $1;
        $current_transcript = $2;
    }
    else{
        my $length_transcript = length;
        if(!exists($max_length_entity{$current_entity}) ||
           $length_transcript > $max_length_entity{$current_entity}){
            $max_length_entity{$current_entity} = $length_transcript;
            $max_transcript{$current_entity} = $current_transcript;
        }
    }
 
}
close($FASTA);
 
foreach my $entity(sort(keys(%max_length_entity))){
    print "Entity : $entity\n";
    print "Longest transcript (length) : $max_transcript{$entity} ($max_length_entity{$entity})\n\n";
}

[Isabella83]

Bonjour et Bonne année !
merci pour votre aide.
Je débute avec la programmation, et j'aimerai comprendre les lignes suivantes du code que vous avez écrit :

Envoyé par Beniou

my $current_entity = q{};
my $current_transcript = q{};

je ne comprends pas le "=q{} "

Merci d'avance pour vos explications

[stoyak]

Si tu peux utiliser un module, je te conseille de faire un petit tour sur le CPAN. Un petit exemple avec Bio::SeqIO. Tu pourras très facilement parser un fichier Fasta et récupérer les informations.

[Beniou]

Bonjour,

Envoyé par Isabella83

je ne comprends pas le "=q{} "

Cela initialise la variable avec une chaîne vide. On peut dire que cela est équivalent à la ligne ci dessous mais la notation est plus propre (à mon goût ):

Code :

my $var = "";

Envoyé par Stoyak

Si tu peux utiliser un module, je te conseille de faire un petit tour sur le CPAN. Un petit exemple avec Bio::SeqIO. Tu pourras très facilement parser un fichier Fasta et récupérer les informations.

J'avais déjà émis l'idée d'utiliser ce module dans une autre discussion mais apparemment, Isabella83 n'a pas choisi cette méthode.

A réfléchir quand même si il y a beaucoup de manipulation de fichiers fasta pour obtenir des informations basiques (longueur séquences, accès facilement à chaque séquence dans un multifasta, lancer des blasts ? etc.)

[Isabella83]

Merci pour vos réponses, je vais regarder en détail Bio::SeqIO