[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.1 - (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 :     use constant BIOINTEL => 2000;
21 :    
22 :     my $STAGE = "sims";
23 :    
24 :     @ARGV == 1 or die "Usage: $0 job-dir\n";
25 :    
26 :     my $jobdir = shift;
27 :    
28 :     -d $jobdir or die "$0: job dir $jobdir does not exist\n";
29 :    
30 :     my $job_id = basename($jobdir);
31 :     my $job = new Job48($job_id);
32 :    
33 :     my $meta = $job->meta;
34 :    
35 :     print "Running job! $jobdir\n";
36 :    
37 :     $job->meta->set_metadata("status.$STAGE", "in_progress");
38 :    
39 :     my $proc = "$jobdir/proc";
40 :     chdir($proc);
41 :    
42 :     #my $fasta = $meta->get_metadata("preprocess.fasta_file");
43 :     my $fasta = "$proc/164.fa";
44 :    
45 :     my $tempf = "/vol/metagenome-48-hour/Data/tblastx-template";
46 :     #my $tempf = "/home/mkubal/TimeLogic_Templates/tblastx-template";
47 :    
48 :     my @databases = ("nr_2007_0208","16s","lsu","gg");
49 :    
50 :     foreach my $db (@databases){
51 :    
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 :     my $job = `/decypher/cli/bin/dc_template -query sequences_to_blast -template $tempf -targ $db`;
67 :     $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 :     }
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