[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.2 - (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 :     print $fhO join("\t",($peg,$func1,$peg2,$func2,$role,$subsys)),"\n";
52 :     }
53 :     }
54 :     }
55 :     }
56 :    
57 :     my $seq = $peg_to_seq->{$peg};
58 :     if ($seq)
59 :     {
60 :     my $pegs2 = $seq_to_pegs->{$seq};
61 :     if ($pegs2)
62 :     {
63 :     foreach my $peg2 (@$pegs2)
64 :     {
65 :     my $func2 = $peg_to_func->{$peg2};
66 :     if ((! $func2) || ($func1 ne $func2))
67 :     {
68 :     if (! ($seen->{"$func1\t$func2"} || $seen->{"$func2\t$func1"}))
69 :     {
70 :     $seen->{"$func1\t$func2"} = 1;
71 : overbeek 1.2 my @anno_pegs = &get_anno_pegs([$peg,$peg2],$anno,$peg_to_seq,$seq_to_pegs);
72 :     print OUT join("\t",($peg,$func1,$peg2,$func2,join(",",@anno_pegs),$role,$subsys)),"\n";
73 : overbeek 1.1 }
74 :     }
75 :     }
76 :     }
77 :     }
78 :     }
79 :     }
80 :     }
81 :    
82 : overbeek 1.2 sub get_anno_pegs {
83 :     my($initial_pegs,$anno,$peg_to_seq,$seq_to_pegs) = @_;
84 :    
85 :     my %anno_pegs;
86 :     foreach my $peg (@$initial_pegs)
87 :     {
88 :     my $same_seq = $seq_to_pegs->{$peg_to_seq->{$peg}};
89 :     my @apegs = grep { $anno->{&SeedUtils::genome_of($_)} } @$same_seq;
90 :     foreach $_ (@apegs) { $anno_pegs{$_} = 1 }
91 :     }
92 :     return keys(%anno_pegs);
93 :     }
94 :    
95 : overbeek 1.1 sub load_funcs {
96 :     my($dataD) = @_;
97 :    
98 :     my $to_func = {};
99 :    
100 :     foreach my $job (`cut -f 3 $dataD/genomes.with.job`)
101 :     {
102 :     chop $job;
103 :     foreach $_ (`cat /vol/rast-prod/jobs/$job/rp/*/proposed*functions`)
104 :     {
105 :     if ($_ =~ /^(\S+)\t(\S.*\S)/)
106 :     {
107 :     my $peg = $1;
108 :     my $func = $2;
109 :     $func =~ s/\s+\#.*$//;
110 :     $to_func->{$peg} = $func;
111 :     }
112 :     }
113 :     }
114 :    
115 :     foreach my $g (`cat $dataD/anno.seed`)
116 :     {
117 :     chop $g;
118 :     foreach $_ (`cat /vol/mirror-seed/Data.mirror/Organisms/$g/assigned_functions`)
119 :     {
120 :     if ($_ =~ /^(\S+)\t(\S[^\t]*\S)/)
121 :     {
122 :     my $peg = $1;
123 :     my $func = $2;
124 :     $func =~ s/\s+\#.*$//;
125 :     $to_func->{$peg} = $func;
126 :     }
127 :     }
128 :     }
129 :     return $to_func;
130 :     }
131 :    
132 :     sub load_seqs {
133 :     my($dataD) = @_;
134 :    
135 :     my $peg_to_seq = {};
136 :     my $seq_to_pegs = {};
137 :    
138 :     foreach my $job (`cut -f 3 $dataD/genomes.with.job`)
139 :     {
140 :     chop $job;
141 :     $/ = "\n>";
142 :     foreach $_ (`cat /vol/rast-prod/jobs/$job/rp/*/Features/peg/fasta`)
143 :     {
144 :     if ($_ =~ /^>?(\S+)[^\n]*\n(.*)/s)
145 :     {
146 :     my $peg = $1;
147 :     my $seq = $2;
148 :     $seq =~ s/\s//gs;
149 :     $peg_to_seq->{$peg} = $seq;
150 :     push(@{$seq_to_pegs->{$seq}},$peg);
151 :     }
152 :     }
153 :     $/ = "\n";
154 :     }
155 :    
156 :     foreach my $g (`cat $dataD/anno.seed`)
157 :     {
158 :     chop $g;
159 :     $/ = "\n>";
160 :     foreach $_ (`cat /vol/mirror-seed/Data.mirror/Organisms/$g/Features/peg/fasta`)
161 :     {
162 :     if ($_ =~ /^>?(\S+)[^\n]*\n(.*)/s)
163 :     {
164 :     my $peg = $1;
165 :     my $seq = $2;
166 :     $seq =~ s/\s//gs;
167 :     $peg_to_seq->{$peg} = $seq;
168 :     push(@{$seq_to_pegs->{$seq}},$peg);
169 :     }
170 :     }
171 :     $/ = "\n";
172 :     }
173 :     return ($peg_to_seq,$seq_to_pegs);
174 :     }
175 :    
176 :     sub load_fams {
177 :     my($dataD) = @_;
178 :    
179 :     my $peg_to_fam = {};
180 :     my $fam_to_pegs = {};
181 :    
182 :     my $fam = 1;
183 :     foreach $_ (`cat $dataD/protein.families`)
184 :     {
185 :     chop;
186 :     my $pegs = [split(/\t/,$_)];
187 :     foreach my $peg (@$pegs)
188 :     {
189 :     $peg_to_fam->{$peg} = $fam;
190 :     $fam_to_pegs->{$fam} = $pegs;
191 :     }
192 :     $fam++;
193 :     }
194 :     return ($peg_to_fam,$fam_to_pegs);
195 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3