Calcul distance genomique/proteique [Résolu]

[shadow19c]

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 :

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 :

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 :

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]

... ai-je bien compris

Code :

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 :

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]

voici le tableau :

Code :

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]

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]

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
#!/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]

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]

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]

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

[shadow19c]

voila je sais que cest tout simple je voudrais

aller dans le premier fichier

Code :

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 :

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]

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

Code :

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]

Envoyé par Jasmine80

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

Code :

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

Code :

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]

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]

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

Code :

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]

et ça va maintenant ton programme?

[shadow19c]

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]

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

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/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]

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

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
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 :

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

[djibril]

Ayant la flemme de tout relire, il faudrait que tu nous explique concrètement ce que tu souhaites faire avec des exemples simples, clairs et concis.

[shadow19c]

Merci, mais j'ai trouve la solution, j'ai refait mon programme et ca marche merci en tout cas

[djibril]

Dans ce cas, merci de taguer le fil en résolu.