#! -w
#
# select-by-gapless-regions  -nogaps foo.nogap
#	[-mismatch m.a2m]	[-err m.err]
#	 < x.a2m > y.a2m
#
# selects those sequences in the a2m file in standard input
#	that have no indels in the regions specified by foo.nogap
#	and outputs them to standard out.
# If the optional -mismatch argument is provided, then rejected
#	sequences are sent there.
# If the optional -err argument is provided, then the reasons for
#	rejection are sent there
#
# The format for the foo.nogap file is a white-space separated
# list of pairs of residues, between which there should be no insertions
# or deletions.
# For example,
#	40 60
#	74 95
# would mean that alignment columns 40, 41, ... ,  60, 74, 75, ... , 95
# must have no deletions, 
# and there can't be insertions after 40, ..., 59, 74, ..., 94.
#
# Comments starting with # are allowed in the input file.



{   use FileHandle;
    STDOUT->autoflush(1);
    STDERR->autoflush(1);
    use English;
    use File::Basename;
    use lib dirname($PROGRAM_NAME);
    use SamT2k;
    

    sub process_command_line;
    sub read_nogap_file($);
    sub check_seq($$);
    sub get_alignment_columns($);    
    # globals
    undef $mis_a2m;
    undef $mis_why;
    undef $nogap_file;
    undef @no_delete;
    undef @no_insert_after;
    
    &process_command_line;
    
    &read_nogap_file($nogap_file);
    # now @nodelete has 1s for alignment columns where deletion is
    #	prohibited, and
    # @no_insert_after has 1s for columns that are not allowed to be
    # follwed by insertions.
    
    open (MIS_A2M, "> $mis_a2m") if (defined $mis_a2m);
    open (MIS_WHY, "> $mis_why") if (defined $mis_why);

    
    while (<STDIN>)
    {   if (/[>]/) 
    	{   &check_seq($id, $seq) if (defined($seq));
	    $id = $_;
	    undef $seq;
	}
	else
	{    $seq .= $_;
	}
    }
    &check_seq($id, $seq) if (defined($seq));
    close MIS_A2M if (defined $mis_a2m);
    close MIS_WHY if (defined $mis_why);
    
}

sub process_command_line
{   local($i) = 0;

    while ($i <= $#ARGV) 
    {	$_ = $ARGV[$i++];
	if (/^-mis/)	{   $mis_a2m  = absolute( $ARGV[$i++]);
			}
	elsif (/^-err/)	{   $mis_why  = absolute( $ARGV[$i++]);
			}
	
	elsif (/^-nogap/)	{   $nogap_file  = absolute( $ARGV[$i++]);
			    (-e $nogap_file) ||
				fatal("$nogap_file doesn't exist");
			}
	else	{ fatal("Invalid option '$ARGV[$i-1]'");
		}
    }

    # Make certain that required parameters were provided.
    defined($nogap_file) 
    	|| fatal("Must specify '-nogaps nogap_file' with $PROGRAM_NAME");

}

# create arrays @nogap_start and @nogap_end
# with the ranges in which indels are prohibited.
sub read_nogap_file($)
{
    my ($nogap_file) = @_;
    open (NOGAP, "<$nogap_file");
    
    undef $start;
    LINE: while(<NOGAP>)
    {   local @words= split;	# split on white space
    	foreach $word (@words)
	{   next LINE if $word =~ /^[#]/;
	    if (!defined($start))
	    {    $start = $word;
	    }
	    else
	    {   for (my $c=$start; $c<$word; $c++)
	        {   $no_delete[$c] = 1;
		    $no_insert_after[$c] = 1;
		}
		$no_delete[$word]=1;
		undef $start;
	    }
	}
    }
    close NOGAP;
    !(defined($start))
      || die "Error: $nogap_file has an odd number of entries, not legal for ranges\n";
    for (my $c=scalar(@no_delete); $c>=0; $c--)
    {  	$no_delete[$c]=0 if (!defined($no_delete[$c]));
   	$no_insert_after[$c]=0 if (!defined($no_insert_after[$c]));
    }
}

# check to see if a sequence has the required gapless regions
# if so, output its id line and the seq to standard out,
# otherwise, output to MIS_A2M (when $mis_a2m defined).
sub check_seq($$)
{
    my ($id, $seq) = @_;
    my @cols = split(//, $seq);
    my $alcol=0;
    my $num_col = scalar(@no_delete);
    for (my $c=0; $c<length($seq); $c++)
    {   my $lett = $cols[$c];
        $alcol++ if ($lett =~ /[-A-Z]/);
	last if ($alcol >= $num_col);	# have checked all important columns
	if ($no_delete[$alcol] && $lett eq "-")
	{ # illegal deletion
	    if (defined $mis_why)
	    {	print MIS_WHY "$id has deletion at $alcol\n";
	    }
	    if (defined $mis_a2m)
	    {    print MIS_A2M $id . $seq;
	    }
	    return;
	}
	if ($no_insert_after[$alcol] && ($lett=~/[a-z]/))
	{ # illegal insertion
	    if (defined $mis_why)
	    {	print MIS_WHY "$id has insertion after $alcol\n";
	    }
	    if (defined $mis_a2m)
	    {    print MIS_A2M $id . $seq;
	    }
	    return;
	}
    }
    # all ok, so print it
    print $id . $seq;
}
