Calcul distance genomique/proteique

**shadow19c** · 07/03/2011, 10h01

Bonjour,
Je suis nouveau en perl car je fais plus de python, mais la je voulais faire un programme de calcul nucleotidique en perl car on m'a dit que ce serait plus simple:
Voila mon soucis:
J'ai un fichier blast:

Code :

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

1
2
3
4
5
 
C169v2-00001	Aster-03513	64	112	1192	1246	
C169v2-00002	Aster-03154	2	83	172	264	
C169v2-00003	Aster-06252	74	116	57	96	
C169v2-00004	Aster-05984	1	218	24	243

Les deux premiers colonne representent mes deux proteines et les 4 autres les positions, exemple C169v2-00001 64 112 et pour Aster-03513 1192 1246
je voudrais a partir des coordonnees genomiques qui sont dans un fichier cds faire la cocordance avec les coordonnees ci dessus.

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
 
# Colonne 1: scaffold
#         2: gene id
#         3: start
#         4: dir
#         5: origine gene model
#         6: proteinId
#         7: ex-nom
#         8: coordonnees genomiques
#         9: structure intron/exon
#        10: segments du CDS couvert par EST
C169-scaffold_1	C169v2-00001	3777	-	ORIGINAL JGI	55057	Genemark1.1_g	3777..3857,4046..4192,4443..4561,4940..5234,5406..5540,5734..5847,6009..6098,6421..6492
C169-scaffold_1	C169v2-00002	12682	-	ORIGINAL JGI	55058	Genemark1.2_g	12682..12691,13195..13445,13694..13711
C169-scaffold_1	C169v2-00003	18095	+	ORIGINAL JGI	31905	fgenesh1_kg.1_#_1_#_4092_1_CBOZ_CBPA	18095..18097,18280..18410,18690..18972
C169-scaffold_1	C169v2-00004	20452	+	ORIGINAL JGI	6968	gw1.1.615.1	20452..20496,20636..20726,20881..21046,21194..21382,21567..21735

par exemple restons sur les proteines CV

Voila l'algorithme que j'ai pense:
Si je prends le premier C169v2-00001 3777..3857,4046..4192,4443..4561,4940..5234,5406..5540,5734..5847,6009..6098,6421..6492
Je dois remplir une table C169v2-00001table=(37777,3778,....22362)
ceci est la meme chose que si je dis d'apres ce que j'ai lu en perl @table=(3777..3857,3777..3857,4046..4192,4443..4561,4940..5234,5406..5540,5734..5847,6009..6098,6421..6492)

Ensuite je fais un

Code :

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

1
2
3
 
for(i=0;i<@tailletable;i+3)
         aa++;

Bon je resume car la je suis pas bon en perl:
si j'ai coco1 10..39,100..129
je dois remplir une table aa--> codon:
1->10
2->13
.
.
.
arrive a 11->39 il y a la virgule(qui correspond a un intron)
12->112(100+12)

a la fin je dois partir au premier fichier et faire la concordance coco1 10->1,129->19 par exple

Est ce que quelqu' un a une idee d'algorithme car la je suis perdu en perl
Le souci c'est que j'ai fait le programme en python mais j'ai un calcul du dico qui marche pas , donc on m'a dit qu'en perl ca serait plus simple.
Merci

**Jasmine80** · 07/03/2011, 12h52

... ai-je bien compris

Code :

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

C169v2-00001	Aster-03513	64	112	1192	1246

blast => donc tu sais que C169v2-00001 du 64ième AA au 112ième s'apparie avec Aster-03513 du 1192ième au 1246ième?

Code :

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

C169-scaffold_1	C169v2-00001	3777	-	ORIGINAL JGI	55057	Genemark1.1_g	3777..3857,4046..4192,4443..4561,4940..5234,5406..5540,5734..5847,6009..6098,6421..6492

fichier cds => tu as la séquence en nucléotides et tu sais que ta protéine correspond à la traduction des nucléotides : 3777..3857,4046..4192,4443..4561,4940..5234,5406..5540,5734..5847,6009..6098,6421..6492

tu veux donc récupérer la séquence cDNA correspondant à ta protéine de la position du 64ième AA au112ième AA?

est-ce bien cela?

**Jasmine80** · 07/03/2011, 13h29

