[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.2 - (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 :     }
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 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3