Construction hash intron

**shadow19c** · 10/03/2011, 17h12

Bonjour,
voila toujours des problemes avec le perl,
mon souci c'est que je dois construire un hash de taille d'intron a partir de fichier cds:

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

Je dois recuperer les coordonnees et construire lhash pour chaque id (CV-****),
pour cela si j'ai une coordonne que voici: a..b,c..d,e..f , les "," correspondent aux introns pour cela donc je dois prendre le couple(b,c) et (d,e) et faire la soustraction du start et end de lintron et calculer par la suite la taille de l'intron.

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
 
#!/usr/local/bin/perl
#use strict;
#use warnings;
 
open (DESC1, "coco_cds");
while(<DESC2>){
	next if /^\#/;
	chomp;
	my @exon=(split/,/,$coord);
	(my @intron)=($coord=~s/(/d+,/d+)/g);
	if ($strand eq '-'){
		reverse(@exon,@intron);
	}
	for(0<$i<întron){
		$exon=$exon[i];
		$intron=$intron[i];
		($s,$e)=(split/,/,$intron);
		$taille=abs($e-$s+1);
		push @cds,($exon);
		$aa=int(@cds/3);
		hash{$id}{$aa}=#taille je ne sais pas comment lui dire
	}
}

Voici mon code du debut mais j'arrive pas a regler le probleme du regex, mais aussi a lui dire que hash{$id}{$aa}=taille de l'intron et pour ca je dois faire pour chaque couple d'intron donc je devrais avoir pour cette exemple deux tailles dintron du couple (b,c) et (d,e).

Est ce que vous pouvez m'aider pour mon code .

Merci

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

bonjour,

je n'ai pas compris la fin, donc il faudra modifier mon code :

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
 
#!/usr/local/bin/perl
 
use strict;
use warnings;
 
 
 
open (DESC1, "coco_cds");
 
