Mercurial > repos > yusuf > dedup_realign_bam

comparison prep_bams @ 0:6b673ffd9e38 default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
initial commit
author Yusuf Ali <ali@yusuf.email>
date Wed, 25 Mar 2015 13:29:45 -0600
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:6b673ffd9e38
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 }