#!/usr/bin/perl -w

#
# Copyright 2011, Ben Langmead <langmea@cs.jhu.edu>
#
# This file is part of Bowtie 2.
#
# Bowtie 2 is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# Bowtie 2 is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with Bowtie 2.  If not, see <http://www.gnu.org/licenses/>.
#

##
# AlignmentCheck.pm
#
# Read in fasta files containing reference sequences that might be
# aligned to, then read in alignment files, checking each alignment to
# be sure it's sane and consistent with the reference sequence it
# aligns to.
#

package AlignmentCheck;
use strict;
use warnings;
use FindBin qw($Bin);
use lib $Bin;
use DNA;
use Data::Dumper;

##
# Parse a fasta file into the %ref hash
#
sub parseFasta($$) {
	my ($fa, $ref) = @_;
	print STDERR "Parsing FASTA file $fa...\n";
	my $fapipe = "$fa";
	$fapipe = "gzip -dc $fa |" if $fa =~ /\.gz$/;
	open(FA, $fapipe) || die "Could not open '$fapipe' for reading";
	my $name = "";
	my $bases = 0;
	while(<FA>) {
		chomp;
		if(/^>/) {
			$name = substr($_, 1);
			$name =~ s/\s.*//;
			print STDERR "  parsing sequence $name...\n";
			$ref->{std}{$name} = "";
		} else {
			$name ne "" || die "sequence before name";
			$bases += length($_);
			$ref->{std}{$name} .= $_;
		}
	}
	print STDERR "Parsed $bases reference bases\n";
	close(FA);
}

##
# Create a new alignment checker
#
sub new {
	my (
		$class,
		$name,   # name of checker
		$fas,    # list of fasta files containing reference sequences
		         # or hash of the references themselves (former if it's
		         # an array ref latter if it's a hash ref)
		$format, # alignment format
		$bisC,   # whether alignment was w/r/t ref w/ all Cs converted to Ys
		$bisCpG, # whether alignment was w/r/t ref w/ all CpGs converted to YpGs
	) = @_;
	(defined($fas)) || die "Must specify non-empty list of fasta files";
	# Set defaults for format, bisC, bisCpG, name
	$format = "bowtie" unless defined($format);
	$bisC   = 0 unless defined($bisC);
	$bisCpG = 0 unless defined($bisCpG);
	$name = "noname" unless defined($name);
	# Parse all the fasta files into the ref hash
	my %ref = ();
	if(ref($fas) eq "HASH") {
		for (keys %$fas) {
			$ref{std}{$_} = $fas->{$_}
		}
	} else {
		ref($fas) eq "ARRAY" || die;
		foreach (@$fas) { parseFasta($_, \%ref); }
	}
	return bless {
		_name       => $name,
		_fas        => $fas,
		_format     => $format,
		_bisC       => $bisC,
		_bisCpG     => $bisCpG,
		_refs       => \%ref,
		_nals       => 0,
		_nedits     => 0
	}, $class;
}
sub name      { return $_[0]->{_name}      }
sub fas       { return $_[0]->{_fas}       }
sub format    { return $_[0]->{_format}    }
sub bisC      { return $_[0]->{_bisC}      }
sub bisCpG    { return $_[0]->{_bisCpG}    }
sub refs      { return $_[0]->{_refs}      }
sub nals      { return $_[0]->{_nals}      }
sub nedits    { return $_[0]->{_nedits}    }
sub nrefs     { return scalar(keys %{$_[0]->{_refs}}); }

