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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3