Calcul taille intergenique

[shadow19c]

Bonjour,
je dois faire un programme qui calcul la taille intergeniques entre deux orthologues.
Pour cela je dois tout d'abord partir d'un fichier de reciprocal best blast hits qui contient le couple de mes orthologues: (ceci est un petit bout du fichier)

Code :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
 
C169v2-04488	Aster-06787
C169v2-04487	Aster-06725
C169v2-04480	Aster-06724
C169v2-04481	Aster-06792
C169v2-03069	Aster-04749
C169v2-03068	Aster-01309
C169v2-03061	Aster-06908
C169v2-03060	Aster-04730
C169v2-03067	Aster-04746
C169v2-03066	Aster-04710
C169v2-03065	Aster-04711
C169v2-03064	Aster-05029
C169v2-03285	Aster-01180

Donc le but est de fixer les C169-***, car je vais dois regarder si deux genes qui se suivent par exemple C169v2-04480 C169v2-04481 n'ont pas de gens orthologues qui se suivent aussi même si ils sont renverses, par exemple si j'avais deux gènes Aster Aster -0002 puis 0001 ca marche .
Donc je voulais savoir je voulais faire un hash qui contiendrait tous les CV et leurs asters puis je pars du premier Cv-** puis je compare si j'ai pas un CV-**** plus grand que un ou plus petit puis je regarde leurs gènes orthologues.

P.S:Je refais l'exemple si je prends le premier C169v2-04488 je lui ajoute +1 je trouve C169v2-04489 je regarde si ce chiffre existe si il existe je vais regarder par la suite si le Aster-06787 correspond soit Aster-06786 ou Aster-06788 si ca ne correspond pas je regarde -1 je fais la meme chose.
Si j'ai fini je passe a l'autre.
merci

