[Bio] / FigKernelScripts / pg_roles_in_some_but_not_X.pl Repository:
ViewVC logotype

Annotation of /FigKernelScripts/pg_roles_in_some_but_not_X.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.8 - (view) (download) (as text)

1 : overbeek 1.1 use strict;
2 :     use Data::Dumper;
3 :     use Getopt::Long;
4 : overbeek 1.2 use SeedEnv;
5 : overbeek 1.1
6 :     my $usage = "usage: pg_in_some_but_not_X -d Data\n";
7 :     my $dataD;
8 :     my $rc = GetOptions('d=s' => \$dataD,);
9 :    
10 :     if ((! $rc) || (! -d $dataD)) { print STDERR $usage; exit }
11 :    
12 : overbeek 1.6 my %anno;
13 :     if (-s "$dataD/anno.seed")
14 :     {
15 :     %anno = map { chop; ($_ => 1) } `cat $dataD/anno.seed`;
16 :     }
17 :     else
18 :     {
19 :     print STDERR "NOTE: you have not included any ASEED genomes in a file $dataD/anno.seed\n";
20 :     }
21 : overbeek 1.2
22 : overbeek 1.1 my$peg_to_func = &load_funcs($dataD);
23 :     my($peg_to_seq,$seq_to_pegs) = &load_seqs($dataD);
24 :     my($peg_to_fam,$fam_to_pegs) = &load_fams($dataD);
25 :    
26 :     open(OUT,">","$dataD/role.inconsistencies") || die "could not open $dataD/role.inconsistencies";
27 :     my %seen;
28 :     foreach my $tuple (map { chop; [split(/\t/,$_)] } `cat $dataD/genomes.with.job.and.genomeID`)
29 :     {
30 :     my($name,undef,$rast_job,$rast_genome) = @$tuple;
31 :     if (open(BINDINGS,"<","/vol/rast-prod/jobs/$rast_job/rp/$rast_genome/Subsystems/bindings"))
32 :     {
33 : overbeek 1.2 &process_genome($peg_to_func,$peg_to_seq,$seq_to_pegs,$peg_to_fam,$fam_to_pegs,\%seen,\*OUT,\*BINDINGS,\%anno);
34 : overbeek 1.1 close(BINDINGS);
35 :     }
36 :     }
37 :     close(OUT);
38 :    
39 :     sub process_genome {
40 : overbeek 1.2 my($peg_to_func,$peg_to_seq,$seq_to_pegs,$peg_to_fam,$fam_to_pegs,$seen,$fhO,$fhB,$anno) = @_;
41 : overbeek 1.1 while (defined($_ = <$fhB>))
42 :     {
43 :     chop;
44 :     my($subsys,$role,$peg) = split(/\t/,$_);
45 :     my $func1 = $peg_to_func->{$peg};
46 :     if ($func1 && (index($role,$func1) >= 0))
47 :     {
48 :     if (my $fam = $peg_to_fam->{$peg})
49 :     {
50 :     my $pegs2 = $fam_to_pegs->{$fam};
51 :     foreach my $peg2 (@$pegs2)
52 :     {
53 :     my $func2 = $peg_to_func->{$peg2};
54 :     if ((! $func2) || ($func1 ne $func2))
55 :     {
56 :     if (! ($seen->{"$func1\t$func2"} || $seen->{"$func2\t$func1"}))
57 :     {
58 :     $seen->{"$func1\t$func2"} = 1;
59 : overbeek 1.3 my @anno_pegs = &get_anno_pegs([$peg,$peg2],$anno,$peg_to_seq,$seq_to_pegs);
60 :     print OUT join("\t",($peg,$func1,$peg2,$func2,join(",",@anno_pegs),$role,$subsys)),"\n";
61 : overbeek 1.1 }
62 :     }
63 :     }
64 :     }
65 :    
66 :     my $seq = $peg_to_seq->{$peg};
67 :     if ($seq)
68 :     {
69 :     my $pegs2 = $seq_to_pegs->{$seq};
70 :     if ($pegs2)
71 :     {
72 :     foreach my $peg2 (@$pegs2)
73 :     {
74 :     my $func2 = $peg_to_func->{$peg2};
75 :     if ((! $func2) || ($func1 ne $func2))
76 :     {
77 :     if (! ($seen->{"$func1\t$func2"} || $seen->{"$func2\t$func1"}))
78 :     {
79 :     $seen->{"$func1\t$func2"} = 1;
80 : overbeek 1.2 my @anno_pegs = &get_anno_pegs([$peg,$peg2],$anno,$peg_to_seq,$seq_to_pegs);
81 :     print OUT join("\t",($peg,$func1,$peg2,$func2,join(",",@anno_pegs),$role,$subsys)),"\n";
82 : overbeek 1.1 }
83 :     }
84 :     }
85 :     }
86 :     }
87 :     }
88 :     }
89 :     }
90 :    
91 : overbeek 1.2 sub get_anno_pegs {
92 :     my($initial_pegs,$anno,$peg_to_seq,$seq_to_pegs) = @_;
93 :    
94 :     my %anno_pegs;
95 :     foreach my $peg (@$initial_pegs)
96 :     {
97 :     my $same_seq = $seq_to_pegs->{$peg_to_seq->{$peg}};
98 :     my @apegs = grep { $anno->{&SeedUtils::genome_of($_)} } @$same_seq;
99 :     foreach $_ (@apegs) { $anno_pegs{$_} = 1 }
100 :     }
101 :     return keys(%anno_pegs);
102 :     }
103 :    
104 : overbeek 1.1 sub load_funcs {
105 :     my($dataD) = @_;
106 :    
107 :     my $to_func = {};
108 :    
109 :     foreach my $job (`cut -f 3 $dataD/genomes.with.job`)
110 :     {
111 :     chop $job;
112 :     foreach $_ (`cat /vol/rast-prod/jobs/$job/rp/*/proposed*functions`)
113 :     {
114 :     if ($_ =~ /^(\S+)\t(\S.*\S)/)
115 :     {
116 :     my $peg = $1;
117 :     my $func = $2;
118 :     $func =~ s/\s+\#.*$//;
119 :     $to_func->{$peg} = $func;
120 :     }
121 :     }
122 :     }
123 : overbeek 1.7 if (-s "$dataD/anno.seed")
124 : overbeek 1.1 {
125 : overbeek 1.7 foreach my $g (`cat $dataD/anno.seed`)
126 : overbeek 1.1 {
127 : overbeek 1.7 chop $g;
128 :     foreach $_ (`cat /vol/mirror-seed/Data.mirror/Organisms/$g/assigned_functions`)
129 : overbeek 1.1 {
130 : overbeek 1.7 if ($_ =~ /^(\S+)\t(\S[^\t]*\S)/)
131 :     {
132 :     my $peg = $1;
133 :     my $func = $2;
134 :     $func =~ s/\s+\#.*$//;
135 :     $to_func->{$peg} = $func;
136 :     }
137 : overbeek 1.1 }
138 :     }
139 :     }
140 :     return $to_func;
141 :     }
142 :    
143 :     sub load_seqs {
144 :     my($dataD) = @_;
145 :    
146 :     my $peg_to_seq = {};
147 :     my $seq_to_pegs = {};
148 :    
149 :     foreach my $job (`cut -f 3 $dataD/genomes.with.job`)
150 :     {
151 :     chop $job;
152 :     $/ = "\n>";
153 :     foreach $_ (`cat /vol/rast-prod/jobs/$job/rp/*/Features/peg/fasta`)
154 :     {
155 : overbeek 1.4 chomp;
156 : overbeek 1.1 if ($_ =~ /^>?(\S+)[^\n]*\n(.*)/s)
157 :     {
158 :     my $peg = $1;
159 :     my $seq = $2;
160 :     $seq =~ s/\s//gs;
161 :     $peg_to_seq->{$peg} = $seq;
162 :     push(@{$seq_to_pegs->{$seq}},$peg);
163 :     }
164 :     }
165 :     $/ = "\n";
166 :     }
167 :    
168 : overbeek 1.8 if (-s "$dataD/anno.seed")
169 : overbeek 1.1 {
170 : overbeek 1.8 foreach my $g (`cat $dataD/anno.seed`)
171 : overbeek 1.1 {
172 : overbeek 1.8 chop $g;
173 :     $/ = "\n>";
174 :     foreach $_ (`cat /vol/mirror-seed/Data.mirror/Organisms/$g/Features/peg/fasta`)
175 : overbeek 1.1 {
176 : overbeek 1.8 chomp;
177 :     if ($_ =~ /^>?(\S+)[^\n]*\n(.*)/s)
178 :     {
179 :     my $peg = $1;
180 :     my $seq = $2;
181 :     $seq =~ s/\s//gs;
182 :     $peg_to_seq->{$peg} = $seq;
183 :     push(@{$seq_to_pegs->{$seq}},$peg);
184 :     }
185 : overbeek 1.1 }
186 : overbeek 1.8 $/ = "\n";
187 : overbeek 1.1 }
188 :     }
189 :     return ($peg_to_seq,$seq_to_pegs);
190 :     }
191 :    
192 :     sub load_fams {
193 :     my($dataD) = @_;
194 :    
195 :     my $peg_to_fam = {};
196 :     my $fam_to_pegs = {};
197 :    
198 :     my $fam = 1;
199 :     foreach $_ (`cat $dataD/protein.families`)
200 :     {
201 :     chop;
202 :     my $pegs = [split(/\t/,$_)];
203 :     foreach my $peg (@$pegs)
204 :     {
205 :     $peg_to_fam->{$peg} = $fam;
206 :     $fam_to_pegs->{$fam} = $pegs;
207 :     }
208 :     $fam++;
209 :     }
210 :     return ($peg_to_fam,$fam_to_pegs);
211 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3