voici le tableau :

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
#!/usr/local/bin/perl
 
 
use strict;
use warnings;
 
 
my @nucleotides_tab = (3777..3857,4046..4192,4443..4561,4940..5234,5406..5540,5734..5847,6009..6098,6421..6492);
 
# AA tous les 3 nucléotides
# AA 1 : 3777
# AA 2 : 3780
# il faut donc prendre une valeur sur 3 
 
 
my $pos  = 1;
for (my  $i = 0; $i <= $#nucleotides_tab; $i+=3){
 
	# affichage des positions des AA
	print 'AA'.$pos.' => '.$nucleotides_tab[$i]."\n";
	$pos++;
}

**shadow19c** · 07/03/2011, 13h52

Bonjour,
oui en fait cest ca, ensuite je dois dire que
C169v2-00001 64->4154 112->4548

car enn fait je dois pas afficher print 'AA'.$pos.' => '.$nucleotides_tab[$i]."\n";

merci je vais d'abord essayer de comprendre comment on fait pour recuperer les valeurs dans un fichier perl pour que j'automatise les systeme merci pour le debut

**Jasmine80** · 07/03/2011, 14h01

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
#!/usr/local/bin/perl
 
 
use strict;
use warnings;
 
 
my @nucleotides_tab = (3777..3857,4046..4192,4443..4561,4940..5234,5406..5540,5734..5847,6009..6098,6421..6492);
 
# AA tous les 3 nucléotides
# AA 1 : 3777
# AA 2 : 3780
# il faut donc prendre une valeur sur 3 stokée dans AA_tab
 
my @AA_tab;
 
for (my  $i = 0; $i <= $#nucleotides_tab; $i+=3){
 
	push (@AA_tab, $nucleotides_tab[$i]);
}
 
 
# comme on commence avec l'indice 0, décalage d'un indice
# AA 64 : indice 63
 
 
my $start = 64;
my $end = 112;
 
print $start." => ".$AA_tab[$start-1]."\n";
print $end." => ".$AA_tab[$end-1]."\n"

**Jasmine80** · 07/03/2011, 14h05

il y a certainement de faire un plus beau code en perl ... mais je ne sais pas comment, il faut attendre l'avis d'autres personnes. Tant que ça marche, c'est le principal mais si on peut l'améliorer tant mieux.

**shadow19c** · 07/03/2011, 14h16

Non c'est tres bien c'est comme ca le plus dur pour moi cest de comprendre comment je vais faire ca pour chaque ligne(cest a dire pour chaque proteine)
car comme je suis nouveau en perl, il faut que je regarde comment on fait pour ouvrir les fichiers, spliter, et garder en memoire.

J'ai fait plus de python, et la ce code est pour beaucoup plus simple.

merci

**Jasmine80** · 07/03/2011, 14h20

bonne continuation et si tu as d'autres questions, n'hésite pas

**shadow19c** · 07/03/2011, 14h35

voila je sais que cest tout simple je voudrais

aller dans le premier fichier

Code :

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

1
2
3
4
5
6
 
 
C169v2-00001	Aster-03513	64	112	1192	1246	
C169v2-00002	Aster-03154	2	83	172	264	
C169v2-00003	Aster-06252	74	116	57	96	
C169v2-00004	Aster-05984	1	218	24	243

aller recuperer C169v2-00001 et 64 112

j'ai fait

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
 
#!/usr/local/bin/perl
 
 
use strict;
use warnings;
 
my @tab;
open (FIC, "test.txt");
my @fic = <FIC> ;
foreach my $line (@fic)
{
       @tab=$line
}
print @tab[0]:

le probleme cest quil maffiche tote la ligne.

**Jasmine80** · 07/03/2011, 14h49

tu dois analyser ligne par ligne avec une expression régulière

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
#!/usr/local/bin/perl
 
 
use strict;
use warnings;
 
 
open my $fic, '<', "test.txt";
 
while (my $line = <$fic>){
 
	if ($line =~ m/([\w-]+)\s+([\w-]+)\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)/){
		my $id1 = $1;
		my $id2 = $2;
		my $start1 = $3;
		my $end1 = $4;
		my $start2 = $5;
		my $end2 = $6;
 
	}
 
}

**Norore** · 07/03/2011, 16h15

Envoyé par Jasmine80

tu dois analyser ligne par ligne avec une expression régulière

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
#!/usr/local/bin/perl
 
 
use strict;
use warnings;
 
 
open my $fic, '<', "test.txt";
 
