ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/scripts/gmap_slice.psx
Revision: 24
Committed: Tue Jul 26 21:46:39 2011 UTC (8 years, 1 month ago) by gpertea
File size: 2750 byte(s)
Log Message:
Line File contents
1 #!/usr/bin/perl
2 use strict;
3 use FindBin;
4
5 umask 0002;
6 #the line below is needed if pvmsx is used
7 # also, the error condition is set only by the presence of $file
8 #$ENV{'PATH'}=$FindBin::Bin.':'.$ENV{'PATH'};
9
10 my $usage=q{
11 gridx/psx slice processing script - never use by itself.
12
13 Usage:
14
15 gridx [gridx_opts] -i <seqfile> gmap_slice.psx <gmapdb_dir>/<gmapdb> [-P] [gmap_opts]
16
17 Example:
18
19 gridx -U -p 10 -n 2000 -i seqs_to_map.fa gmap_slice.psx /szfs/szannotation/zebrafish/gmapdb/Zv6chr -P -k 0.95 -n 3
20
21 Options:
22 if -P is used, pmap will be run instead of gmap
23 any other options except -P will be passed directly to gmap
24 (defaults are -B 2 -f 2 for batch, GFF gene format output)
25 };
26
27 #so for pvmsx to consider the task was successful, $file must be deleted!
28 #==============
29 # 1 is the name of the fasta sequence input file
30 # 2 is the # of sequences in ${1} should = 1 for this script
31 # 3 is the slice no. being processed by sx
32 # 4 is 0 if not the last file, 1 if the last file
33 # 5 is the # of sequences skipped initially
34 # 6 is the # of sequences to be processed (-1 = ALL)
35 # 7 user parameter
36 # 1 2 3 4 5 6
37 my ($file, $numpass, $slice_num, $last, $skipped, $total, $gmapdbpath, @otheropts)=@ARGV;
38
39
40 my $dbdir=$gmapdbpath;
41 my ($gmapdb) = ($gmapdbpath=~m/\/(\w+)$/);
42 $dbdir=~s/\/(\w+)$//;
43 die "$usage\n" unless $gmapdb && -d $dbdir ;
44
45 my $log_file='log_std';
46 my $err_file='err_log';
47 open(STDERR, '>>'.$err_file);
48 open(STDOUT, '>>'.$log_file);
49 my @gmapopts=grep( !/^\-?P$/, @otheropts );
50 my $protsearch = $#gmapopts < $#otheropts;
51
52
53
54 # my $masklc = $flags=~/M/;
55 # $masklc=0 if $protsearch;
56 #
57 # #my $toskip=($file =~ m/_\@(\d+)_v\d+\.\d+/) ? $1 : $skipped+$numpass*($slice_num-1);
58 # my $fsrch=($masklc)? $file.'_msk.fa' : $file;
59 # my $ftrimcoords=$file.'_msk.trim';
60 # if ($masklc) {
61 # my $prepcmd="mdust -v32 -mN < $file | trimpoly -n20.00 -o $fsrch > $ftrimcoords";
62 # &runCmd($prepcmd);
63 # }
64 my $fsrch=$file;
65
66 my $cmd=$protsearch?'pmap':'gmap';
67 my $gmap_res=$file.".$cmd.gff3";
68 my $gmapopts=join(' ',@gmapopts);
69 my $defopts=($gmapopts=~m/\-B\s*\d/) ? '' : '-B 2';
70 $defopts.=($gmapopts=~m/\-f\s*\d/) ? '' : ' -f 2';
71 $defopts.=' -Y';
72 $cmd.=" -D $dbdir -d $gmapdb $defopts $gmapopts $fsrch > $gmap_res";
73 my $slno=sprintf("slice:%09d",$slice_num);
74 print STDERR ">>$slno: $cmd\n";
75 &runCmd($cmd, $gmap_res);
76
77 print STDERR "<<$slno: done.\n";
78
79 unlink($file);
80 unlink($fsrch);
81 exit 0;
82
83 sub runCmd {
84 my ($docmd, @todel) = @_;
85 my $errmsg = `($docmd) 2>&1`;
86 if ($? || ($errmsg=~/ERROR/si) || ($errmsg=~/Segmentation/si) || ($errmsg=~/Failed/s) || $errmsg=~/Invalid/s) {
87 print STDERR "!Error at:\n$docmd\n";
88 print STDERR "$errmsg\n";
89 foreach (@todel) {
90 unlink($_);
91 }
92 exit(1);
93 }
94 }

Properties

Name Value
svn:executable *