[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.8 - (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.8 my $job = `/decypher/cli/bin/dc_template -query sequences_to_blast -template $tempf -targ $db -priority 9`;
67 : mkubal 1.1 $job =~ s/^OK1\s+//;
68 :     chomp($job);
69 :     push(@jobs,$job);
70 :     }
71 :    
72 :     unlink "sequences_to_blast";
73 :    
74 :     #check status of jobs,
75 :    
76 :     my $jobs_not_done = 1;
77 :     my %completed_jobs;
78 :     while ($jobs_not_done){
79 :     open(RUN,"/decypher/cli/bin/dc_qstatus |");
80 :     my $record_jobs = 0;
81 :     while($_ = <RUN>){
82 :     if($_ =~/SOC/){
83 :     if($record_jobs){
84 :     chomp($_);
85 :     $completed_jobs{$_} = 1;
86 :     }
87 :     }
88 :     if($_ =~/Completed/){
89 :     $record_jobs = 1;
90 :     }
91 :     if($_ =~/Running/){
92 :     $record_jobs = 0;
93 :     }
94 :    
95 :     }
96 :     close(RUN);
97 :    
98 :     $jobs_not_done = 0;
99 :     foreach my $job (@jobs){
100 :     if(!$completed_jobs{$job}){
101 :     $jobs_not_done = 1;
102 :     }
103 :     }
104 : mkubal 1.3 sleep(60);
105 : mkubal 1.1 }
106 :    
107 :     #once all jobs finish, parse results and write to $jobdir/proc
108 :    
109 :     my $prefix;
110 :     my $raw_sims;
111 :     my $expanded_sims;
112 :    
113 :     if($db eq "nr_2007_0208"){
114 :     $prefix = "seed";
115 :     $expanded_sims = "sims.".$prefix.".raw.exp";
116 :     open(EXP,">$expanded_sims");
117 :     }
118 :     else{$prefix = $db;}
119 :    
120 :     $raw_sims = "sims.".$prefix.".raw";
121 :     open(RAW,">$raw_sims");
122 :    
123 :     #use sims hash to avoid duplicate entries
124 :     my %sims;
125 :     foreach my $job (@jobs){
126 :     open(IN,"/decypher/output/$job.out");
127 :     while($_ = <IN>){
128 :     #trim header
129 :     if($_ !~/QUERY/){
130 :     if(! $sims{$_}){
131 :     my @fields =split("\t",$_);
132 :     my $primary = $fields[1];
133 :     if(! $primary){next}
134 :     print RAW "$_";
135 :     if($db eq "nr_2007_0208"){
136 :     if($primary =~/fig/){
137 :     print EXP "$_";
138 :     }
139 :     else{
140 :     my @other_ids = $fig->mapped_prot_ids($primary);
141 :     foreach my $set (@other_ids){
142 :     my($id,$len) = @$set;
143 :     if($id !~/xxx/){
144 :     $fields[1] = $id;
145 :     my $temp_string = join("\t",@fields);
146 :     print EXP "$temp_string";
147 :     }
148 :     }
149 :     }
150 :     }
151 :     $sims{$_} = 1;
152 :     }
153 :     }
154 :     }
155 :     close(IN);
156 :     }
157 : mkubal 1.4
158 : mkubal 1.7 }
159 : mkubal 1.4
160 : mkubal 1.7 $job->meta->set_metadata("status.$STAGE", "complete");
161 :     $job->meta->set_metadata("$STAGE.running", "no");

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3