[Bio] / FortyEight / rp_compute_peer_sims.pl Repository:
ViewVC logotype

Annotation of /FortyEight/rp_compute_peer_sims.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : olson 1.1 #
2 :     # Compute the all-to-all sims between a set of organism directories.
3 :     #
4 :     # Sims are stored in the RAST organism directory, in a directory named
5 :     # by the target genome ID:
6 :     # jobdir/rp/XXXXXX.YY/sims/QQQQQ.RR
7 :     #
8 :     # Within that directory the forward sims (from XXXXXX.YY => QQQQQ.RR) are in a file sims,
9 :     # with a btree index sims.index.
10 :     # Expanded sims are in sims_exp with an index in sims_exp.indx.
11 :     #
12 :     # Flipped sims not stored; they will be found in the other organism's job dir.
13 :     #
14 :    
15 :     use strict;
16 :     use Data::Dumper;
17 :     use Carp;
18 :     use DB_File;
19 :     use strict;
20 :     use FIG;
21 :     use FIG_Config;
22 :     use File::Basename;
23 :     use GenomeMeta;
24 :     use FileHandle;
25 :     use Sim;
26 :    
27 :     my $meta;
28 :    
29 :     @ARGV > 1 or die "Usage: $0 work-dir job-id job-id ...\n";
30 :    
31 :     my $workdir = shift;
32 :     my @job_ids = @ARGV;
33 :    
34 :     if (! -d $workdir)
35 :     {
36 :     mkdir($workdir) or &fatal("mkdir $workdir failed: $!\n");
37 :     }
38 :    
39 :     my @blast_args = qw(-p blastp -m 8 -e 1.0e-5);
40 :    
41 :     #
42 :     # Map the job ids to org dirs, and collect a fasta file.
43 :     #
44 :    
45 :     my $info = new FileHandle "$workdir/INFO", "w";
46 :    
47 :     my @org_dirs;
48 :     my @fasta_files;
49 :     my %genome_to_dir;
50 :    
51 :     for my $job_id (@job_ids)
52 :     {
53 :     my $org_dir;
54 :     my $genome;
55 :     if ($job_id =~ /^\d+$/)
56 :     {
57 :     my $job_dir = "$FIG_Config::rast_jobs/$job_id";
58 :     $genome = &FIG::file_head("$job_dir/GENOME_ID", 1);
59 :     chomp $genome;
60 :     $genome or &fatal("cannot find genome id for job $job_dir\n");
61 :    
62 :     $org_dir = "$job_dir/rp/$genome";
63 :     }
64 :     elsif ($job_id =~ m,^/,)
65 :     {
66 :     # Handle the case where we pass job directories or org dirs.
67 :    
68 :     if (-f "$job_id/GENOME_ID" and -d "$job_id/rp")
69 :     {
70 :     my $job_dir = $job_id;
71 :     $genome = &FIG::file_head("$job_dir/GENOME_ID", 1);
72 :     chomp $genome;
73 :     $genome or &fatal("cannot find genome id for job $job_dir\n");
74 :    
75 :     $org_dir = "$job_dir/rp/$genome";
76 :     }
77 :     elsif (-f "$job_id/GENOME" and -d "$job_id/Features" and $job_id =~ m,/^(\d+\.\d+)$,)
78 :     {
79 :     $org_dir = $job_id;
80 :     $genome = $1;
81 :     }
82 :     else
83 :     {
84 :     &fatal("Unknown jobid $job_id\n");
85 :     }
86 :     }
87 :     else
88 :     {
89 :     &fatal("Unknown jobid $job_id\n");
90 :     }
91 :     -d $org_dir or &fatal("org_dir $org_dir does not exist\n");
92 :    
93 :     my $fasta_file = "$org_dir/Features/peg/fasta";
94 :     if (! -f $fasta_file)
95 :     {
96 :     warn "No fasta found for $org_dir\n";
97 :     next;
98 :     }
99 :    
100 :     push(@fasta_files, $fasta_file);
101 :    
102 :     print $info join("\t", $genome, $org_dir, $fasta_file), "\n";
103 :     push(@org_dirs, $org_dir);
104 :     $genome_to_dir{$genome} = $org_dir;
105 :     }
106 :    
107 :     #
108 :     # Now write the combined fasta.
109 :     #
110 :    
111 :     my $all_fasta = "$workdir/fasta";
112 :     my $fasta_fh = new FileHandle $all_fasta, "w";
113 :    
114 :     #goto xx;
115 :     for my $file (@fasta_files)
116 :     {
117 :     open(my $fh, "<", $file) or &fatal("cannot open $file: $!");
118 :     my $buf;
119 :     while (read($fh, $buf, 4096))
120 :     {
121 :     print $fasta_fh $buf;
122 :     }
123 :     close($fh);
124 :     }
125 :     close($fasta_fh);
126 :    
127 :     &run("$FIG_Config::ext_bin/formatdb", "-p", "t", "-i", $all_fasta);
128 :    
129 :     &run("$FIG_Config::ext_bin/blastall", @blast_args, "-i", $all_fasta, "-d", $all_fasta,
130 :     "-o", "$workdir/sims.raw");
131 :    
132 : olson 1.2 &run("$FIG_Config::bin/reformat_sims $all_fasta < $workdir/sims.raw > $workdir/sims.final");
133 : olson 1.1
134 :     #
135 :     # Read the generated sims (id1, id2, vals).
136 :     # A sim for id1, id2 is written to the file
137 :     # orgdir(genome-of(id1))/sims/genome-of(id2)
138 :     #
139 :     xx:
140 :     my $cur_genome;
141 :     open(S, "<", "$workdir/sims.final") or &fatal("Cannot open $workdir/sims.final: $!");
142 :     my($id1, $id2, $g1, $g2);
143 :     my($last_g1, $last_g2);
144 :     my %fhh;
145 :     my @new_files;
146 :     while (<S>)
147 :     {
148 :     if (($id1, $g1, $id2, $g2) = /^(fig\|(\d+\.\d+)\S+)\t(fig\|(\d+\.\d+)\S+)/)
149 :     {
150 :     my $fh = $fhh{$g1, $g2};
151 :    
152 :     if (!$fh)
153 :     {
154 :     my $sdir = "$genome_to_dir{$g1}/sims";
155 :     -d $sdir or mkdir $sdir or &fatal("cannot mkdir $sdir: $!");
156 :    
157 :     open($fh, ">", "$sdir/$g2") or &fatal("cannot write $sdir/$g2: $!");
158 :     push(@new_files, "$sdir/$g2");
159 :     $fhh{$g1, $g2} = $fh;
160 :     }
161 :     print $fh $_;
162 :     }
163 :     }
164 :     close(S);
165 :     map { close($_) } values(%fhh);
166 :    
167 :     for my $nf (@new_files)
168 :     {
169 :     print "Index $nf\n";
170 :     &index_sims($nf, "$nf.index");
171 :     }
172 :    
173 :     exit(0);
174 :    
175 :    
176 :     #
177 :     # Use the C index_sims_file app to create a berkeley db index
178 :     # of the sims file.
179 :     #
180 :    
181 :     sub index_sims
182 :     {
183 :     my($sims, $index_file) = @_;
184 :    
185 :     my $path = &FIG::find_fig_executable("index_sims_file");
186 :    
187 :     open(IDX, "$path 0 < $sims |") or &fatal("Cannot open index_sims_file pipe: $!\n");
188 :    
189 :     my %index;
190 :     my $tied = tie %index, 'DB_File', $index_file, O_RDWR | O_CREAT, 0666, $DB_BTREE;
191 :    
192 :     $tied or &fatal("Creation of hash $index_file failed: $!\n");
193 :    
194 :     while (<IDX>)
195 :     {
196 :     chomp;
197 :     my($peg, undef, $seek, $len) = split(/\t/);
198 :    
199 :     $index{$peg} = "$seek,$len";
200 :     }
201 :     close(IDX);
202 :    
203 :     $tied->sync();
204 :     untie %index;
205 :     }
206 :    
207 :     sub fatal
208 :     {
209 :     my($msg) = @_;
210 :    
211 :     if ($meta)
212 :     {
213 :     $meta->add_log_entry($0, ['fatal error', $msg]);
214 :     $meta->set_metadata("scenario.running", "no");
215 :     $meta->set_metadata("status.scenario", "error");
216 :     }
217 :    
218 :     croak "$0: $msg";
219 :     }
220 :    
221 :     sub run
222 :     {
223 :     my(@cmd) = @_;
224 :     print "Run @cmd\n";
225 :     my $rc = system(@cmd);
226 :     if ($rc != 0)
227 :     {
228 :     &fatal("error $rc running @cmd\n");
229 :     }
230 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3