#!/usr/common/usg/languages/perl/5.16.0/bin/perl -w
#
# merauder.pl by Jarrod Chapman <jchapman@lbl.gov> Tue Jun  2 07:48:59 PDT 2009
# Copyright 2009 Jarrod Chapman. All rights reserved.
#

use Getopt::Std;

my %opts = ();
my $validLine = getopts('c:b:s:m:D:S:G:PVR:FQ:ZL:NA', \%opts);
my @required = ("c","b","s","m",);
my $nRequired = 0;
map {$nRequired += exists($opts{$_})} @required;
$validLine &= (($nRequired == @required) && (exists($opts{"S"}) || exists($opts{"G"}) ) );
if ($validLine != 1) {
    print "Usage: ./merauder.pl <-b bmaDataFile> <-s scaffoldReportFile> <-m merSize> <-c contigFastaFile> (<-S seqFileGlob> || <-G gapDataFile>) <<-D minDepth>> <<-P(olymorphicMode?)>> <<-V(erbose?)>>  <<-R repeatCopyCount>> <<-F(astQ?)>> <<-Q qualOffset>> <<-Z(gZipped?)>> <<-L libID>> <<-N(ibbleMode?)>> <<-A(ggressiveClosures?)>>\n";
    exit;
}

my $date = `date`;
chomp $date;
print STDERR "$date\nmerauder4";
while (my ($opt,$val) = each(%opts)) {
    print STDERR " -$opt=$val";
}
print STDERR "\n";

my $bmaDataFile = $opts{"b"};
my $scaffReportFile = $opts{"s"};
my $merSize = $opts{"m"};
my $contigFastaFile = $opts{"c"};
my $minDepth = 2;
if (exists($opts{"D"})) {
    $minDepth = $opts{"D"};
}
my $gapDataFile = undef;
if (exists($opts{"G"})) {
    $gapDataFile = $opts{"G"};
}
my $seqFileGlob = undef;
if (exists($opts{"S"})) {
    $seqFileGlob = $opts{"S"};
}
my $fastQ = 0;
if (exists($opts{"F"})) {
    $fastQ = 1;
}

my $gzipped = 0;
if (exists($opts{"Z"})) {
    $gzipped = 1;
}

my $nibbleMode = 0;
if (exists($opts{"N"})) {
    $nibbleMode = 1;
}

my $aggression = 0;
if (exists($opts{"A"})) {
    $aggression = 1;
}

my $keepLib = undef;
if (exists($opts{"L"})) {
    $keepLib = $opts{"L"};
}

my $qualOffset = 64;
if (exists($opts{"Q"})) {
    $qualOffset = $opts{"Q"};
}

my $polymorphicMode = 0;
if (exists($opts{"P"})) {
    $polymorphicMode = 1;
}

my $verboseMode = 0;
if (exists($opts{"V"})) {
    $verboseMode = 1;
}

my $excludeRepeats = 0;
if (exists($opts{"R"})) {
    $excludeRepeats = $opts{"R"};
}

my $depthInfoAvailable = 0;
my %scaffoldDepth = ();
my %scaffoldLength = ();

my $gapClosingMode = 0;
my %contigsOfInterest = ();
my %scaffoldsOfInterest = ();
if (defined($gapDataFile)) {
    $gapClosingMode = 1;
}

# If gapDataFile has already been generated, decide which contigs/scaffolds to store in memory

if ($gapClosingMode) {

    open (G,$gapDataFile) || die "Couldn't open $gapDataFile\n";
    while (my $line = <G>) {
	chomp $line;

	my ($scaffold,$contig1,$primer1,$contig2,$primer2,$gapSizeEstimate,$gapUncertainty,@fillerSeq) = split(/\t/,$line);
	my ($c1ID) = $contig1 =~ /^(.+)\.[35]$/;
	my ($c2ID) = $contig2 =~ /^(.+)\.[35]$/;
	$contigsOfInterest{$c1ID} = 1;
	$contigsOfInterest{$c2ID} = 1;
	$scaffoldsOfInterest{$scaffold} = 1;
    }
    close G;
}

# Read in contig sequence

print STDERR "Reading $contigFastaFile...\n";

my %contigSequence = ();
my $currentEntry = undef;
my @contigNames = ();
my $seq = "";
open (F,$contigFastaFile) || die "Couldn't open $contigFastaFile\n";
while (my $line = <F>) {
    chomp $line;
    if ($line =~ /^>/) {
	if (defined($currentEntry)) {
	    if (($gapClosingMode==0) || (exists($contigsOfInterest{$currentEntry}))) { 
		$contigSequence{$currentEntry} = $seq;
	    }
	}
	$seq = "";
	($currentEntry) = $line =~ />(.+)/;
	if (($gapClosingMode==0) || (exists($contigsOfInterest{$currentEntry}))) { 
	    push(@contigNames,$currentEntry);
	}
    } else {
	$seq .= $line;
    }
}
close F;

if (defined($currentEntry)) {
    if (($gapClosingMode==0) || (exists($contigsOfInterest{$currentEntry}))) { 
	$contigSequence{$currentEntry} = $seq;
    }
}

print STDERR "Done.\n";

# Read in scaffold report

my @scaffoldIDs = ();
my %scaffoldInfo = ();
my %contigScaffoldMap = ();
print STDERR "Reading $scaffReportFile...\n";
open (S,$scaffReportFile) || die "Couldn't open $scaffReportFile\n";
my $contigOfInterest = 0;
while (my $line = <S>) {
    chomp $line;
    my @cols = split(/\t/,$line);
    my $scaffID = $cols[0];
    
    unless (exists ($scaffoldInfo{$scaffID})) {
        push(@scaffoldIDs,$scaffID);
        $scaffoldInfo{$scaffID} = [];
	$scaffoldDepth{$scaffID} = [0,0];
	$scaffoldLength{$scaffID} = 0;
    }
    if ($cols[1] =~ /^CONTIG/) {
	$contigOfInterest = 0;

	my ($cStrand,$cName) = $cols[2] =~ /^([\+\-])(.+)$/;
#	if (($gapClosingMode==0) || (exists($scaffoldsOfInterest{$scaffID}))) { 
	if (($gapClosingMode==0) || (exists($contigsOfInterest{$cName}))) { 

	    push(@{$scaffoldInfo{$scaffID}},"$cols[2],$cols[3],$cols[4]");
	    $contigOfInterest = 1;

	    if (exists($contigScaffoldMap{$cName})) {
		die "Current version only supports unique contig->scaffold mappings ($cName)\n";
	    } else {
		$contigScaffoldMap{$cName} = "$cStrand,$scaffID,$cols[3],$cols[4]";
	    }
	}

	$scaffoldLength{$scaffID} = $cols[4];

	if ($#cols > 4) {
	    $depthInfoAvailable = 1;
	    my $depth = $cols[5];
	    my $weight = $cols[4]-$cols[3]+1;
	    $scaffoldDepth{$scaffID}->[0] += $depth*$weight;
	    $scaffoldDepth{$scaffID}->[1] += $weight;
	} elsif ($excludeRepeats) {

	    warn "Warning: Depth information not available. Repeat exclusion turned off.\n";
	    $excludeRepeats = 0;

	}

    } elsif ($cols[1] =~ /^GAP/) {

#	if (($gapClosingMode==0) || (exists($scaffoldsOfInterest{$scaffID}))) { 
	if (($gapClosingMode==0) || ($contigOfInterest == 1)) { 
	    push(@{$scaffoldInfo{$scaffID}},"$cols[2],$cols[3]");	
	}
    }
}
close S;

my $peakDepth = 0;
if ($depthInfoAvailable) {
    my %depthHist = ();
    foreach my $scaffID (@scaffoldIDs) {
	my $scaffWeight = $scaffoldDepth{$scaffID}->[1];
	my $meanDepth = $scaffoldDepth{$scaffID}->[0] / $scaffWeight;
	$scaffoldDepth{$scaffID} = $meanDepth;
	my $depthBin = sprintf("%.0f",$meanDepth);
	$depthHist{$depthBin} += $scaffWeight;
    }
    my $maxWeight = 0;
    while (my ($depthBin,$binWeight) = each(%depthHist)) {
	if ($binWeight > $maxWeight) {
	    $peakDepth = $depthBin;
	    $maxWeight = $binWeight;
	}
    }
    print STDERR "Modal scaffold depth $peakDepth\n";
    foreach my $scaffID (@scaffoldIDs) {
	$scaffoldDepth{$scaffID} /= $peakDepth;
    }
}

