problème de comparaison d'association de deux données entre deux fichiers

**benlh** · 18/02/2009, 17h30

Bonjour,
J'essaie désespérément de vérifier si une association de deux chaines de caractères est correcte par rapport à une association originelle. Je m'explique :
j'avais à l'origine un fichier (tad delimited txt) avec plusieurs colonnes de données dont deux m'intéressaient. Une colonne "oligo" et une colonne "contig" qui définissaient une association d'un oligo à son contig. Pour vérifier la spécificité de cette association, j'ai extrait les listes d'oligos et de contigs pour ensuite générer deux banques de séquences à comparer par blastn. Je récupère alors un fichier csv contenant les meilleurs hits pour chaque oligo. Dans ce fichier j'ai bien une colonne "oligo", une colonne "contig" et d'autres colonnes sur les qualités d'alignements. Je voudrais maintenant vérifier si ces "nouvelles" associations oligo/contig sont les mêmes qu'au départ.

Pour cela j'ai voulu faire une comparaison de clés de table de hachage pour générer un fichier identique au fichier csv de résultats des blast mais avec une colonne de plus qui indique une bonne (OK) ou une mauvaise (NA) correspondance oligo/contig par rapport aux données d'origine.
Finalement mon idée ne marche pas, ça a l'air de déconner au niveau du dernier test car je me retrouve avec le bon fichier de sortie mais avec que des OK ou que des NA.

Si quelqu'un peu m'aider je lui en serai reconnaissant.
Je joins mon script ci dessous :

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
use strict;
 
my ($lineB,$contig,$oligo,$lineA,$query_name,$query_accession,$query_description,$query_length,$hit_name,$hit_accession,$hit_description,$hit_length,$score,$evalue,$hsp_number,$hsp_query_start,$hsp_query_end,$hsp_hit_start,$hsp_hit_end,$hsp_percent_identity,$direction);
my (@elmtsA,@elmtsB);
my %contig2oligo;
 
open OLIGO, "< oligo_otherV6.txt" or
	die "Unable to open input file #2 : $!\n";
while ($lineB=<OLIGO>) {
	chomp $lineB;
	if ($lineB =~ /^#/) {#Header
	    next;
  }
	@elmtsB=split "\t", $lineB;
	$contig=$elmtsB[0];
	$oligo=$elmtsB[4];
	if (defined $contig2oligo{$contig}) {
		$contig2oligo{$contig}=$contig2oligo{$contig}."\t$oligo";
		#$agilent2Oligos{$agilent}.="\t$oligo";
	} else {
		$contig2oligo{$contig}=$oligo;
	}
}
close(OLIGO);
 
open OUT, "> oligo_contig_perfect.txt" or
	die "Unable to open oligo_contig_perfect file : $!\n";
 
open CONTIG, "< blastn-oligo_restrict_otherV6-contig_restrict_otherV6-1e-3-1e-3-1,1.csv" or
	die "Unable to open input file #1 : $!\n";
my $noLine=0;
while ($lineA=<CONTIG>) {
	$noLine++;
	chomp $lineA;
	if ($lineA =~/^#/ || $noLine==1) {#Header
	    print OUT "query_name\tquery_accession\tquery_description\tquery_length\thit_name\thit_accession\thit_description\thit_length\tscore\tevalue\thsp_number\thsp_query_start\thsp_query_end\thsp_hit_start\thsp_hit_end\thsp_percent_identity\tdirection\tmatch\n";
 
	    next;
  }
	@elmtsA=split ";", $lineA;
	$query_name=$elmtsA[0];
  $query_accession=$elmtsA[1];
  $query_description=$elmtsA[2];
  $query_length=$elmtsA[3];
  $hit_name=$elmtsA[4];
  $hit_accession=$elmtsA[5];
  $hit_description=$elmtsA[6];
  $hit_length=$elmtsA[7];
  $score=$elmtsA[8];
  $evalue=$elmtsA[9];
  $hsp_number=$elmtsA[10];
  $hsp_query_start=$elmtsA[11];
  $hsp_query_end=$elmtsA[12];
  $hsp_hit_start=$elmtsA[13];
  $hsp_hit_end=$elmtsA[14];
  $hsp_percent_identity=$elmtsA[15];
  $direction=$elmtsA[16];
 
 
	if (defined $contig2oligo{$hit_name}) {
 
		$contig2oligo{$hit_name}=$contig2oligo{$hit_name}."\t$query_name";
		#$agilent2Oligos{$agilent}.="\t$oligo";
	} else {
		$contig2oligo{$hit_name}=$query_name;
	}	
  if ($contig2oligo{$contig}==$contig2oligo{$hit_name}){
    print OUT "$query_name\t$query_accession\t$query_description\t$query_length\t$hit_name\t$hit_accession\t$hit_description\t$hit_length\t$score\t$evalue\t$hsp_number\t$hsp_query_start\t$hsp_query_end\t$hsp_hit_start\t$hsp_hit_end\t$hsp_percent_identity\t$direction\tOK\n";
	} else {
      print OUT "$query_name\t$query_accession\t$query_description\t$query_length\t$hit_name\t$hit_accession\t$hit_description\t$hit_length\t$score\t$evalue\t$hsp_number\t$hsp_query_start\t$hsp_query_end\t$hsp_hit_start\t$hsp_hit_end\t$hsp_percent_identity\t$direction\tNA\n";
  }
}  
close(CONTIG);
close(OUT);

Merci par avance.

Benoit

**benlh** · 19/02/2009, 10h10

je précise, j'ai un switch de résultats entre OK et NA quand je change l'opérateur d'égalité à la fin du script :
== que des OK
eq que des NA

**benlh** · 19/02/2009, 14h43

ça y est mon test fonctionne mais j'ai un nouveau problème qui est que cette méthode me génère énormément de doublons. En gros, je passe de 2000 lignes de données à plus de 6 000 000...
voilà le nouveau script si quelqu'un à une idée je suis preneur:

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
use strict;
 
my ($lineB,$contig,$oligo,$lineA,$query_name,$query_accession,$query_description,$query_length,$hit_name,$hit_accession,$hit_description,$hit_length,$score,$evalue,$hsp_number,$hsp_query_start,$hsp_query_end,$hsp_hit_start,$hsp_hit_end,$hsp_percent_identity,$direction);
my (@elmtsA,@elmtsB);
my %contig2oligo;
 
