view release on metacpan or  search on metacpan

scripts/bam_split.pl  view on Meta::CPAN

use File::Basename;
use Path::Class;
use Bio::ViennaNGS::Bam qw(split_bam);
use Bio::ViennaNGS::Util qw(bed_or_bam2bw mkdircheck);

#^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^#
#^^^^^^^^^^ Variables ^^^^^^^^^^^#
#^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^#
my ($bam_p,$bam_n,$bed_p,$bed_n,$size_p,$size_n,$basename,$bamdir,$bamext,$lf);
my ($rev,$wantuniq,$wantbed,$wantnorm,$bw) = (0)x5;
my $logext = ".bam_split.log";
my $cs_in = "-";
my $bam_in = "-";
my $outdir = ".";
my $visdir = "vis";
my $scale   = 1000000;
my @result  = ();

#^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^#
#^^^^^^^^^^^^^^ Main ^^^^^^^^^^^^^#
#^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^#

Getopt::Long::config('no_ignore_case');
pod2usage(-verbose => 1) unless GetOptions("bam=s"      => \$bam_in,
					   "bed"        => sub{$wantbed = 1},
					   "bw"         => sub{$bw = 1},
					   "bwdir=s"    => \$visdir,
					   "cs=s"       => \$cs_in,
					   "norm"       => sub{$wantnorm = 1},
					   "o|out=s"    => \$outdir,
					   "r|reverse"  => sub{$rev = 1},
					   "scale=s"    => \$scale,
					   "uniq"       => sub{$wantuniq = 1},
					   "l|log=s"    => \$logext,
					   "man"        => sub{pod2usage(-verbose => 2)},
					   "help|h"     => sub{pod2usage(1)}
					  );

unless ($bam_in =~ /^\// || $bam_in =~ /\.\//){$bam_in = "./".$bam_in;}
unless (-f $bam_in){
  warn "Could not find input file $bam_in given via --bam option";
  pod2usage(-verbose => 0);
}

if ($bw == 1) {
  unless (-f $cs_in){
    warn "Could not find input file $cs_in given via --cs option";
    pod2usage(-verbose => 0);
  }
  unless ($wantbed == 1){$wantbed = 1;}
}

#TODO check if we are allowed to write to $outdir
unless (-d $outdir){mkdircheck($outdir)}

($basename,$bamdir,$bamext) = fileparse($bam_in,qr/\..*/);

$lf = file($outdir,$basename.$logext);

@result = split_bam($bam_in,$rev,$wantuniq,$wantbed,$outdir,$lf);
$bam_p  = $result[0]; # BAM file containing fragments of [+] strand
$bam_n  = $result[1]; # BAM file containing fragments of [-] strand
$size_p = $result[2]; # of alignments on [+] strand
$size_n = $result[3]; # of alignments on [-] srand
$bed_p  = $result[4]; # BED file containing fragments of [+] strand
$bed_n  = $result[5]; # BED file containing fragments of [-] strand

if ($bw == 1) {
  my $od = dir($outdir,$visdir);
  mkdircheck($od->stringify());
  bed_or_bam2bw("bed",$bed_p,$cs_in,"+",$od,$wantnorm,$size_p,$scale,$lf);
  bed_or_bam2bw("bed",$bed_n,$cs_in,"-",$od,$wantnorm,$size_n,$scale,$lf);
}


__END__


=head1 NAME

bam_split.pl - Split a BAM file by strands

=head1 SYNOPSIS

bam_split.pl [--bam I<FILE>] [options]

=head1 DESCRIPTION

Split a BAM file by strands and create two new BAM file: One
containing all reads that map to the positive strand and another one
with all reads mapped to the negative strand. Optionally filter unique
alignments by inspecting NH:i SAM attribute.

Optionally create bedGraph and (stranded |normalized) bigWig coverage
for UCSC visualization

=head1 OPTIONS

=over

=item B<--bam>

Input file in BAM format

=item B<--bed>

Create a BED6 file for each split BAM file

=item B<--bw>

Create BedGraph and bigWig coverage files for e.g. genome browser
visualization.

=item B<--bwdir>

Directory name for resulting bigWig files. This directory is created
as subdirectory of the output directory. Default is 'vis'.

=item B<--cs>

Chromosome sizes file (required if B<--bw> is given).

=item B<--norm>

Normalize resulting bigWig files