##
# Given a sequence that represents the read oriented s.t. the Watson-
# upstream end is on the left, and given a set of edits to apply to the
# read, also oriented assuming that Watson-upstream is on the left,
# return the string corresponding to the read mutated to match the
# reference.
#
my $nedits = 0;
sub applyEdits($$$) {
	my ($seq, $edits, $line) = @_;
	my $rfseq = $seq;
	my $lpos = length($seq)+1;
	$nedits += scalar(@$edits);
	foreach (reverse @$edits) {
		next unless defined($_);
		#print STDERR "Applying edit at $_->{pos}\n";
		# Apply the edit
		$_->{pos} <= $lpos || die "Edit position $_->{pos} was greater than previous $lpos";
		if($_->{qchr} eq "-") {
			# Insert
			$_->{pos} <= length($rfseq) || die "Read gap pos $_->{pos} was not <= string len ".length($rfseq)."\n$line";
			substr($rfseq, $_->{pos}, 0) = $_->{chr};
		} elsif($_->{chr} eq "-") {
			# Deletion
			$_->{pos} < length($rfseq) || die "Ref gap pos $_->{pos} was not < string len ".length($rfseq)."\n$line";
			my $dc = substr($rfseq, $_->{pos}, 1);
			$dc eq $_->{qchr} ||
				die "Edit: $_->{pos}:$_->{chr}>$_->{qchr} but ref char was $dc".
				    "\n$rfseq\n$line";
			substr($rfseq, $_->{pos}, 1) = "";
		} else {
			# Mismatch
			$_->{pos} < length($rfseq) || die "Mismatch pos $_->{pos} was not < string len ".length($rfseq)."\n$line";
			substr($rfseq, $_->{pos}, 1) eq $_->{qchr} ||
				die "Edit: $_->{pos}:$_->{chr}>$_->{qchr}\n$rfseq\n$line";
			substr($rfseq, $_->{pos}, 1) = $_->{chr};
		}
	}
	return $rfseq;
}

##
# Given a list of Bowtie edits, invert them by reversing the list and
# changing the poss to be with respect to the other end of the read.
#
sub invertEdits($$) {
	my ($edits, $len) = @_;
	@$edits = reverse @$edits;
	for (@$edits) {
		next unless defined($_);
		defined($_->{qchr}) || die;
		if($_->{qchr} eq "-") {
			$_->{pos} = $len - $_->{pos};
			length($_->{chr}) >= 1 || die;
			$_->{chr} = reverse $_->{chr};
		} else {
			$_->{pos} = $len - $_->{pos} - 1;
			length($_->{chr}) == 1 || die;
		}
	}
}

##
# Given an edit string, parses it into a list of hashes and returns
# the list.
#
sub parseEdits($) {
	my $editstr = shift;
	return undef if (!defined($editstr) || $editstr eq "-" || $editstr eq "");
	my @edits = ();
	# For each edit
	for (split(/,/, $editstr)) {
		# Parse pos
		my ($pos, $ed) = split(/:/);
		defined($ed) || die;
		# Parse chr, qchr
		my ($chr, $qchr) = split(/>/, $ed);
		push @edits, { pos => $pos, chr => $chr, qchr => $qchr };
	}
	return \@edits;
}

##
# Given an edit string, possibly with 4 semicolon-delimited fields,
# parse it into a set of 4 lists of hashes and return the set as an
# array ref.
#
sub parseAllEdits($) {
	my $editstr = shift;
	return [ undef, undef, undef, undef ] if ($editstr eq "-" || $editstr eq "");
	my @editls = split(/;/, $editstr, -1);
	if(scalar(@editls) > 1) {
		scalar(@editls) == 4 || die;
		return [
			parseEdits($editls[0]),
			parseEdits($editls[1]),
			parseEdits($editls[2]),
			parseEdits($editls[3]) ];
	} else {
		scalar(@editls) == 1 || die;
		return [
			parseEdits($editls[0]),
			undef,
			undef,
			undef ];
	}
}