[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
#!/usr/bin/perl
 
use strict;
use warnings;
 
 
my %hash = (
'C169v2-03069'	=>	'Aster-04749',
'C169v2-03068'	=>	'Aster-01309',
'C169v2-03061'	=>	'Aster-06908',
'C169v2-03060'	=>	'Aster-04730',
'C169v2-03067'	=>	'Aster-04746',
'C169v2-03066'	=>	'C169v2-03065',
'Aster-04711'	=>	'Aster-04712',
 
);
 
while  ( my ($k, $v) = each %hash){
 
	my ($gene_k, $num_k) = $k =~ m/^(\w+)-(\d+)$/;
	my ($gene_v, $num_v) = $v =~ m/^(\w+)-(\d+)$/;
 
	if ((abs ($num_k - $num_v) == 1) && ($gene_k eq $gene_v) ){
 
		print "$k\t$v\n";
 
	}
}

[shadow19c]

Merci pour ta reponse toujours aussi parfaite, je suis entrain de chercher la complexite des sens des genes c'est pas drole merci en tout cas.
Par ce que le but c'est que je dois verfier dans le fichier cds le sens des genes + ou - je dois verfier la compatiblite des SA*SB

[Jasmine80]

je suis contente d'aider ... pour une fois que j'ai bien compris le problème

[shadow19c]

Par contre je voulais savoir , le code me fait il toutes les possibilites.
Cest a dire en ajoutant +1 verification et le -1.
Je pense que ta pas compris y a til un moyen de tenvoyer le fichier reciprocal best blast hit
merci

[Jasmine80]

oui, grâce à la valeur absolue (abs)

n'oublie pas d'indiquer quand ce sera résolu

[shadow19c]

Ok je garde , car ca se peut que j'aurais des questions a poser concernant la suite de mon programme.

PS: rien ne s'affiche si je fais:

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
 
 
#!/usr/local/bin/perl
open IN, "intron2.txt";
while(<IN>){
	my	($A,$B)=(split/\t/);	
	 my %hash;
	  $hash{$A}=$B;	
 
}
close IN;
 
 
 
 
 
while  ( my ($k, $v) = each %hash){
 
	my ($gene_k, $num_k) = $k =~ m/^(\w+)-(\d+)$/;
	my ($gene_v, $num_v) = $v =~ m/^(\w+)-(\d+)$/;
 
	if ((abs ($num_k - $num_v) == 1) && ($gene_k eq $gene_v) ){
 
		print "$k\t$v\n";
 
	}
}

[Jasmine80]

avec ma hash, ça fonctionne

Code :

1
2
3
4
5
6
7
8
9
10
my %hash = (
'C169v2-03069'	=>	'Aster-04749',
'C169v2-03068'	=>	'Aster-01309',
'C169v2-03061'	=>	'Aster-06908',
'C169v2-03060'	=>	'Aster-04730',
'C169v2-03067'	=>	'Aster-04746',
'C169v2-03066'	=>	'C169v2-03065',
'Aster-04711'	=>	'Aster-04712',
 
);

essaie de voir ce qui ne correspond pas

[shadow19c]

Oui, c'est normal que ca marche car vous avez fait:

'C169v2-03066' => 'C169v2-03065',
'Aster-04711' => 'Aster-04712',

Or le but c'est partir d'un fichier

Code :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
 
C169v2-04488	Aster-06787
C169v2-04487	Aster-06725
C169v2-04480	Aster-06724
C169v2-04481	Aster-06792
C169v2-03069	Aster-04749
C169v2-03068	Aster-01309
C169v2-03061	Aster-06908
C169v2-03060	Aster-04730
C169v2-03067	Aster-04746
C169v2-03066	Aster-04710
C169v2-03065	Aster-04711
C169v2-03064	Aster-05029
C169v2-03285	Aster-01180
C169v2-03283	Aster-04356
C169v2-03282	Aster-01211
C169v2-03281	Aster-05299
C169v2-03280	Aster-01473

Comme ca je remplis mon hash avec cv et aster et je verifie

Merci

[Jasmine80]

C'est bien, maintenant que tu as trouvé l'erreur, tu n'as plus qu'à la corriger.

[shadow19c]

J'ai essaye de modifier le script mais je ne vois pas d'erreur , on rentre bien les deux nombres a achaque fois et on fait le calcul.
Ca doit etre surement dans le test if ((abs ($num_k - $num_v) == 1)

[Jasmine80]

Tu as donc modifié :

Code :

if ((abs ($num_k - $num_v) == 1) && ($gene_k eq $gene_v) )

en

Code :

if (abs ($num_k - $num_v) == 1)

Si cela ne fonctionne pas, c'est un problème au niveau des expressions régulières (le ^ ou $ en trop)

Code :

1
2
	my ($gene_k, $num_k) = $k =~ m/^(\w+)-(\d+)$/;
	my ($gene_v, $num_v) = $v =~ m/^(\w+)-(\d+)$/;

Fais un copier-coller d'une partie de ton fichier

[shadow19c]

Le copier coller je vous le met sur le site?

[shadow19c]

Par ce que moi je pars d'un fichier qui a cette allure la:

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
 
C169v2-04488	Aster-06787
C169v2-04487	Aster-06725
C169v2-04480	Aster-06724
C169v2-04481	Aster-06792
C169v2-03069	Aster-04749
C169v2-03068	Aster-01309
C169v2-03061	Aster-06908
C169v2-03060	Aster-04730
C169v2-03067	Aster-04746
C169v2-03066	Aster-04710
C169v2-03065	Aster-04711
C169v2-03064	Aster-05029
C169v2-03285	Aster-01180
C169v2-03283	Aster-04356
C169v2-03282	Aster-01211
C169v2-03281	Aster-05299
C169v2-03280	Aster-01473
C169v2-03289	Aster-06954
C169v2-07315	Aster-00866
C169v2-07314	Aster-00182
C169v2-02221	Aster-04882
C169v2-02882	Aster-02471
C169v2-02883	Aster-02472
C169v2-02880	Aster-02469
C169v2-02881	Aster-02470
C169v2-01756	Aster-06737
C169v2-01755	Aster-04564
C169v2-01753	Aster-01750
C169v2-01752	Aster-02418
C169v2-07312	Aster-01105
C169v2-01759	Aster-02413
C169v2-01758	Aster-04902
C169v2-03931	Aster-01485
C169v2-03930	Aster-03765
C169v2-03933	Aster-04375
C169v2-03938	Aster-04358
C169v2-01287	Aster-00793
C169v2-01286	Aster-00794

Par exemple un couple qui fonctionne bien c'est :
C169v2-02882 Aster-02471
C169v2-02883 Aster-02472
C169v2-01287 Aster-00793
C169v2-01286 Aster-00794
Ceci est un exemple car les genes se suivent bien dans le fichier.

[Jasmine80]

Envoyé par shadow19c

Par exemple un couple qui fonctionne bien c'est :
C169v2-02882 Aster-02471
C169v2-02883 Aster-02472

Ceci est un exemple car les genes se suivent bien dans le fichier.

... donc, dans ce fichier, certains couples sont bien trouvés par le programme mais d'autres non? Peux-tu me donner un exemple d'un couple qui devrait être trouvé mais qui ne l'est pas?

Le principe est donc de vérifier que les nombres associés, (02471 et 02472) on une différence de 1 ... est-ce bien ça?

[shadow19c]

Non , le programme me sort aucun couple, en fait c'est moi qui suit parti chercher,
en fait Ce que moi j'avais en tete comme je ne sais pas si tous les genes se suivent dans le fichier c'est de faire une verfication a chaque fois au cas ou si ils etaient nommes ensuite.
En effet, pour faire le calcu de la taille intergeniques jai besoin de bloc de genes Aster et CV.
Donc par exemple jai un gene A de CV suivi de d'un gene B de Cv je dois avoir en face un gene de aster A suivi de gene B de aster.

[shadow19c]

Je vais essayer de faire un schema de ce que je dois construire pour que ca soit plus clair

[shadow19c]

J'ai mis en piece jointe le schema

merci

P.S: EN fait si j'ai resolu le probleme des sens ce qui va me permettre a chaque fois d'aller chercher les fins et Start de chaque gene et pouvoir faire le calcul de la taille de la region intergenique.

Voici un exemple des deux fichiers Cds des genes Aster et CV:

Code :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
 
# 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
C169-scaffold_1	C169v2-00005	21893	-	ORIGINAL JGI	34662	fgenesh1_pm.1_#_3	21893..22051,22201..22293,22590..22747,22864..23102,23192..23364
C169-scaffold_1	C169v2-00006	24500	+	ORIGINAL JGI	38811	fgenesh1_pg.1_#_4	24500..24674,25216..25503,25746..25907,26012..26196,26516..26752,26951..27159,27373..27550,27728..27961,28207..28371,28619..28890,29131..29389
C169-scaffold_1	C169v2-00007	30817	-	ORIGINAL JGI	26526	estExt_Genewise1Plus.C_10009	30817..30879,30998..31038,31263..31332,31515..31717,31942..32029,32208..32219
C169-scaffold_1	C169v2-00008	32689	+	ORIGINAL JGI	11063	e_gw1.1.493.1	32689..32773,32890..32972,33092..33167,33269..33438,33586..33720,34083..34176,34283..34353,34483..34693,34821..34879
C169-scaffold_1	C169v2-00009	35390	-	ORIGINAL JGI	11236	e_gw1.1.376.1	35390..35542,35740..35841,36023..36207,36571..36760,37035..37229
C169-scaffold_1	C169v2-00010	39592	+	ORIGINAL JGI	64051	estExt_Genemark1.C_10009	39592..39597,39755..39901,40021..40124,40250..40385,40560..40688,40859..41028,41177..41294
C169-scaffold_1	C169v2-00011	43327	-	ORIGINAL JGI	64052	estExt_Genemark1.C_10010	43327..43404,43570..43731,43980..44174,44358..44658,44831..44940,45022..45225

Code :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
 
# 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
scaffold_00001	Aster-00001	17594	-	ORIGINAL JGI	14863	gw1.00001.44.1	17594..17852,17997..18023,18357..18572,18689..18759,18842..18967
scaffold_00001	Aster-00002	19928	-	ORIGINAL JGI	21312	e_gw1.00001.25.1	19928..20186,20453..20616,20744..20792,20913..21081,21507..21714
scaffold_00001	Aster-00003	31765	+	ORIGINAL JGI	14113	gw1.00001.37.1	31765..31855,32261..32514,32748..32831,32917..32998,33154..33362,33646..33753
scaffold_00001	Aster-00004	44240	-	ORIGINAL JGI	36055	fgenesh1_pg.00001_#_3	44240..44511,44793..44943,45050..45133,45256..45318,45686..45721,45782..45904,46025..46342,46459..46642,46822..46871
scaffold_00001	Aster-00005	48441	+	ORIGINAL JGI	21325	e_gw1.00001.70.1	48441..48583,48791..48920
scaffold_00001	Aster-00006	55067	-	ORIGINAL JGI	21277	e_gw1.00001.42.1	55067..55339,55579..55734,56180..56338,56530..56657,56838..57039,57120..57170
scaffold_00001	Aster-00007	58248	+	ORIGINAL JGI	36057	fgenesh1_pg.00001_#_5	58248..58402,58562..58728,58853..58974,59081..59191,59517..59666,59867..60088,60293..60521,60838..61073
scaffold_00001	Aster-00008	61579	-	ORIGINAL JGI	33873	fgenesh1_pm.00001_#_1	61579..61967,62122..62293,62396..62484,62589..62681,63116..63272,63649..63882,63915..63942,64316..64491,64815..64889,65118..65392,65635..65917,66265..66381
scaffold_00001	Aster-00009	69555	+	ORIGINAL JGI	21315	e_gw1.00001.28.1	69555..69742,69982..70097,70278..70425,70804..70940,71137..71262,71511..71674,71886..72040,72347..72515,72776..72901
scaffold_00001	Aster-00010	73153	-	ORIGINAL JGI	36060	fgenesh1_pg.00001_#_8	73153..73374,73652..73873,73989..74084
scaffold_00001	Aster-00011	74471	+	ORIGINAL JGI	36061	fgenesh1_pg.00001_#_9	74471..74583,74663..74699,75513..75688,75829..75952

[Jasmine80]

J'essaie de t'aider mais c'es compliqué


D'abord, tu ne peux pas utiliser un hash car il ne garde pas l'ordre des entrées, tu dois utiliser un array

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
open my $fh, '<', 'test.txt' or die;
 
# liste des positions du gène C
my @C_array;
 
# liste des positions du gène Aster
my @Aster_array;
 
while ( my $l = <$fh>){
 
	if ($l =~ m/C169v2-(\d+)\tAster-(\d+)/){
		push (@C_array, $1);
		push (@Aster_array, $2);
	}
}
 
for (my $i =0; $i<$#C_array; $i++){
 
	# on prend les gènes C deux par deux (j'ai gardé ta nomination A et B))
 
	my $A = $C_array[$i];
	my $B = $C_array[$i+1];
 
	# on vérifie que la différence
	# soit au plus de 1
	# On n'en tient pas compte pour Aster
	if (abs ($A - $B) == 1){
 
		# on va dans le fichier CDS
	}
}

je ne comprends pas le SA*SB=1 ... le sens de A doit être le même que celui de B

Peux-tu me donner un fichier CDS
A quoi servent les données sur Aster?

[shadow19c]

Oui, je reprend en fait le fichier reciprocal best blast hits cest un fcihier contenant les genes orthologues de deux genes CV-*** et Aster-***.
Ainsi,la premiere partie du programme doit construire les blocs de genes qui se suivent.
Ensuuite, je vais dans le fichier cds (qui contiennent) les distances des genes et leur sens.
Donc pour chaue bloc par exemple si j'ai C169v2-02882 strand +(A) et C169v2-02883 -(B) SA*SB =-1

Aster-02471(A) + SA*SB=-1
Aster-02472(B) -

Ceci n'est qu 'un exemple,a la fin je dois faire le calcul SA*SB des deux blocs de genes si ils sont egaux je calcul la taille intergenique.

Sauf pour ce cas C169v2-02882 strand +(A) C169v2-02883 -(B) SA*SB =-1
Aster-02471(A) - Aster-02472(B) + SA*SB=-1

On a bien les deux strand -1 sauf que cest faux cette lectutre a enlever aussi

merci

Je vous envoi quelque chose?