view release on metacpan or  search on metacpan

lib/Bio/ViennaNGS/Bam.pm  view on Meta::CPAN

# -*-CPerl-*-
# Last changed Time-stamp: <2017-06-10 18:20:31 michl>

package Bio::ViennaNGS::Bam;

use strict;
use warnings;
use Exporter;
use Bio::ViennaNGS;
use Bio::DB::Sam 1.37;
use Data::Dumper;
use File::Basename qw(fileparse);
use File::Temp qw(tempfile);
use Path::Class;
use Carp;
use version; our $VERSION = version->declare("$Bio::ViennaNGS::VERSION");

our @ISA = qw(Exporter);
our @EXPORT = ();

our @EXPORT_OK = qw ( split_bam uniquify_bam uniquify_bam2 );

#^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^#
#^^^^^^^^^^^ Subroutines ^^^^^^^^^^#
#^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^#

sub split_bam {
  my %data = ();
  my @NHval = ();
  my @processed_files = ();
  my ($verbose,$nh_warning_issued) = (0)x2;
  my ($bamfile,$reverse,$want_uniq,$want_bed,$dest_dir,$log) = @_;
  my ($bam,$sam,$bn,$path,$ext,$header,$flag,$NH,$eff_strand,$tmp);
  my ($bam_pos,$bam_neg,$tmp_bam_pos,$tmp_bam_neg,$bamname_pos,$bamname_neg);
  my ($bed_pos,$bed_neg,$bedname_pos,$bedname_neg);
  my ($seq_id,$start,$stop,$strand,$target_names,$id,$score);
  my $this_function = (caller(0))[3];
  my %count_entries = (
		       total     => 0,
		       uniq      => 0,
		       pos       => 0,
		       neg       => 0,
		       skip      => 0,
		       cur       => 0,
		       mult      => 0,
		       se_alis   => 0,
		       pe_alis   => 0,
		       flag      => 0,
		       unmapped  => 0,
		      );
  $data{count} = \%count_entries;
  $data{flag} = ();
  $data{nh_issues} = 0; # will be set to 1 if NH attribute is missing
  $nh_warning_issued = 0;

  croak "ERROR [$this_function] $bamfile does not exist\n"
    unless (-e $bamfile);
  croak "ERROR [$this_function] $dest_dir does not exist\n"
    unless (-d $dest_dir);

  (undef,$tmp_bam_pos) = tempfile('BAM_POS_XXXXXXX',UNLINK=>0);
  (undef,$tmp_bam_neg) = tempfile('BAM_NEG_XXXXXXX',UNLINK=>0);

  $bam = Bio::DB::Bam->open($bamfile, "r");
  $header = $bam->header;
  $target_names = $header->target_name;

  ($bn,$path,$ext) = fileparse($bamfile, qr /\..*/);
  unless ($dest_dir =~ /\/$/){$dest_dir .= "/";}
  $bamname_pos = $dest_dir.$bn.".pos".$ext;
  $bamname_neg = $dest_dir.$bn.".neg".$ext;
  $bam_pos = Bio::DB::Bam->open($tmp_bam_pos,'w')
    or croak "ERROR [$this_function] Could not open bam_pos file for writing: $!";
  $bam_neg = Bio::DB::Bam->open($tmp_bam_neg,'w')
    or croak "ERROR [$this_function] Could not open bam_neg file for writing: $!";

  if ($want_bed == 1){
     $bedname_pos = $dest_dir.$bn.".pos.bed";
     $bedname_neg = $dest_dir.$bn.".neg.bed";
     open($bed_pos, ">", $bedname_pos); open($bed_neg, ">", $bedname_neg);
  }

  $bam_pos->header_write($header);$bam_neg->header_write($header);
  if($reverse == 1) {  # switch +/- strand mapping
    $tmp = $bam_pos;$bam_pos = $bam_neg;$bam_neg = $tmp;
    $tmp = $bed_pos;$bed_pos = $bed_neg;$bed_neg = $tmp;
  }

  while (my $read= $bam->read1() ) {
    @NHval = ();
    $data{count}{total}++;

    # collect statistics of SAM flags
    $flag = $read->flag;
    unless (exists $data{flag}{$flag}){
      $data{flag}{$flag} = 0;
    }

lib/Bio/ViennaNGS/Bam.pm  view on Meta::CPAN

    printf STDERR "%20d unmapped reads\n",$data{count}{unmapped};
    printf STDERR "%20d total alignment/read count\n",$data{count}{cur};
    printf STDERR "ERROR: pos alignments + neg alignments != total alignments\n";
    print Dumper(\%data);
    croak $!;
  }
  foreach (keys %{$data{flag}}){
    $data{count}{flag} += $data{flag}{$_};
  }
  unless ($data{count}{flag} == $data{count}{cur} + $data{count}{unmapped}){
    printf STDERR "%20d alignments considered\n",$data{count}{cur};
    printf STDERR "%20d unmapped reads\n",$data{count}{unmapped};
    printf STDERR "%20d alignments found in flag statistics\n",$data{count}{flag};
    printf STDERR "ERROR: alignments considered + unmapped reads != alignments from flag stat\n";
    print Dumper(\%data);
    croak $!;
  }

  # logging output
  open(LOG, ">", $log) or croak $!;
  printf LOG "# bam_split.pl log for file \'$bamfile\'\n";
  printf LOG "#-----------------------------------------------------------------\n";
  printf LOG "%20d total alignments (unique/multi/unmapped)\n",$data{count}{total};
  printf LOG "%20d unique-mappers   (%7.2f%% of total) ",
    $data{count}{uniq},eval(100*$data{count}{uniq}/$data{count}{total});
  if ($data{nh_issues}){printf LOG " *** NH attribute issues ***";}
  printf LOG "\n";
  printf LOG "%20d multi-mappers    (%7.2f%% of total) ",
    $data{count}{mult},eval(100*$data{count}{mult}/$data{count}{total});
  if ($data{nh_issues}){printf LOG " *** NH attribute issues ***";}
  printf LOG "\n";
  printf LOG "%20d skipped\n", $data{count}{skip};
  printf LOG "%20d alignments considered\n", $data{count}{cur};
  printf LOG "%20d paired-end\n", $data{count}{pe_alis};
  printf LOG "%20s single-end\n", $data{count}{se_alis};
  if($data{count}{cur}>0){
    printf LOG "%20d fragments on [+] strand  (%7.2f%% of considered)\n",
      $data{count}{pos},eval(100*$data{count}{pos}/$data{count}{cur});
    printf LOG "%20d fragments on [-] strand  (%7.2f%% of considered)\n",
      $data{count}{neg},eval(100*$data{count}{neg}/$data{count}{cur});
  }
  else{
     printf LOG "%20d fragments on [+] strand\n",$data{count}{pos};
     printf LOG "%20d fragments on [-] strand\n",$data{count}{neg};
  }
  printf LOG "%20d unmapped\n", $data{count}{unmapped};
  printf LOG "#-----------------------------------------------------------------\n";
  printf LOG "Dumper output:\n". Dumper(\%data);
  close(LOG);
  return @processed_files;
}

