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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3