print STDERR "Done.\n";

# Build scaffold data structures

my %gapInfo = ();
my @gapOrder = ();
my @gapBoundaries = ();
my %scaffGapRange = ();
my $nGaps = 0;

foreach my $scaffID (@scaffoldIDs) {
    my @scaffSegments = @{$scaffoldInfo{$scaffID}};

    my $nScaffSegments = scalar(@scaffSegments);

#    print STDERR "scaffID: $scaffID: nScaffSegments = $nScaffSegments\n";

    if ($nScaffSegments > 1) {

	$scaffGapRange{$scaffID} = "$nGaps.";

	for (my $cIndex = 0; $cIndex < $nScaffSegments-2; $cIndex += 2) {
	    my $contigInfo1 = $scaffSegments[$cIndex];
	    my ($cID1,$sStart1,$sEnd1) = split(/\,/,$contigInfo1);
	    my ($cStrand1,$cName1) = $cID1 =~ /^([\+\-])(.+)$/;
	    my $cSeq1 = $contigSequence{$cName1};
	    my $contig1 = $cName1;
#	    print STDERR "[$cIndex][$contigInfo1][$cName1][$cSeq1]\n";
	    if ($cStrand1 eq "-") {
		my $rc = reverse($cSeq1);
		$rc =~ tr/acgtACGT/tgcaTGCA/;
		$cSeq1 = $rc;
		$contig1 .= ".5";
	    } else {
		$contig1 .= ".3";
	    }
#	    unless (defined($cSeq1)) {
#		print "cSeq1 undefined for $cName1\n";
#		die;
#	    }
	    my $primer1 = substr($cSeq1,-$merSize);
	    
	    my $gapInfo = $scaffSegments[$cIndex+1];
	    my ($gapSize,$gapUncertainty) = split(/\,/,$gapInfo);

	    my $contigInfo2 = $scaffSegments[$cIndex+2];
	    my ($cID2,$sStart2,$sEnd2) = split(/\,/,$contigInfo2);
	    my ($cStrand2,$cName2) = $cID2 =~ /^([\+\-])(.+)$/;
	    my $cSeq2 = $contigSequence{$cName2};
	    my $contig2 = $cName2;
	    if ($cStrand2 eq "-") {
		my $rc = reverse($cSeq2);
		$rc =~ tr/acgtACGT/tgcaTGCA/;
		$cSeq2 = $rc;
		$contig2 .= ".3";
	    } else {
		$contig2 .= ".5";
	    }
	    my $primer2 = substr($cSeq2,0,$merSize);

	    $gapInfo{$contig1} = [$nGaps,$scaffID,$contig1,$primer1,$gapSize,$contig2,$primer2,$gapUncertainty];
	    $gapInfo{$contig2} = [$nGaps,$scaffID,$contig1,$primer1,$gapSize,$contig2,$primer2,$gapUncertainty];
	    push(@gapOrder,$contig1);

	    push(@gapBoundaries,"$sEnd1:$sStart2");
	    $nGaps++;
	}

	$scaffGapRange{$scaffID} .= "$nGaps";

    }
}

print STDERR "Done. Found $nGaps (potentially) closable gaps.\n";

my $nMappedReads = 0;

my %projectors = ();
my %fillers = ();


my %qBins = ();
$qBins{"A0"} = 0;
$qBins{"A1"} = 1;
$qBins{"A2"} = 2;
$qBins{"A3"} = 3;
$qBins{"C0"} = 4;
$qBins{"C1"} = 5;
$qBins{"C2"} = 6;
$qBins{"C3"} = 7;
$qBins{"G0"} = 8;
$qBins{"G1"} = 9;
$qBins{"G2"} = 10;
$qBins{"G3"} = 11;
$qBins{"T0"} = 12;
$qBins{"T1"} = 13;
$qBins{"T2"} = 14;
$qBins{"T3"} = 15;


# Read in bmaData for each blastMapAnalysis

my $nBlastMaps = 0;
my %mapInfo = ();
my @libs = ();
my @bmaFiles = ();
my $maxInsertSize = 0;
open (B,$bmaDataFile) || die "Couldn't open $bmaDataFile\n";
while (my $line = <B>) {
    chomp $line;
    my ($lib,$insertSize,$stdDev,$bmaFile,$gapModelFile) = split(/\t/,$line);

    if (defined($keepLib)) {
	unless ($lib eq $keepLib) {
	    next;
	}
    }

    if ($insertSize > $maxInsertSize) {
	$maxInsertSize = $insertSize;
    }

    $mapInfo{$lib} = [$insertSize,$stdDev,$bmaFile,$gapModelFile];
    push(@libs,$lib);
    push(@bmaFiles,$bmaFile);
    $nBlastMaps++;
}
close B;

# If gapDataFile has already been generated, close gaps

