-
Notifications
You must be signed in to change notification settings - Fork 58
/
make_stampy_jobs.pl
56 lines (50 loc) · 1.91 KB
/
make_stampy_jobs.pl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
#!/usr/bin/perl -w
use strict;
use Env qw(HOME);
use File::Spec;
use Getopt::Long;
my $jobdir = File::Spec->catdir($HOME,'jobs');
my $reference_genome;
my ($outdir,$genome_pref);
my $force = 0;
$outdir = ".";
my $solexaqual = 0;
my $tmpdir = "$HOME/bigdata/tmp/stampy";
#my $otheropts = '--substitutionrate=0.05';
my $table = 'table.tab';
GetOptions(
'r|ref:s' => \$reference_genome,
'o|outdir:s' => \$outdir,
'p|prefix:s' => \$genome_pref,
'sq|solexaqual!' => \$solexaqual,
't|tmp:s' => \$tmpdir,
'j|jobdir:s' => \$jobdir,
'f|force!' => \$force);
$outdir = File::Spec->rel2abs($outdir);
if( ! defined $reference_genome ) {
die "must provide reference genome with -r or --ref\n";
}
$reference_genome = File::Spec->rel2abs($reference_genome);
my @fqfiles = @ARGV;
my $opts = $solexaqual ? "--solexa " : "";
for my $fqfile ( @fqfiles ) {
$fqfile = File::Spec->rel2abs($fqfile);
my (undef,$bdir,$fname) = File::Spec->splitpath($fqfile);
next unless ( $fname =~ /(\S+)\.fq$/ || $fname =~ /(\S+)\.fq\.gz$/);
my $base = $1;
$base =~ s/\.trim//;
$base =~ s/\.filter//;
my $smbase = $base;
$smbase =~ s/\.run\d+//;
# $bdir =~ s/\/$//;
open(my $jobfh => ">$jobdir/$base.$genome_pref.stampy.sh") || die $!;
print $jobfh "#!/bin/bash\n",
"#PBS -N $base.stampy -l nodes=1:ppn=4 -j oe \n",
"#PBS -d $bdir\n",
"module load stajichlab\n", "module load stampy\n";
if( $force ) {
unlink("$outdir/$base.$genome_pref.sam");
}
print $jobfh "if [ ! -f $fqfile.recaldata ]; then\n stampy.py -g $reference_genome -h $reference_genome -R $fqfile \nfi\n";
print $jobfh "if [ ! -f $outdir/$base.$genome_pref.sam ]; then\n stampy.py $opts -g $reference_genome -h $reference_genome --readgroup=ID:$base,SM:$smbase --baq --bwaoptions=\"-q10 $reference_genome\" --bwatmpdir=$tmpdir -M $fqfile -f sam -o $outdir/$base.$genome_pref.sam\nfi\n";
}