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

Annotation of /FigKernelScripts/make_sim_part_data.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 ########################################################################
2 :    
3 :     use DB_File;
4 :    
5 :     use FIG;
6 :     my $fig = new FIG;
7 :    
8 :     use strict;
9 :     my $usage = "usage: make_sim_part_data Dir < stdout.from.partition_prots";
10 :    
11 :     my $dir;
12 :     ($dir = shift @ARGV)
13 :     || die $usage;
14 :    
15 :     if (-d $dir) { die "$dir already exists" }
16 :    
17 :     mkdir($dir,0777) || die "could not make $dir";
18 :    
19 :     my $partD = "$dir/Partitions";
20 :     my $simF = "$dir/SimilarityPartitions";
21 :     my $relF = "$dir/ToSimilarityPartition";
22 :    
23 :     my %prot_hash;
24 :     &build_expand_hash($fig,\%prot_hash);
25 :    
26 :     open(SIMPART,">$simF")
27 :     || die "could not open $simF";
28 :     open(REL,">$relF")
29 :     || die "could not open $relF";
30 :    
31 :     mkdir("$partD",0777) || die "could not make $partD";
32 :    
33 :     my $n = 1;
34 :    
35 :     my $x = <STDIN>;
36 :     while (defined($x) && ($x =~ /^IN\t(\S+)\t(\S+)/))
37 :     {
38 :     my $currset = $1;
39 :     my @setI = ();
40 :     while (defined($x) && ($x =~ /^IN\t(\S+)\t(\S+)/) && ($1 eq $currset))
41 :     {
42 :     push(@setI,$2);
43 :     $x = <STDIN>;
44 :     }
45 :    
46 :     my $bug;
47 :     my @setS = ();
48 :     while (defined($x) && ($x =~ /^SET\t(\S+)\t(\S+)/) && ($1 eq $currset))
49 :     {
50 :     # if ($1 eq 'fig|314291.3.peg.2784') { $bug = 1 }
51 :     push(@setS,$2);
52 :     $x = <STDIN>;
53 :     }
54 :     &process_set($fig,\*SIMPART,\*REL,$partD,\@setI,\@setS,\$n,\%prot_hash,$bug);
55 :     }
56 :     close(SIMPART);
57 :     close(REL);
58 :    
59 :     sub process_set {
60 :     my($fig,$entity_fh,$rel_fh,$partD,$setI,$setS,$nP,$prot_hash,$bug) = @_;
61 :     # if ($bug) { print STDERR &Dumper($setI,$setS) }
62 :    
63 :     my @reduced = &reduce($setS,$prot_hash);
64 :     # if ($bug) { print STDERR &Dumper('reduced1',\@reduced); }
65 :     if (@reduced > 1)
66 :     {
67 :     my $md5;
68 :     foreach $md5 (@reduced)
69 :     {
70 :     print $entity_fh join("\t",($$nP,$md5)),"\n";
71 :     }
72 :     my @reduced_in = &reduce($setI,$prot_hash);
73 :     # if ($bug) { print STDERR &Dumper('reduced_in',\@reduced_in); }
74 :     foreach $md5 (@reduced_in)
75 :     {
76 :     print $rel_fh join("\t",($$nP,$md5)),"\n";
77 :     }
78 :     my $dir1 = $$nP % 1000;
79 :     &FIG::verify_dir("$partD/$dir1");
80 :     open(FASTA,">$partD/$dir1/$$nP") || die "could not make $partD/$dir1/$$nP";
81 :     foreach $md5 (@reduced)
82 :     {
83 :     my @pegs = $fig->pegs_with_md5($md5);
84 :     # if ($bug) { print STDERR &Dumper([$md5,\@pegs]); }
85 :     if (@pegs == 0)
86 :     {
87 :     print STDERR "could not handle hash=$md5\n";
88 :     }
89 :     else
90 :     {
91 :     my($i,$x);
92 :     for ($i=0; ($i < @pegs) && (! ($x = $fig->get_translation($pegs[$i]))); $i++) {}
93 :     if ($i < @pegs)
94 :     {
95 :     print FASTA ">$md5\n$x\n";
96 :     }
97 :     }
98 :     }
99 :     close(FASTA);
100 :    
101 :     $$nP++;
102 :     }
103 :     # if ($bug) { die "aborted" }
104 :     }
105 :    
106 :     sub reduce {
107 :     my($set,$prot_hash) = @_;
108 :    
109 :     my @setE = ();
110 :     my $x;
111 :     foreach $x (@$set)
112 :     {
113 :     my $toL;
114 :     if ($toL = $prot_hash->{$x})
115 :     {
116 :     push(@setE,split(/,/,$toL));
117 :     }
118 :     elsif ($x =~ /^fig\|/)
119 :     {
120 :     push(@setE,$x);
121 :     }
122 :     }
123 :     # print STDERR &Dumper(['expanded',\@setE]);
124 :     my %reduced = map { $_ => 1 } @setE;
125 :     my @reduced_set = sort { &FIG::by_fig_id($a,$b) } keys(%reduced);
126 :    
127 :     my %md5s;
128 :     foreach my $peg (@reduced_set)
129 :     {
130 :     my $md5 = $fig->md5_of_peg($peg);
131 :     if ($md5)
132 :     {
133 :     push(@{$md5s{$md5}},$peg);
134 :     }
135 :     }
136 :     my @md5_set = sort keys(%md5s);
137 :     return @md5_set;
138 :     }
139 :    
140 :     sub build_expand_hash {
141 :     my($fig,$prot_hash) = @_;
142 :    
143 :     # my $prot_hash_tie = tie %$prot_hash, 'DB_File','make_sim_part_prot_hash.db',O_RDWR,0666,$DB_BTREE; return;
144 :     my $prot_hash_tie = tie %$prot_hash, 'DB_File','make_sim_part_prot_hash.db',O_CREAT,0666,$DB_BTREE;
145 :     $prot_hash_tie || die "tie failed";
146 :    
147 :     open(SYMS,"<$FIG_Config::global/peg.synonyms")
148 :     || die "could not open peg.synonyms";
149 :     while (defined($_ = <SYMS>))
150 :     {
151 :     chop;
152 :     my($head,$rest) = split(/\t/,$_);
153 :     if ($rest =~ /fig\|/)
154 :     {
155 :     my @fig = map { ($_ =~ /(fig\|\d+\.\d+\.peg\.\d+)/) ? $1 : () } split(/;/,$rest);
156 :     @fig = grep { $fig->is_complete(&FIG::genome_of($_)) } @fig;
157 :     if (@fig > 0)
158 :     {
159 :     $head =~ s/,\d+$//;
160 :     if ($head =~ /^fig\|/) { push(@fig,$head) }
161 :     $prot_hash->{$head} = join(",",@fig);
162 :     }
163 :     }
164 :     }
165 :     close(SYMS);
166 :     # print STDERR "prot_hash has been constructed\n";
167 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3