if (defined($gapDataFile)) {

    my $success = 0;
    my $nSuccess = 0;
    my $nFailure = 0;
    my $maxGapZ = 3;
    my $gapUncertaintyWiggle = 5.5;

    my $maxFillReads = 5000;

    open (G,$gapDataFile) || die "Couldn't open $gapDataFile\n";
    while (my $line = <G>) {
	chomp $line;

	my $reportLine = "";
	my $reportStatus = "";
	my $reportNote = "";

	$success = 0;
	my ($scaffold,$contig1,$primer1,$contig2,$primer2,$gapSizeEstimate,$gapUncertainty,@fillerSeq) = split(/\t/,$line);
	$reportLine = "REPORT:\t$scaffold\t$contig1\t$primer1\t$contig2\t$primer2\t$gapSizeEstimate\t$gapUncertainty";

	my $nFillReads = scalar(@fillerSeq);
	my $gapInfo = "[$nFillReads\t$primer1\t$gapSizeEstimate\t$primer2]";

	if ($nFillReads > $maxFillReads) {
	    print STDERR "*************\nGap excluded - excessive filler reads. (nFillReads = $nFillReads) $scaffold:$contig1<-[$gapSizeEstimate +/- $gapUncertainty]->$contig2\n";
	    $reportStatus = "FAILED\tExcessReads=$nFillReads\n";
	    print STDERR "$reportLine\t$reportStatus";
	    next;
	}

	if ($excludeRepeats) {
	    my $scaffDepth = $scaffoldDepth{$scaffold};
	    if ($scaffDepth > $excludeRepeats) {
		print STDERR "*************\nRepeat scaffold gap excluded. (depth = $scaffDepth) $scaffold:$contig1<-[$gapSizeEstimate +/- $gapUncertainty]->$contig2\n";
		$reportStatus = "FAILED\tscaffDepth=$scaffDepth\n";
		print STDERR "$reportLine\t$reportStatus";
		next;
	    }
	}


	print STDERR "*************\nAttempt to close $scaffold:$contig1<-[$gapSizeEstimate +/- $gapUncertainty]->$contig2\n";

	$gapUncertainty = $maxGapZ*$gapUncertainty + $gapUncertaintyWiggle;

	my ($c1ID) = $contig1 =~ /^(.+)\.[35]$/;
	unless (exists($contigSequence{$c1ID})) {
	    die "[$c1ID] no seq\n";
	}
	unless (exists($contigScaffoldMap{$c1ID})) {
	    die "[$c1ID] no scaffMapInfo\n";
	}
	my $c1Seq = $contigSequence{$c1ID};
	my $c1Info = $contigScaffoldMap{$c1ID};
	my ($c1Strand) = $c1Info =~ /^([\+\-])/;
	if ($c1Strand eq "-") {
	    $c1Seq = reverse($c1Seq);
	    $c1Seq =~ tr/acgtACGT/tgcaTGCA/;
	}
	
	my ($c2ID) = $contig2 =~ /^(.+)\.[35]$/; 
	unless (exists($contigSequence{$c2ID})) {
	    die "[$c2ID] no seq\n";
	}
	unless (exists($contigScaffoldMap{$c2ID})) {
	    die "[$c2ID] no scaffMapInfo\n";
	}
	my $c2Seq = $contigSequence{$c2ID};
	my $c2Info = $contigScaffoldMap{$c2ID};
	my ($c2Strand) = $c2Info =~ /^([\+\-])/;
	if ($c2Strand eq "-") {
	    $c2Seq = reverse($c2Seq);
	    $c2Seq =~ tr/acgtACGT/tgcaTGCA/;
	}

	unless (($c1Seq =~ /$primer1$/) && ($c2Seq =~ /^$primer2/)) {
	    die "Contig sequence doesn't match primer ($contig1:$primer1 || $contig2:$primer2)\n";
	}

	my $closure = "";
	
	#  Attempt to cross the gap using splinting reads

	my $spanClosure = span($primer1,$primer2,\@fillerSeq);
	my $spanCheck = 0;

	if ($spanClosure) {
	    my $spanCheck = checkClosure($spanClosure,$primer1,$primer2,$gapSizeEstimate,$gapUncertainty);
	    $reportNote .= "spanClosure=$spanClosure;spanCheck=$spanCheck;";
	    unless ($spanCheck) {
		$closure = $spanClosure;
	    }
	}

	#  If splints fail try a mer-walk
	my $bridgeClosureRight = 0;
	my $rightWalk = "";
	my $rightFail = "";
	my $rightCheck = 0;

	my $bridgeClosureLeft = 0;
	my $leftWalk = "";
	my $leftFail = "";
	my $leftCheck = 0;

	my $patchCheck = 0;
	my $patchClosure = 0;

	my $nibbleClosure = 0;

	unless ($closure) {
	    my $maxLeftWalk = "";
	    my $maxRightWalk = bridgeIter($c1Seq,$c2Seq,\@fillerSeq,"right");
	    print STDERR "MAXRIGHT [$maxRightWalk]\n";
	    ($rightWalk,$rightFail) = $maxRightWalk =~ /^([ACGT]*)(.*)$/;

	    unless ($rightFail) {   # Rightward walk succeeded
		$bridgeClosureRight = $rightWalk;
		$rightCheck = checkClosure($rightWalk,$primer1,$primer2,$gapSizeEstimate,$gapUncertainty);
		$reportNote .= "rightClosure=$rightWalk;rightCheck=$rightCheck;";
		unless ($rightCheck) {
		    $closure = $bridgeClosureRight;
		}
	    }

	    unless ($closure) {
		$maxLeftWalk = bridgeIter($c1Seq,$c2Seq,\@fillerSeq,"left");
		print STDERR "MAXLEFT [$maxLeftWalk]\n";
		($leftWalk,$leftFail) = $maxLeftWalk =~ /^([ACGT]*)(.*)$/;
		$leftWalk = reverse($leftWalk);
		$leftWalk =~ tr/ACGT/TGCA/;
		$leftFail =~ tr/ACGT/TGCA/;

		unless ($leftFail) {   # Leftward walk succeeded
		    $bridgeClosureLeft = $leftWalk;
		    $leftCheck = checkClosure($leftWalk,$primer1,$primer2,$gapSizeEstimate,$gapUncertainty);
		    $reportNote .= "leftClosure=$leftWalk;leftCheck=$leftCheck;";
		    unless ($leftCheck) {
			$closure = $bridgeClosureLeft;
		    }
		}
	    }
	}

	# If walks failed to close, try to patch between left and right walk
	unless ($closure || $bridgeClosureRight || $bridgeClosureLeft) {
	    $patchClosure = patch($rightWalk,$leftWalk,$gapSizeEstimate,$gapUncertainty);
	    if ($patchClosure) {
		$patchCheck = checkClosure($patchClosure,$primer1,$primer2,$gapSizeEstimate,$gapUncertainty);
		$reportNote .= "patchClosure=$patchClosure;patchCheck=$patchCheck;";
		unless ($patchCheck) {
		    $closure = $patchClosure;
		}
	    }
	}
	
	if ($nibbleMode && (!$closure)) {

	    $nibbleClosure = nibblePatch($rightWalk,$leftWalk,$gapSizeEstimate);
	    if ($nibbleClosure) {
		my $nibbleCheck = checkClosure($nibbleClosure,$primer1,$primer2,$gapSizeEstimate,$gapUncertainty);
		$reportNote .= "nibbleClosure=$nibbleClosure;nibbleCheck=$nibbleCheck;";
		unless ($nibbleCheck) {
		    $closure = $nibbleClosure;
		}
	    }

	}

	if ($closure) {

	    my $closedGapSize = length($closure) - 2*$merSize;
	    if ($nibbleClosure) { 
		$nFailure++;
		print STDERR "$gapInfo failed to close ($spanCheck;$rightCheck:$rightFail;$leftCheck:$leftFail) (gap nibbled)\n";
		$reportStatus = "FAILED\t$reportNote\n";
	    } else {
		$success = 1;
		$nSuccess++;
		print STDERR "$gapInfo successfully closed : $closedGapSize $closure\n";
		$reportStatus = "SUCCESS\t$reportNote\n";
	    }

#  Lower case as little of the primer sequences as possible:

	    my $minGapMask = 2*5;
	    if ($closure =~ /$primer1.{$minGapMask,}$primer2/) {
		my ($p1,$g,$p2) = $closure =~ /($primer1)(.+)($primer2)/;
		$g = lc($g);
		$closure = $p1.$g.$p2;
	    } else {
		my $l = length($closure);
		if ($l % 2 == 1) {
		    $minGapMask++;
		}
		my $d = ($l-$minGapMask)/2;
		my ($p1,$g,$p2) = $closure =~ /^(.{$d})(.{$minGapMask})(.{$d})/;
		$g = lc($g);
		$closure = $p1.$g.$p2;
	    }
	    
	    print "$scaffold\t$contig1\t$primer1\t$contig2\t$primer2\t$closure\n";
	    
	} else {
	    $nFailure++;
	    print STDERR "$gapInfo failed to close ($spanCheck;$rightCheck:$rightFail;$leftCheck:$leftFail)\n";
	    $reportStatus = "FAILED\t$reportNote\n";
	}

	unless($success) {
	    print STDERR "OPEN:$line\n";
	}

	print STDERR "$reportLine\t$reportStatus";


    }

    close G;

    print STDERR "Done.  Successfully closed $nSuccess gaps. ($nFailure failed to close)\n";

    exit;
}

# If gapDataFile hasn't been provided, generate one from bma information