while (my $line = <DESC1>){
 
	next if ($line =~ m/^\#/);
	chomp $line;
 
	my ( $id, $strand, $coord) = $line =~ m/([\w-]+)\t[\w-]+\t\d+\t([-+])\tORIGINAL JGI\t\d+\t[\#\w\.]+\t([,\d\.]+)/; 
	my @exon = split (/,/, $coord);
 
	my @intron = ($coord =~ m/(\d+,\d+)/g);
 
	if ($strand eq '-'){
		@exon = reverse(@exon);
		@intron = reverse(@intron);
	}
 
	my $cds;
	my %hash; 
 
# 3777..3857,4046..4192,4443..4561,4940..5234,5406..5540,5734..5847,6009..6098,6421..6492
#          324
 
	for my $i (0..$#intron){
 
		my ($s,$e)=(split/,/,$intron[$i]);
		my $taille=abs($e-$s+1);
 
		my ($s2,$e2)=(split/\.\./,$exon[$i]);
		my $taille2 =abs($e2-$s2+1);
 
		$cds += $taille2;
		my $aa=int($cds/3);
 
		$hash{$id}{$aa}=$taille;
 
	}
 
}

**shadow19c** · 11/03/2011, 17h01

Bonjour,
voila le but de mon long programme cest de chercher la taille des introns orthologues:
Pour cel a je pars de deux fichiers Cds Aster et CV(comme d'hab lol), ensuite je construis un hash, apres cela je lis mon fichier de reciprocal best blast hits qui va me permettre davoir le couple Cv-***, Aster_****.
Ensuite je dois construire un pipeline qui va me recuperer pour chaque CV et aster de couple me creer un fichier fasta donc je vais faire un alignement blast de deux sequences.
A partir de ce fichier d'aligenement je vais parcourir les deux sequences jusqu ce que j'obtiens le bon intron calcule lors du 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
91
92
93
94
95
96
97
98
99
100
 
#!/usr/local/bin/perl
 
use strict;
use warnings;
 
 
 
open (DESC1, "coco_cds");
 
while (my $line = <DESC1>){
 
	next if ($line =~ m/^\#/);
	chomp $line;
 
	my ( $id, $strand, $coord) = $line =~ m/([\w-]+)\t[\w-]+\t\d+\t([-+])\tORIGINAL JGI\t\d+\t[\#\w\.]+\t([,\d\.]+)/; 
	my @exon = split (/,/, $coord);
 
	my @intron = ($coord =~ m/(\d+,\d+)/g);
 
	if ($strand eq '-'){
		@exon = reverse(@exon);
		@intron = reverse(@intron);
	}
 
	my $cds;
	my %hash; 
 
# 3777..3857,4046..4192,4443..4561,4940..5234,5406..5540,5734..5847,6009..6098,6421..6492
#          324
 
	for my $i (0..$#intron){
 
		my ($s,$e)=(split/,/,$intron[$i]);
		my $taille=abs($e-$s+1);
 
		my ($s2,$e2)=(split/\.\./,$exon[$i]);
		my $taille2 =abs($e2-$s2+1);
 
		$cds += $taille2;
		my $aa=int($cds/3);
		 $hash{$id}{$aa}=$taille."\n";
 
 
	}
 
}
 
open (DESC2, "aster_cds");
 
while (my $line = <DESC2>){
 
	next if ($line =~ m/^\#/);
	chomp $line;
 
	my ( $id, $strand, $coord) = $line =~ m/([\w-]+)\t[\w-]+\t\d+\t([-+])\tORIGINAL JGI\t\d+\t[\#\w\.]+\t([,\d\.]+)/; 
	my @exon = split (/,/, $coord);
 
	my @intron = ($coord =~ m/(\d+,\d+)/g);
 
	if ($strand eq '-'){
		@exon = reverse(@exon);
		@intron = reverse(@intron);
	}
 
	my $cds;
	my %hash; 
 
# 3777..3857,4046..4192,4443..4561,4940..5234,5406..5540,5734..5847,6009..6098,6421..6492
#          324
 
	for my $i (0..$#intron){
 
		my ($s,$e)=(split/,/,$intron[$i]);
		my $taille=abs($e-$s+1);
 
		my ($s2,$e2)=(split/\.\./,$exon[$i]);
		my $taille2 =abs($e2-$s2+1);
 
		$cds += $taille2;
		my $aa=int($cds/3);
 
			$hash{$id}{$aa}=$taille."\n";
 
 
	}
 
}
 
open (DESC3, "R_B_B_H_A_C.txt");
while (my $ligne1 = <DESC3>){
 
	my ($id1,$id2)=(split/\t/,$ligne1);
 
	system("formatdb -o -pT -i DB.pep");
	system("fastacmd -d DB.pep -s $id1>i");
	system("fastacmd -d DB.pep -s $id2>j");
	system("bl2seq -i i -j j -p blastp -o out.bl -e 1.e-5");
 
}

Voci mon code du debut de pipeline, mon soucis cest que j'arrive pas lui dire de faire le blast2seq que pour le premier couple ensuite je vais faire la recherche de parsing puis des que cest termine de recommence.

Merci

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

pourquoi avoir ajouter un \n?

Code :

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

$hash{$id}{$aa}=$taille."\n";

pour ton blast tu peux utiliser le module Bio::Tools::Run::StandAloneBlast

Code :

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

1
2
3
4
5
6
7
	my @params = ('program' => 'blastn');
	my $factory = Bio::Tools::Run::StandAloneBlast->new(@params);
	$factory->outfile($out_file);
	$factory->executable("bl2seq", "C:/BLAST/bin/bl2seq.exe");
	my $input = Bio::Seq->new(-id=>$id1 ,-seq=> $seq1);
	my $input2 = Bio::Seq->new(-id=>$id2 ,-seq=> $seq2);
	my $bl2seq_report = $factory->bl2seq($input, $input2);

**shadow19c** · 14/03/2011, 10h39

Bonjour,
le "\n" cetait juste pour moi quand je voulais verfier, mais je l'ai enleve, par contre il y a une erreur a cette ligne:

Code :

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

1
2
 
my ( $id, $strand, $coord) = $line =~ m/([\w-]+)\t[\w-]+\t\d+\t([-+])\tORIGINAL JGI\t\d+\t[\#\w\.]+\t([,\d\.]+)/;

Car le id doit correspondre a Aster-** ou CV-** et non a scalfold

J'ai corrige en faisant

Code :

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

1
2
 
my ( $id, $strand, $coord) = $line =~ m/[\w-]+\t([\w-]+)\t\d+\t([-+])\tORIGINAL JGI\t\d+\t[\#\w\.]+\t([,\d\.]+)/;

Merci, en tout cas , par contre je ne peux utiliser Bioseq car ce n'est pas installe sur les machnes et je nai pas les droits donc je lance des commandes systemes , mais ca marche comme ca j'apprend perl en meme temps.

**Jasmine80** · 14/03/2011, 10h55

oui, au lieu de cette regexp, tu peux utiliser split sur \t comme l'avait fait remarquer Norore précédemment, ça a l'avantage d'être simple et général

**shadow19c** · 14/03/2011, 11h12

Oui, aussi, en tout cas maintenant je vois qu'il y a plus de soultions, et j'arrive a code pas mal de chose avec perl en une semaine,
par contre j'ai un souci que voici:

Le but de ce programme c'est de chercher les tailles des introns orthologues.
Pour cela j'ai construit les hashs, pour qu'ils puissent pour chaque couple Aster-** Cv-** calculer la taille des introns par rapport a un aa donne
d'ou le

Code :

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

1
2
3
my $aa=int($cds/3);
 
		$hash{$id}{$aa}=$taille;

Ensuite, je pars d'un fichier reciprocal best blast hits,je recupere les id ensuite,pour chaque id j'aligne les deux sequences grace a blast2seq, ensuite je dois faire un double pointeur qui va essayer de cherche les aa bien aligne entre les deux sequences, en faisant attention au gap( car il ne peut avoir de gap dans un intron), si il ya alignement je donne la taille des introns calcule auparavant dans le hash.

J'ai deux exemples de fichiers du debut du programme pour le hash:

Code :

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

1
2
 
C169-scaffold_3	C169v2-01776	1181194	+	ORIGINAL JGI	12907	e_gw1.3.640.1	1181194..1181244,1181375..1181484,1181561..1181783,1182045..1182257,1182511..1182669

Code :

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

1
2
 
scaffold_00064	Aster-03003	463467	-	ORIGINAL JGI	34851	fgenesh1_pm.00064_#_22	463467..463832,464369..464591,464724..464833,464904..464954

Ce couple C169v2-01776 Aster-03003 permet de faire le blas2seq que j'obtiens:

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
 
Query= lcl|C169v2-01776 jgi|Coc_C169_1|12907|e_gw1.3.640.1
         (251 letters)
 
 
 
>lcl|Aster-03003 geneId:"fgenesh1_pm.00064_#_22" protId:"34851"
           coord:[scaffold_00064;strand="-";exons="463467..463832,
           464369..464591,464724..464833,464904..464954"]
          Length = 249
 
 Score =  305 bits (781), Expect = 6e-88
 Identities = 149/241 (61%), Positives = 179/241 (74%)
 
Query: 1   MRRRPGIQGLQRTVQARDQYKELGKNVAETKLEQMRAQMASFKEHLEEFALKHRDDIRKD 60
           MRRRPGIQGLQ     R  YK LG  V +TKLE ++AQMA FK  LEEFA+KHRDDIR+D
Sbjct: 1   MRRRPGIQGLQLDAATRSSYKVLGNQVQQTKLETIKAQMAVFKRSLEEFAIKHRDDIRQD 60
 
Query: 61  PVFRAQFHAMCANIGVDPLASNKGVWAQVLGFGDFYYELGVQIVEACLASRSLNGGLMDM 120
           P FRAQFH MCAN+GVDPLASNKG++AQ+LG GDFYY+LGVQI+EACLA+R+ NGG++++
Sbjct: 61  PTFRAQFHVMCANVGVDPLASNKGMFAQLLGIGDFYYQLGVQIIEACLATRTHNGGMIEL 120
 
Query: 121 QSLMRYVARRRGSKADPVTEDDVLRAIDKLQXXXXXXXXXXXXDRRLVRSVPGELNTDKN 180
             L R V RRRG  ADP++EDD++RAIDKL+               LVRSVPGELNTDKN
Sbjct: 121 SLLTRLVQRRRGEAADPISEDDIVRAIDKLKKLEGGYNIIKLGGHTLVRSVPGELNTDKN 180
 
Query: 181 QALLLAQGRGHITKRQLTEAYKWEEARVVETLWSLLKEGIAMADDQGPDGQRLYWFPCLE 240
             L LAQ +G ++K Q+    +W E R  ET+ SLLKEG AM DD   DG RLYWFPCL+
Sbjct: 181 AVLELAQQQGFVSKPQILRDLRWAEPRATETVTSLLKEGFAMLDDGALDGVRLYWFPCLD 240
 
Query: 241 S 241
           +
Sbjct: 241 A 241
 
 
Lambda     K      H
   0.320    0.135    0.401 
 
Gapped
Lambda     K      H
   0.267   0.0410    0.140 
 
 
Matrix: BLOSUM62
Gap Penalties: Existence: 11, Extension: 1
Number of Sequences: 1
Number of Hits to DB: 214
Number of extensions: 4
Number of successful extensions: 1
Number of sequences better than 1.0e-05: 1
Number of HSP's better than  0.0 without gapping: 1
Number of HSP's gapped: 1
Number of HSP's successfully gapped: 1
Length of query: 251
Length of database: 249
Length adjustment: 24
Effective length of query: 227
Effective length of database: 225
Effective search space:    51075
Effective search space used:    51075
Neighboring words threshold: 11
Window for multiple hits: 40
X1: 16 ( 7.4 bits)
X2: 38 (14.6 bits)
X3: 64 (24.7 bits)
S1: 41 (21.8 bits)
S2: 72 (32.3 bits)

Ce que je dois faire cest parser et recuperer

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
 
Query: 1   MRRRPGIQGLQRTVQARDQYKELGKNVAETKLEQMRAQMASFKEHLEEFALKHRDDIRKD 60
           MRRRPGIQGLQ     R  YK LG  V +TKLE ++AQMA FK  LEEFA+KHRDDIR+D
Sbjct: 1   MRRRPGIQGLQLDAATRSSYKVLGNQVQQTKLETIKAQMAVFKRSLEEFAIKHRDDIRQD 60
 
Query: 61  PVFRAQFHAMCANIGVDPLASNKGVWAQVLGFGDFYYELGVQIVEACLASRSLNGGLMDM 120
           P FRAQFH MCAN+GVDPLASNKG++AQ+LG GDFYY+LGVQI+EACLA+R+ NGG++++
Sbjct: 61  PTFRAQFHVMCANVGVDPLASNKGMFAQLLGIGDFYYQLGVQIIEACLATRTHNGGMIEL 120
 
Query: 121 QSLMRYVARRRGSKADPVTEDDVLRAIDKLQXXXXXXXXXXXXDRRLVRSVPGELNTDKN 180
             L R V RRRG  ADP++EDD++RAIDKL+               LVRSVPGELNTDKN
Sbjct: 121 SLLTRLVQRRRGEAADPISEDDIVRAIDKLKKLEGGYNIIKLGGHTLVRSVPGELNTDKN 180
 
Query: 181 QALLLAQGRGHITKRQLTEAYKWEEARVVETLWSLLKEGIAMADDQGPDGQRLYWFPCLE 240
             L LAQ +G ++K Q+    +W E R  ET+ SLLKEG AM DD   DG RLYWFPCL+
Sbjct: 181 AVLELAQQQGFVSKPQILRDLRWAEPRATETVTSLLKEGFAMLDDGALDGVRLYWFPCLD 240
 
Query: 241 S 241
           +
Sbjct: 241 A 241

En commencant dabord par la premiere ligne d'aligenement:

Code :

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

1
2
3
4
 
Query: 1   MRRRPGIQGLQRTVQARDQYKELGKNVAETKLEQMRAQMASFKEHLEEFALKHRDDIRKD 60
           MRRRPGIQGLQ     R  YK LG  V +TKLE ++AQMA FK  LEEFA+KHRDDIR+D
Sbjct: 1   MRRRPGIQGLQLDAATRSSYKVLGNQVQQTKLETIKAQMAVFKRSLEEFAIKHRDDIRQD 60

mettre CV index =1 @= MRRRPGIQGLQRTVQARDQYKELGKNVAETKLEQMRAQMASFKEHLEEFALKHRDDIRKD
Aster index =1 @=MRRRPGIQGLQLDAATRSSYKVLGNQVQQTKLETIKAQMAVFKRSLEEFAIKHRDDIRQD

Voila, si je trouve le bon aligenement cest bon, sinon je continue pour lautre ligne, le souci pour moi cest l'aligenement car comment je fais pour lire en meme temps les deux lignes , et aussi si il y a un gap pour la premiere je saute mais je dois continuer pour l'autre ligne.

Merci

PS:

Voici mon code du debut:

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
#!/usr/local/bin/perl
 
use strict;
use warnings;
 
 
 
open (DESC1, "coco_cds");
 
while (my $line = <DESC1>){
 
	next if ($line =~ m/^\#/);
	chomp $line;
 
	my ( $id, $strand, $coord) = $line =~ m/[\w-]+\t([\w-]+)\t\d+\t([-+])\tORIGINAL JGI\t\d+\t[\#\w\.]+\t([,\d\.]+)/;
 
	my @exon = split (/,/, $coord);
 
	my @intron = ($coord =~ m/(\d+,\d+)/g);
 
	if ($strand eq '-'){
		@exon = reverse(@exon);
		@intron = reverse(@intron);
	}
 
	my $cds;
	my %hash; 
 
# 3777..3857,4046..4192,4443..4561,4940..5234,5406..5540,5734..5847,6009..6098,6421..6492
#          324
 
	for my $i (0..$#intron){
 
		my ($s,$e)=(split/,/,$intron[$i]);
		my $taille=abs($e-$s+1);
 
		my ($s2,$e2)=(split/\.\./,$exon[$i]);
		my $taille2 =abs($e2-$s2+1);
 
		$cds += $taille2;
		my $aa=int($cds/3);
		$hash{$id}{$aa}=$taille;
 
 
	}
 
}
 
open (DESC2, "aster_cds");
 
while (my $line = <DESC2>){
 
	next if ($line =~ m/^\#/);
	chomp $line;
 
	my ( $id, $strand, $coord) = $line =~ m/[\w-]+\t([\w-]+)\t\d+\t([-+])\tORIGINAL JGI\t\d+\t[\#\w\.]+\t([,\d\.]+)/; 
 
	my @exon = split (/,/, $coord);
 
	my @intron = ($coord =~ m/(\d+,\d+)/g);
 
	if ($strand eq '-'){
		@exon = reverse(@exon);
		@intron = reverse(@intron);
	}
 
	my $cds;
	my %hash; 
 
# 3777..3857,4046..4192,4443..4561,4940..5234,5406..5540,5734..5847,6009..6098,6421..6492
#          324
 
	for my $i (0..$#intron){
 
		my ($s,$e)=(split/,/,$intron[$i]);
		my $taille=abs($e-$s+1);
 
		my ($s2,$e2)=(split/\.\./,$exon[$i]);
		my $taille2 =abs($e2-$s2+1);
 
		$cds += $taille2;
		 my $aa=int($cds/3);
 
		$hash{$id}{$aa}=$taille;
 
 
	}
 
}
 
open (DESC3, "reciprocal.txt");
while (my $ligne1 = <DESC3>){
 
	my ($id1,$id2)=(split/\t/,$ligne1);
 
	system("formatdb -o -pT -i DB.pep");
	system("fastacmd -d DB.pep -s $id1>cv");
	system("fastacmd -d DB.pep -s $id2>ast");
	system("bl2seq -i cv -j ast -p blastp -o out.bl -e 1.e-5");
 
}

**Jasmine80** · 14/03/2011, 11h19

Envoyé par shadow19c

mettre CV index =1 @= MRRRPGIQGLQRTVQARDQYKELGKNVAETKLEQMRAQMASFKEHLEEFALKHRDDIRKD
Aster index =1 @=MRRRPGIQGLQLDAATRSSYKVLGNQVQQTKLETIKAQMAVFKRSLEEFAIKHRDDIRQD

Voila, si je trouve le bon aligenement cest bon, sinon je continue pour lautre ligne, le souci pour moi cest l'aligenement car comment je fais pour lire en meme temps les deux lignes , et aussi si il y a un gap pour la premiere je saute mais je dois continuer pour l'autre ligne.

[/code]

... je ne comprends pas

c'est dommage que tu ne puisses pas installer le module Boulder::Blast qui permet de lire des fichiers blast

pourquoi ne pas récupérer les 2 séquences en concaténant toutes les lignes puis les analyser dans un seconde temps?

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

Oui, je pense que j'ai mal explique:
Lors du hash, jai construit pour chaque aa l'intron et ceci pour chaque couple d'orthologue.

Donc maintenant, je dois parcourir pour chaque alignement, donc je vais chercher le bon aligenement existant construis dans le hash, et ainsi je construit plusieurs hsp d'orthologues.
C'est pour cela que les index du debut sont importants pour me situe a quel aa, mais aussi en faisant attention au gap signale par XXXXXXXXXXXX, voudra me dire que pour lautre sequence il nexiste pas une correspondance.

**Jasmine80** · 14/03/2011, 11h46

pourquoi ne pas prendre la séquence consensuelle entre la query et la subjct et vérifier qu'elle ne contienne pas d'espace?

**shadow19c** · 14/03/2011, 11h58

C'est a dire,oui aussi je peux faire comme cela, pour une premiere partie , mais ensuite il faudra deux lectures en meme temps pour query et subject.

**Jasmine80** · 14/03/2011, 12h58

voici comment récupérer les séquences mais ensuite que faut-il faire?

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
#!/usr/local/bin/perl
 
use strict;
use warnings;
 
 
 
open my $blast_file, '<', "blast.txt" or die;
 
my $query_seq = '';
my $sbjct_seq = '';
 
while (my $line = <$blast_file>){
 
	if ($line =~ m/Query: \d+ ([A-Z]+)/){
		$query_seq .= $1
	}
	elsif ($line =~ m/Sbjct: \d+ ([A-Z]+)/){
		$sbjct_seq .= $1
	}
}
 
 
print "$query_seq\n$sbjct_seq\n";

**shadow19c** · 14/03/2011, 13h22

Merci, je reexplique, en fait lors de la construction du hash j'aio construis pour chaque couple d'intron , calcule la taille , ainsi que le positionnement de l'aa:

Code :

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

1
2
3
 
my $aa=int($cds/3);
		 $hash{$id}{$aa}=$taille."\n";

ainsi, le but c'est d'aller chercher si elles existent bien dans l'alignement faites blast2seq et si il y a correspondance entre les deux sortir la taille des introns du couple , en sachant qu'il peut y en avoir plusieurs, puisque pour chaque aa on a correpondance de la taille.

Exemple:
si dans ma collection de hash j'ai un aa 125 pour Cv et aster 138:
je recupere lindex des deux entrees ensuite je commence a lire chaque aa et si par hasard aa 125 et aa 138 sont alignes jecris Cv aa 125 taille intron et Aster aa 138 taille intron et je continu pour le reste de l'alignement car il peut y avoir d'autres correspondance.En faisant attention au gap.

Merci

**Jasmine80** · 14/03/2011, 13h32

recoucou,

si j'ai bien compris tu compares des espèces différentes ( par exemple C169v2-01776 et Aster-03003) pour lesquelles tu veux comparer les introns

ce que je ne comprends pas c'est le lien entre les introns et les aa étant donnés que les introns ne sont pas traduits

désolée mais ce n'est pas claire pour moi

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

Re re coucou lol, oui plus serieux en fait Je compare les deux genomes car on se rend compte que il y a plus de genes chez Cv que Aster mais le genome est beaucoup plus grand chez Aster que chez Cv, ceci est du aux introns et donc c'est pour cela que je fais cet etude.

Ceci va pouvoir me dire les tailles des introns orthologues, ainsi je pourrai tracer un plot en voyant que la taille des introns sont superieurs que chez Cv.

**shadow19c** · 14/03/2011, 13h58

Re, ce que j'avais penser c'etait de recuperer l'index cest a dire le chiffre a lentree de chaque segment d'alignement, ensuite mettre chaque qery et subject dans une table pour mieux les lires??
Mais aussi je pense quil faut que je fasse bloc par bloc , ainsi je sais a quel endroit je dois commence.
Mon siouci etant a chaque fois au debut de dire a la boucle de commencer par que on commence la lecture par la premiere lettre et quil commence par une position donne par l'index;

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
 
#!/usr/local/bin/perl
 
use strict;
use warnings;
 
open my $blast_file, '<', "out.bl" or die;
 
my @seq1=(24,48,54,92);
 
my @seq2=(64,137,14);
 
 
while (<$blast_file>){
 
	if ($_ =~ /^Query:\s(\d+)\s+(.+)\s+\d+$/){
#	if (/Query: \d+ ([A-Z]+)/){
		#print $1." : ".$2."\n";
#		$query_seq .= $1;
#		 my @tab1 = split(//,$query_seq);
 
	}
	elsif ($_ =~ /^Sbjct:\s(\d+)\s+(.+)\s+\d+$/){
		#print $1." : ".$2."\n";
		#$sbjct_seq .= $1;
		 my @tab2 = split(//,$2);
 
			foreach my $i(@tab2) {
 				foreach my $j(@seq2){
					if $i==$j;
}
			}
 
	}
 
}

**shadow19c** · 16/03/2011, 09h48

Bon la je suis perdu, j'ai fait un petit code pour m'aider mais je n'y arrive pas
voila je pars avec une seule donne de position:
outblast file

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
 
Query= lcl|C169v2-06337 jgi|Coc_C169_1|67012|estExt_Genemark1.C_120368
         (728 letters)
 
 
 
>lcl|Aster-03202 geneId:"fgenesh1_pg.00067_#_97" protId:"38715"
           coord:[scaffold_00067;strand="+";exons="1001085..1001276
           ,1001343..1001407,1001644..1001734,1001900..1002040,
           1002292..1002453,1002520..1002649,1002929..1003095,
           1003825..1003897,1004217..1004404,1005008..1005196,
           1005549..1005743,1006263..1006406,1006686..1006781,
           1007364..1007516,1008285..1008624,1009037..1009798,
           1010020..1010603,1011111..1011497"]
          Length = 1352
 
 Score = 62.8 bits (151), Expect = 1e-13
 Identities = 59/248 (23%), Positives = 101/248 (40%), Gaps = 20/248 (8%)
 
Query: 160 DSVIAAMSKRGQDRLHHVIRILAEKGHVKELATAFYGNYTVQDLLDATHXXXXXXXXXXX 219
           +  I A+ ++   +L  V+ IL E G V EL T ++GNY++QDL++A             
Sbjct: 264 EQAINALQRKSYQKLTRVVDILIESGKVPELVTHYFGNYSIQDLMEACFKLRKSFEKQGR 323
 
Query: 220 XXXXXSQDIERMLGFSDGQDSLGRMARGIAGAMPEGAMHKQGIYVVLKLLEVGDCADILP 279
                 Q        +   D  GRM + +   +PE   H  G +V+ KLL+     +++ 
Sbjct: 324 KADDAVQYF-----LNPDTDPFGRMFQEVLKTVPESIKHTAGAFVMQKLLQWLTPQEVVT 378
 
Query: 280 LASKLFPGGPKLTTDLHSPEAWRGVKIMTGLVDAMSELAHQGGHMSAAAMGLLDSLCGLY 339
           L      G   +T  +      +G+  M  L       + Q     AA    L ++C  Y
Sbjct: 379 L------GKVSVTHVMDLANDTKGLHFMLKLSQISLTNSEQDYVKDAAQQ--LHAICHQY 430
 
Query: 340 LKDEQAVIKLAQHKNY-GPMLVDAASAGLPRPK----CKQVKSWTNQ--PSPAFGASLVW 392
           L         +   ++ G MLV+A +  +P P        +  W+++    P    S  W
Sbjct: 431 LSSITGSALPSWATSFAGQMLVEAIAGSMPLPSSCTLANNIAGWSSKLLQLPLKDKSTPW 490
 
Query: 393 KILQRLDD 400
            +L++L D
Sbjct: 491 LMLEQLYD 498
 
 
 
 Score = 51.2 bits (121), Expect = 3e-10
 Identities = 34/144 (23%), Positives = 51/144 (35%), Gaps = 13/144 (9%)
 
Query: 10  GTADAGIKCPFCRQFVEGYFPINDGAGKGNXXXXXXXXXXXXXXXXXXXXXXXXXXXXHD 69
           G A +G++CP+CR  V GY P++                                   H 
Sbjct: 63  GQASSGVRCPWCRALVGGYEPMDPSHANDELREANRSARLSALRDIAPMQRQDPALRAH- 121
 
Query: 70  VLPPS-----DWECAKCQNINFSARSKCNKCGQPGPKGAPAIIHGAPSGDVLKCTDKELR 124
             PP      DW C +C+  N   R+ C +C    P             + +  T +++ 
Sbjct: 122 -APPDHRGADDWTCPRCRFDNKPNRATCRECKAKKPDAV------RKDQNYMAATAEDVY 174
 
Query: 125 DCALEKMHPNFQQAFQEAGTPAQT 148
               +K+HP    AFQ AG P  T
Sbjct: 175 TWCTQKLHPALDGAFQAAGVPITT 198
 
 
 
 Score = 40.0 bits (92), Expect = 8e-07
 Identities = 21/86 (24%), Positives = 47/86 (54%), Gaps = 3/86 (3%)
 
Query: 393 KILQRLDDEDENSWVKSIVNELFKCLGSIQSQLKALDMLSDAL-CLKGVPDHAVASHMQM 451
           K++ RLD+++E  WV+ I+ +L   +  +     A D++ + L C+  + D ++ ++ ++
Sbjct: 553 KLILRLDEQEEKEWVQHIMEDLTSSMQKLFKVTHARDVVLEGLTCVHSIDDDSLNTYFEV 612
 
Query: 452 LHDVL--GDARAEEMLKNVERQRGAK 475
           +   L  GD R + ++  +   R AK
Sbjct: 613 ITRQLGKGDMRWQSIMVKIRAARDAK 638
 
 
Lambda     K      H
   0.318    0.135    0.417 
 
Gapped
Lambda     K      H
   0.267   0.0410    0.140 
 
 
Matrix: BLOSUM62
Gap Penalties: Existence: 11, Extension: 1
Number of Sequences: 1
Number of Hits to DB: 2709
Number of extensions: 169
Number of successful extensions: 6
Number of sequences better than 1.0e-05: 1
Number of HSP's better than  0.0 without gapping: 1
Number of HSP's gapped: 4
Number of HSP's successfully gapped: 3
Length of query: 728
Length of database: 1352
Length adjustment: 44
Effective length of query: 684
Effective length of database: 1308
Effective search space:   894672
Effective search space used:   894672

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
 
#!/usr/bin/perl
 
#Variables
my @positions=(24,48,54,92,137,235,275,324); #Liste des positions que tu souhaites verifier
my @positions2=(26,65,145,189);
my $query_seq='';
my $sbjct_seq='';
my ($depart_bloc,$fin_bloc)=(0,0);
my ($query_trouve,$sbjct_trouve,$numero_bloc)=(0,0,1);
 
open (BLAST_FILE,"<out.bl") or die ("Impossible de lire out.bl\n");
while(<BLAST_FILE>)
{
	if(/Query: (\d+)\s+(.+)\s+(\d+)/)	#Ajout de \s+ pour gerer s'il y a plus d'un espacement
	{
		$depart_bloc=$1;  #Index de depart du bloc
		$query_seq=$2;    #Sequence Query
		$fin_bloc=$3;     #Index de fin du bloc
		$query_trouve=1;  #Booleen pour indiquer qu'une sequence a ete trouvee
	}
	elsif(/Sbjct: \d+\s+(.+)/)
	{
		$sbjct_seq=$1;   #Sequence Sbjct
		$sbjct_trouve=1; #Booleen pour indiquer qu'une sequence a ete trouvee		
	}
        #Si les deux sequences "Query" et "Sbjct" ont ete trouvees, on a un nouveau bloc:
	if(($query_trouve==1)&&($sbjct_trouve==1))
	{
                #Separation des caracteres des sequences dans un tableau
		my @tab1 = split(//,$query_seq);
		my @tab2 = split(//,$sbjct_seq);
 
                #Creation des 2 tableaux contenant chacun "$depart_bloc" elements
		my @query=(1..$depart_bloc);
		my @sbjct=(1..$depart_bloc);
                #Ajout de la sequence a la fin de ces tableaux 
                #(ainsi les indices du tableau correspondent aux positions)
		push(@query,@tab1);    
		push(@sbjct,@tab2);
 
 
 
		for (my $i=$depart_bloc,my $j=$depart_bloc;$i<@query,$j<@sbjct;$i++,$j++){
		my $c=$i;
		my $d=$j;
		foreach my $a(@positions){
			if($c eq $a){my $k=1};
			else {my $k=0}
		}
		foreach my $b(@positions2){		
			if($d eq $b){my $q=1};
			else {my $q=0};
		}
		if(($k||$q)==0){next};
		if(($k&&$q)==1){print $d,$c)}}

En fait j'arrive pas a lui dire si j'ai deux couples de positions daller verfier si ils sont allignes , en plus la cest que j'ai fait cest la comparaison de lettre mais cest pas ca ce que je dois faire, cest plus verfier pour chaque position si ils sont alignes.

Exemple:si jai pour la sequence une my @positions=(24,48,54,92,137,235,275,324);
et sequence deux my @positions=(26,65,145,189);

Je dois par exemple si je suis a la position 24 de la premiere sequence je verifie a quelle position je suis dans l'autre sequence si par exemple je suis a une position qui n'est pas dans la collection my @positions=(26,65,145,189); ben je passe a la position suivante pour les deux sequences jusqua ce que , si je suis a la position 65 seq 1 et qu'en face je sois a la position 145 , les deux lettres sont en face la cest bon je renvoi la taille des introns , ainsi de suite.
Je pense qu'il faut que je fasse deux foreach pour chaque table de position?
merci

PS: il faut que je fasse attention aux gaps, si il y en a je dois sauter

**Jasmine80** · 17/03/2011, 10h59

Code :

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

1
2
3
4
5
6
7
                #Creation des 2 tableaux contenant chacun "$depart_bloc" elements
		my @query=(1..$depart_bloc);
		my @sbjct=(1..$depart_bloc);
                #Ajout de la sequence a la fin de ces tableaux 
                #(ainsi les indices du tableau correspondent aux positions)
		push(@query,@tab1);    
		push(@sbjct,@tab2);

... ainsi les indices du tableau correspondent aux positions ... = > pas vraiment

dans ton tableau tu auras des chiffres suivis des acides aminés 1, 2, 3, ... $depart_bloc, M, R, R, ...

**Jasmine80** · 17/03/2011, 11h04

il vaut mieux faire quelque chose comme ça :

Code :

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

1
2
3
4
5
6
		# tableau : clé : AA   valeur : position
		my %positions_query;
		foreach my $aa ( @tab1){
			$positions_query{$aa} = $depart_bloc;
			$depart_bloc ++;
		}

**Jasmine80** · 17/03/2011, 11h14

Voici un début de code, pour la suite je ne sais pas ce qu'il faut faire

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
 
#!/usr/bin/perl
 
use strict;
use warnings;
 
 
#Variables
my @positions=(24,48,54,92,137,235,275,324); #Liste des positions que tu souhaites verifier
my @positions2=(26,65,145,189);
my $query_seq;
my $sbjct_seq;
my ($depart_bloc,$fin_bloc)=(0,0);
my $numero_bloc = 1;
 
open (BLAST_FILE,"<blast.txt") or die ("Impossible de lire blast.txt\n");
while(<BLAST_FILE>)
{
	if(/Query: (\d+)\s+(.+)\s+(\d+)/)	#Ajout de \s+ pour gerer s'il y a plus d'un espacement
	{
		$depart_bloc=$1;  #Index de depart du bloc
		$query_seq=$2;    #Sequence Query
		$fin_bloc=$3;     #Index de fin du bloc
 
	}
	elsif(/Sbjct: \d+\s+(.+)/)
	{
		$sbjct_seq=$1;   #Sequence Sbjct	
	}
 
        #Si les deux sequences "Query" et "Sbjct" ont ete trouvees, on a un nouveau bloc:
	if((defined $query_seq) && (defined $sbjct_seq) && ($query_seq =~ m/\w/) && ($sbjct_seq =~ m/\w/))
	{
                #Separation des caracteres des sequences dans un tableau
		my @tab1 = split(//,$query_seq);
		my @tab2 = split(//,$sbjct_seq);
 
 
		my $r = $depart_bloc;
 
		# tableau : clé : AA   valeur : position
		my %positions_query;
		foreach my $aa ( @tab1){
			if ($aa ne 'X'){
				$positions_query{$aa} = $r;
			}
				$r ++;
		}
 
 
 
		my $s = $depart_bloc;
 
		# tableau : clé : AA   valeur : position
		my %positions_sbjct;
		foreach my $aa ( @tab2){
			if ($aa ne 'X'){
				$positions_sbjct{$aa} = $s;
			}
				$s ++;
		}
 
 
		# on vide les variables
		$query_seq = '';
		$sbjct_seq = '';
 
 
	}
 
}

Construction hash intron

Bioinformatique Perl

Discussions similaires

Partager

Partager