open OUT, "> oligo_contig_perfect.txt" or
	die "Unable to open oligo_contig_perfect file : $!\n";
 
open OLIGO, "< oligo_otherV6.txt" or
	die "Unable to open input file #2 : $!\n";
while ($lineB=<OLIGO>) {
	chomp $lineB;
	if ($lineB =~ /^#/) {#Header
	    next;
  }
	@elmtsB=split "\t", $lineB;
	$contig=$elmtsB[0];
	$oligo=$elmtsB[4];
	if (defined $contig2oligo{$contig}) {
		$contig2oligo{$contig}=$contig2oligo{$contig}."\t$oligo";
		#$agilent2Oligos{$agilent}.="\t$oligo";
	} else {
		$contig2oligo{$contig}=$oligo;
	}
 
open CONTIG, "< blastn-oligo_restrict_otherV6-contig_restrict_otherV6-1e-3-1e-3-1,1.csv" or
	die "Unable to open input file #1 : $!\n";
my $noLine=0;
while ($lineA=<CONTIG>) {
	$noLine++;
	chomp $lineA;
	if ($lineA =~/^#/ || $noLine==1) {#Header
	    ##print OUT "query_name\tquery_accession\tquery_description\tquery_length\thit_name\thit_accession\thit_description\thit_length\tscore\tevalue\thsp_number\thsp_query_start\thsp_query_end\thsp_hit_start\thsp_hit_end\thsp_percent_identity\tdirection\tmatch\n";
 
	    next;
  }
	@elmtsA=split ";", $lineA;
	$query_name=$elmtsA[0];
  $query_accession=$elmtsA[1];
  $query_description=$elmtsA[2];
  $query_length=$elmtsA[3];
  $hit_name=$elmtsA[4];
  $hit_accession=$elmtsA[5];
  $hit_description=$elmtsA[6];
  $hit_length=$elmtsA[7];
  $score=$elmtsA[8];
  $evalue=$elmtsA[9];
  $hsp_number=$elmtsA[10];
  $hsp_query_start=$elmtsA[11];
  $hsp_query_end=$elmtsA[12];
  $hsp_hit_start=$elmtsA[13];
  $hsp_hit_end=$elmtsA[14];
  $hsp_percent_identity=$elmtsA[15];
  $direction=$elmtsA[16];
 
 
	if (defined $contig2oligo{$hit_name}) {
 
		$contig2oligo{$hit_name}=$contig2oligo{$hit_name}."\t$query_name";
	} else {
		$contig2oligo{$hit_name}=$query_name;
	}	
 
  if ($contig2oligo{$contig} eq $contig2oligo{$hit_name}){
   print OUT "$query_name\t$hit_name\tOK\n";
	} else {
      print OUT "$query_name\t$hit_name\tNA\n";
  }
}
 
}
close(OLIGO);
close(CONTIG);
close(OUT);

**benlh** · 19/02/2009, 15h54

finalement, fausse alerte, ça me génère des faux positifs...
l'approche de comparaison des clés n'est peut être pas la bonne solution mais pour l'instant je n'ai pas d'autres idées. Pour l'instant je tâtonne en perl, ça ne fait qu'un mois que je m'y suis mis. Si quelqu'un pouvait m'aiguiller un peu ça m'aiderait.
Merci d'avance.

**Jasmine80** · 20/02/2009, 15h18

... sans commentaires dans ton script, c'est difficile de comprendre ce que tu veux faire. Pourrais-tu en rajouter, en expliquant bien ce que contient chaque fichier ainsi que ce que tu veux faire. Utilise les balises CODE afin de rendre ton script lisible.

Pour vérifier la spécificité de cette association, j'ai extrait les listes d'oligos et de contigs pour ensuite générer deux banques de séquences à comparer par blastn.

Fais-tu un blast en ligne (module RemoteBlast) ou en local (module StandaloneBlast)?

Qu'entends-tu par générer 2 banques de séquences? La base de données de blast est-elle générée à partir de tes séquences ou est-elle en ligne? Je suppose que c'est en local, que tu blastes une à une tes séquences d'oligo versus une base de données blast crée à partir de tes séquences de contig.

Puis tu récupères donc 1 fichier de résultat par séquence d'oligo? (nb il existe le module BoulderBlast permettant de lire un fichier output.blast sans devoir te casser la tête)

Je récupère alors un fichier csv contenant les meilleurs hits pour chaque oligo. Dans ce fichier j'ai bien une colonne "oligo", une colonne "contig" et d'autres colonnes sur les qualités d'alignements. Je voudrais maintenant vérifier si ces "nouvelles" associations oligo/contig sont les mêmes qu'au départ.

Pourquoi utiliser un fichier csv? Pourrais-tu montrer un exemple de lignes de ce fichier?

Ce que je ferais.
1) créer un hash $hash{oligo} = contig; # à condition que oligo soit unique
2) faire un blast et récupérer le résultat dans un fichier output.blast
3) utiliser RemoteBlast afin de lire le résultat facilement et comparer avec ce que tu as dans $hash{oligo}

**benlh** · 23/02/2009, 12h19

désolé pour le manque de clarté de mon script, je débute n'ai pas l'habitude de mettre des commentaires mais vais remédier à la situation et détailler un plus ce que j'essaie de faire.

Pour commencer, voilà en résumé mon projet :
je veux modifier une micropuce à oligos de 15K construite à partir d'EST. Cette modification est réalisé car nous possédons au labo des banques spécifiques de notre étude et voulons enrichir la puce avec ces mêmes séquences. Le problème c'est que pour des raisons techniques, on doit enlever de la puce autant d'oligos qu'on en rajoute. Autre difficulté, c'est que nos banques sont constituées de contig d'EST. De ce fait sur la puce initiale un contig peu être représenté par 0, 1 ou plusieurs EST sur lesquelles ont été déssinés les oligos.

Ma démarche a donc été la suivante :

associer chaque contig à son/ses oligo(s) par la présence d'EST composant le contig sur la puce en s'appuyant sur les correspondance ID_genbank/ID_oligo du fichier de design de la puce initiale
vérification des associations contig/oligo et de leur spécificité