sub uniquify_bam {
  my %data = ();
  my ($bamfile,$dest,$log) = @_;
  my ($bam, $bn,$path,$ext,$read,$header);
  my ($tmp_uniq,$tmp_mult,$fn_uniq,$fn_mult,$bam_uniq,$bam_mult,$NH);
  my ($count_all,$count_uniq,$count_mult,$unmapped,$nh_warning_issued) = (0)x5;
  my @processed_files = ();
  my @NHval = ();
  my $this_function = (caller(0))[3];
  my %count_entries = (
		       total     => 0,
		       uniq      => 0,
		       mult      => 0,
		       unmapped  => 0,
		       nonhtag   => 0,
		      );
  $data{count} = \%count_entries;
  $data{nh_issues} = 0; # will be set to 1 if NH attribute is missing

  croak "ERROR [$this_function] Cannot find $bamfile\n"
    unless (-e $bamfile);
  croak "ERROR [$this_function] $dest does not exist\n"
    unless (-d $dest);

  ($bn,$path,$ext) = fileparse($bamfile, qr /\..*/);

  (undef,$tmp_uniq) = tempfile('BAM_UNIQ_XXXXXXX',UNLINK=>0);
  (undef,$tmp_mult) = tempfile('BAM_MULT_XXXXXXX',UNLINK=>0);

  $bam     = Bio::DB::Bam->open($bamfile, "r");
  $fn_uniq = file($dest,$bn.".uniq".$ext);
  $fn_mult = file($dest,$bn.".mult".$ext);
  $header  = $bam->header; # TODO: modify header, leave traces ...

  $bam_uniq = Bio::DB::Bam->open($tmp_uniq,'w')
    or croak "ERROR [$this_function] Cannot open temp file for writing: $!";
  $bam_mult = Bio::DB::Bam->open($tmp_mult,'w')
    or croak "ERROR [$this_function] Cannot open temp file for writing: $!";
  $bam_uniq->header_write($header);
  $bam_mult->header_write($header);

  while ($read = $bam->read1() ) {
    $data{count}{total}++;

    # skip unmapped reads
    if ( $read->get_tag_values('UNMAPPED') ){
      $data{count}{unmapped}++;
      next;
    }

    # check if NH (the SAM tag used to indicate multiple mappings) is set
    if ($read->has_tag("NH")) {
      @NHval = $read->get_tag_values("NH");
      $NH = $NHval[0];
      if ($NH == 1) {
	$bam_uniq->write1($read);
	$data{count}{uniq}++;
      }
      else {
	$bam_mult->write1($read);
	$data{count}{mult}++;
      }
    }
    else{ # no NH tag found
      $data{nh_issues} = 1; # set this once and for all
      unless ($nh_warning_issued == 1){
	$data{count}{nonhtag}++;
	carp "ERROR [$this_function] Read ".$read->query->name.
	  " does not have NH attribute\nCannot continue ...";
	$nh_warning_issued = 1;
      }
    }
  }

  unless ($data{count}{uniq} + $data{count}{mult} ==
	  $data{count}{total} - $data{count}{unmapped} - $data{count}{nonhtag}){
    printf STDERR "%20d unique alignments\n",$data{count}{uniq};
    printf STDERR "%20d multiple alignments\n",$data{count}{mult};
    printf STDERR "%20d total alignments\n",$data{count}{total};
    printf STDERR "%20d unmapped reads\n",$data{count}{unmapped};
    croak "ERROR [$this_function] Read counts don't match\n";
  }

  rename ($tmp_uniq, $fn_uniq);
  rename ($tmp_mult, $fn_mult);
  push (@processed_files, ($fn_uniq,$fn_mult));

  if (defined $log){
    my $lf = file($dest,$log);
    open(LOG, ">", $lf) or croak $!;
    printf LOG "%15d reads total\n%15d unique alignments\n%15d multiple alignments\n%15d unmapped\n",
      $data{count}{total},$data{count}{uniq},$data{count}{mult},$data{count}{unmapped};
    close(LOG);
  }
}