##
# Given array refs for two lists of edits, one corresponding to the
# nucleotide edit list and the other corresponding to the resolved
# ambiguous base list, eliminate any edits that appear in both lists.
# Really this shouldn't happen, but I observe that merman does report
# an edit in both categories if the reference base is being resolved to
# an incompatible nucleotide, e.g. R>C.
#
sub removeDups($$) {
	my ($ned, $aed) = @_;
	return unless (defined($ned) && defined($aed));
	for my $i (0..scalar(@$ned)-1) {
		next unless defined($ned->[$i]);
		for my $j (0..scalar(@$aed)-1) {
			next unless defined($aed->[$j]);
			if($ned->[$i]->{qchr} eq $aed->[$j]->{qchr} &&
			   $ned->[$i]->{chr}  eq $aed->[$j]->{chr} &&
			   $ned->[$i]->{pos}  eq $aed->[$j]->{pos})
			{
				#print STDERR "  Eliminated a duplicate edit\n";
				$aed->[$j] = undef;
			}
		}
	}
}

##
# Take all the references in %ref and make both Watson and Crick
# versions where the sequence is in-silico bisulfite treated.
#
sub bisulfiteC($) {
	my $ref = shift;
	for(keys %{$ref->{std}}) {
		$ref->{bisc_fw}{$_} = $ref->{std}{$_};
		$ref->{bisc_fw}{$_} = s/C/Y/g;
		$ref->{bisc_rc}{$_} = $ref->{std}{$_};
		$ref->{bisc_rc}{$_} = s/G/R/g;
	}
}

##
# Take all the references in %ref and make both Watson and Crick
# versions where the sequence is in-silico bisulfite treated.
#
sub bisulfiteCpG($) {
	my $ref = shift;
	for(keys %{$ref->{std}}) {
		$ref->{biscpg_fw}{$_} = $ref->{std}{$_};
		$ref->{biscpg_fw}{$_} =~ s/CG/YG/g;
		$ref->{biscpg_fw}{$_} =~ s/C/T/g;
		$ref->{biscpg_rc}{$_} = $ref->{std}{$_};
		$ref->{biscpg_rc}{$_} =~ s/CG/CR/g;
		$ref->{biscpg_rc}{$_} =~ s/G/A/g;
	}
}

##
# Given a Bowtie orientation string, return true iff the 5' end of the
# read is at the left of the printed sequence.
#
sub fivePrimeLeft($) {
	return ($_[0] eq "+" || $_[0] eq "W" || $_[0] eq "CR");
}

##
# Given the orientation of the read and the state of the global
# bisulfite variables, determine which version of the reference to
# compare against.
#
sub calcRefType {
	my ($self, $orient) = @_;
	if($self->bisC || $self->bisCpG) {
		if($orient eq "W" || $orient eq "WR" || $orient eq "+") {
			return $self->bisC ? "bisc_fw" : "biscpg_fw";
		} else {
			$orient eq "C" || $orient eq "CR" || $orient eq "-" || die;
			return $self->bisC ? "bisc_rc" : "biscpg_rc";
		}
	} else {
		return "std";
	}
}

##
# Parse a CIGAR string into parallel arrays of CIGAR operations (M, I, D)
#
sub cigarParse($$$) {
	my ($cigar, $ops, $runs) = @_;
	my $i = 0;
	while($i < length($cigar)) {
		substr($cigar, $i) =~ /^([0-9]+)/;
		defined($1) || die "Could not parse number at pos $i: '$cigar'";
		$i += length($1);
		$i < length($cigar) || die;
		push @$runs, $1;
		my $op = substr($cigar, $i, 1);
		defined($op) || die "Could not parse operation at pos $i: '$cigar'";
		push @$ops, $op;
		$i++;
	}
}

##
# Trim a read sequence according to the soft clipping in the CIGAR string.
#
sub cigarTrim($$) {
	my ($seq, $cigar) = @_;
	my @ops = ();
	my @runs = ();
	cigarParse($cigar, \@ops, \@runs);
	my ($trimup, $trimdn) = (0, 0);
	if($ops[0] eq 'S') {
		$runs[0] < length($seq) || die "Soft clipped entire alignment!";
		$seq = substr($seq, $runs[0]);
		$trimup = $runs[0];
	}
	if(scalar(@ops) > 1 && $ops[-1] eq 'S') {
		$runs[-1] < length($seq) || die "Soft clipped entire alignment!";
		$seq = substr($seq, 0, -$runs[-1]);
		$trimdn = $runs[-1];
	}
	return ($seq, $trimup, $trimdn);
}