for (my $mapIndex = 0; $mapIndex < $nBlastMaps; $mapIndex++) {
    my $libID = $libs[$mapIndex];
    my $libSize = $mapInfo{$libID}->[0];
    my $libStdDev = $mapInfo{$libID}->[1];
    my $blastMapFile = $bmaFiles[$mapIndex];
    
    print STDERR "Reading $blastMapFile (lib $libID $libSize +/- $libStdDev)\n";
    
    open (B,$blastMapFile) || die "Couldn't open $blastMapFile.\n";
    while (my  $line = <B>) {
	chomp $line;
	my @cols = split(/\t/,$line);
	
	my @contigEnds = (0,0,0,0);
	my @aligns = (0,0,0,0);
	my @gapChoices = (0,0,0,0);

	my $read1gap1 = 0;
	my $read1gap2 = 0;
	my $read2gap1 = 0;
	my $read2gap2 = 0;

	my $pairStatus = $cols[0];
	my $linkType1 = 0;
	my $linkType2 = 0;
	if ($pairStatus eq "SINGLE") {
	    $linkType1 = $cols[1];
	} else {
	    ($linkType1,$linkType2) = $cols[1] =~ /^(.+)\.(.+)$/;
	}

#	print "{$pairStatus}{$linkType1}{$linkType2}\n";

	$contigEnds[0] = $cols[2];
	$aligns[0] = $cols[3];
	if ($linkType1 eq "SPLINT") {
	    $contigEnds[1] = $cols[4];
	    $aligns[1] = $cols[5];
	}
	if ($linkType2) {   # a PAIR
	    if ($linkType1 ne "SPLINT") {
		$contigEnds[2] = $cols[4];
		$aligns[2] = $cols[5];
		if ($linkType2 eq "SPLINT") {
		    $contigEnds[3] = $cols[6];
		    $aligns[3] = $cols[7];
		}
	    } else {
		$contigEnds[2] = $cols[6];
		$aligns[2] = $cols[7];
		if ($linkType2 eq "SPLINT") {
		    $contigEnds[3] = $cols[8];
		    $aligns[3] = $cols[9];
		}
	    }
	}

	for (my $i=0;$i<4;$i++) {
	    if ($aligns[$i]) {
		$gapChoices[$i] = assignGap($aligns[$i]);
	    }
	}

	# assignGap returns "[+/-]$readID:$flankingGapID:$scaffold" or "ANCHOR:[+/-]$readID:$scaffold" 
	# chooseGap returns index (0 or 1) of the gap choice in the larger scaffold

	my $read1GapChoiceIndex = chooseGap($gapChoices[0],$gapChoices[1]);
	my $read2GapChoiceIndex = 2 + chooseGap($gapChoices[2],$gapChoices[3]);
	
	my $read1GapChoice = $gapChoices[$read1GapChoiceIndex];
	my $read2GapChoice = $gapChoices[$read2GapChoiceIndex];

	if ( (! $read1GapChoice) && (! $read2GapChoice) ) { 
	    next;
	} elsif ($read1GapChoice && (! $read2GapChoice) ) { 

	    unless ($read1GapChoice =~ /^ANCHOR/) {
		my ($str1,$r1,$gid1,$scaff1) = $read1GapChoice =~ /^([\+\-])(.+):(\d+):([^:]+)$/;
		$fillers{$r1} = $str1 . $gid1;
	    }
	    # only one read is placed, try to project the mate
	    my $align = $aligns[$read1GapChoiceIndex];
	    my $contigEnd = $contigEnds[$read1GapChoiceIndex];
	    $contigEnd = fixEnd($align,$contigEnd);

	    my $projectPairLocation = projectPair($mapIndex,$contigEnd,$align);
	    if ($projectPairLocation) {
		my ($readName,$readEnd,$scaffID,$bestGapID) = $projectPairLocation =~ /^(.+):(\d):([^:]+):(\d+)$/;
		$projectors{$readName} = "$readEnd:$contigEnd:$bestGapID";
	    }

	} elsif ( (! $read1GapChoice) && $read2GapChoice ) { 

	    unless ($read2GapChoice =~ /^ANCHOR/) {
		my ($str2,$r2,$gid2,$scaff2) = $read2GapChoice =~ /^([\+\-])(.+):(\d+):([^:]+)$/;
		$fillers{$r2} = $str2 . $gid2;
	    }
	    # only one read is placed, try to project the mate
	    my $align = $aligns[$read2GapChoiceIndex];
	    my $contigEnd = $contigEnds[$read2GapChoiceIndex];
	    $contigEnd = fixEnd($align,$contigEnd);

	    my $projectPairLocation = projectPair($mapIndex,$contigEnd,$align);
	    if ($projectPairLocation) {
		my ($readName,$readEnd,$scaffID,$bestGapID) = $projectPairLocation =~ /^(.+):(\d):([^:]+):(\d+)$/;
		$projectors{$readName} = "$readEnd:$contigEnd:$bestGapID";
	    }

	} else {  # Both reads have a value returned from assignGap

	    my $read1 = undef;
	    my $strand1 = undef;
	    my $scaff1 = undef;
	    my $gap1 = undef;

	    if ($read1GapChoice =~ /^ANCHOR/) {
		($strand1,$read1,$scaff1) = $read1GapChoice =~ /^ANCHOR:([\+\-])(.+):([^:]+)$/;
	    } else {
		($strand1,$read1,$gap1,$scaff1) = $read1GapChoice =~ /^([\+\-])(.+):(\d+):([^:]+)$/;
	    }

	    my $read2 = undef;
	    my $strand2 = undef;
	    my $scaff2 = undef;
	    my $gap2 = undef;

	    if ($read2GapChoice =~ /^ANCHOR/) {
		($strand2,$read2,$scaff2) = $read2GapChoice =~ /^ANCHOR:([\+\-])(.+):([^:]+)$/;
	    } else {
		($strand2,$read2,$gap2,$scaff2) = $read2GapChoice =~ /^([\+\-])(.+):(\d+):([^:]+)$/;
	    }

	    if ($scaff1 eq $scaff2) {
		unless ($read1GapChoice =~ /^ANCHOR/) {
		    $fillers{$read1} = $strand1 . $gap1;
		}
		unless ($read2GapChoice =~ /^ANCHOR/) {
		    $fillers{$read2} = $strand2 . $gap2;
		}

	    } else {

		my $scaffLen1 = $scaffoldLength{$scaff1}; 
		my $align1 = $aligns[$read1GapChoiceIndex];
		my $contigEnd1 = $contigEnds[$read1GapChoiceIndex];
		$contigEnd1 = fixEnd($align1,$contigEnd1);
		my $projectPairLocation1 = projectPair($mapIndex,$contigEnd1,$align1);

		my $scaffLen2 = $scaffoldLength{$scaff2}; 
		my $align2 = $aligns[$read2GapChoiceIndex];
		my $contigEnd2 = $contigEnds[$read2GapChoiceIndex];
		$contigEnd2 = fixEnd($align2,$contigEnd2);
		my $projectPairLocation2 = projectPair($mapIndex,$contigEnd2,$align2);

		if ($projectPairLocation1 && ($scaffLen1 > $scaffLen2)) {
		    my ($readName,$readEnd,$scaffID,$bestGapID) = $projectPairLocation1 =~ /^(.+):(\d):([^:]+):(\d+)$/;
		    $projectors{$readName} = "$readEnd:$contigEnd1:$bestGapID";		    
		    unless ($read1GapChoice =~ /^ANCHOR/) {
			$fillers{$read1} = $strand1 . $gap1;
		    }
		} elsif ($projectPairLocation2 && ($scaffLen2 > $scaffLen1)) {
		    my ($readName,$readEnd,$scaffID,$bestGapID) = $projectPairLocation2 =~ /^(.+):(\d):([^:]+):(\d+)$/;
		    $projectors{$readName} = "$readEnd:$contigEnd2:$bestGapID";		    
		    unless ($read2GapChoice =~ /^ANCHOR/) {
			$fillers{$read2} = $strand2 . $gap2;
		    }
		} else {

		    unless ($read1GapChoice =~ /^ANCHOR/) {
			$fillers{$read1} = $strand1 . $gap1;
		    }
		    unless ($read2GapChoice =~ /^ANCHOR/) {
			$fillers{$read2} = $strand2 . $gap2;
		    }
		}
	    }
	}
    }
    close B; 
}

my $nFillers = keys(%fillers);
my $nProjectors = keys(%projectors);

print STDERR "Done. Found $nFillers (potential) gap fillers; $nProjectors (potential) gap projectors.\n";

my %fillSeq = ();
my @allReadSequence = ();
my @allReadNames = ();
my $readIndex = 0;

$nFillers = 0;

my @seqFiles = glob($seqFileGlob);

foreach my $seqFile (@seqFiles) {

    print STDERR "Reading $seqFile...\n";

    if ($gzipped) {
        open (S,"gunzip -c $seqFile |") || die "Couldn't open pipe to $seqFile\n";
    } else {
        open (S,$seqFile) || die "Couldn't open $seqFile\n";
    }

    while (my $line = <S>) {
	chomp $line;

	my $readName = undef;
	my $nts = undef;
	my $quals = undef;

	if ($fastQ) {
            chomp $line;
            unless ($line =~ /^\@/) {
                die "Invalid fastQ $line\n";
            }
	    my $name = undef;
	    my $end = undef;
           # Latest Illumina header format
            if ($line =~ /^\@\S+\s+[12]:[YN]:\d+:/) {
                ($name,$end) = $line =~ /^\@(\S+)\s+([12]):[YN]:/;
            } else {
                # older format
		($name,$end) = $line =~ /^\@(\S+)\/(\d)/;
            }

	    unless (defined($name) && defined($end)) {
		die "Unexpected read name format in $line\n";
	    }
            $readName = $name ."/$end";
            $line = <S>;
            chomp $line;
            $nts = $line;
            $line = <S>;
            unless ($line =~ /^\+/) {
                die "Invalid fastQ $line\n";
            }
            $line = <S>;
            chomp $line;
            $quals = $line;
	} else {
	    my @cols = split(/\:/,$line);
	    $readName = "$cols[0]:$cols[1]:$cols[2]:$cols[3]:$cols[4]";
	    $nts = $cols[5];
	    $quals = $cols[6];
	}

	my ($pairName,$pairEnd) = $readName =~ /(.+)\/(\d)/;

	if (exists($projectors{$pairName})) {
	    my $pInfo = $projectors{$pairName};
	    my ($placedEnd,$contig,$projectedGapID) = $pInfo =~ /^(\d):(.+):(\d+)$/;
	    if ($pairEnd ne $placedEnd) {
	    
		$nFillers++;	    
		unless (exists($gapInfo{$contig})) {
		    print STDERR "Error: projected read ($pairName,$pairEnd) associated with ungapped end: $contig\n";
		    die;
		}
	    
		my ($flankingGapID,$scaffold,$contig1,$primer1,$gapSize,$contig2,$primer2,$gapUncertainty) = @{$gapInfo{$contig}};
		if ($contig eq $contig1) {
		    $nts = reverse($nts);
		    $nts =~ tr/ACGTacgt/TGCAtgca/;
		    $quals = reverse($quals);
		}
	    
		my $gapID = $projectedGapID;

		if (exists($fillSeq{$gapID})) {
		    $fillSeq{$gapID} .= ":$readIndex";
		} else {
		    $fillSeq{$gapID} = "$readIndex";
		}

		push(@allReadSequence,"$nts:$quals");
		push(@allReadNames,$readName);
		$readIndex++;

	    }

	}

	if (exists($fillers{$readName})) {

	    $nFillers++;
	    my ($strand,$gapID) = $fillers{$readName} =~ /^([\+\-])(\d+)$/;
	
	    if ( $strand eq "-" ) {
		$nts = reverse($nts);
		$nts =~ tr/ACGTacgt/TGCAtgca/;
		$quals = reverse($quals);
	    }
	
	    if (exists($fillSeq{$gapID})) {
		$fillSeq{$gapID} .= ":$readIndex";
	    } else {
		$fillSeq{$gapID} = "$readIndex";
	    }

	    push(@allReadSequence,"$nts:$quals");
	    push(@allReadNames,$readName);
	    $readIndex++;

	}
    }
    close S;
}
print STDERR "Done. Confirmed $nFillers (potential) gap-filling reads.\n";