sub uniquify_bam2 {
  my %data = ();
  my ($bamfile,$dest,$log) = @_;
  my ($bam, $bn,$path,$ext,$read,$header,$ali);
  my ($tmp_uniq,$tmp_mult,$fn_uniq,$fn_mult,$bam_uniq,$bam_mult,$NH);
  my ($count_all,$count_uniq,$count_mult,$unmapped,$nh_warning_issued) = (0)x5;
  my $pushme = 0;
  my $allgomulti = 0;
  my @processed_files = ();
  my @NHval = ();
  my $lastqname = undef;
  my @band = ();
  my $this_function = (caller(0))[3];
  my %count_entries = (
		       total     => 0,
		       uniq      => 0,
		       mult      => 0,
		       unmapped  => 0,
		       nonhtag   => 0,
		      );
  $data{count} = \%count_entries;
  $data{nh_issues} = 0; # will be set to 1 if NH attribute is missing

  croak "ERROR [$this_function] Cannot find $bamfile\n"
    unless (-e $bamfile);
  croak "ERROR [$this_function] $dest does not exist\n"
    unless (-d $dest);

  ($bn,$path,$ext) = fileparse($bamfile, qr /\..*/);

  (undef,$tmp_uniq) = tempfile('BAM_UNIQ_XXXXXXX',UNLINK=>0);
  (undef,$tmp_mult) = tempfile('BAM_MULT_XXXXXXX',UNLINK=>0);

  $bam     = Bio::DB::Bam->open($bamfile, "r");
  $fn_uniq = file($dest,$bn.".uniq.band.".$ext);
  $fn_mult = file($dest,$bn.".mult.band.".$ext);
  $header  = $bam->header; # TODO: modify header, leave traces ...

  $bam_uniq = Bio::DB::Bam->open($tmp_uniq,'w')
    or croak "ERROR [$this_function] Cannot open temp file for writing: $!";
  $bam_mult = Bio::DB::Bam->open($tmp_mult,'w')
    or croak "ERROR [$this_function] Cannot open temp file for writing: $!";
  $bam_uniq->header_write($header);
  $bam_mult->header_write($header);

  while ($read = $bam->read1() ) {
    $data{count}{total}++;

    # skip unmapped reads
    if ( $read->get_tag_values('UNMAPPED') ){
      $data{count}{unmapped}++;
      next;
    }

    my $qname = $read->qname;
    $qname =~ s/\/[12]$//;

    if (!defined($lastqname) || $qname eq $lastqname ){
      #print "$qname\n";
      push @band, $read;
    }
    else {
      # band is now completely read, processing it ...
      foreach $ali (@band){
	# check if NH (the SAM attribute used to indicate multiple mappings) is set
	if ($ali->has_tag("NH")) {
	  @NHval = $ali->get_tag_values("NH");
	  $NH = $NHval[0];
	  if ($NH > 1) {
	    $allgomulti = 1;
	    $data{count}{mult}++;
	  }
	  else {