j'en suis donc à la 2nd étape.
J'ai en ma possession la liste des association ID_oligo/ID_contig (faites par le biais des ID_genbank). J'ai extrait les liste des ID_oligo et des ID_contig séparément pour pouvoir générer deux banques distinctes "oligos" et "contigs" en récupérant les séquences associées à chaque ID_oligo et ID_contig.
Ensuite, j'ai blasté ("dseablast" en local) la banque oligos contre la banque contigs afin de vérifier la spécificité des oligos. Le dseablast est un blast installé en local à l'INRA qui réalise l'ensemble des alignement, fait un tri selon un cutoff sur la e-value et un parse des hits significatifs. Je me suis arranger pour paramétrer des conditions relativement stringentes. Au final, les résultats significatifs des hits sont consignés dans un fichier csv contenant entre autre l'ID_oligo, l'ID_contig et les renseignements relatifs à l'alignement.
Ce que je voudrait faire maintenant c'est vérifier si les correspondances ID_oligo/ID_contig issues du blast sont les mêmes que celles générées lors de la première étape d'association conyig/oligo.
D'où mon script que je vais maintenant détailler :

première étape : lecture et enregistrement des associations originales ID_oligo/ID_contig

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
 
use strict;
#déclaration des variables à extraire des deux fichiers
my ($lineB,$contig,$oligo,$lineA,$query_name,$query_accession,$query_description,$query_length,$hit_name,$hit_accession,$hit_description,$hit_length,$score,$evalue,$hsp_number,$hsp_query_start,$hsp_query_end,$hsp_hit_start,$hsp_hit_end,$hsp_percent_identity,$direction);
#déclaration de tableaux pour spliter les lignes de données
my (@elmtsA,@elmtsB);
#déclaration de la table de hachage pour consigner les associations ID_oligo/ID_contig avaec en clé les ID_contig
my %contig2oligo;
 
#création et ouverture du fichier pour consigner les résultats
open OUT, "> oligo_contig_perfect.txt" or
    die "Unable to open oligo_contig_perfect file : $!\n";
 
#ouverture du fichier d'associations ID_oligo/ID_contig par le biais des ID_genbank
open OLIGO, "< oligo_otherV6.txt" or
    die "Unable to open input file #2 : $!\n";
#extraction des données ligne par ligne en sautant le header
while ($lineB=<OLIGO>) {
    chomp $lineB;
    if ($lineB =~ /^#/) {#Header
       next;
}
#splitage de chaque ligne dans le tableau elmtsB et mise en mémoire uniquement des ID_oligo et ID_contig dans les variables $coligo et $contig
    @elmtsB=split "\t", $lineB;
    $contig=$elmtsB[0];
    $oligo=$elmtsB[4];
#mémorisation dans la table de hachage des associassions ID_contig/ID_oligo sous les clés ID_contig
    if (defined $contig2oligo{$contig}) {
#plusieurs oligos peuvent être associés au même contig
       $contig2oligo{$contig}=$contig2oligo{$contig}."\t$oligo";
    } else {
       $contig2oligo{$contig}=$oligo;
}

seconde étape : lecture et enregistrement des associations ID_oligo/ID_contig après blast

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
 
#ouverture du fichier de résultats des blasts
open CONTIG, "< blastn-oligo_restrict_otherV6-contig_restrict_otherV6-1e-3-1e-3-1,1.csv" or
    die "Unable to open input file #1 : $!\n";
#extraction des données ligne par ligne en sautant le header et en créant le header du fichier OUT
my $noLine=0;
while ($lineA=<CONTIG>) {
    $noLine++;
    chomp $lineA;
    if ($lineA =~/^#/ || $noLine==1) {#Header
       print OUT "query_name\tquery_accession\tquery_description\tquery_length\thit_name\thit_accession\thit_description\thit_length\tscore\tevalue\thsp_number\thsp_query_start\thsp_query_end\thsp_hit_start\thsp_hit_end\thsp_percent_identity\tdirection\tmatch\n";
 
    next;
}
#splitage des lignes et extraction des données à mémoriser
    @elmtsA=split ";", $lineA;
    $query_name=$elmtsA[0]; #=ID_oligo
    $query_accession=$elmtsA[1];#=ID_oligo
    $query_description=$elmtsA[2];
    $query_length=$elmtsA[3];
    $hit_name=$elmtsA[4];#=ID_contig
    $hit_accession=$elmtsA[5];#=ID_contig
    $hit_description=$elmtsA[6];
    $hit_length=$elmtsA[7];
    $score=$elmtsA[8];
    $evalue=$elmtsA[9];
    $hsp_number=$elmtsA[10];
    $hsp_query_start=$elmtsA[11];
    $hsp_query_end=$elmtsA[12];
    $hsp_hit_start=$elmtsA[13];
    $hsp_hit_end=$elmtsA[14];
    $hsp_percent_identity=$elmtsA[15];
    $direction=$elmtsA[16];
 
#enregistrement des associations ID_oligo($query_name)/ID_contig ($hit_name) dans la table de hachage sous la clé ID_contig ($hit_name)
    if (defined $contig2oligo{$hit_name}) {
#plusieurs ID_oligos peuvent être associés à un même oligo
       $contig2oligo{$hit_name}=$contig2oligo{$hit_name}."\t$query_name";
    } else {
       $contig2oligo{$hit_name}=$query_name;
}

et enfin le test pour savoir si les associations avant blast et après blast sont les mêmes

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
 
if ($contig2oligo{$contig} eq $contig2oligo{$hit_name}){
    print OUT "$query_name\t$query_accession\t$query_description\t$query_length\t$hit_name\t$hit_accession\t$hit_description\t$hit_length\t$score\t$evalue\t$hsp_number\t$hsp_query_start\t$hsp_query_end\t$hsp_hit_start\t$hsp_hit_end\t$hsp_percent_identity\t$direction\tOK\n";
    } else {
       print OUT "$query_name\t$query_accession\t$query_description\t$query_length\t$hit_name\t$hit_accession\t$hit_description\t$hit_length\t$score\t$evalue\t$hsp_number\t$hsp_query_start\t$hsp_query_end\t$hsp_hit_start\t$hsp_hit_end\t$hsp_percent_identity\t$direction\tNA\n";
    }
}
 
}
close(OLIGO);
close(CONTIG);
close(OUT);

voilà, j'espère que ma démarche est plus claire maintenant.
Peut être n'est ce pas la méthode la plus simple. Merci d'avance pour vos suggestions et votre aide.

**Jasmine80** · 23/02/2009, 13h40

Oui, cela est beaucoup plus clair.

Au lieu de déclarer toutes tes variables au début du script, tu peux le faire au fur et à mesure en ce qui concerne les lignes