for (my $g=0;$g < $nGaps; $g++) {
    my $leadContig = $gapOrder[$g];
    my ($gapID,$scaffold,$contig1,$primer1,$gapSizeEstimate,$contig2,$primer2,$gapUncertainty) = @{$gapInfo{$leadContig}};

    my ($c1ID) = $contig1 =~ /^(.+)\.[35]$/;
    unless (exists($contigSequence{$c1ID})) {
	die "[$c1ID] no seq\n";
    }
    unless (exists($contigScaffoldMap{$c1ID})) {
	die "[$c1ID] no scaffMapInfo\n";
    }
    my $c1Seq = $contigSequence{$c1ID};
    my $c1Info = $contigScaffoldMap{$c1ID};
    my ($c1Strand) = $c1Info =~ /^([\+\-])/;
    if ($c1Strand eq "-") {
	$c1Seq = reverse($c1Seq);
	$c1Seq =~ tr/acgtACGT/tgcaTGCA/;
    }

    my ($c2ID) = $contig2 =~ /^(.+)\.[35]$/; 
    unless (exists($contigSequence{$c2ID})) {
	die "[$c2ID] no seq\n";
    }
    unless (exists($contigScaffoldMap{$c2ID})) {
	die "[$c2ID] no scaffMapInfo\n";
    }
    my $c2Seq = $contigSequence{$c2ID};
    my $c2Info = $contigScaffoldMap{$c2ID};
    my ($c2Strand) = $c2Info =~ /^([\+\-])/;
    if ($c2Strand eq "-") {
	$c2Seq = reverse($c2Seq);
	$c2Seq =~ tr/acgtACGT/tgcaTGCA/;
    }

    unless (($c1Seq =~ /$primer1$/) && ($c2Seq =~ /^$primer2/)) {
	die "Contig sequence doesn't match primer ($contig1:$primer1 || $contig2:$primer2)\n";
    }

    if ($fillSeq{$gapID}) {

	my $fillerIndex = $fillSeq{$gapID};
	my @fillerIndex = split(/\:/,$fillerIndex);
	my $fillerSeq = join("\t",@allReadSequence[@fillerIndex]);

	print "$scaffold\t$contig1\t$primer1\t$contig2\t$primer2\t$gapSizeEstimate\t$gapUncertainty\t$fillerSeq\n";
    }
}

#
# SUBROUTINES
#

sub printFasta {
    my ($name,$seq) = @_;
    my $seqLen = length($seq);

    my $bpl = 50;
    my $nLines = sprintf("%d",$seqLen/$bpl);
    if ($seqLen%$bpl != 0) {
        $nLines++;
    }

    print ">$name\n";
    for (my $j = 0;$j<$nLines;$j++) {
        my $seqLine = substr($seq,$j*$bpl,$bpl);
        my $text = $seqLine;
        $text.= "\n";
        print $text;
    }
}

sub span {
    my ($primer1,$primer2,$seqRef) = @_;
    
    my @reads = @{$seqRef};
    my $nReads = scalar(@reads);

    if ($verboseMode) {
	print STDERR "Attempting span: $primer1 -> $primer2 ($nReads read(s) available)\n";
    }

    my %pureSpans = ();
    foreach my $read (@reads) {
	my ($nts,$quals) = $read =~ /^([^:]+)\:(.+)$/;
	if (($nts =~ /$primer1/) && ($nts =~ /$primer2/)) {
	    my ($tmpSpan) = $nts =~ /($primer1.*)$/;
	    my ($span) = $tmpSpan =~ /^(.*$primer2)/; 
	
	    if ($span) {
		$pureSpans{$span}++;
	    }
	}
    }
    my @spans = keys(%pureSpans);
    my $nSpans = scalar(@spans);

    if ($verboseMode) {
	print STDERR "$nSpans distinct spanning sequence(s) found\n";
    }

    if (($nSpans == 1) && ($pureSpans{$spans[0]} > 1)) {
	print STDERR "Unique spanning sequence found: $pureSpans{$spans[0]}/$nReads reads span the gap.\n";
	return $spans[0];

    } elsif ($polymorphicMode) {

	my $maxSpanIndex = 0;
	my $maxSpanFreq = 0;
	for (my $i=0; $i<$nSpans; $i++) {
	    my $spanFreq = $pureSpans{$spans[$i]};
	    if ($spanFreq > $maxSpanFreq) {
		$maxSpanIndex = $i;
		$maxSpanFreq = $spanFreq;
	    }
	}

	if ($maxSpanFreq > 1) {
	    print STDERR "Maximum frequency spanning sequence found: $pureSpans{$spans[$maxSpanIndex]}/$nReads reads span the gap.\n";
	    return $spans[$maxSpanIndex];
	} else {
	    if ($verboseMode) {
		print STDERR "No spanning sequence with frequency > 1 found\n";
	    } 
	    return 0;
	}

    } else {
	if ($verboseMode) {
	    print STDERR "No unique spanning sequence with frequency > 1 found\n";
	} 
	return 0;
    }
}


# Iterative bridging allows variable k-mer size