##
# Parse a CIGAR string into a string of operators.  Operators are expanded into
# runs where appropriate.  = and X are collapsed into M.
#
sub parse_cigar($) {
	my ($cigar) = @_;
	my $ret = "";
	my $i = 0;
	my ($rdlen, $rflen) = (0, 0);
	while($i < length($cigar)) {
		substr($cigar, $i) =~ /^([0-9]+)/;
		defined($1) || die "Could not parse number at pos $i: '$cigar'";
		my $runlen = $1;
		$i += length($1);
		$i < length($cigar) || die;
		my $op = substr($cigar, $i, 1);
		defined($op) || die "Could not parse operation at pos $i: '$cigar'";
		if($op eq "X" || $op eq "=") {
			$op = "M";
		}
		$rdlen += $runlen if $op ne "D";
		$rflen += $runlen if $op ne "I";
		$ret .= ($op x $runlen);
		$i++;
	}
	return ($ret, $rdlen, $rflen);
}

##
# Parse an MD:Z string into a string with length equal to query length.  Each
# position contains either a space, if the read matches the reference at that
# position, or a character, if the reference contains a character that doesn't
# match its opposite in the alignment.  In the latter case, the character in
# the string is the reference character.
#
sub parse_md($) {
	my ($md) = @_;
	my $i = 0;
	my $ret = "";
	while($i < length($md)) {
		# Starts with a number?
		my $ch = substr($md, $i, 1);
		if($ch =~ /[0-9]/) {
			# Parse the number off the beginning
			substr($md, $i) =~ /^([0-9]+)/;
			defined($1) || die "Could not parse number at pos $i: '$md'";
			my $runlen = $1;
			$ret .= (" " x $runlen) if $runlen > 0;
			$i += length($runlen);
		} elsif($ch eq "^") {
			# Read gap
			$i++;
			substr($md, $i) =~ /^([A-Za-z]+)/;
			defined($1) || die "Could not parse read gap at pos $i: '$md'";
			my $chrs = $1;
			$i += length($chrs);
			$ret .= $chrs;
		} else {
			# DNA character
			$ch =~ /[ACGTN.]/i || die "Bad char '$ch' at pos $i: '$md'";
			$ret .= $ch;
			$i++;
		}
	}
	return $ret;
}

