[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.4 - (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 :     my $tempf = "/vol/metagenome-48-hour/Data/tblastx-template";
44 :     #my $tempf = "/home/mkubal/TimeLogic_Templates/tblastx-template";
45 :    
46 :     my @databases = ("nr_2007_0208","16s","lsu","gg");
47 :    
48 :     foreach my $db (@databases){
49 :    
50 :     my $fa=raelib->read_fasta("$fasta"); # this is a generic method of reading a fasta file into a ref to a hash
51 :     my @keys=keys %$fa;
52 :    
53 :     my @jobs;
54 :     while (@keys)
55 :     {
56 :     # remove however many sequences we want to submit
57 :     my @subset=splice(@keys, 0, BIOINTEL);
58 :    
59 :     # write them to a temporary file
60 :     open(OUT, ">sequences_to_blast") || die "Can't write to sequences_to_blast";
61 :     map {print OUT ">$_\n",$fa->{$_},"\n"} @subset;
62 :     close OUT;
63 :    
64 :     my $job = `/decypher/cli/bin/dc_template -query sequences_to_blast -template $tempf -targ $db`;
65 :     $job =~ s/^OK1\s+//;
66 :     chomp($job);
67 :     push(@jobs,$job);
68 :     }
69 :    
70 :     unlink "sequences_to_blast";
71 :    
72 :     #check status of jobs,
73 :    
74 :     my $jobs_not_done = 1;
75 :     my %completed_jobs;
76 :     while ($jobs_not_done){
77 :     open(RUN,"/decypher/cli/bin/dc_qstatus |");
78 :     my $record_jobs = 0;
79 :     while($_ = <RUN>){
80 :     if($_ =~/SOC/){
81 :     if($record_jobs){
82 :     chomp($_);
83 :     $completed_jobs{$_} = 1;
84 :     }
85 :     }
86 :     if($_ =~/Completed/){
87 :     $record_jobs = 1;
88 :     }
89 :     if($_ =~/Running/){
90 :     $record_jobs = 0;
91 :     }
92 :    
93 :     }
94 :     close(RUN);
95 :    
96 :     $jobs_not_done = 0;
97 :     foreach my $job (@jobs){
98 :     if(!$completed_jobs{$job}){
99 :     $jobs_not_done = 1;
100 :     }
101 :     }
102 : mkubal 1.3 sleep(60);
103 : mkubal 1.1 }
104 :    
105 :     #once all jobs finish, parse results and write to $jobdir/proc
106 :    
107 :     my $prefix;
108 :     my $raw_sims;
109 :     my $expanded_sims;
110 :    
111 :     if($db eq "nr_2007_0208"){
112 :     $prefix = "seed";
113 :     $expanded_sims = "sims.".$prefix.".raw.exp";
114 :     open(EXP,">$expanded_sims");
115 :     }
116 :     else{$prefix = $db;}
117 :    
118 :     $raw_sims = "sims.".$prefix.".raw";
119 :     open(RAW,">$raw_sims");
120 :    
121 :     #use sims hash to avoid duplicate entries
122 :     my %sims;
123 :     foreach my $job (@jobs){
124 :     open(IN,"/decypher/output/$job.out");
125 :     while($_ = <IN>){
126 :     #trim header
127 :     if($_ !~/QUERY/){
128 :     if(! $sims{$_}){
129 :     my @fields =split("\t",$_);
130 :     my $primary = $fields[1];
131 :     if(! $primary){next}
132 :     print RAW "$_";
133 :     if($db eq "nr_2007_0208"){
134 :     if($primary =~/fig/){
135 :     print EXP "$_";
136 :     }
137 :     else{
138 :     my @other_ids = $fig->mapped_prot_ids($primary);
139 :     foreach my $set (@other_ids){
140 :     my($id,$len) = @$set;
141 :     if($id !~/xxx/){
142 :     $fields[1] = $id;
143 :     my $temp_string = join("\t",@fields);
144 :     print EXP "$temp_string";
145 :     }
146 :     }
147 :     }
148 :     }
149 :     $sims{$_} = 1;
150 :     }
151 :     }
152 :     }
153 :     close(IN);
154 :     }
155 : mkubal 1.4
156 :     $job->meta->set_metadata("status.$STAGE", "complete");
157 :     $job->meta->set_metadata("$STAGE.running", "no");
158 :    
159 : mkubal 1.1 }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3