sub bridgeIter {
    my ($c1seq,$c2seq,$seqRef,$direction) = @_;
    
    my $minQ = 19;
    my $minMerLen = 13;
    my $maxReadLength = 0;

    my @reads = @{$seqRef};
    my $nReads = scalar(@reads);

    if ($verboseMode) {
	print STDERR "Attempting $direction bridge: ($nReads read(s) available)\n";
    }

    my $merLen = $merSize;

    my $iterating = 1;
    my $downshift = 0;
    my $upshift = 0;

    if ($direction eq "left") {
	$c1seq = reverse($c1seq);
	$c1seq =~ tr/ACGT/TGCA/;
	$c2seq = reverse($c2seq);
	$c2seq =~ tr/ACGT/TGCA/;
	($c1seq,$c2seq) = ($c2seq,$c1seq);

	for (my $r=0; $r < $nReads; $r++) {
	    my $read = $reads[$r];
	    my ($nts,$quals) = $read =~ /^([^:]+)\:(.+)$/;
	    $quals = reverse($quals);
	    $nts = reverse($nts);
	    $nts =~ tr/ACGT/TGCA/;
	    $reads[$r] = "$nts:$quals";
	}

    }

    my $maxWalk = undef;
    my $maxWalkLen = 0;
    my ($truePrimer1) = $c1seq =~ /(.{$merSize})$/;
    my ($truePrimer2) = $c2seq =~ /^(.{$merSize})/;

    while ($iterating) {

	my ($primer1) = $c1seq =~ /(.{$merLen})$/;
	my ($primer2) = $c2seq =~ /^(.{$merLen})/;

	# Allow k-mers less than merSize to be used 

	if ($merLen < $merSize) {
	    my $theRest = $merSize-$merLen;
	    ($primer1) = $c1seq =~ /(.{$merLen}).{$theRest}$/;
	    ($primer2) = $c2seq =~ /^.{$theRest}(.{$merLen})/;
	}

	unless ($primer1 && $primer2) {
	    if ($verboseMode) {
		print STDERR "Unable to find k=$merLen seeds\n";
	    }
	    if (defined($maxWalk)) {
		return $maxWalk;
	    } else {
		return "X";
	    }
	}

	my %mers = ();
	my %fullInfo = ();

	foreach my $read (@reads) {

	    my ($nts,$quals) = $read =~ /^([^:]+)\:(.+)$/;
	    my $readLength = length($nts);
	    if ($readLength > $maxReadLength) {
		$maxReadLength = $readLength;
	    }
	    if ($merLen >= $readLength) {
		next;
	    }

	    my $merPlus = $merLen+1;
	    while ($nts =~ /([ACGT]{$merPlus})/g) {
		my $merExt = $1;

		my ($mer,$extension) = $merExt =~ /^([^:]+)(.)$/;
		my $q = substr($quals,pos($nts)-1,1);
		my $Q = ord($q)-$qualOffset;

		my $qBin = sprintf("%d",$Q/10);
		if ($qBin > 3) {
		    $qBin = 3;
		}
		$qBin = $extension.$qBin;
		unless (exists($fullInfo{$mer})) {
		    $fullInfo{$mer} = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
		}
		$fullInfo{$mer}->[$qBins{$qBin}]++;

		pos($nts) -= $merLen;
	    }
	}

# Full analysis of quality/frequency profile

	while (my ($mer,$extMatrix) = each(%fullInfo)) {
	    my @extensions = @{$extMatrix};

	    my %ext = ();
	    $ext{"A"} = [@extensions[0..3]];
	    $ext{"C"} = [@extensions[4..7]];
	    $ext{"G"} = [@extensions[8..11]];
	    $ext{"T"} = [@extensions[12..15]];
	
	    my %nExt = ();
	    my %nHiQ = ();
	    my %ratings = ();
	    my $nTotal = 0;
	    my @bases = ("A","C","G","T");
	    
	    foreach my $base (@bases) {
#    Ignore q<10 bases
#		$nExt{$base} = $ext{$base}->[0]+$ext{$base}->[1]+$ext{$base}->[2]+$ext{$base}->[3];

		$nExt{$base} = $ext{$base}->[1]+$ext{$base}->[2]+$ext{$base}->[3];
		$nHiQ{$base} = $ext{$base}->[2]+$ext{$base}->[3];
		$nTotal += $nExt{$base};
		$ratings{$base} = categorizeExtension(@{$ext{$base}});		
	    }

	    my @sortedBases = sort {($ratings{$b} <=> $ratings{$a}) || 
					($nHiQ{$b} <=> $nHiQ{$a}) || 
					($nExt{$b} <=> $nExt{$a})} (@bases);

	    my @sortedRatings = @ratings{@sortedBases};

# Rules for choosing next base
# ratings:
# 0 = No votes
# 1 = One vote
# 2 = nVotes < minViable
# 3 = minDepth > nVotes >= minViable, nHiQ < minViable
# 4 = minDepth > nVotes >= minViable ; nHiQ >= minViable
# 5 = nVotes >= minDepth ; nHiQ < minViable
# 6 = nVotes >= minDepth ; minViable <= nHiQ < minDepth 
# 7 = nHiQ >= minDepth 

	    my $topRating = $sortedRatings[0];
	    my $runnerUp = $sortedRatings[1];
	    my $topRatedBase = $sortedBases[0];

	    if ($topRating < 3) {         # must have at least minViable bases
		$mers{$mer} = "X";
	    } elsif ($topRating == 3) {    # must be uncontested   
		if ($runnerUp == 0) {
		    $mers{$mer} = $topRatedBase;
		} else {
		    $mers{$mer} = "X";
		}
	    } elsif ($topRating < 6) {
		if ($runnerUp < 3) {
		    $mers{$mer} = $topRatedBase;
		} else {
		    $mers{$mer} = "X";
		}   
	    } elsif ($topRating == 6) {  # viable and fair hiQ support
		if ($runnerUp < 4) {
		    $mers{$mer} = $topRatedBase;
		} else {
		    $mers{$mer} = "X";
		}
	    } else {                     # strongest rating trumps
		if ($runnerUp < 7) {       
		    $mers{$mer} = $topRatedBase;
		} else {
		    my $fork = "F";
		    for (my $b = 0; $b < 4; $b++) {
			my $base = $sortedBases[$b];
			my $rating = $sortedRatings[$b];
			my $n = $nExt{$base};
			if ($rating == 7) {
			    $fork .= "$base$n";
			} else {
			    last;
			}
		    }
		    $mers{$mer} = $fork;
		}
	    }
	}

	my $walk = $primer1;
	unless (defined($maxWalk)) {
	    $maxWalk = $walk;
	    $maxWalkLen = length($maxWalk);
	}
	my $step = $primer1;
	my %loopCheck = ();
	my $success = 0;
	my $nForks = 0;
	my $nSteps = 0;
	while () {
	    if (exists($loopCheck{$step})) {
		$walk .= "R";
		last;
	    } else {
		$loopCheck{$step} = 1;
	    }

	    if (exists($mers{$step})) {
	    
		my $next = $mers{$step};

		$nSteps++;
		if ($verboseMode) {
		    print STDERR "$nSteps : $step->$next\t";
		    print STDERR "[@{$fullInfo{$step}}[0..3]][@{$fullInfo{$step}}[4..7]][@{$fullInfo{$step}}[8..11]][@{$fullInfo{$step}}[12..15]]\n";
		}

		if ($next =~ /^([ACGT])$/) {
		    $step = substr($step,1);
		    $step .= $next;
		    $walk .= $next;

		    if ($walk =~ /$truePrimer2$/) {
			$success = 1;
			last;
		    }


		} elsif (($polymorphicMode == 1) && ($nForks == 0) && ($next =~ /^F/)) {

		    #  Biallelic positions only (for now) maximum vote path is taken
		    if ($next =~ /^F[ACGT]\d+[ACGT]\d+$/) {
			my ($o1,$v1,$o2,$v2) = $next =~ /^F([ACGT])(\d+)([ACGT])(\d+)$/;
			$next = ($v1 > $v2) ? $o1 : $o2;
			if ($verboseMode) {
			    print STDERR "Polymorphic conditions met .. attempting max frequency resolution.\n";
			}
		    } else {
			$walk .= $next;
			last;
		    }

		    $step = substr($step,1);
		    $step .= $next;
		    $walk .= $next;
		    $nForks++;

		    if ($walk =~ /$truePrimer2$/) {
			$success = 1;
			last;
		    }

#		    if ($step eq $primer2) {
#			$success = 1;
#			last;
#		    }

		} else {
		    $walk .= $next;
		    last;
		}
	    } else {
		$walk .= "X";
		last;
	    }
	}

	# Trim off extra lead bases if upshifted
	my $additionalBases = $merLen-$merSize;
	if ($additionalBases > 0) {
	    my $reducedWalk = $walk;
	    ($walk) = $reducedWalk =~ /^.{$additionalBases}(.+)$/;
	}

	my ($walkBases) = $walk =~ /^([ACGT]*)/;
	my $walkLen = length($walkBases);

	if ($success || ($walkLen > $maxWalkLen)) {
	    $maxWalk = $walk;
	    $maxWalkLen = $walkLen;
	}
	
	if (($walk =~ /F/) || ($walk =~ /R/)) {
	    $merLen += 2;
	    $upshift = 1;

	    if (($downshift==1) || ($merLen >= $maxReadLength)) {
		return($maxWalk);
	    }
	    if ($verboseMode) {
		print STDERR "Degeneracy encountered; upshifting (k->$merLen)\n";
	    }
	} elsif ($walk =~ /X/) {
	    $merLen -= 2;
	    $downshift = 1;

	    if (($upshift==1) || ($merLen < $minMerLen)) {
		return($maxWalk);
	    }
	    if ($verboseMode) {
		print STDERR "Termination encountered; downshifting (k->$merLen)\n";
	    }
	} else {
	    return($maxWalk);
	}
    }
}