##
# Given a read sequence (with characters in upstream-to-downstream order with
# respect to the reference - NOT necessarily 5'-to-3') and a CIGAR string and
# an MD:Z string, build the alignment strings.  The alignment strings will only
# contain the portion of the read that aligned.  Any portions that were either
# hard-trimmed or soft-trimmed are trimmed from this function's result.
#
# For now, I'm assuming that the MD:Z string only describes aligned characters,
# i.e. *after* trimming.
#
sub _read_md_cigar_to_al($$$) {
	my ($seq, $md, $cigar) = @_;
	my $alread = "";
	my $alref = "";
	$cigar ne "*" || die "CIGAR string was star!";
	$seq ne "" || die "Empty sequence given to _read_md_cigar_to_al";
	my $parsed_md  = parse_md($md);
	my ($parsed_cig, $cigar_rdlen, $cigar_rflen) = parse_cigar($cigar);
	my ($rdoff, $mdoff) = (0, 0);
	my ($htriml, $striml, $htrimr, $strimr) = (0, 0, 0, 0);
	my $nonsh = 0; # have I seen a non-S, non-H CIGAR op?
	my $nonh = 0;  # have I seen a non-H CIGAR op?
	for(my $i = 0; $i < length($parsed_cig); $i++) {
		my $cigop = substr($parsed_cig, $i, 1);
		$nonh++ if $cigop ne "H";
		$nonsh++ if ($cigop ne "H" && $cigop ne "S");
		if($cigop eq "S") {
			if($nonsh) {
				$strimr++;
			} else {
				$striml++;
			}
			$rdoff++;
			next;
		}
		if($cigop eq "H") {
			if($nonh) {
				$htrimr++;
			} else {
				$htriml++;
			}
			next;
		}
		$cigop = "M" if $cigop eq "=" || $cigop eq "X";
		if($cigop eq "P") {
			# Padding
			$alread .= "-";
			$alref .= "-";
		} elsif($cigop eq "M") {
			my $rdc = substr($seq, $rdoff, 1);
			$mdoff < length($parsed_md) ||
				die "Bad mdoff ($mdoff)\nlength(parsed_md)=".length($parsed_md)."\nseq:\n$seq\ncigar:\n$cigar\nmd:\n$md\nparsed md:\n$parsed_md";
			my $rfc = substr($parsed_md, $mdoff, 1);
			$rfc = $rdc if $rfc eq " ";
			$alread .= $rdc;
			$alref .= $rfc;
			$rdoff++;
			$mdoff++;
		} elsif($cigop eq "D") {
			# Read gap
			#  Read: AAA-AAA
			#   Ref: AAAAAAA
			my $rfc = substr($parsed_md, $mdoff, 1);
			$rfc ne " " ||
				die "Must have a ref character opposite a gap in the read:\n".
				    "cig: $parsed_cig ($i)\nmd:  $parsed_md ($mdoff)\n";
			$alread .= "-";
			$alref .= $rfc;
			$mdoff++;
		} else {
			# Reference gap
			#  Read: AAAAAAA
			#   Ref: AAA-AAA
			$cigop eq "I" || die "Unsupported cigop: $cigop in cigar $cigar";
			my $rdc = substr($seq, $rdoff, 1);
			$alread .= $rdc;
			$alref .= "-";
			$rdoff++;
			# $mdoff SHOULD NOT be incremented here
		}
		$rdoff <= length($seq) ||
			die "Bad rdoff:$rdoff for seq '$seq' cigop=$cigop\nseq: $seq\ncigar=$cigar\nmd=$md";
	}
	return ($alread, $alref, $htriml, $striml, $htrimr, $strimr);
}