while (my $line = <$fic>){
 
	if ($line =~ m/([\w-]+)\s+([\w-]+)\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)/){
		my $id1 = $1;
		my $id2 = $2;
		my $start1 = $3;
		my $end1 = $4;
		my $start2 = $5;
		my $end2 = $6;
 
	}
 
}

Le code est exact mais il y a encore plus direct :

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
#!/usr/bin/perl
use strict;
use warnings;
 
# ouvre le fichier test.txt si existant sinon arrête le programme en retournant le type d'erreur
open (FIC, '<', "test.txt") or ( warn "Can't open file test.txt\n$!\n" and die);
 
while (<FIC>){
        # retire le \n à la fin de chaque ligne
	chomp;
        # découpe la chaîne par les \t et récupére chaque valeur dans une variable
        my ($id1, $id2, $start1, $end1, $start2, $end2) = split(/\t/, $_);
        # affichage ligne par ligne des infos souhaitées
        print $id1."\t".$start2."\t".$end1."\n";
}
# fermeture du fichier
close FIC;

Jasmine80, pourquoi faire une expression régulière lorsqu'un simple split() suffit pour ce genre de traitement ? Y-a-t-il une raison particulière ?

**Jasmine80** · 08/03/2011, 09h08

Envoyé par Norore

Jasmine80, pourquoi faire une expression régulière lorsqu'un simple split() suffit pour ce genre de traitement ? Y-a-t-il une raison particulière ?

