[Bio] / FortyEightMeta / mg_compute_sims_on_timelogic.pl Repository:
ViewVC logotype

Annotation of /FortyEightMeta/mg_compute_sims_on_timelogic.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : mkubal 1.1 use FIG;
2 :     use strict;
3 :     use raelib;
4 :     use warnings;
5 :     use FIG_Config;
6 :     use File::Basename;
7 :     use GenomeMeta;
8 :    
9 :     use DBMaster;
10 :     use Job48;
11 : mkubal 1.4
12 :     my $STAGE = "sims";
13 : mkubal 1.1
14 :     my $fig = new FIG;
15 :    
16 :     # this defines the number of sequences for each job.
17 :     use constant BIOBLASTA => 10000;
18 :     use constant BIOINTEL => 2000;
19 :    
20 :     my $STAGE = "sims";
21 :    
22 :     @ARGV == 1 or die "Usage: $0 job-dir\n";
23 :    
24 :     my $jobdir = shift;
25 :    
26 :     -d $jobdir or die "$0: job dir $jobdir does not exist\n";
27 :    
28 :     my $job_id = basename($jobdir);
29 :     my $job = new Job48($job_id);
30 :    
31 :     my $meta = $job->meta;
32 :    
33 :     print "Running job! $jobdir\n";
34 :    
35 :     $job->meta->set_metadata("status.$STAGE", "in_progress");
36 :    
37 :     my $proc = "$jobdir/proc";
38 :     chdir($proc);
39 :    
40 : mkubal 1.2 my $fasta = $meta->get_metadata("preprocess.fasta_file");
41 :    
42 : mkubal 1.1
43 : mkubal 1.5 my $tempf;
44 : mkubal 1.1
45 : mkubal 1.6 my @databases = ("nr_2007_0208","16s","lsu","gg","ssu");
46 : mkubal 1.1
47 :     foreach my $db (@databases){
48 : mkubal 1.5
49 : mkubal 1.8 if($db eq "nr_2007_0208"){ $tempf = "/vol/metagenome-48-hour/Data/tblastx-template-nr";}
50 :     else{$tempf = "/vol/metagenome-48-hour/Data/tblastn-template-others";}
51 : mkubal 1.1
52 :     my $fa=raelib->read_fasta("$fasta"); # this is a generic method of reading a fasta file into a ref to a hash
53 :     my @keys=keys %$fa;
54 :    
55 :     my @jobs;
56 :     while (@keys)
57 :     {
58 :     # remove however many sequences we want to submit
59 :     my @subset=splice(@keys, 0, BIOINTEL);
60 :    
61 :     # write them to a temporary file
62 :     open(OUT, ">sequences_to_blast") || die "Can't write to sequences_to_blast";
63 :     map {print OUT ">$_\n",$fa->{$_},"\n"} @subset;
64 :     close OUT;
65 :    
66 : mkubal 1.9 my $tl_job = `/decypher/cli/bin/dc_template -query sequences_to_blast -template $tempf -targ $db -priority 5`;
67 :     $tl_job =~ s/^OK1\s+//;
68 :     chomp($tl_job);
69 :     push(@jobs,$tl_job);
70 :     $job->meta->set_metadata("$tl_job", "submitted");
71 : mkubal 1.1 }
72 :    
73 :     unlink "sequences_to_blast";
74 :    
75 :     #check status of jobs,
76 :    
77 :     my $jobs_not_done = 1;
78 :     my %completed_jobs;
79 :     while ($jobs_not_done){
80 :     open(RUN,"/decypher/cli/bin/dc_qstatus |");
81 :     my $record_jobs = 0;
82 :     while($_ = <RUN>){
83 :     if($_ =~/SOC/){
84 :     if($record_jobs){
85 :     chomp($_);
86 :     $completed_jobs{$_} = 1;
87 :     }
88 :     }
89 :     if($_ =~/Completed/){
90 :     $record_jobs = 1;
91 : mkubal 1.9
92 : mkubal 1.1 }
93 :     if($_ =~/Running/){
94 :     $record_jobs = 0;
95 :     }
96 :    
97 :     }
98 :     close(RUN);
99 :    
100 :     $jobs_not_done = 0;
101 : mkubal 1.9 foreach my $tl_job (@jobs){
102 :     if(!$completed_jobs{$tl_job}){
103 : mkubal 1.1 $jobs_not_done = 1;
104 :     }
105 : mkubal 1.9 else{
106 :     $job->meta->set_metadata("$tl_job", "completed");
107 :     }
108 : mkubal 1.1 }
109 : mkubal 1.3 sleep(60);
110 : mkubal 1.1 }
111 :    
112 :     #once all jobs finish, parse results and write to $jobdir/proc
113 :    
114 :     my $prefix;
115 :     my $raw_sims;
116 :     my $expanded_sims;
117 :    
118 :     if($db eq "nr_2007_0208"){
119 :     $prefix = "seed";
120 :     $expanded_sims = "sims.".$prefix.".raw.exp";
121 :     open(EXP,">$expanded_sims");
122 :     }
123 :     else{$prefix = $db;}
124 :    
125 :     $raw_sims = "sims.".$prefix.".raw";
126 :     open(RAW,">$raw_sims");
127 :    
128 :     #use sims hash to avoid duplicate entries
129 :     my %sims;
130 :     foreach my $job (@jobs){
131 :     open(IN,"/decypher/output/$job.out");
132 :     while($_ = <IN>){
133 :     #trim header
134 :     if($_ !~/QUERY/){
135 :     if(! $sims{$_}){
136 :     my @fields =split("\t",$_);
137 :     my $primary = $fields[1];
138 :     if(! $primary){next}
139 :     print RAW "$_";
140 :     if($db eq "nr_2007_0208"){
141 :     if($primary =~/fig/){
142 :     print EXP "$_";
143 :     }
144 :     else{
145 :     my @other_ids = $fig->mapped_prot_ids($primary);
146 :     foreach my $set (@other_ids){
147 :     my($id,$len) = @$set;
148 :     if($id !~/xxx/){
149 :     $fields[1] = $id;
150 :     my $temp_string = join("\t",@fields);
151 :     print EXP "$temp_string";
152 :     }
153 :     }
154 :     }
155 :     }
156 :     $sims{$_} = 1;
157 :     }
158 :     }
159 :     }
160 :     close(IN);
161 :     }
162 : mkubal 1.4
163 : mkubal 1.7 }
164 : mkubal 1.4
165 : mkubal 1.7 $job->meta->set_metadata("status.$STAGE", "complete");
166 :     $job->meta->set_metadata("$STAGE.running", "no");

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3