#!/usr/bin/env perl

=head1 NAME

getReadsInUnique.pl

=head1 SYNOPSIS
  
  getReadsInUnique.pl 
                      -iRead  <454Contigs.ace.1.pairs>    
                      -cName  <contig name>
                      -cPos   <left|right>
	              -bStart <unique boundry start>
	              -bEnd   <unique boundry end>
                      -oRead  <path to read pair output file>
                      -help 
  

=head1 DESCRIPTION

This script takes read pair file, contig name, position ( left or right of gap ),
unique boundary start, unique boundary end, path to output file of read ids
and an optional file name to output a read pair subset file of read pair lines
in the target region.

All reads in boundaries specified are examined.  Unpaired reads are output 
if a read or portion of the read falls within the boundaries.  Additionally,
paired reads are output the read or portion of the read falls within the 
boundary and it is oriented in such a way that its mate lies in the 
direction of the gap.

=head1 VERSION

$Revision: 1.6 $

$Date: 2009-08-26 17:18:15 $

=head1 HISTORY

=over

=item *

Brian Foster 2008/10/30 Creation

=back

=cut

# Includes

  
use strict;
use Carp;
use vars qw( $RealBin $optiRead $optcName $optcPos $optbStart $optbEnd $optHelp $optoUniq $optoRead);
use FindBin qw($RealBin);

use Getopt::Long;
use Pod::Usage;
use File::Path;
use File::Basename;


#============================================================================#
# INPUT VALIDATION
#============================================================================#

if( !GetOptions(
		"iRead=s" =>\$optiRead,
		"cName=s" =>\$optcName,
		"cPos=s"  =>\$optcPos,
		"bStart=s"=>\$optbStart,
		"bEnd=s" =>\$optbEnd,
		"oRead=s" =>\$optoRead,
		"h"=>\$optHelp)
) {
    printHelp();
}

pod2usage(-verbose=>2) and exit 1 if defined $optHelp;

if( ! defined $optiRead || ! defined $optcName 
 || ! defined $optcPos || ! defined $optbStart|| ! defined $optbEnd){
    printHelp();
}
if (! -e $optiRead){
    print STDERR "Cannot find your read pair file $optiRead\n";
    printHelp();    
}
$optcPos = lc($optcPos);
if ($optcPos ne "left" && $optcPos ne "right"){
    print STDERR "Position $optcPos must be left or right\n";
    print Help();
}
if ( ! -d dirname($optoRead)){
    mkpath (dirname($optoRead));
}
if ( ! -w dirname($optoRead)){
    print STDERR dirname($optoRead), "Not writable\n";
    printHelp();
}


#============================================================================#
# MAIN
#============================================================================#
# get read information from pair file
my (@F,$template);
my (@result);
open (IN , $optiRead) or confess "can't open $optiRead\n";
while(<IN>){
    @F = (split/\t/,$_);

    #reject wrong contig
    next if ($F[4] ne $optcName); 
    
    # reject outside optb range with overlaps allowed
    next if (($F[1] < $optbStart && $F[2] < $optbStart ) 
	  || ($F[1] > $optbEnd && $F[2]   > $optbEnd   )); 

    #reformat ids
    ($template = $F[0]) =~ s/^(.*?)(?:\_|\.).*$/$1/;
    
    #consider pairs
    if ($F[6] =~/(?:M|P)/){
	my $accept = 0;
	# accept if facing into gap while in left contig
	$accept = 1 if ( $F[3]  eq '+' && $optcPos eq 'left');
	
	# accept if facing into gap while in right contig
	$accept = 1 if ( $F[3] eq '-' && $optcPos eq 'right');

	next unless ($accept);
    }
    push (@result,$_);
}
close(IN);
if (@result){
    open(OUT, ">$optoRead") or confess "can't open $optoRead for writing\n";
    print OUT join("",@result);
    close(OUT);
}

#============================================================================#
sub printHelp {
    
    my $optVerbose = shift || 1;
    pod2usage(-verbose=>$optVerbose);
    exit 1;
}

#============================================================================#

exit 0;