non, tu as raison, je n'y avais pas pensé mais il faut être certain qu'il n'a pas d'autres lignes que celles citées plus haut (pas d'entête au fichier).

**shadow19c** · 08/03/2011, 09h13

Non il n'y a rien pour ce fichier la mais pour le deuxieme il ya l'en tete a enleve:

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
 
# Colonne 1: scaffold
#         2: gene id
#         3: start
#         4: dir
#         5: origine gene model
#         6: proteinId
#         7: ex-nom
#         8: coordonnees genomiques
#         9: structure intron/exon
#        10: segments du CDS couvert par EST
C169-scaffold_1	C169v2-00001	3777	-	ORIGINAL JGI	55057	Genemark1.1_g	3777..3857,4046..4192,4443..4561,4940..5234,5406..5540,5734..5847,6009..6098,6421..6492
C169-scaffold_1	C169v2-00002	12682	-	ORIGINAL JGI	55058	Genemark1.2_g	12682..12691,13195..13445,13694..13711
C169-scaffold_1	C169v2-00003	18095	+	ORIGINAL JGI	31905	fgenesh1_kg.1_#_1_#_4092_1_CBOZ_CBPA	18095..18097,18280..18410,18690..18972
C169-scaffold_1	C169v2-00004	20452	+	ORIGINAL JGI	6968	gw1.1.615.1	20452..20496,20636..20726,20881..21046,21194..21382,21567..21735

**Jasmine80** · 08/03/2011, 09h14

et ça va maintenant ton programme?

**shadow19c** · 08/03/2011, 09h17

Ben j'ai pas continuer depuis hier car je dois modifier un premier programme python que je galere un peu ensuite je vais attaquer celui la car j'ai fait que ouvrir le premier fichier et recuperer l'ID et les positions.

Merci

PS: Je sais que je pose pas au bon endroit ma question , mais sur un forum on m'ai aide pour un souci de recuperation de donnees mais par SQL.
Est ce qu'il faut forcement que j'ai sgbd pour faire ca?

Par ce que j'ai fait un programme en python qui me fait un reciprocal blast best hits, comme pour le premier fichier, sauf que j'ai oublie de recuperer les positions pour chaque couple, et la je suis bloque car je ne sais ou mettre mon split sans modifier mes tests de comparaison.

**shadow19c** · 08/03/2011, 13h23

Voila mon code mais j'ai un souci dans le calcul car ca ne m'affiche pas la correspondance.

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
 
#!/usr/local/bin/perl
#use strict;
#use warnings;
 
open (DESC1, "fic1.txt") ;
open (DESC2, "fic2");
while ($ligne = <DESC1>){
  $pos=1;
  @table = split(/	/g,$ligne);
 
	while ($ligne1 = <DESC2>){
		 @table1 = split(/	/g,$ligne1);
		 if($table1[0] eq "#") {next ;}
		 elsif($table1[1] eq $table[0]){
 
 
		 		if($table1[3] eq "+"){
 
		 				@tableau=$table1[7];
		 			}
 
		 			for($i=0; $i<@tableau; $i+=3){
		 				push (@AA_tab, $tableau[$i]);
		 				print " ".$tableau[$i];
						$pos++;
		 			}
		 		}
 
 
		 		if($table1[3] eq "-"){
		 			for($i=@table1-1; $i>7; $i--){
		 				$tableau[$i-7]=$table1[$i];
		 			}
 
		 			for($i=0; $i<@tableau; $i+=3){
		 				push (@AA_tab, $tableau[$i]);
						$pos++;
		 			}
		 		}
 
 
 
	}
	print $table[0]." => ".$table[2]." => ".$AA_tab[$table[2]-1]."\n";
	print $table[0]." => ".$table[3]." => ".$AA_tab[$table[3]-1]."\n";
 
 
}

PS: Je fais le test du - et + car cest pour le sens de lecture.

Par contre je voulais vous poser une question , car le but de cette methode est l'etude de la variation genomique entre deux especes.
Donc a partir de ce parseur,je part a partir du blast proteique en connaissant pour la proteine Cv et Aster leur positionnement proteique.
Et donc je vais dans le fichier cds(genomique) faire le positionnement de notre proteine , a partir de la faire la difference entre la fin et le debut de chaue gene.EN faisant attention aux introns car c'est les introns qui font la variation des genes.
Donc est ce que le code du debut en mettant les my @nucleotides_tab = (3777..3857,4046..4192,4443..4561,4940..5234,5406..5540,5734..5847,6009..6098,6421..6492); comme ceci prend t il en compte les introns symbolyse par les , .

Merci

**shadow19c** · 09/03/2011, 13h54

Voici un nouveau code mais jai un souci danns le renversement et ca me prend les bonnes valeurs

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
 
 
#!/usr/local/bin/perl
#use strict;
#use warnings;
 
open (DESC1, "test.txt") ;
 
open (DESC2, "coco_cds");
@toto=<DESC2>;
while ($ligne = <DESC1>){
  $pos=1;
  @table = split(/	/g,$ligne);
	#print(" deb $table[2]  fin $table[3]\n");
	foreach my $ligne1 (@toto){
	#while ($ligne1 = <DESC2>){
 
 
 
		 @table1 = split(/	/g,$ligne1);
	         #print("tab de 0: $table1[0]\n");
 
		if(($table1[0] =~ m/#/) ){
			#print("bonjour\n"); 
			next ;
		}
 
		elsif($table1[1] eq $table[0]){
 
		 		if($table1[3] eq "+"){
						@tableau=eval($table1[7]);
						#print("tab de 3 ++ $table1[3]=========> eval:  $tableau[0];\n");
 
 
					for($i=0; $i<= $#tableau; $i+=3){
						push (@AA_tab, $tableau[$i]);
						#print " ".$tableau[$i];
						$pos++;
					}
				}
 
				elsif($table1[3] eq "-"){
 
					@tab=eval($table1[7]);
					#print("tab de 3 -- $table1[3]=========> eval:  $tab[0];\n");
					#print("tab de 3 -- $table1[3]=========> eval:  $tab[-1];\n");
					for($i=@tab; $i>=0; $i--){
						#print("inverse\n");
						$tableau[$i]=$tab[$i];
					}
 
					for($i=0; $i< $#tableau; $i+=3){
						push (@AA_tab, $tableau[$i]);
						print 'AA'.$pos.' => '.$tableau[$i]."\n";
						$pos++;
					}
				}
 
		 }next;
 
 
 
	}
 
 
	print $table[0]." => ".$table[2]." => ".$AA_tab[$table[2]-1]."\n";
	print $table[0]." => ".$table[3]." => ".$AA_tab[$table[3]-1]."\n";
}

]

Par ce que voici ce quil me fait, il ne caclu que pour C169v2-00001, et je pense que l'inverse n'est pas bon, et il met les valeurs de calcul aux autres CV-****

Code :

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

1
2
3
4
5
6
7
8
9
10
11
 
 
 deb 2  fin 83
C169v2-00002 => 2 => 3780
C169v2-00002 => 83 => 4461
 deb 74  fin 116
C169v2-00003 => 74 => 4184
C169v2-00003 => 116 => 4560
 deb 1  fin 218
C169v2-00004 => 1 => 3777
C169v2-00004 => 218 => 5415

Merci