|
0
|
1 #!/usr/bin/env perl
|
|
|
2
|
|
|
3 use strict;
|
|
|
4 use warnings;
|
|
|
5 use File::Basename;
|
|
|
6
|
|
|
7 # configuration file stuff
|
|
|
8 my %config;
|
|
|
9 my $dirname = dirname(__FILE__);
|
|
|
10 my $tool_dir = shift @ARGV;
|
|
|
11
|
|
|
12 if(not -e "$tool_dir/dedup_realign_bam.loc"){
|
|
|
13 system("cp $dirname/tool-data/dedup_realign_bam.loc $tool_dir/dedup_realign_bam.loc");
|
|
|
14 }
|
|
|
15 open CONFIG, '<', "$tool_dir/dedup_realign_bam.loc";
|
|
|
16 while(<CONFIG>){
|
|
|
17 next if $_ =~ /^#/;
|
|
|
18 (my $key, my $value) = split(/\s+/, $_);
|
|
|
19 $config{$key} = $value;
|
|
|
20 }
|
|
|
21 close CONFIG;
|
|
|
22
|
|
|
23 @ARGV > 4 or die "Usage: $0 <output log> <reference.fa> <temporary directory> <deduped_realigned.bam> <input1.bam> [input2.bam] ...\n";
|
|
|
24
|
|
|
25 my $log = shift @ARGV;
|
|
|
26 my $ref = $config{"dbs_dir"} . shift @ARGV;
|
|
|
27 my $tmpdir = shift @ARGV;
|
|
|
28 my $outfile = shift @ARGV;
|
|
|
29 $SIG{__DIE__} = $SIG{INT} = $SIG{HUP} = $SIG{TERM} = $SIG{ABRT} = \&cleanup;
|
|
|
30 if(@ARGV > 1){
|
|
|
31 system ("merge_bam_headers $log $tmpdir/tmp$$.headers.sam ".join(" ", @ARGV))>>8 and die "merge_bam_headers failed with exit status ", ($?>>8), ", please check the log\n";
|
|
|
32 system("samtools merge -h $tmpdir/tmp$$.headers.sam $tmpdir/tmp$$.bam ".join(" ", @ARGV)." 2>> $log")>>8 and die "Samtools merge failed with exit status ", ($?>>8), ", please check the log\n";
|
|
|
33 system("samtools flagstat $tmpdir/tmp$$.bam > $outfile.before_dedup.flagstat.txt");
|
|
|
34 system("samtools rmdup $tmpdir/tmp$$.bam $tmpdir/tmp$$.rmdup.bam 2>> $log")>>8 and die "Samtools rmdup failed with exit status ", ($?>>8), ", please check the log\n";
|
|
|
35 unlink "$tmpdir/tmp$$.bam";
|
|
|
36 }
|
|
|
37 else{
|
|
|
38 system("samtools flagstat $ARGV[0] > $outfile.before_dedup.flagstat.txt");
|
|
|
39 system("samtools rmdup $ARGV[0] $tmpdir/tmp$$.rmdup.bam 2>> $log")>>8 and die "Samtools rmdup failed with exit status ", ($?>>8), ", please check the log\n";
|
|
|
40 }
|
|
|
41 die "Samtools rmdup did not generate the expected output file" if not -e "$tmpdir/tmp$$.rmdup.bam";
|
|
|
42 die "Samtools generated a blank output file" if -z "$tmpdir/tmp$$.rmdup.bam";
|
|
|
43 system "samtools index $tmpdir/tmp$$.rmdup.bam";
|
|
|
44 system "java -Xmx4G -jar $dirname/GenomeAnalysisTK.jar -I $tmpdir/tmp$$.rmdup.bam -R $ref -T RealignerTargetCreator -o $tmpdir/tmp$$.rmdup.gatk_realigner.intervals 2>> $log";
|
|
|
45 die "GATK did not generate the expected intervals file" if not -e "$tmpdir/tmp$$.rmdup.gatk_realigner.intervals";
|
|
|
46 system "java -Xmx4G -jar $dirname/GenomeAnalysisTK.jar -I $tmpdir/tmp$$.rmdup.bam -R $ref -T IndelRealigner -targetIntervals $tmpdir/tmp$$.rmdup.gatk_realigner.intervals -o $outfile 2>> $log";
|
|
|
47 die "GATK did not generate the expected realigned BAM output file" if not -e $outfile;
|
|
|
48 die "GATK generated a blank output file" if -z $outfile;
|
|
|
49 my $outfile_bai = $outfile;
|
|
|
50 $outfile_bai =~ s/\.bam$/.bai/;
|
|
|
51 if(not -e $outfile_bai){
|
|
|
52 if(not -e "$outfile.bai"){
|
|
|
53 system "samtools index $outfile";
|
|
|
54 }
|
|
|
55 system "cp $outfile.bai $outfile_bai";
|
|
|
56 }
|
|
|
57 else{
|
|
|
58 system "cp $outfile_bai $outfile.bai"; # some tools expect name.bai, some expect name.bam.bai, so provide both
|
|
|
59 }
|
|
|
60 &cleanup;
|
|
|
61
|
|
|
62 sub cleanup{
|
|
|
63 print STDERR @_;
|
|
|
64 unlink "$tmpdir/$tmpdir/tmp$$.headers.sam";
|
|
|
65 unlink "$tmpdir/tmp$$.bam";
|
|
|
66 unlink "$tmpdir/tmp$$.rmdup.bam";
|
|
|
67 unlink "$tmpdir/tmp$$.rmdup.bam.bai";
|
|
|
68 unlink "$tmpdir/tmp$$.rmdup.gatk_realigner.intervals";
|
|
|
69 }
|