##
# Parse a line from a Bowtie alignment file and check that the
# alignment is sane and consistent with the reference.
#
sub parseBowtieLines {
	my ($self, $lines) = @_;
	for my $line (@$lines) {
		chomp($line);
		my ($rdname, $orient, $refname, $off, $seq, $qual, $oms, $editstr,
		    $flags) = split(/\t/, $line, -1);
		next if $refname eq "*";
		$flags =~ /XC:([^,\s]+)/;
		my $cigar = $1;
		defined($cigar) ||
			die "Could not parse CIGAR string from flags: '$flags'";
		defined($editstr) || die "Failed to get 8 tokens from line:\n$_";
		$off == int($off) || die "Offset field (col 4) must be an integer:\n$_";
		$oms == int($oms) || die "OMS field (col 7) must be an integer:\n$_";
		my $reftype = $self->calcRefType($orient);
		defined($self->refs->{$reftype}{$refname}) ||
			die "No such refname as $refname for reftype $reftype:\n".
			Dumper($self->refs->{$reftype});
		my $edits4 = parseAllEdits($editstr);
		my ($ned, $aed) = ($edits4->[0], $edits4->[1]);
		removeDups($ned, $aed);
		my $fpl = fivePrimeLeft($orient);
		# Trim seq according to CIGAR string
		my $rfseq = $seq;
		my ($trimup, $trimdn);
		($rfseq, $trimup, $trimdn) = cigarTrim($rfseq, $cigar);
		invertEdits($ned, length($rfseq)) unless ($fpl || !defined($ned));
		invertEdits($aed, length($rfseq)) unless ($fpl || !defined($aed));
		$rfseq = applyEdits($rfseq, $ned, $line) if defined($ned);
		$rfseq = applyEdits($rfseq, $aed, $line) if defined($aed);
		# Check if our alignment falls off the end of the reference, in
		# which case we need to pad the reference string with Ns
		my $exoff = $off;
		my $padleft = "";
		my $exlen = length($rfseq);
		my $tlen = length($self->refs->{$reftype}{$refname});
		if($exoff < 0) {
			# Alignment hangs off LHS; pad it
			my $npad = -$exoff;
			$padleft = "N" x $npad;
			$exlen += $exoff;
			$exlen >= 0 ||
				die "Read was entirely off the LHS of the reference\n".
					"Referemce len=$tlen\n".
					"Alignment referemce len=$tlen\n".
					"$line\n";
			$exoff = 0;
		}
		my $padright = "";
		my $roverhang = $off + length($rfseq) - $tlen;
		if($roverhang > 0) {
			$padright = "N" x $roverhang;
			$exlen -= $roverhang;
			$exlen >= 0 ||
				die "Read was entirely off the RHS of the reference\n".
					"Referemce len=$tlen\n".
					"Alignment referemce len=$tlen\n".
					"$line\n";
		}
		my $refsub = substr($self->refs->{$reftype}{$refname}, $exoff, $exlen);
		length($refsub) == $exlen ||
			die "Tried to extract ref substring of length $exlen, got ".
			    "\"$refsub\" from \"".$self->refs->{$reftype}{$refname}."\"".
				"\n$line\n".
				"\noff=$off, rfseq=$rfseq\n";
		$refsub = DNA::iupacSubN($refsub);
		my $trueRfseq = $padleft . $refsub . $padright;
		length($trueRfseq) == length($rfseq) ||
			die "Different lengths for edited read and ref:\n".
			"       Read: $seq\n".
			"Edited read: $rfseq\n".
			"        Ref: $trueRfseq\n";
		$rfseq eq $trueRfseq ||
			die "Did not match:\n".
			"       Read: $seq\n".
			"Edited read: $rfseq\n".
			"        Ref: $trueRfseq\n";
		$self->{_nals}++;
	}
}