sub categorizeExtension {
    my (@ext) = @_;

    my $minViable = 3;
    if ($minViable > $minDepth) {
	warn "Warning: in categorizeExtension minViable reset to match minDepth ($minDepth)\n";
	$minViable = $minDepth;
    }

# 0 = No votes
# 1 = One vote
# 2 = nVotes < minViable
# 3 = minDepth > nVotes >= minViable, nHiQ < minViable
# 4 = minDepth > nVotes >= minViable ; nHiQ >= minViable
# 5 = nVotes >= minDepth ; nHiQ < minViable
# 6 = nVotes >= minDepth ; minViable < nHiQ < minDepth 
# 7 = nHiQ >= minDepth 

#    Ignore q<10 bases
#    my $n = $ext[0]+$ext[1]+$ext[2]+$ext[3];
    my $n = $ext[1]+$ext[2]+$ext[3];
    my $nHiQ = $ext[2]+$ext[3];
    my $nLoQ = $n-$nHiQ;

    my $category = undef;
    if ($n == 0) {
	$category = 0;
    } elsif ($n == 1) {
	$category = 1;
    } elsif ($n < $minViable) {
	$category = 2;
    } else {
	if (($n < $minDepth) || ($n == $minViable)) {
	    if ($nHiQ < $minViable) {
		$category = 3;
	    } else {
		$category = 4;
	    }
	} else {
	    if ($nHiQ < $minViable) {
		$category = 5;
	    } elsif ($nHiQ < $minDepth) {
		$category = 6;
	    } else {
		$category = 7;
	    }
	}
    }

    unless (defined($category)) {
	die "Undefined extension category for [@ext]\n";
    }
    return $category;
};

sub checkClosure {
    my ($closureToCheck,$primer1,$primer2,$gapSizeEstimate,$gapUncertainty) = @_;

    my $badClosure = 0;

    if (! (($closureToCheck =~ /^$primer1/) && ($closureToCheck =~ /$primer2$/))) {
	$badClosure += 1;
	if ($verboseMode) {
	    print STDERR "closure [$closureToCheck] rejected because it disagrees with primers\n";
	}
    }
    
    my $closedGapSize = length($closureToCheck) - 2*$merSize;
    my $gapDiff = $closedGapSize - $gapSizeEstimate;
    
    if (abs($gapDiff) > $gapUncertainty) {
	$badClosure += 2;
	if ($verboseMode) {
	    print STDERR "closure [$closureToCheck] rejected due to gap estimate differential: |$gapDiff| > $gapUncertainty\n";
	}
    }

    # In aggressive mode ignore size differences between closure and estimate
    if ($aggression && ($badClosure == 2) && ($gapSizeEstimate < 2*$maxInsertSize))
    {
	$badClosure = 0;
	if ($verboseMode) {
	    print STDERR "closure [$closureToCheck] allowed via AGRESSIVE mode\n";
	}
    }

    return $badClosure;

}

sub nibblePatch {
    my ($rightWalk,$leftWalk,$gapSizeEstimate) = @_;

    my $rightLen = length($rightWalk);
    my $leftLen = length($leftWalk);

    my $minGapSize = 10;
    my $gapPad = "N"x$minGapSize;

    my $idealLength = 2*$merSize+$gapSizeEstimate;

    if ($rightLen + $leftLen + $minGapSize <=  $idealLength) {

	my $extraPadLength = $idealLength - ($rightLen+$leftLen+$minGapSize);
	if ($extraPadLength > 0) {
	    $gapPad .= "N"x$extraPadLength;
	}

	return "$rightWalk$gapPad$leftWalk";

    } else {

	my $extraLength = ($rightLen+$leftLen+$minGapSize)-$idealLength;
	if ($rightLen > $leftLen) {
	    my $lengthDiff = $rightLen-$leftLen;
	    if ($lengthDiff >= $extraLength) {
		my ($trimmedRightWalk) = $rightWalk =~ /^(.+).{$extraLength}$/;
		$rightWalk = $trimmedRightWalk;

		return "$rightWalk$gapPad$leftWalk";

	    } else {
		my ($trimmedRightWalk) = $rightWalk =~ /^(.+).{$lengthDiff}$/;
		$rightWalk = $trimmedRightWalk;
	    }
	} elsif ($leftLen > $rightLen) {
	    my $lengthDiff = $leftLen-$rightLen;
	    if ($lengthDiff >= $extraLength) {
		my ($trimmedLeftWalk) = $leftWalk =~ /^.{$extraLength}(.+)$/;
		$leftWalk = $trimmedLeftWalk;

		return "$rightWalk$gapPad$leftWalk";

	    } else {
		my ($trimmedLeftWalk) = $leftWalk =~ /^.{$lengthDiff}(.+)$/;
		$leftWalk = $trimmedLeftWalk;
	    }
	}

	$rightLen = length($rightWalk);
	$leftLen = length($leftWalk);

	my $rightTurn = 1;
	while ($rightLen + $leftLen + $minGapSize >  $idealLength) {
	    if ($rightTurn) {
		chop $rightWalk;
		$rightLen--;
		$rightTurn = 0;
	    } else {
		$leftWalk = substr($leftWalk,1);
		$leftLen--;
		$rightTurn = 1;
	    }
	}

	return "$rightWalk$gapPad$leftWalk";
    }
}




sub patch {
    my ($rightWalk,$leftWalk,$gapSizeEstimate,$gapUncertainty) = @_;

    my $minAcceptableOverlap = 10;

    my $rightLen = length($rightWalk);
    my $leftLen = length($leftWalk);

    my $minOffset = sprintf("%.0f",-$gapUncertainty);
    my $maxOffset = sprintf("%.0f",$gapUncertainty);

    my $r0 = 1;
    my $r1 = $rightLen;
    my $l1 = 2*$merSize+$gapSizeEstimate;
    my $l0 = $l1-$leftLen+1;

    if ($verboseMode) {
	print STDERR "Attempting patch [$rightLen][$leftLen] (gap: $gapSizeEstimate +/- $gapUncertainty)\n";
    }
	
    my %overlaps = ();
    for (my $o = $minOffset; $o <= $maxOffset; $o++) {

	my $testRight = $rightWalk;
	my $testLeft = $leftWalk;

	my $l0shifted = $l0 + $o;
	my $l1shifted = $l1 + $o;
	
	my $chopLeft = $r0 - $l0;
	my $chopRight = $r1 - $l1;

	if ($chopLeft > 0) {
	    my ($chop,$keep) = $testLeft =~ /^(.{$chopLeft})(.*)$/;
	    $testLeft = $keep;
	}
	my $tlLen = length($testLeft);

	if ($chopRight > 0) {
	    my ($keep,$chop) = $testRight =~ /^(.*)(.{$chopRight})$/;
	    $testRight = $keep;
	}
	my $trLen = length($testRight);

	my $overlap = $trLen+$tlLen-(2*$merSize+$gapSizeEstimate+$o);
	unless (($overlap >= $minAcceptableOverlap) && ($overlap <= $trLen) && ($overlap <= $tlLen)) {
	    next;
	}

	my ($p1Prefix,$p1Suffix) = $testRight =~ /^(.*)(.{$overlap})$/;
	my ($p2Prefix,$p2Suffix) = $testLeft =~ /^(.{$overlap})(.*)$/;
	if ($p1Suffix eq $p2Prefix) {
	    $overlaps{$o} = $testRight.$p2Suffix;
	}
    }

    my $bestGuessSeq = "";
    my $bestGuessDelta = $gapUncertainty+1;
    my $nBestGuesses = 0;
    my $nGuesses = 0;
    while (my ($offset,$guessSeq) = each(%overlaps)) {
	my $delta = abs($offset);
	$nGuesses++;
	if ($delta < $bestGuessDelta) {
	    $bestGuessSeq = $guessSeq;
	    $bestGuessDelta = $delta;
	    $nBestGuesses = 1;
	} elsif ($delta == $bestGuessDelta) {
	    if (length($guessSeq) < length($bestGuessSeq)) {
		$bestGuessSeq = $guessSeq;
		$bestGuessDelta = $delta;
	    }
	    $nBestGuesses++;
	}
    }

    if ($nGuesses) {
	print STDERR "$nGuesses potential patches identified.  $nBestGuesses best guess(es) differ from gap estimate by $bestGuessDelta\n";
	return $bestGuessSeq;
    } else {
	print STDERR "No valid patches found.\n";
	return 0;
    }
}

# Given a bma alignment, assign the read to a gapID with orientation