Code :

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

1
2
3
4
#déclaration des variables à extraire des deux fichiers
my ($lineB,$contig,$oligo,$lineA,$query_name,$query_accession,$query_description,$query_length,$hit_name,$hit_accession,$hit_description,$hit_length,$score,$evalue,$hsp_number,$hsp_query_start,$hsp_query_end,$hsp_hit_start,$hsp_hit_end,$hsp_percent_identity,$direction);
while ($lineB=<OLIGO>) {
while ($lineA=<CONTIG>) {

Tu peux utiliser la même variable $line

Code :

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

1
2
while (my $line=<OLIGO>) {
while (my $line=<CONTIG>) {

inutile de créer la variable $noLine, cette information est contenue dans $.

Code :

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

1
2
if ($lineA =~/^#/ || $noLine==1)
if ($lineA =~/^#/ || $. == 1)

Le problème :

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
# récupération des données du fichier 1
while ($lineB=<OLIGO>) {
    if (defined $contig2oligo{$contig}) {
    #plusieurs oligos peuvent être associés au même contig
       $contig2oligo{$contig}=$contig2oligo{$contig}."\t$oligo";
    } else {
       $contig2oligo{$contig}=$oligo;
}
 
# récupération des données du fichier 2
while ($lineA=<CONTIG>) {
 
    #enregistrement des associations ID_oligo($query_name)/ID_contig ($hit_name) dans la table de hachage sous la clé ID_contig ($hit_name)
    if (defined $contig2oligo{$hit_name}) {
    #plusieurs ID_oligos peuvent être associés à un même oligo
       $contig2oligo{$hit_name}=$contig2oligo{$hit_name}."\t$query_name";
    } else {
       $contig2oligo{$hit_name}=$query_name;
}
 
 
 
if ($contig2oligo{$contig} eq $contig2oligo{$hit_name}){
    print OUT "$query_name\t$query_accession\t$query_description\t$query_length\t$hit_name\t$hit_accession\t$hit_description\t$hit_length\t$score\t$evalue\t$hsp_number\t$hsp_query_start\t$hsp_query_end\t$hsp_hit_start\t$hsp_hit_end\t$hsp_percent_identity\t$direction\tOK\n";
    } else {
       print OUT "$query_name\t$query_accession\t$query_description\t$query_length\t$hit_name\t$hit_accession\t$hit_description\t$hit_length\t$score\t$evalue\t$hsp_number\t$hsp_query_start\t$hsp_query_end\t$hsp_hit_start\t$hsp_hit_end\t$hsp_percent_identity\t$direction\tNA\n";
    }
}

$contig contient la dernière valeur trouvée dans le fichier et idem pour $hit_name. Et surtout idem pour toutes les valeurs de ce que tu veux écrire dans print out!

Un hash ne peut pas avoir 2 fois la même clé. Si lors de la lecture du premier fichier tu as trouvé une valeur pour $contig2oligo{$contig} ... si tu retrouve la même valeur de clé lors de la lecture du second fichier $contig2oligo{$hit_name} ... les valeurs seront ajoutées à $contig2oligo{$contig}

solution : créer 2 hashes et garder @elmtsA comme valeur.

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
my %oligo_tab;
# récupération des données du fichier 1
while ( my $line=<OLIGO>) {
    #enregistrement des associations ID_oligo($query_name)/ID_contig ($hit_name) dans la table de hachage sous la clé ID_contig ($hit_name)
    if (defined $oligo_tab{$hit_name}) {
    #plusieurs ID_oligos peuvent être associés à un même oligo
       $oligo_tab{$hit_name} .= "\t$query_name";
    } else {
       $oligo_tab{$hit_name} = $query_name;
}
 
 
my %contig_tab; 
# récupération des données du fichier 2
while ($line=<CONTIG>) {
 
    chomp $line;
    if ($line =~/^#/ || $.==1) {#Header
       print OUT "query_name\tquery_accession\tquery_description\tquery_length\thit_name\thit_accession\thit_description\thit_length\tscore\tevalue\thsp_number\thsp_query_start\thsp_query_end\thsp_hit_start\thsp_hit_end\thsp_percent_identity\tdirection\tmatch\n";
 
    next;
    }
 
    # clé = it_name valeur = liste des informations
    @elmtsA=split ";", $lineA;
    # $hit_name=$elmtsA[4];
    $contig2oligo{$elmtsA[4]} = @elmtsA;
}
 
# comparaison des clés
foreach my $key (keys %contig_tab){
 
	if  ( exists $contig_tab{$key}){
 
	}
}

**benlh** · 23/02/2009, 13h59

Certaines subtilités m'échappent encore mais je comprends mieux maintenant. Je vais essayer de mettre cela en application et merci beaucoup.

**Jasmine80** · 23/02/2009, 14h02

Code :

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

1
2
3
4
5
6
7
# comparaison des clés
foreach my $key (keys %contig_tab){
 
	if  ( exists $contig_tab{$key}){
 
	}
}

comparer tes valeurs et voir si elles sont égales va poser problème étant donné qu'elles sont concaténées via des tabulations dans les valeurs de tes hashes et que leur ordre n'est probablement pas le même. Il faudrait untiliser un hash de hah ou un hash de array. Je vais manger, puis je continue de regarder à ton script.

**Jasmine80** · 23/02/2009, 15h02

Code :

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

1
2
    $query_name=$elmtsA[0]; #=ID_oligo
    $hit_name=$elmtsA[4];#=ID_contig

Donc, cela vient d'un blast de l'oligo sur une DB contenant les différents contig. Pour un même $query_name, tu as donc plusieurs $hit_name ou bien n'as-tu garder que le meilleur hit dans ton fichier des résultats?

Pour un même oligo tu as donc plusieurs contig et pour un même contig, tu as également plusieurs oligo associés ... si tu as blasté tes 2 banques l'une contre l'autre.

**Jasmine80** · 23/02/2009, 15h57

Voici le script mis au propre. Vérifie qu'il fasse bien ce que tu veux.

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
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
 
use strict;
use warnings;
 
 
# module permettant de comparer des listes
# http://search.cpan.org/~jkeenan/List-Compare-0.37/lib/List/Compare.pm
use List::Compare;
 
 
 
#---------------------------------------------------------# 
# PARTIE 1 RECUPERATION DES DONNEES DU FICHIER OLIGO
# ET STOCKAGE DES DONNES DANS %oligo_tab;
#---------------------------------------------------------# 
 
#déclaration de la table de hachage pour le fichier OLIGO
my %oligo_tab;
 
#ouverture du fichier d'associations ID_oligo/ID_contig par le biais des ID_genbank
open OLIGO, "< oligo_otherV6.txt" or
    die "Unable to open input file #2 : $!\n";
 
#extraction des données ligne par ligne en sautant le header
while ($line = <OLIGO>) {
 
    if ($line =~ /^#/) {#Header
       next;
    }
 
    chomp $line; # le chomp mis après le next évite d'utiliser cette commande inutilement
 
    #splitage de chaque ligne dans le tableau elmtsB et mise en mémoire uniquement des ID_oligo et ID_contig dans les variables $coligo et $contig
    my @elmts = split "\t", $line;
 
    my $contig = $elmtsB[0];
    my $oligo = $elmtsB[4];
 
    #mémorisation dans la table de hachage des associassions ID_contig/ID_oligo sous les clés ID_contig
    push (@{$oligo_tab{$contig}}, $oligo);
} 
close(OLIGO); 
 
 
 
 
#---------------------------------------------------------# 
# PARTIE 2 RECUPERATION DES DONNEES DU FICHIER CONTIG
# ET STOCKAGE DES DONNES DANS %contig_tab;
#---------------------------------------------------------# 
 
#déclaration de la table de hachage pour le fichier CONTIG
my %contig_tab;
 
#ouverture du fichier de résultats des blasts
open CONTIG, "< blastn-oligo_restrict_otherV6-contig_restrict_otherV6-1e-3-1e-3-1,1.csv" or
    die "Unable to open input file #1 : $!\n";
#extraction des données ligne par ligne en sautant le header et en créant le header du fichier OUT
 
while ($line = <CONTIG>) {
 
	if ($line =~ /^#/ || $. == 1) {#Header
		next;
	}
 
	chomp $line;
 
	#splitage des lignes et extraction des données à mémoriser
	my @elmts = split ";", $line;
 
 
	$query_name=$elmtsA[0];
	$hit_name=$elmtsA[4];
 
	#enregistrement des associations ID_oligo($query_name)/ID_contig ($hit_name) dans la table de hachage sous la clé ID_contig ($hit_name)
	push (@{$contig_tab}{$hit_name}, $query_name);
 
	# on garde en mémoire les valeurs que l'on écrira dans le fichier de sortie OUT
	#	       ID_contig   $ID_oligo
	${$contig_tab}{$hit_name}{$query_name} = @elmts;
 
}
close(CONTIG); 
 
 
 
 
#---------------------------------------------------------# 
# PARTIE 3 COMPARAISON DES VALEURS RECUPEREES
# ET ECRITURE DANS LE FICHIER OUT
#---------------------------------------------------------# 
 
 
#création et ouverture du fichier pour consigner les résultats
open OUT, "> oligo_contig_perfect.txt" or
    die "Unable to open oligo_contig_perfect file : $!\n";
 
# entêtes des différentes colonnes
print OUT "query_name\tquery_accession\tquery_description\tquery_length\thit_name\thit_accession\thit_description\thit_length\tscore\tevalue\thsp_number\thsp_query_start\thsp_query_end\thsp_hit_start\thsp_hit_end\thsp_percent_identity\tdirection\tmatch\n";
 
 
# comparaison des clés
foreach my $contig (keys %oligo_tab){
 
	if  ( exists $contig_tab{$contig}){
 
		# on regarde si la liste des oligos qui se sont blastés à ce contig (fichier CONTIG)
		# est la même que celle qui ont été trouvées dans le fichier OLIGO.		
		my $lc = List::Compare->new(\@{$oligo_tab{$contig}}, \@{$contig_tab}{$hit_name});
 
		# @LorRonly est la liste des oligos n'apparaissant que dans un seul des 2 array
		# il doit donc être nul afin de conclure que les listes des oligo associés à un 
		# même contig dans les 2 fichiers d'entrées sont les mêmes
		my @LorRonly = $lc->get_symmetric_difference;
 
		my $switch = 'NA';
 
		if ( @LorRonly == 0 ){
			$switch = 'OK';
		}
 
		foreach my $query_name (keys %{${$contig_tab}{$hit_name}}){
			print OUT join ("\t", @{${$contig_tab}{$hit_name}{$query_name}})." $switch\n";
		}
	}
}
close(OUT);

PS : télécharge le module List::Compare

Le problème de vouloir que la liste des oligos associés aux contigs soit exactement la même dans les deux fichiers risque de poser problème si les paramètres de blast ont été mal réglés. Si Blast pour un contig trouve ne fusse qu'un hit de plus ou de moins que ce qui se trouve dans le fichier OLIGO ... cela sera donc considéré comme une différence d'oligos associés à ce contig.

**benlh** · 23/02/2009, 16h37

Merci beaucoup.
A vrai dire je me sens un peu dépassé en voyant ton script qui à l'air si simple par rapport à ce que j'avais commencé à faire. Par contre, là où ça coince c'est que je suis sur une machine windows. Pour mes "bricolages" perl et bioinformatique, je me connecte alors à un serveur qui est en distribution linux via PuTTY , du coup je ne peux pas "installer" le module List::Compare. Par contre je peux essayer d'installer sur mon poste une distribution windows de perl. Dans ce cas est-ce que je pourrais installer le module List::Compare et cela fonctionnera-t-il ?

**benlh** · 23/02/2009, 17h31

bon, j'ai installé ActivePerl ainsi que le module list::compare sur mon poste.
par contre lorsque je vérifie le script avant de l'exécuter (perl -cw "nom_du_script") j'ai plusieurs erreurs dont la plupart est réglé (déclaration de variables en début de script). Par contre là où je m'interroge, c'est qu'il réclame la déclaration de $contig_tab dans les lignes suivantes :

Code :

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

1
2
3
4
5
6
7
8
9
 
#ligne 75
push (@{$contig_tab}{$hit_name}, $query_name);
#ligne 79
${$contig_tab}{$hit_name}{$query_name} = @elmts;
my $lc = List::Compare->new(\@{$oligo_tab{$contig}}, \@{$contig_tab}{$hit_name});
#lignes 121 et 122
foreach my $query_name (keys %{${$contig_tab}{$hit_name}}){
			print OUT join ("\t", @{${$contig_tab}{$hit_name}{$query_name}})." $switch\n";

enfin, j'ai ce message d'erreur pour la ligne 75 :
Type of arg 1 to push must be array (not hash slice) at conca_blastoligo_targete
dcontigs_D.pl line 75, near "$query_name)"

J'avoue que là je suis pommé. Peux tu m'expliquer un peu tout ça ?
merci d'avance.

ps: en parallèle je tente de poursuivre le script d'origine

**Jasmine80** · 24/02/2009, 08h47

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
#!/usr/local/bin/perl
 
use strict;
use warnings;
 
 
# structure de %contig_tab
# hash de array
# clé : hit name     array correspondante : liste des query name associés à ce hit name
 
my %contig_tab = (
	hit_name1 => ["query_name1", "query_name2", "query_name3"],
        hit_name2 => ["query_name4", "query_name5", "query_name6"],
	hit_name3 => ["query_name7", "query_name8", "query_name9"],
);
 
 
=h
Au lieu de :
#ligne 75
push (@{$contig_tab}{$hit_name}, $query_name);
=cut
 
# ajout d'un query name à la liste
my $query_name = "query_name10";
push (@{$contig_tab{hit_name3}}, $query_name);
 
 
# affichage de l'array associé à hit_name3
map{print "$_\n"} @{$contig_tab{hit_name3}};
 
=h
Au lieu de :
#ligne 79
${$contig_tab}{$hit_name}{$query_name} = @elmts;
my $lc = List::Compare->new(\@{$oligo_tab{$contig}}, \@{$contig_tab}{$hit_name});
=cut
 
# tableau qui contiendra les informations du blast
# hash de hash
# clés : hit name  et query_name
my %association_tab;
 
# array avec les informations du blast
my @elmts = ('query_accession', '$query_description', 'query_length', 'hit_name', 'hit_accession', 'hit_description', 'hit_length');
 
# utilisation d'une référence de l'array signe \ placé devant @elmts
$association_tab{hit_name3}{query_name10} = \@elmts;
 
=h
Au lieu de :
#lignes 121 et 122
foreach my $query_name (keys %{${$contig_tab}{$hit_name}}){
			print OUT join ("\t", @{${$contig_tab}{$hit_name}{$query_name}})." $switch\n";
=cut
 
foreach my $hit_name (keys %contig_tab){ 
 
	# les $query_name sont contenus dans un array et non dans un hash
	foreach my $query_name (@{$contig_tab{$hit_name}}){
 
		print "HIT NAME $hit_name : ".join ("\t", @{$association_tab{hit_name3}{query_name10}})."\n";
	}
}

J'avoue que là je suis pommé. Peux tu m'expliquer un peu tout ça ?
merci d'avance.

Dans ce script, je rajoute $query_name dans l'array de la clé hit_name3 du tableau. Est-ce plus clair ainsi? Exécute ce script, cela t'aidera à comprendre documentation

Les erreurs sont normales, je m'étais trompées dans mes crochets. J'ai corrigé mes erreurs dans le script ci-joint. Essaie de comprendre comment cela fonctionne. L'ennui est que j'écris le code sans pouvoir l'exécuter. Si tu me donnes quelques de tes 2 fichiers, en remplaçant les données si elles sont confidentielles, je pourrai faire tourner mon script et plus facilement le corriger.

ps: en parallèle je tente de poursuivre le script d'origine

Dans ton script d'origine, j'ai l'impression que
1) tu lis le fichier OLIGO
2) à chaque ligne tu ouvres et lis le fichier CONTIG
... essaies de ne lire qu'une seul fois chaque fichier en gardant l'ensemble des informations dans des tableaux.
... ton script ayant été posté en 3 morceaux et tes accolades n'étant pas alignées, il est difficile de vraiment savoir si la lecture de ton second fichier est imbriquée dans celle du premier ou pas.

N'hésite pas si tu as d'autres questions.

**benlh** · 24/02/2009, 10h17

Merci, j'y vois un peu plus clair.
C'est juste que j'ai encore du mal à me représenter les structures des différents "objets" et les contextes d'utilisation. Je découvre également petit à petit des nouvelles fonctions (comme push) et des raccourcis d'écriture et surtout, ma logique n'est pas forcément très bonne...

Je vais essayer de te donner des "templates" pour tester le script. Les données ne sont pas confidentielles mais il y a une masse importante d'informations et je voudrait essayer de te fournir un échantillon représentatif des différents cas qui peuvent être rencontrés.

J'essaie de faire ça rapidement et vais re-publier mon script en un seul morceau sous une forme qui s'approche de ce que tu me fournis car c'est vrai que c'est nettement plus facile à déchiffrer.

**benlh** · 24/02/2009, 11h47

Voilà des templates simplifiés des deux fichiers. Les noms sont les mêmes avec en préfixe "Copie de ". Par contre j'ai été obligé de changer l'extension .csv du fichier blast en .txt.

Je continue de remanier mon script pour le rendre plus clair et le colle sur le forum dès que possible (désolé que cela prenne du temps mais je partage mon temps entre la paillasse et l'analyse bioinfo...).

**benlh** · 24/02/2009, 12h22

voilà les deux premières parties du script pour récupérer les données. Là où je bloque encore c'est sur la comparaison des données des deux hash.

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
 
use strict;
use warnings;
 
#---------------------------------------------------------# 
# DECLARATION DES VARIABLES
#---------------------------------------------------------# 
 
my ($lineB,$contig,$oligoA,$oligoB,$lineA,$query_name,$hit_name);
my (@elmtsA,@elmtsB);
my %contigA2oligoA; #hash liant ID_oligo/ID_contig avant blast
my %contigB2oligoB; #hash liant ID_oligo/ID_contig après blast
 
 
#création et ouverture du fichier résultats
open OUT, "> oligo_contig_perfect.txt" or
	die "Unable to open oligo_contig_perfect file : $!\n";
 
#---------------------------------------------------------# 
# PARTIE 1 RECUPERATION DES DONNEES DU FICHIER OLIGO
# ET STOCKAGE DES DONNES DANS %contigA2oligoA;
#---------------------------------------------------------# 
 
#ouverture du fichier d'associations ID_oligo/ID_contig par le biais des ID_genbank
open ASSO, "< oligo_otherV6.txt" or
	die "Unable to open input file #2 : $!\n";
 
	# saut du header
  while ($lineA=<ASSO>) {
	if ($lineA =~ /^#/) {#Header
    next;
  }
 
  #extraction des données ligne par ligne
  chomp $lineA;
 
  #splitage de chaque ligne dans le tableau elmtsA et mise en mémoire uniquement des ID_oligo et ID_contig dans les variables $oligoA et $contig
	@elmtsA=split "\t", $lineA;
	$contig=$elmtsA[0];
	$oligoA=$elmtsA[4];
 
	#mémorisation dans la table de hachage des associassions ID_contig/ID_oligo sous les clés ID_contig
  if (defined $contigA2oligoA{$contig}) {
		$contigA2oligoA{$contig}=$contigA2oligoA{$contig}."\t$oligoA";
	} else {
		$contigA2oligoA{$contig}=$oligoA;
	}
}
close (ASSO);
 
 
#---------------------------------------------------------# 
# PARTIE 2 RECUPERATION DES DONNEES DU FICHIER CONTIG
# ET STOCKAGE DES DONNES DANS %contig_tab;
#---------------------------------------------------------#
 
#ouverture du fichier de résultats des blasts                                                
open ASSO, "< blastn-oligo_restrict_otherV6-contig_restrict_otherV6-1e-3-1e-3-1,1.csv" or
	die "Unable to open input file #1 : $!\n";
 
# saut du header ou d'une ligne vide
## my $noLine=0; ## inutile de créer la variable $noLine, cette information est contenue dans $. 
while ($lineB=<ASSO>) {
	##$noLine++;
	chomp $lineB;
	##if ($lineB =~/^#/ || $noLine==1) {#Header
	if ($lineB =~/^#/ || $. == 1)
      print OUT "$lineB;match\n"; #création du header pour le fichier résulats
      next;
  }
 
  #splitage des lignes et extraction des données à mémoriser
	@elmtsB=split ";", $lineB;
	$query_name=$elmtsB[0]; # =ID_oligo
	$hit_name=$elmtsB[4]; # =ID_contig
 
  #enregistrement des associations ID_oligo($query_name)/ID_contig ($hit_name) dans le hash sous la clé ID_contig ($hit_name)
	if (defined $contigB2oligoB{$hit_name}) {
		$contigB2oligoB{$hit_name}=$contigB2oligoB{$hit_name}."\t$query_name";
	} else {
		$contigB2oligoB{$hit_name}=$query_name;
  }
}
close (ASSO);
 
 
#---------------------------------------------------------# 
# PARTIE 3 COMPARAISON DES VALEURS RECUPEREES
# ET ECRITURE DANS LE FICHIER OUT
#---------------------------------------------------------#

**Zwiter** · 24/02/2009, 15h36

voici une solution, ou du moins une idée de solution :
1- fait un hash de hash avec ton contig et ton oligo a partir du 1er fichier, et aau fure et mesure que tu lis ton 2eme fichier,
2- tu verifies la presence ou non de ces paires. Au passage, tu marque quand c'est trouver.
3- Puis tu parcours une derniere fois a la recherche des non marquées
Autrement dit (attention, il peut y avoir des fautes de codes, j'ai pas executé) :

Code :

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

1
2
#mémorisation dans la table de hachage des associassions ID_contig/ID_oligo sous les clés ID_contig
$contigA2oligoA{$contig}{$oligoA} = 0;

Code :

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

1
2
3
4
5
6
7
8
9
#enregistrement des associations ID_oligo($query_name)/ID_contig ($hit_name) dans le hash sous la clé ID_contig ($hit_name)
if (exists $contigA2oligoA{$hit_name}{$query_name}) {
	print "$query_name/$hit_name : présent!!\n";
	$contigA2oligoA{$hit_name}{$query_name} = 1; #switch a 1 : il est présent
	} 
else {
	print "$query_name/$hit_name : non présent!!\n";
	#action dans ce cas genre print OUT si c'est les non présent que tu veux recup
	}

Code :

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

1
2
3
4
5
6
7
8
for my $i (keys %contigA2oligoA) {
	for my $j (keys %{$contigA2oligoA{$i}}) {
		if (!$contigA2oligoA{$i}{$j}) { #switch toujours a 0, donc pas retrouvé
			print "$j/$i : non présent!!\n";
			#action dans ce cas genre print OUT si c'est les non présent que tu veux recup
			}
		}
	}

Z.

**Jasmine80** · 24/02/2009, 17h05

Voici la dernière version du script.

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
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
 
 
 
use strict;
use warnings;
 
 
# module permettant de comparer des listes
# http://search.cpan.org/~jkeenan/List-Compare-0.37/lib/List/Compare.pm
use List::Compare;
 
 
 
#---------------------------------------------------------# 
# PARTIE 1 RECUPERATION DES DONNEES DU FICHIER OLIGO
# ET STOCKAGE DES DONNES DANS %oligo_tab;
#---------------------------------------------------------# 
 
#déclaration de la table de hachage pour le fichier OLIGO
my %oligo_tab;
 
#ouverture du fichier d'associations ID_oligo/ID_contig par le biais des ID_genbank
open OLIGO, "< oligo_otherV6.txt" or
    die "Unable to open input file #2 : $!\n";
 
#extraction des données ligne par ligne en sautant le header
while ( my $line = <OLIGO>) {
 
    if ($line =~ /^#/) {#Header
       next;
    }
 
    chomp $line; 
 
    #splitage de chaque ligne dans le tableau elmtsB et mise en mémoire uniquement des ID_oligo et ID_contig dans les variables $coligo et $contig
    my @elmts = split "\t", $line;
 
    my $contig = $elmts[0];
    my $oligo = $elmts[4];
 
    #mémorisation dans la table de hachage des associassions ID_contig/ID_oligo sous les clés ID_contig
    push (@{$oligo_tab{$contig}}, $oligo);
} 
close(OLIGO); 
 
 
 
 
#---------------------------------------------------------# 
# PARTIE 2 RECUPERATION DES DONNEES DU FICHIER CONTIG
# ET STOCKAGE DES DONNES DANS %contig_tab;
#---------------------------------------------------------# 
 
#déclaration de la table de hachage pour le fichier CONTIG
my %contig_tab;
 
# tableau qui contiendra les informations du blast
# hash de hash
# clés : hit name  et query_name
my %association_tab;
 
 
#ouverture du fichier de résultats des blasts
open CONTIG, "< blastn-oligo_restrict_otherV6-contig_restrict_otherV6-1e-3-1e-3-1,1.txt" or
    die "Unable to open input file #1 : $!\n";
#extraction des données ligne par ligne en sautant le header et en créant le header du fichier OUT
 
while ( my $line = <CONTIG>) {
 
	if ($line =~ /^#/ || $. == 1) {#Header
		next;
	}
 
	chomp $line;
 
	#splitage des lignes et extraction des données à mémoriser
	my @elmts = split ";", $line;
 
 
	my $query_name=$elmts[0];
	my $hit_name=$elmts[4];
 
	#enregistrement des associations ID_oligo($query_name)/ID_contig ($hit_name) dans la table de hachage sous la clé ID_contig ($hit_name)
	push (@{$contig_tab{$hit_name}}, $query_name);
 
	# on garde en mémoire les valeurs que l'on écrira dans le fichier de sortie OUT
	#	       ID_contig   $ID_oligo
	# utilisation d'une référence de l'array signe \ placé devant @elmts
	$association_tab{$hit_name}{$query_name} = \@elmts;
 
}
close(CONTIG); 
 
 
 
 
#---------------------------------------------------------# 
# PARTIE 3 COMPARAISON DES VALEURS RECUPEREES
# ET ECRITURE DANS LE FICHIER OUT
#---------------------------------------------------------# 
 
 
#création et ouverture du fichier pour consigner les résultats
open OUT, "> oligo_contig_perfect.txt" or
    die "Unable to open oligo_contig_perfect file : $!\n";
 
# entêtes des différentes colonnes
print OUT "query_name\tquery_accession\tquery_description\tquery_length\thit_name\thit_accession\thit_description\thit_length\tscore\tevalue\thsp_number\thsp_query_start\thsp_query_end\thsp_hit_start\thsp_hit_end\thsp_percent_identity\tdirection\tmatch\n";
 
 
# comparaison des clés
foreach my $contig (keys %oligo_tab){
 
	if  ( exists $contig_tab{$contig}){
 
		# on regarde si la liste des oligos qui se sont blastés à ce contig (fichier CONTIG)
		# est la même que celle qui ont été trouvées dans le fichier OLIGO.		
		my $lc = List::Compare->new(\@{$oligo_tab{$contig}}, \@{$contig_tab{$contig}});
 
		# @LorRonly est la liste des oligos n'apparaissant que dans un seul des 2 array
		# il doit donc être nul afin de conclure que les listes des oligo associés à un 
		# même contig dans les 2 fichiers d'entrées sont les mêmes
		my @LorRonly = $lc->get_symmetric_difference;
 
		my $switch = 'NA';
 
		if ( @LorRonly == 0 ){
			$switch = 'OK';
		}
 
		foreach my $query_name (@{$contig_tab{$contig}}){
			print OUT join ("\t", @{$association_tab{$contig}{$query_name}})." $switch\n";
		}
	}
}
close(OUT);

Et voici le contenu du fichier de sortie (j'ai gardé les données de @elmts dans l'ordre initial)

query_name query_accession query_description query_length hit_name hit_accession hit_description hit_length score evalue hsp_number hsp_query_start hsp_query_end hsp_hit_start hsp_hit_end hsp_percent_identity direction match
A_70_P000311 A_70_P000311 60 EE752570.p.ov.6 EE752570.p.ov.6 50 84 6.00E-18 1 1 49 2325 2374 98 -1 NA
A_70_P000311_bis A_70_P000311_bis 60 EE752570.p.ov.6 EE752570.p.ov.6 50 84 6.00E-18 1 1 49 2325 2374 98 -1 NA
A_70_P000001 A_70_P000001 60 EE747509.p.ov.6 EE747509.p.ov.6 60 119 1.00E-28 1 1 60 1435 1494 100 1 OK
A_70_P000331 A_70_P000331 60 EE775091.p.ov.6 EE775091.p.ov.6 60 119 1.00E-28 1 1 60 1305 1364 100 1 OK
A_70_P000035 A_70_P000035 60 EE751947.p.ov.6 EE751947.p.ov.6 60 119 1.00E-28 1 1 60 447 506 100 1 OK
A_70_P000181 A_70_P000181 60 EE752292.p.ov.6 EE752292.p.ov.6 61 98 4.00E-22 1 1 60 1343 1403 96.72131148 1 OK
A_70_P000121 A_70_P000121 60 CD286222.p.ov.6 CD286222.p.ov.6 60 119 1.00E-28 1 1 60 2002 2061 100 1 OK
A_70_P000041 A_70_P000041 60 EE751953.p.ov.6 EE751953.p.ov.6 60 119 1.00E-28 1 1 60 181 240 100 -1 OK
A_70_P000116 A_70_P000116 60 EE752170.p.ov.6 EE752170.p.ov.6 60 119 1.00E-28 1 1 60 1650 1709 100 1 OK
A_70_P000216 A_70_P000216 60 EE752407.p.ov.6 EE752407.p.ov.6 56 80 9.00E-17 1 1 54 747 802 96.42857143 1 OK

**benlh** · 24/02/2009, 18h25

Je n'arrivais pas à faire fonctionner ton script avec les dernières modification mais maintenant je comprends mieux, j'avais gardé les $hit_name et $query_name dans la troisième partie au lieu de prendre $contig et $oligo. Du coup j'avais une erreur du type "$query_name non définie dans le hash".
J'ai en plus modifié le print final pour avoir les OK et les NA dans la colonne "match" : à la place de la concaténation j'ai mis une tabulation.

Code :

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

1
2
 
print OUT join ("\t", @{$association_tab{$contig}{$query_name}})."\t$switch\n";

Sinon, j'ai essayé différents cas de figures ce script fonctionne relativement bien.
juste un défaut : si une association ID_contig/ID_oligo existe dans le fichier résultat du blast avec un ID_contig non présente dans le fichier d'origine des associations ID_contig/ID_oligo, celle-ci n'apparait plus dans plus le fichier OUT. Je pense qu'il faudrait ajouter une catégorie KO qui rende compte de ce cas de figure même si a priori il ne devrait pas exister à cause de la construction des banques de séquences restreintes.

Merci beaucoup Jasmine. Je vais essayer de faire cette troisième catégorie demain.

Merci également à Zwiter pour ses suggestions, je suis en train de regarder ça aussi.

Merci à tous.

problème de comparaison d'association de deux données entre deux fichiers

Bioinformatique Perl

Discussions similaires

Partager

Partager