##
# Parse a line from a SAM alignment file and check that the
# alignment is sane and consistent with the reference.
#
sub parseSamLines {
	my ($self, $lines) = @_;
	my ($lastseq, $lastqual) = ("", "");
	my $idx = 0;
	for my $line (@$lines) {
		$idx++;
		print STDERR "Processing line...\n";
		chomp($line);
		next if $line =~ /^\@/;
		my @toks = split(/\t/, $line, -1);
		my (
			$qname, #1
			$flag,  #2
			$rname, #3
			$pos,   #4
			$mapq,  #5
			$cigar, #6
			$rnext, #7
			$pnext, #8
			$tlen,  #9
			$seq,   #10
			$qual) = @toks;
		defined($qual) || die "Not enough SAM tokens:\n$line\n";
		my @opt_flags_list = @toks[11..$#toks];
		my %opt_flags = ();
		next if $cigar eq "*"; # Failed to align
		# Get the read sequence & qualities from previous record if necessary
		if($seq eq "*") {
			$lastseq ne "" || die "Line #$idx:\n$line";
			$seq = $lastseq;
			$qual = $lastqual;
		} else {
			$lastseq = $seq;
			$lastqual = $qual;
		}
		$seq ne "*" || die;
		my ($parsed_cigar, $rdlen_cigar, $rflen_cigar) = parse_cigar($cigar);
		length($seq) == $rdlen_cigar ||
			die "Sequence length and parsed cigar string length ($rdlen_cigar) mismatch:\n".
			    "$seq\n$parsed_cigar\nLine:\n$line";
		# Stick optional flags into a hash
		for my $fl (@opt_flags_list) {
			my @fs = split(/:/, $fl, -1);
			scalar(@fs) > 2 || die "Bad optional flag: $fl\n$line\n";
			$opt_flags{$fs[0]}{type} = $fs[1];
			$opt_flags{$fs[0]}{value} = join(":", @fs[2..$#fs]);
		}
		defined($opt_flags{"MD"}) || die "No MD:Z flag:\n$line\n";
		$opt_flags{"MD"}{type} eq "Z" || die "Bad type for MD:Z flag\n$line\n";
		my $md = $opt_flags{"MD"}{value};
		$pos   == int($pos)   || die "POS field (col 4) must be an int:\n$line\n";
		$pnext == int($pnext) || die "PNEXT field (col 8) must be an int:\n$line\n";
		$tlen  == int($tlen)  || die "TLEN field (col 9) must be an int:\n$line\n";
		$mapq  == int($mapq)  || die "MAPQ field (col 5) must be an int:\n$line\n";
		# TODO: deal with bisulfite strands??
		my $fw = (($flag & 0x10) == 0);
		my $orient = $fw ? "+" : "-";
		my $reftype = $self->calcRefType($orient);
		defined($self->refs->{$reftype}{$rname}) ||
			die "No such refname as $rname for reftype $reftype:\n$line\n".
			Dumper($self->refs->{$reftype});
		my $exoff = $pos-1; # expected 0-based reference offset
		my ($alread, $alref, $htriml, $striml, $htrimr, $strimr) =
			_read_md_cigar_to_al($seq, $md, $cigar);
		print STDERR "$alread\n$alref\n";
		my $rfseq = $alref;
		$rfseq =~ s/[ -]//g; # remove spaces & gaps
		my $exlen = length($rfseq);
		my $refsub = substr($self->refs->{$reftype}{$rname}, $exoff, $exlen);
		length($refsub) == $exlen ||
			die "Tried to extract ref substring of length $exlen from:\n".
			    $self->refs->{$reftype}{$rname}.
				"\ngot string of length ".length($refsub).":\n".
				$refsub.
				"\nfrom:\n".
				$line.
				"\nexlen is the length of:\n$rfseq\npos=$pos, rfseq=$rfseq\n";
		$refsub = DNA::iupacSubN($refsub);
		my $trueRfseq = $refsub;
		length($trueRfseq) == length($rfseq) ||
			die "Different lengths for edited read and ref:\n".
			"       Read: $seq\n".
			"Edited read: $rfseq\n".
			"        Ref: $trueRfseq\n";
		$rfseq eq $trueRfseq ||
			die "Did not match:\n".
			"       Read: $seq\n".
			"Edited read: $rfseq\n".
			"        Ref: $trueRfseq\n".
			"$line";
		$self->{_nals}++;
	}
}

##
# Parse lines from a Bowtie alignment file and check that the
# alignments are sane and consistent with the reference.
#
sub parseBowtie {
	my ($self, $fh) = @_;
	while(<$fh>) {
		$self->parseBowtieLines([$_]);
	}
}

##
# Parse lines from a SAM alignment file and check that alignments are
# sane and consistent with the reference.
#
sub parseSam {
	my ($self, $fh) = @_;
	my @lines = ();
	while(<$fh>) { push @lines, $_; }
	$self->parseSamLines(\@lines);
}

##
# Parse lines from an alignment file of the type given by self->format
#
sub parseLines {
	my ($self, $lines) = @_;
	if($self->format eq "bowtie") {
		$self->parseBowtieLines($lines);
	} else {
		$self->format eq "sam" || die;
		$self->parseSamLines($lines);
	}
}

##
# Parse lines from an alignment file of the type given by self->format
#
sub parse {
	my ($self, $fh) = @_;
	if($self->format eq "bowtie") {
		$self->parseBowtie($fh);
	} else {
		$self->format eq "sam" || die;
		$self->parseSam($fh);
	}
}

##
# Print summary of how many alignments and edits were checked.
#
sub printSummary {
	my $self = shift;
	print STDERR "--- Summary ---\n";
	print STDERR "Read ".scalar(keys %{$self->refs})." reference strings\n";
	print STDERR "Checked $self->{_nals} alignments, $self->{_nedits} edits\n";
	print STDERR "---------------\n";
	print STDERR "PASSED\n";
}

##
# Check the given batch of alignments.  We check that they're
# internally consistent in some basic ways, and we check that the
# sequence and edits are consistent with the reference.
#
# The $als argument is either a list of (possibly compressed) filenames
# of files containing alignments, or a list of alignment strings.  If
# the former, $filenames is non-zero.
#
sub checkAlignments {
	my ($self, $als, $filenames) = @_;
	if($self->bisC) {
		print STDERR "Generating all-C bisulfite-treated references\n";
		bisulfiteC($self->refs);
	}
	if($self->bisCpG) {
		print STDERR "Generating all-CpG bisulfite-treated references\n";
		bisulfiteCpG($self->refs);
	}
	if($filenames) {
		foreach (@$als) {
			my $alnpipe = $_;
			print STDERR "Processing alignment file '$_'\n";
			$alnpipe = "gzip -dc $_ |" if ($_ =~ /\.gz$/);
			my $alnfh = undef;
			open($alnfh, $alnpipe) || die "Could not open '$alnpipe' for reading";
			$self->parse($alnfh);
			close($alnfh);
		}
	} else {
		$self->parseLines($als);
	}
}

##
# Check simple alignments
#
sub test1 {
	my $ac = AlignmentCheck->new(
		"AlignmentCheck.pm test1 checker",
		{ "r1" => "TTGTTCGT" },
		"bowtie",
		0,
		0
	);
	$ac->checkAlignments([
		"0\t+\tr1\t1\tTGTTCGT\tIIIIIII\t40\t-",
		"1\t+\tr1\t0\tTTGTTCG\tIIIIIII\t40\t-",
		"2\t+\tr1\t2\tGTTCGTA\tIIIIIII\t40\t6:N>A",
		"3\t+\tr1\t-1\tATTGTTC\tIIIIIII\t40\t0:N>A"], 0);
	return 1;
}

##
# Check simple alignments from files
#
sub test2 {
	open(TMP, ">/tmp/.AlignmentCheck.pm.fa") || die;
	print TMP ">r1\n";
	print TMP "TTGTTCGT\n";
	close(TMP);
	my $ac = AlignmentCheck->new(
		"AlignmentCheck.pm test1 checker",
		[ "/tmp/.AlignmentCheck.pm.fa" ],
		"bowtie",
		0,
		0
	);
	$ac->checkAlignments([
		"0\t+\tr1\t1\tTGTTCGT\tIIIIIII\t40\t-",
		"1\t+\tr1\t0\tTTGTTCG\tIIIIIII\t40\t-",
		"2\t+\tr1\t2\tGTTCGTA\tIIIIIII\t40\t6:N>A",
		"3\t+\tr1\t-1\tATTGTTC\tIIIIIII\t40\t0:N>A"], 0);
	return 1;
}

if($0 =~ /[^0-9a-zA-Z_]?AlignmentCheck\.pm$/) {
	my @fas = ();
	my @als = ();
	my $format = "sam";
	my $bisC = 0;
	my $bisCpG = 0;
	my $test = 0;
	
	use Getopt::Long;
	GetOptions (
		"test"                  => \$test,
		"fasta|ref=s"           => \@fas,
		"als|aln=s"             => \@als,
		"format=s"              => \$format,
		"bis-C|bisulfite-C"     => \$bisC,
		"bis-CpG|bisulfite-CpG" => \$bisCpG,
		"bowtie"                => sub {$format = "bowtie"},
		"sam"                   => sub {$format = "sam"}) || die;

	if($test) {
		use Test;
		print "Running unit tests\n";
		# Run unit tests
		Test::shouldSucceed("test1", \&test1);
		Test::shouldSucceed("test2", \&test2);
		exit 0;
	}

	my $ac = AlignmentCheck->new(
		"AlignmentCheck.pm checker",
		\@fas,
		$format,
		$bisC,
		$bisCpG);
	$ac->checkAlignments(\@als, 1);
}

1;