sub assignGap {
    my ($alignment) = @_;

    unless ($alignment) {
	return 0;
    }

    my ($align) = $alignment =~ /^\[(.+)\]$/;

    my ($rStat,$r,$r0,$r1,$rL,$c,$c0,$c1,$cL,$strand) = split(/\s+/,$align);
    my ($readName,$readEnd) = $r =~ /^(.+)\/(\d)$/;

    my ($r5stat,$r3stat,$inOrOut) = $rStat =~ /^([A-Z]+)\.([A-Z]+)\.([A-Z]+)$/;

    unless (($r5stat eq "GAP") || ($r3stat eq "GAP")) {
	if (($r5stat ne "TRUNC") && ($r3stat ne "TRUNC")) {

	    # Adding null check for contigs not included in the assembly 1/9/2013
	    unless (exists($contigScaffoldMap{$c})) {
		return 0;
	    }
	    my $contigScaffInfo = $contigScaffoldMap{$c};
	    my ($sStrand,$sID,$sStart,$sEnd) = split(/\,/,$contigScaffInfo);
	    my $cStrand = ($strand eq "Plus") ? "+" : "-";
	    my $rStrand = ($cStrand eq $sStrand) ? "+" : "-";
	    return "ANCHOR:$rStrand$r:$sID";  
	} else {
	    return 0;
	}
    }


    if ($r5stat eq "GAP") {    #  read 5'-end gaps will be preferentially used 
	my $contig = $c;
	if ($strand eq "Plus") {
	    $contig .= ".5";
	} else {
	    $contig .= ".3";
	}
	
	if (exists($gapInfo{$contig})) {
	    my ($flankingGapID,$scaffold,$contig1,$primer1,$gapSize,$contig2,$primer2,$gapUncertainty) = @{$gapInfo{$contig}};
	    if ($contig eq $contig1) {
		return "-$r:$flankingGapID:$scaffold";
	    } elsif ($contig eq $contig2) {
		return "+$r:$flankingGapID:$scaffold";
	    } else {
		print STDERR "Error: $contig does not match gap $flankingGapID contigs ($contig1,$contig2)\n";
		die;
	    }
	}
    }

    if ($r3stat eq "GAP") {    
	my $contig = $c;
	if ($strand eq "Plus") {
	    $contig .= ".3";
	} else {
	    $contig .= ".5";
	}
	
	if (exists($gapInfo{$contig})) {
	    my ($flankingGapID,$scaffold,$contig1,$primer1,$gapSize,$contig2,$primer2,$gapUncertainty) = @{$gapInfo{$contig}};
	    if ($contig eq $contig1) {
		return "+$r:$flankingGapID:$scaffold";
	    } elsif ($contig eq $contig2) {
		return "-$r:$flankingGapID:$scaffold";
	    } else {
		print STDERR "Error: $contig does not match gap $flankingGapID contigs ($contig1,$contig2)\n";
		die;
	    }
	}
    }

    return 0;   
}

# If two gap placement options are available for a read, choose the one in the larger scaffold
# Returns 0 for option1, 1 for option2;

sub chooseGap {
    my ($option1,$option2) = @_;

    unless ($option1 || $option2) {
	return 0;
    }

    if ($option1) {
	if ($option2) {

	    if ($option1 eq $option2) {
		return 0;
	    }

	    my ($strand1,$read1,$gap1,$scaffold1) = $option1 =~ /^([\+\-])(.+):(\d+):([^:]+)$/;
	    my ($strand2,$read2,$gap2,$scaffold2) = $option2 =~ /^([\+\-])(.+):(\d+):([^:]+)$/;

	    unless ($read1 eq $read2) {
		print STDERR "Error in chooseGap: $read1 ne $read2\n";
		die;
	    }

	    my $scaffLen1 = $scaffoldLength{$scaffold1};
	    my $scaffLen2 = $scaffoldLength{$scaffold2};

	    if ($scaffLen2 > $scaffLen1) {
		return 1;
	    } else {
		return 0;
	    }

	} else {
	    return 0;
	}
    } else {
	return 1;
    }
}

# Code for projecting pair into the best fit gap:
# assumes projected read is same length as projecting read

sub projectPair {
    my ($mapIndex,$contigEnd,$alignment) = @_;
    my ($align) = $alignment =~ /^\[(.+)\]$/;

    my ($contig,$end) = $contigEnd =~ /^(.+)\.([35])$/;

    unless (exists($contigScaffoldMap{$contig})) {
	warn "No contig->scaffold information available for $contig.\n";
	return 0;
    }

    unless (exists($gapInfo{$contigEnd})) {
#	warn "No gapInfo available for $contigEnd.\n";
	return 0;
    }

    my ($rStat,$r,$r0,$r1,$rL,$c,$c0,$c1,$cL,$strand) = split(/\s+/,$align);
    my ($readName,$readEnd) = $r =~ /^(.+)\/(\d)$/;

    my $libID = $libs[$mapIndex];
    my $libSize = $mapInfo{$libID}->[0];
    my $libStdDev = $mapInfo{$libID}->[1];

    my $projectedLoc = sprintf("%.0f",$libSize - $rL/2. - $r0 + 1);
    if ($strand eq "Plus") {
	$projectedLoc = $c0+$projectedLoc;
    } else {
	$projectedLoc = $c1-$projectedLoc;
    }

    my $contigScaffInfo = $contigScaffoldMap{$contig};
    my ($sStrand,$sID,$sStart,$sEnd) = split(/\,/,$contigScaffInfo);

    my $scaffProjectedLoc = ($sStrand eq "+") ? ($sStart + $projectedLoc - 1) : ($sEnd - $projectedLoc + 1);
    
    unless(exists($scaffGapRange{$sID})) {
	warn "No gaps associated with scaffold $sID.\n";
	return 0;
    }
    
    my $flankingGapInfo = $gapInfo{$contigEnd};
    my $flankingGapID = $flankingGapInfo->[0];
    
    my $gapRange = $scaffGapRange{$sID};
    my ($g0,$g1) = $gapRange =~ /^(\d+)\.(\d+)$/;
    
# Gap choice optimizer
    my $bestGapZ = undef;
    my $bestGap = $flankingGapID;
    for (my $g = $flankingGapID; $g < $g1; $g++) {
	my $gapBoundaries = $gapBoundaries[$g];
	my ($leftContigEnd,$rightContigStart) = $gapBoundaries =~ /^(\d+):(\d+)$/;
	my $zL = abs(($scaffProjectedLoc-$leftContigEnd)/$libStdDev);
	my $zR = abs(($scaffProjectedLoc-$rightContigStart)/$libStdDev);
	my $z = ($zL < $zR) ? $zL : $zR;

	if (!defined($bestGapZ)) {
	    $bestGapZ = $z;
	    $bestGap = $g;
	} elsif ($z < $bestGapZ) {
	    $bestGapZ = $z;
	    $bestGap = $g;
	} else {
	    last;
	}
    }
    for (my $g = $flankingGapID-1; $g >= $g0; $g--) {
	my $gapBoundaries = $gapBoundaries[$g];
	my ($leftContigEnd,$rightContigStart) = $gapBoundaries =~ /^(\d+):(\d+)$/;
	my $zL = abs(($scaffProjectedLoc-$leftContigEnd)/$libStdDev);
	my $zR = abs(($scaffProjectedLoc-$rightContigStart)/$libStdDev);
	my $z = ($zL < $zR) ? $zL : $zR;

	if (!defined($bestGapZ)) {
	    $bestGapZ = $z;
	    $bestGap = $g;
	} elsif ($z < $bestGapZ) {
	    $bestGapZ = $z;
	    $bestGap = $g;
	} else {
	    last;
	}
    }

    my ($scaffoldLength) = $scaffoldLength{$sID};

    my $zOffScaffLeftEnd = abs(($scaffProjectedLoc-0)/$libStdDev);
    my $zOffScaffRightEnd = abs(($scaffProjectedLoc-$scaffoldLength)/$libStdDev);

    if ($zOffScaffLeftEnd < $bestGapZ) {
#	warn "Projected off left scaffold end\n";
	return 0;
    } elsif ($zOffScaffRightEnd < $bestGapZ) {
#	warn "Projected off right scaffold end\n";
	return 0;
    }

    return "$readName:$readEnd:$sID:$bestGap";

}    


sub fixEnd {
    my ($alignment,$contigEnd) = @_;

    my ($align) = $alignment =~ /^\[(.+)\]$/;

    my ($rStat,$r,$r0,$r1,$rL,$c,$c0,$c1,$cL,$strand) = split(/\s+/,$align);
    my ($contig,$end) = $contigEnd =~ /^(.+)\.([35])$/;

    my ($r5stat,$r3stat,$inOrOut) = $rStat =~ /^([A-Z]+)\.([A-Z]+)\.([A-Z]+)$/;
    if (($r5stat eq "GAP") && ($inOrOut eq "OUT")) {
	$end = ($end == 3 ? 5 : 3);
	$contigEnd = "$contig.$end";
    }

    return $contigEnd;

}
