[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.1 - (view) (download) (as text)

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3