# -*-Perl-*-

##############################################
# Copyright (C) 1995,1996 Simon Dear
# Genome Research Limited, Cambridge, United Kingdom
# Copyright (C) 1994, 1995, 1996 Gabor T. Marth
# Washington University, St. Louis
# All Rights Reserved.
##############################################

##################################################################################
# LIBRARY OF CAF-RELATED PARSING SUBROUTINES
# AUTHORS: Simon Dear
#          Gabor T. Marth
##################################################################################

package library_parse;

use strict;

#############################################################
# Methods to be exported
#############################################################

require Exporter;

@library_parse::ISA = qw (Exporter);
@library_parse::EXPORT = qw (findFile readAceFile2 writeAceFile2 writeAceFileSpec theContigs2 theAssembly2 
			     trashRead2 browse_ace
			     get_config get_problem_types get_subclone_types get_read_types get_parameters
			     get_calling_arguments
			     get_problems get_ordered_reads parse_caf_problemtypes
			     get_caf_read 
			     get_caf_contig 
			     write_read_into_caf 
			     write_contig_into_caf 
			     get_caf_order 
			     get_region_list
			     parse_crossmatch_output 
);


#############################################################
# Methods to be imported
#############################################################

use lib (".", split(':', $ENV{CAF_LIB}));

use Objects;
use library_math;
use library_convolution;

##################################################################################
# SUBROUTINES
##################################################################################

#############################################################
# Find a file in the path
#############################################################

sub findFile {

    # Input
    my ($file,@path) = @_;

    # Internal
    my ($dir);

    foreach $dir ( @path ) {
	return "$dir/$file" if (-f "$dir/$file");
    }
    return $file;
}

##################################################################################
# ACE FILE PARSING ROUTINES
# CREATED BY: Simon Dear
##################################################################################

#############################################################
# $ace = &readAceFile2(HANDLE)
# Read in and parse an ace file from a HANDLE.
# Returns ref to structure
#############################################################

sub readAceFile2 {
    my ($file) = @_ ;
    my ($ace) = {};
    my ($class, $colon, $name);

    # Set input record separator to paragraph mode
    $/ = "";

    while (<$file>) {

        s/\/\/.*// ;			# ace file comments

	if ( ($class, $colon, $name) = /^(\S+)(\s+:)?\s+(\S+)/ ) {
	    s/^.*\n//;
	    s/^\s+//;                   # leading white space
	    if ($class eq "DNA") {
		$ace->{"DNA"}->{$name} = $_;
		$ace->{"DNA"}->{$name} =~ s/\s//g;
            } elsif ($class eq "BaseQuality") {
		$ace->{"BaseQuality"}->{$name} = [ split(/\s+/, $_) ];
            } elsif ($class eq "BasePosition") {
		$ace->{"BasePosition"}->{$name} = [ split(/\s+/, $_) ];
            } elsif ($class eq "DistributionIndex") {
		$ace->{"DistributionIndex"}->{$name} = [ split(/\s+/, $_) ];
            } elsif ($class eq "DistributionValue") {
		$ace->{"DistributionValue"}->{$name} = [ split(/\s+/, $_) ];
	    } else {
		if (! defined $ace->{$class}->{$name}) {
		    $ace->{$class}->{$name} = {};
		}
		foreach ( split(/\n/,$_) ) {
		    s/\\\n\n// ;		# continuation lines
		    my ($tag,@fields) = &parseAceLine($_);
		    push(@{$ace->{$class}->{$name}->{$tag}},[@fields]);
		}
	    }
	}
	if (eof($file)) { last ; }
    }

    # Reinstate input field separator
    $/ = "\n";

    $ace;
}


###################################################
# Parses a single line in ACE format
###################################################

sub parseAceLine {
    my ($line) = @_;
    my (@fields) = ();
    $line =~ /\s*/;     # Consume opening space
    while ($line = $') {
        if ($line =~ /^"((([^"\\])|(\\.))*)"\s*/ || $line =~ /^(\S+)\s*/ || $line =~ /.*/) { #"
            push(@fields,$1);
        }
    }
    @fields;
}



#############################################################
# &writeAceFile2(HANDLE,$ace)
# Writes out an ace file to a HANDLE.
#############################################################

sub writeAceFile2 {
    my ($file, $ace) = @_ ;
    my ($class, $len, $name, $tag, $i);

    foreach $class ( keys %{$ace} ) {

	foreach $name ( keys %{$ace->{$class}} ) {

	    print $file "$class $name\n";
	    if ($class eq "DNA") {
		$len = length ($ace->{$class}->{$name});
		for ($i = 0 ; $i < $len ; $i += 60) {
		    print $file substr ($ace->{$class}->{$name}, $i, 60), "\n" ;
		}
	    } 
	    elsif ($class eq "BaseQuality") {
		print $file "@{$ace->{$class}->{$name}}\n";
	    }
	    elsif ($class eq "DistributionIndex") {
		print $file "@{$ace->{$class}->{$name}}\n";
	    } 
	    elsif ($class eq "DistributionValue") {
		print $file "@{$ace->{$class}->{$name}}\n";
	    } 
	    else {
		foreach $tag ( keys %{$ace->{$class}->{$name}} ) {
		    for ($i = 0; $i <= $#{$ace->{$class}->{$name}->{$tag}}; $i++) {
			print $file &createAceLine($tag,@{$ace->{$class}->{$name}->{$tag}->[$i]}),"\n";
		    }
		}
	    }
	    print $file "\n";
	}

    }
}

#############################################################
# &writeAceFileSpec(HANDLE,$ace)
# Writes out an ace object to a HANDLE.
# If name specified, only writes that object
#############################################################

sub writeAceFileSpec {

    my ($file, $ace, $name) = @_ ;
    my ($class, $len, $tag, $i);

    foreach $class ( keys %{$ace} ) {

	if ( $name ) {
	    if ( $ace->{$class}->{$name} ) {

		print $file "$class $name\n";
		if ($class eq "DNA") {
		    $len = length ($ace->{$class}->{$name});
		    for ($i = 0 ; $i < $len ; $i += 50) {
			print $file substr ($ace->{$class}->{$name}, $i, 50), "\n" ;
		    }
		} 
		elsif ($class eq "BaseQuality") {
		    print $file "@{$ace->{$class}->{$name}}\n";
		}
		elsif ($class eq "BasePosition") {
		    print $file "@{$ace->{$class}->{$name}}\n";
		}
		elsif ($class eq "DistributionIndex") {
		    print $file "@{$ace->{$class}->{$name}}\n";
		} 
		elsif ($class eq "DistributionValue") {
		    print $file "@{$ace->{$class}->{$name}}\n";
		} 
		else {
		    foreach $tag ( keys %{$ace->{$class}->{$name}} ) {
			for ($i = 0; $i <= $#{$ace->{$class}->{$name}->{$tag}}; $i++) {
			    print $file &createAceLine($tag,@{$ace->{$class}->{$name}->{$tag}->[$i]}),"\n";
			}
		    }
		}
		print $file "\n";
	    }
	}
	else {
	    foreach $name ( keys %{$ace->{$class}} ) {

		print $file "$class $name\n";
		if ($class eq "DNA") {
		    $len = length ($ace->{$class}->{$name});
		    for ($i = 0 ; $i < $len ; $i += 60) {
			print $file substr ($ace->{$class}->{$name}, $i, 50), "\n" ;
		    }
		} 
		elsif ($class eq "BaseQuality") {
		    print $file "@{$ace->{$class}->{$name}}\n";
		}
		elsif ($class eq "BasePosition") {
		    print $file "@{$ace->{$class}->{$name}}\n";
		}
		elsif ($class eq "DistributionIndex") {
		    print $file "@{$ace->{$class}->{$name}}\n";
		} 
		elsif ($class eq "DistributionValue") {
		    print $file "@{$ace->{$class}->{$name}}\n";
		} 
		else {
		    foreach $tag ( keys %{$ace->{$class}->{$name}} ) {
			for ($i = 0; $i <= $#{$ace->{$class}->{$name}->{$tag}}; $i++) {
			    print $file &createAceLine($tag,@{$ace->{$class}->{$name}->{$tag}->[$i]}),"\n";
			}
		    }
		}
		print $file "\n";
	    }
	}

    }
}

###################################################
# Composes a printable ACE format line
###################################################

sub createAceLine {
    my (@nice) = ();
    foreach ( @_ ) {
	if (/\s/) {
	    push(@nice,"\"$_\"");
        } else {
	    push(@nice,$_);
        }
    }

    join(" ",@nice);
}

#############################################################
# &browse(*var,$level)
# Recurses through a structure, printing out all substructures as it goes.
#############################################################

sub browse_ace {
    my ($x,$level) = @_;
    my ($type) = ref($x);
    my ($i);


    print " " x $level;
    if (defined $type && $type ne "") {
	print "TYPE=$type";
    } else {
	print "$x\n"; return;
    }

    if ($type eq "REF") {
	print "\n";
	&browse_ace(${$x},$level+4);
    } elsif ($type eq "SCALAR") {
	print "${$x}\n";
    } elsif ($type eq "ARRAY") {
	print "\n";
	for($i=0;$i<=$#{$x};$i++) {
	    &browse_ace($x->[$i],$level+4);
	}
    } elsif ($type eq "HASH") {
	print "\n";
	foreach $i ( keys %$x ) {
	    print " " x ($level+4), $i,"\n";
	    &browse_ace($x->{$i},$level+8);
	}
    } elsif ($type eq "CODE") {
	print "\n";
    } elsif ($type eq "GLOB") {
	print "\n";
    } else {
	$x->browse;
    }

}

############################################################
# @contigs == &theContigs2($ace);
# Returns a list of all Sequences that are contigs.
############################################################

sub theContigs2 {
    my ($ace) = @_;

    grep(defined $ace->{"Sequence"}->{$_}->{"Is_contig"}, keys %{$ace->{"Sequence"}});
}



############################################################
# @reads == &theReads2($ace);
# Returns a list of all Sequences that are reads.
############################################################

sub theReads2 {
    my ($ace) = @_;

    grep(defined $ace->{"Sequence"}->{$_}->{"Is_read"}, keys %{$ace->{"Sequence"}});
}

############################################################
# @assembly == &theAssembly2($ace);
# Returns a list of all Sequences that are reads.
############################################################

sub theAssembly2 {
    my ($ace) = @_;

    grep(defined $ace->{"Sequence"}->{$_}->{"Is_assembly"}, keys %{$ace->{"Sequence"}});
}



############################################################
# &trashRead2($ace,$read);
# Removes read for data structures.
############################################################

sub trashRead2 {
    #
    # Delete a reading from relevant structures
    #
    my ( $ace, $read ) = @_;
    undef $ace->{"Sequence"}->{$read};
    undef $ace->{"DNA"}->{$read};
}


############################################################
# Read list file of regions
############################################################

sub get_region_list {

    # Input
    my ($file) = @_;

    # Internal
    my ($line, $contig, $from, $to);
    my ($region, @REGIONS);

    # Output
    my ($regions);

    while ( $line = <$file> ) {
	chomp $line;
	($contig, $from, $to) = split(' ', $line);
	if ( defined $contig && defined $from && defined $to ) {
	    $region = new Region;
	    $region->add_field('CONTIGNAME', $contig);
	    $region->add_field('LEFT', $from);
	    $region->add_field('RIGHT', $to);
	    
	    push(@REGIONS, $region);
	}
    }

    # Return output
    if ( $#REGIONS >= 0 ) {
	\@REGIONS;
    }
    else {
	undef;
    }
}


############################################################
# Read configuration file into caf hash structure
############################################################

sub get_config {

    # Input
    my ($filename) = @_;

    # Output
    my ($caf);

    # Internal
    my ($filehandle);

    $filename = &findFile($filename,".",@INC);

    open(FILE, $filename) || die "Could not open configuration file $filename: $!\n";
    $filehandle = "FILE";
    $caf = &readAceFile2(*FILE);
    close(FILE);
    
    # Return output
    $caf;
}

############################################################
# Obtain internal representation of problem types
############################################################

sub get_problem_types {

    # Input
    my ($file) = @_;

    # Output
    my ($caf);

    # Read CAF format file into CAF hash
    $caf = &get_config($file);

    $caf;
}

############################################################
# Obtain internal representation of subclone types
############################################################

sub get_subclone_types {

    # Input
    my ($file) = @_;

    # Internal
    my ($caf);

    # Output
    my ($SUBCLONETYPES_POINTER, $SUBCLONETYPETABLE_POINTER);

    # Read CAF format file into CAF hash
    $caf = &get_config($file);


    # Build SubcloneType objects
    ($SUBCLONETYPES_POINTER, $SUBCLONETYPETABLE_POINTER) = &parse_subclonetypes_to_objects($caf);

    # Return output
    ($SUBCLONETYPES_POINTER, $SUBCLONETYPETABLE_POINTER);
}

sub parse_subclonetypes_to_objects {

    # Input
    my ($caf) = @_;

    # Internal
    my ($name, $line, @lines);
    my ($subclonetype, $subclonetypename, @SUBCLONETYPENAMES);
    my (@sorted_indices, @indices, @values, $i, $dist);
    my ($left, $right, $distr);

    # Output
    my ($SUBCLONETYPES_POINTER, $SUBCLONETYPETABLE_POINTER);

    # Obtain a list of all subclone types
    @SUBCLONETYPENAMES = keys %{$caf->{'SubcloneType'}};

    # Read all fields of all subclone types
    foreach $subclonetypename (@SUBCLONETYPENAMES) {

	$subclonetype = new SubcloneType;

	$subclonetype->{'NAME'} = $subclonetypename;

	# Subclone length distribution
	@indices = ();
	@indices = @{$caf->{'DistributionIndex'}->{$subclonetypename}};
	@values = ();
	@values = @{$caf->{'DistributionValue'}->{$subclonetypename}};

	@sorted_indices = sort { $a <=> $b } @indices;
	$left  = $sorted_indices[0];
	$right = $sorted_indices[$#sorted_indices];

	if ( $#indices >= 0 && $#values >= 0 && $#indices == $#values ) {

	    # Create new Distribution object
	    $distr = new Distribution;

	    $distr->add_field('INDICES', \@indices);
	    $distr->add_field('VALUES', \@values);
	    $distr->add_field('LEFT', $left);
	    $distr->add_field('RIGHT', $right);
	}

	$subclonetype->add_field('LENGTH_DISTRIBUTION', $distr);

	# Enter $subclonetype into list of all subclonetypes
	push(@{$SUBCLONETYPES_POINTER}, $subclonetype);

	# Enter $subclonetype into subclonetype recognition table
	$SUBCLONETYPETABLE_POINTER->{$subclonetype->{'NAME'}} = $subclonetype;
    }

    # Return pointers to SUBCLONETYPES list and SUBCLONETYPETABLE hash
    ($SUBCLONETYPES_POINTER, $SUBCLONETYPETABLE_POINTER);
}

############################################################
# Obtain internal representation of read types
############################################################

sub get_read_types {

    # Input
    my ($file) = @_;

    # Internal
    my ($caf);

    # Output
    my ($READTYPES_POINTER, $READTYPETABLE_POINTER);

    # Read CAF format file into $caf hash
    $caf = &get_config($file);

    # Build ReadType objects
    ($READTYPES_POINTER, $READTYPETABLE_POINTER) = &parse_readtypes_to_objects($caf);

    # Return output
    ($READTYPES_POINTER, $READTYPETABLE_POINTER);
}

sub parse_readtypes_to_objects {

    # Input
    my ($caf) = @_;

    # Internal
    my ($name, $line, @lines);
    my ($readtype, $readtypename, @READTYPENAMES);
    my (@sorted_indices, @indices, @values, $i);
    my ($left, $right, $distr);


    # Output
    my ($READTYPES_POINTER, $READTYPETABLE_POINTER);

    # Obtain a list of all read types
    @READTYPENAMES = keys %{$caf->{'ReadType'}};

    # Determine all fileds of all directed read types
    foreach $readtypename (@READTYPENAMES) {

	$readtype = new ReadType;

	$readtype->{'NAME'} = $readtypename;

	# Readtype length distribution
	@indices = ();
	@indices = @{$caf->{'DistributionIndex'}->{$readtypename}};
	@values = ();
	@values = @{$caf->{'DistributionValue'}->{$readtypename}};

        # Create new Distribution object
	$distr = new Distribution;

	@sorted_indices = sort { $a <=> $b } @indices;
	$left  = $sorted_indices[0];
	$right = $sorted_indices[$#sorted_indices];

	if ( $#indices >= 0 && $#values >= 0 && $#indices == $#values ) {

	    # Create new Distribution object
	    $distr = new Distribution;

	    $distr->add_field('INDICES', \@indices);
	    $distr->add_field('VALUES', \@values);
	    $distr->add_field('LEFT', $left);
	    $distr->add_field('RIGHT', $right);
	}

	$readtype->add_field('LENGTH_DISTRIBUTION', $distr);

	for ( $i = 0; $i <= $#indices; $i++ ) {
	    $readtype->{'LENGTH_DISTRIBUTION'}->{$indices[$i]} = $values[$i];
	}

	# Cost
	@lines = @{$caf->{'ReadType'}->{$readtypename}->{'Cost'}};
	foreach $line (@lines) {
	    $name = $line->[0];
	    $readtype->{'COST'} = $name;
	}

	# Name of machine
	@lines = @{$caf->{'ReadType'}->{$readtypename}->{'Machine'}};
	foreach $line (@lines) {
	    $name = $line->[0];
	    $readtype->{'MACHINE'} = $name;
	}

	# Name of dye
	@lines = @{$caf->{'ReadType'}->{$readtypename}->{'Dye'}};
	foreach $line (@lines) {
	    $name = $line->[0];
	    $readtype->{'DYE'} = $name;
	}

	# Name of chemistry
	@lines = @{$caf->{'ReadType'}->{$readtypename}->{'Chemistry'}};
	foreach $line (@lines) {
	    $name = $line->[0];
	    $readtype->{'CHEMISTRY'} = $name;
	}

	# Name of enzyme
	@lines = @{$caf->{'ReadType'}->{$readtypename}->{'Enzyme'}};
	foreach $line (@lines) {
	    $name = $line->[0];
	    $readtype->{'ENZYME'} = $name;
	}

	# Name of gel
	@lines = @{$caf->{'ReadType'}->{$readtypename}->{'Gel'}};
	foreach $line (@lines) {
	    $name = $line->[0];
	    $readtype->{'GEL_LENGTH'} = $name;
	}

	# Name of primer
	@lines = @{$caf->{'ReadType'}->{$readtypename}->{'Primer'}};
	foreach $line (@lines) {
	    $name = $line->[0];
	    $readtype->{'PRIMER'} = $name;
	}

	# Name of strand
	@lines = @{$caf->{'ReadType'}->{$readtypename}->{'Strand'}};
	foreach $line (@lines) {
	    $name = $line->[0];
	    $readtype->{'STRAND'} = $name;
	}

	# Runtime
	@lines = @{$caf->{'ReadType'}->{$readtypename}->{'Runtime'}};
	foreach $line (@lines) {
	    $name = $line->[0];
	    $readtype->{'RUNTIME'} = $name;
	}

	# Designator
	@lines = @{$caf->{'ReadType'}->{$readtypename}->{'Designator'}};
	foreach $line (@lines) {
	    $name = $line->[0];
	    $readtype->{'DESIGNATOR'} = $name;
	}

	# Enter $readtype into list of all readtypes
	push(@{$READTYPES_POINTER}, $readtype);

	# Enter $readtype into readtype recognition table
	$READTYPETABLE_POINTER->{$readtype->{'NAME'}} = $readtype;
    }

    # Return pointers to READTYPES list and READTYPETABLE hash
    ($READTYPES_POINTER, $READTYPETABLE_POINTER);
}

############################################################
# Obtain internal representation of parameters
############################################################

sub get_parameters {

    # Input
    my ($struct, $file) = @_;

    # Internal
    my ($caf, $key, $hash);

    # Output
    my ($errors, $out);

    # Read CAF format file into $caf hash if filename is defined
    if ( defined($file) ) {
	$caf = &get_config($file);
    }

    foreach $key (keys %{$caf->{'Parameter'}}) {
	$hash->{$key} = $caf->{'Parameter'}->{$key}->{'Value'}->[0]->[0];
    }

    ($errors, $out) = &parse_input_hash($struct, $hash);

    # Return hash
    ($errors, $out);
}

############################################################
# Parse calling arguments
############################################################

sub get_calling_arguments {

    # Input
    my ($struct, @arglist) = @_;

    # Internal
    my ($caf, $key, $hash, $item);

    # Output
    my ($errors, $out);

    foreach $item (@arglist) {

	# Create list of arguments for each key on command line
	if ( $item =~ /^-/ ) {            # If item is a key
	    $key = $item;
	    $hash->{$key} = 1;
	}
	else {
	    $hash->{$key} = $item;
	}
    }

    ($errors, $out) = &parse_input_hash($struct, $hash);

    # Return hash
    ($errors, $out);
}

sub parse_input_hash {

    # Input
    my ($struct, $hash) = @_;

    # Internal
    my ($key, $validelement, $min, $max);

    # Output
    my ($errors, $out);

    foreach $key (keys %{$struct}) {

	# Check if value is within  range
	if ( ($struct->{$key}->{'VALUE'} eq 'int') || ($struct->{$key}->{'VALUE'} eq 'string') ) {
	    if ( defined $struct->{$key}->{'RANGE'} ) {
		foreach $validelement (@{$struct->{$key}->{'RANGE'}}) {
		    if ( ( defined $hash->{$key} ) && $hash->{$key} == $validelement ) {
			$out->{$key} = $hash->{$key};
		    }
		}
	    }
	    else {
		$out->{$key} = $hash->{$key};
	    }
	}
	elsif ( $struct->{$key}->{'VALUE'} eq 'switch' ) {
	    if ( defined $hash->{$key} ) {
		$out->{$key} = 1;
	    }
	}
	elsif ( $struct->{$key}->{'VALUE'} eq 'float' ) {
	    if ( defined $struct->{$key}->{'RANGE'} ) {
		$min = $struct->{$key}->{'RANGE'}->[0];
		$max = $struct->{$key}->{'RANGE'}->[1];
		if ( ( defined $hash->{$key} ) && ($hash->{$key} >= $min) && ($hash->{$key} <= $max) ) {
		    $out->{$key} = $hash->{$key};
		}
	    }
	    else {
		$out->{$key} = $hash->{$key};
	    }
	}
		
	# If parameter not defined
	if ( ! defined $out->{$key} ) {

	    # If mandatory, report error
	    if ( $struct->{$key}->{'TYPE'} eq 'mandatory' ) {
		push(@{$errors}, "Mandatory argument not defined: $key ...\n");
	    }
	    # If optional, use default
	    elsif ( $struct->{$key}->{'TYPE'} eq 'optional' ) {
		$out->{$key} = $struct->{$key}->{'DEFAULT'};
	    }
	}
    }

    # Return output
    ($errors, $out);
}

############################################################
# Parse problem resolution usage file
############################################################

sub parse_caf_problemtypes {

    # Input
    my ($caf) = @_;

    # Internal
    my ($problemtypename, @problemtypenames);

    # Output
    my ($problemtypetable);

    # Obtain a list of all problem types
    @problemtypenames = keys %{$caf->{'Resolution'}};

    # Enter $problem into subclonetype recognition table
    foreach $problemtypename (@problemtypenames) {
	$problemtypetable->{$problemtypename} = 1;
    }

    # Return output
    $problemtypetable;
}

############################################################
# Obtain internal representation of problems
############################################################

sub get_problems {

    # Input
    my ($file) = @_;

    # Internal
    my ($caf);


    # Read CAF format file into CAF hash
    $caf = &get_config($file);

    $caf;
}


############################################################
# Obtain internal representation of selected directed reads
############################################################

sub get_ordered_reads {

    # Input
    my ($file) = @_;

    # Internal
    my ($caf);

    # Output
    my ($ORDEREDREADS_POINTER, $ORDEREDREADTABLE_POINTER);

    # Read CAF format file into $caf hash
    $caf = &get_config($file);

    # Build ReadType objects
    ($ORDEREDREADS_POINTER, $ORDEREDREADTABLE_POINTER) = &parse_orderedreads_to_objects($caf);

    # Return output
    ($ORDEREDREADS_POINTER, $ORDEREDREADTABLE_POINTER);
}

sub parse_orderedreads_to_objects {

    # Input
    my ($caf) = @_;

    # Internal
    my ($orderedread);
    my ($name, $line, @lines);
    my ($item, @NAMES);

    # Output
    my ($ORDEREDREADS_POINTER, $ORDEREDREADTABLE_POINTER);

    # Obtain a list of all read types
    @NAMES = keys %{$caf->{'Sequence'}};

    # Determine all fields of all directed reads
    foreach $item (@NAMES) {

	$orderedread = new DirectedRead;

	$orderedread->{'NAME'} = $item;

	# Readtype
	@lines = @{$caf->{'Sequence'}->{$item}->{'ReadType'}};
	foreach $line (@lines) {
	    $name = $line->[0];
	    $orderedread->{'READTYPE'} = $name;
	}

	# Template
	@lines = @{$caf->{'Sequence'}->{$item}->{'Template'}};
	foreach $line (@lines) {
	    $name = $line->[0];
	    $orderedread->{'TEMPLATE'} = $name;
	}

	# Template
	@lines = @{$caf->{'Sequence'}->{$item}->{'Clone'}};
	foreach $line (@lines) {
	    $name = $line->[0];
	    $orderedread->{'CLONE'} = $name;
	}


	# Enter $readtype into list of all readtypes
	push(@{$ORDEREDREADS_POINTER}, $orderedread);

	# Enter $readtype into readtype recognition table
	$ORDEREDREADTABLE_POINTER->{$orderedread->{'NAME'}} = $orderedread;
    }

    # Return pointers to READTYPES list and READTYPETABLE hash
    ($ORDEREDREADS_POINTER, $ORDEREDREADTABLE_POINTER);
}

############################################################
# Read CAF format read sequence entry
############################################################

sub get_caf_read {

    # Input
    my ($file, $readname, $sequence, $dna, $basequality, $baseposition, $rules) = @_;

    # Internal
    my (%obj_caf);
    my ($line, $new_paragraph);
    my ($class, $colon, $name);
    my ($complete, $sequence_ok, $dna_ok, $basequality_ok, $baseposition_ok);
    my ($obj);
    my (@fields, $field, $value);
    my ($relleft, $relright, $rellength, $reldir);
    my ($absleft, $absright, $abslength, $absdir);
    my ($stolen_flag, $stolen_comment, $stolen_dna);
    my ($DNA, $BASEQUALITY, $BASEPOSITION);

    # Output
    my ($read);

    # Translation between CAF and object field names
    %obj_caf = ( 
		'Name' => 'NAME',
		'State' => 'STATE',
		'DNA' => 'DNA',
		'BaseQuality' => 'BASEQUALITY',
		'BasePosition' => 'BASEPOSITION',
		'SCF_File' => 'TRACEFILE',
		'Phd_version' => 'PHDVERSION',
		'Template' => 'TEMPLATE',
		'Sequencing_vector' => 'TEMPLATETYPE',
		'Primer' => 'PRIMER',
		'Dye' => 'DYE',
		'Strand' => 'STRAND',
		'Clone' => 'CLONE',
		'Seq_vec' => 'SEQVECTOR',
		'Clone_vec' => 'CLONEVECTOR',
		'Clipping' => 'CLIPPING',
		'Align_to_SCF' => 'TRACEALIGNMENT',
		'Tag' => 'TAGS',
		'Assembled_from' => 'BASESEGMENTS',
		'Golden_path' => 'TILES'
		);
    
    # Define new Read object
    $read = new Read;

    # Add name to object
    $read->add_field('NAME', $readname);

    # Persist until contig is complete or file ends 
    $new_paragraph = 1;
    until ( $complete || eof($file)) {

	$line = <$file>;
	chomp $line;

	# Discard CAF file comments
        $line =~ s/\/\/.*//;

	# If line is empty, discard it, set new paragraph flag
	if ( $line eq "" ) {
	    $new_paragraph = 1;
	    next;
	}

	
	# If first line of paragraph is recognized as valid
	if ( $new_paragraph && ( ($class, $colon, $name) = $line =~ /^(\S+)(\s+:)?\s+(\S+)/ ) ) {

	    ###############################################################
	    # Branch on class
	    ###############################################################
	    # Sequence class
	    if ( $sequence && ($class eq 'Sequence') && ($name eq $readname) ) {

		# Register presence of Sequence paragraph
		$sequence_ok = 1;

		# Read and parse remaining lines of the paragraph
		while ( ($line = <$file>) && ($line ne "\n") ) {
		    chomp $line;

		    # Parse ace format attribute line
		    @fields = &parseAceLine($line);
		    $field = $fields[0];

		    # Special treatment : State
		    if ( $field =~ /Padded/ || $field =~ /Unpadded/ ) {
			$value = $field;
			$field = 'State';
			if ( defined $field && defined $value ) {
			    $read->add_field($obj_caf{$field}, $value);
			}
		    }

		    # Only continue if field is wanted
		    if ( (! defined $rules) || $rules->{$obj_caf{$field}} ) {

			# Special treatment :  Tag
			if ( $field eq 'Tag' ) {
			
			    shift(@fields);

			    # Create new Tag object
			    $obj = new Tag;
			    $obj->add_field('METHOD', $fields[0]);
			    $obj->add_field('LEFT', $fields[1]);
			    $obj->add_field('RIGHT', $fields[2]);
			    $obj->add_field('COMMENT', $fields[3]);
			
			    $value = $obj;
			}
			# Special treatment :  Stolen tag
			elsif ( $field eq 'Stolen' ) {
			    
			    shift(@fields);
			    
			    $stolen_flag = 1;
			    $stolen_comment = $fields[0];
			}
			# Special treatment: Vector or Clipping
			elsif ( $field eq 'Seq_vec' || $field eq 'Clone_vec' || $field eq 'Clipping' ) {
			
			    shift(@fields);
			    
			    # Create new Region object
			    $obj = new Region;
			    $obj->add_field('METHOD', $fields[0]);
			    $obj->add_field('LEFT', $fields[1]);
			    $obj->add_field('RIGHT', $fields[2]);
			    $obj->add_field('TYPE', $fields[3]);
			    
			    $value = $obj;
			}
			# Special treatment: Trace alignment
			elsif ( $field eq 'Align_to_SCF' ) {
			    
			    shift(@fields);
			    
			    # Calculate derived characteristics of a basesegment relative to contig
			    ($relleft, $relright, $rellength, $reldir) =
				&start_end_to_left_right_length_direction($fields[0], $fields[1]);
			    
			    # Calculate derived characteristics of a basesegment relative to assembled read
			    ($absleft, $absright, $abslength, $absdir) =
				&start_end_to_left_right_length_direction($fields[2], $fields[3]);
			    
			    # Create new BaseSegment object
			    $obj = new TraceAlignment;
			    
			    $obj->add_field('RELSTART', $fields[0]);
			    $obj->add_field('RELEND', $fields[1]);
			    $obj->add_field('RELLEFTEND', $relleft);
			    $obj->add_field('RELRIGHTEND', $relright);
			    $obj->add_field('RELLENGTH', $rellength);
			    $obj->add_field('RELDIRECTION', $reldir);
			    
			    $obj->add_field('ABSSTART', $fields[2]);
			    $obj->add_field('ABSEND', $fields[3]);
			    $obj->add_field('ABSLEFTEND', $absleft);
			    $obj->add_field('ABSRIGHTEND', $absright);
			    $obj->add_field('ABSLENGTH', $abslength);
			    $obj->add_field('ABSDIRECTION', $absdir);
			    
			    $value = $obj;
			}
			else {
			    shift(@fields);			
			    $value = $fields[0];
			}
		    
			# Register $value in $field of $read if rules are not defined or rules permit
			if ( ((! defined $rules) || $rules->{$obj_caf{$field}}) && defined $field && defined $value ) {
			    $read->add_field($obj_caf{$field}, $value);
			}
		    }
		}
	    }
	
	    # DNA class
	    elsif ($dna && ($class eq 'DNA') && ($name eq $readname) ) {

		# Register presence of DNA paragraph
		$dna_ok = 1;

		# If needed, parse it
		if ( (! defined $rules) || $rules->{$obj_caf{'DNA'}} ) {

		    # Initialise DNA
		    undef $DNA;
		    
		    # Read and parse remaining lines of the paragraph
		    while ( ($line = <$file>) && ($line ne "\n") ) {
			chomp $line;
			
			# Discard white space
			$line =~ s/\s//g;
			
			# Concatenate to DNA
			$DNA .= $line;		
		    } 
		    
		    # Register DNA sequence
		    $read->add_field($obj_caf{'DNA'}, $DNA);
		}
		# If not needed, just go to end of paragraph
		else {
		    while ( ($line = <$file>) && ($line ne "\n") ) {
		    } 
		}
		$new_paragraph = 1;
            } 
	    
	    # BaseQuality class
	    elsif ($basequality && ($class eq 'BaseQuality') && ($name eq $readname) ) {

		# Register presence of BaseQuality  paragraph
		$basequality_ok = 1;

		# If needed, parse it
		if ( (! defined $rules) || $rules->{$obj_caf{'BaseQuality'}} ) {

		    # Read and parse remaining lines of the paragraph
		    while ( ($line = <$file>) && ($line ne "\n") ) {
			chomp $line;	
	    
			$BASEQUALITY = [split(/\s+/, $line)];

			# Register BaseQuality sequence if rules are not defined or rules permit
			$read->add_field($obj_caf{'BaseQuality'}, $BASEQUALITY);

			undef $BASEQUALITY;
		    } 
		}

		# If not needed, just go to end of paragraph
		else {
		    while ( ($line = <$file>) && ($line ne "\n") ) {
		    } 
		}
		$new_paragraph = 1;
	    } 
	    
	    # BasePosition class
	    elsif ($baseposition && ($class eq 'BasePosition') && ($name eq $readname) ) {

		# Register presence of BasePosition  paragraph
		$baseposition_ok = 1;

		# If needed, parse it
		if ( (! defined $rules) || $rules->{$obj_caf{'BasePosition'}} ) {

		    # Read and parse remaining lines of the paragraph
		    while ( ($line = <$file>) && ($line ne "\n") ) {
			chomp $line;	
	    
			$BASEPOSITION = [split(/\s+/, $line)];

			# Register BaseQuality sequence if rules are not defined or rules permit
			$read->add_field($obj_caf{'BasePosition'}, $BASEPOSITION);

			undef $BASEPOSITION;
		    } 
		}

		# If not needed, just go to end of paragraph
		else {
		    while ( ($line = <$file>) && ($line ne "\n") ) {
		    } 
		}
		$new_paragraph = 1;
	    } 
	    else {
		undef $new_paragraph;
	    }
	}

	# Check if reading is complete
	$complete = (! $sequence || $sequence_ok) && (! $dna || $dna_ok) && ( ! $basequality || $basequality_ok) && ( ! $baseposition || $baseposition_ok);
    }

    # Return $read if complete
    if ( $complete ) {

	# If there was a stolen key, create Tag object
	if ( defined $stolen_flag ) {

	    $stolen_dna = $read->get_field('DNA');

	    # Create new Tag object
	    $obj = new Tag;
	    $obj->add_field('METHOD', 'STOL');
	    $obj->add_field('COMMENT', $stolen_comment);
	    
	    $value = $obj;
	    $obj->add_field('LEFT', '1');
	    $obj->add_field('RIGHT', length($stolen_dna) );

	    $field = 'Tag';
	    $value = $obj;

	    # Register $value in $field of $read if rules are not defined or rules permit
	    if ( ((! defined $rules) || $rules->{$obj_caf{$field}}) && defined $field && defined $value ) {
		$read->add_field($obj_caf{$field}, $value);
	    }
	}

	$read;
    }
    else {
	undef;
    }
}

############################################################
# Read CAF format contig sequence entry
############################################################

sub get_caf_contig {

    # Input
    my ($file, $contigname, $sequence, $dna, $basequality, $rules) = @_;

    # Internal
    my (%obj_caf);
    my ($line, $new_paragraph);
    my ($class, $colon, $name);
    my ($complete, $sequence_ok, $dna_ok, $basequality_ok);
    my ($obj);
    my ($ar, %ASSEMBLEDREAD);
    my (@fields, $field, $value);
    my ($bs, $bsname);
    my ($tile);
    my ($bsrelleft, $bsrelright, $bsrellength, $bsreldir, $bsrelstart, $bsrelend);
    my ($bsabsleft, $bsabsright, $bsabslength, $bsabsdir, $bsabsstart, $bsabsend);
    my ($arrelleft, $arrelright, $arrellength, $arreldir, $arrelstart, $arrelend);
    my ($arabsleft, $arabsright, $arabslength, $arabsdir, $arabsstart, $arabsend);
    my ($firstbasesegment);
    my ($DNA, $BASEQUALITY);

    my(@LINES);

    # Output
    my ($contig);

    # Translation between CAF and object field names
    %obj_caf = ( 
		'Name' => 'NAME',
		'State' => 'STATE',
		'DNA' => 'DNA',
		'BaseQuality' => 'BASEQUALITY',
		'Tag' => 'TAGS',
		'Assembled_from' => 'BASESEGMENTS',
		'Golden_path' => 'TILES'
		);
    
    # Define new Contig object
    $contig = new Contig;

    # Add name to object
    $contig->add_field('NAME', $contigname);

    # Set firstbasesegment flag
    $firstbasesegment = 1;

    # Persist until contig is complete or file ends 
    $new_paragraph = 1;
    until ( $complete || eof($file)) {

	$line = <$file>;
	chomp $line;

	# Discard CAF file comments
        $line =~ s/\/\/.*//;

	# If line is empty, discard it, set new paragraph flag
	if ( $line eq "" ) {
	    $new_paragraph = 1;
	    next;
	}
	
	# If first line of paragraph is recognized as valid
	if ( $new_paragraph && ( ($class, $colon, $name) = $line =~ /^(\S+)(\s+:)?\s+(\S+)/ ) ) {

	    ###############################################################
	    # Branch on class
	    ###############################################################
	    # Sequence class
	    if ( $sequence && ($class eq 'Sequence') && ($name eq $contigname) ) {

		# Register presence of Sequence paragraph
		$sequence_ok = 1;

		# Read and parse remaining lines of the paragraph
		while ( ($line = <$file>) && ($line ne "\n") ) {
		    chomp $line;

		    # Parse ace format attribute line
		    @fields = &parseAceLine($line);
		    $field = $fields[0];

		    # Special treatment : State
		    if ( $field =~ /Padded/ || $field =~ /Unpadded/ ) {
			$value = $field;
			$field = 'State';
			if ( defined $field && defined $value ) {
			    $contig->add_field($obj_caf{$field}, $value);
			}
		    }

		    # Only continue if field is wanted
		    if ( (! defined $rules) || $rules->{$obj_caf{$field}} ) {

			# Special treatment: Golden path
			if ( $field eq 'Golden_path' ) {
			    
			    shift(@fields);
			    
			    $tile = [ [$fields[0], $fields[1],$fields[2] ] ];
			    $contig->add_field('TILES', $tile);
			    
			    undef $value;
			}
			# Special treatment: BaseSegment
			elsif ( $field eq 'Assembled_from' ) {
			    
			    shift(@fields);
			    
			    # Get basesegment coordinates
			    $bsname = $fields[0];
			    $bsrelstart = $fields[1];
			    $bsrelend = $fields[2];
			    $bsabsstart = $fields[3];
			    $bsabsend = $fields[4];
			    
			    # Calculate derived characteristics of a basesegment relative to contig
			    ($bsrelleft, $bsrelright, $bsrellength, $bsreldir) =
				&start_end_to_left_right_length_direction($bsrelstart, $bsrelend);
			    
			    # Calculate derived characteristics of a basesegment relative to assembled read
			    ($bsabsleft, $bsabsright, $bsabslength, $bsabsdir) =
				&start_end_to_left_right_length_direction($bsabsstart, $bsabsend);
			    

			    
			    # If basesegment new, create new AssembledRead
			    if ( ! defined $ASSEMBLEDREAD{$bsname} ) {
				
				# Define new AssembledRead object
				$ar = new AssembledRead;	
				$ar->add_field('NAME', $bsname);
				$ASSEMBLEDREAD{$bsname} = $ar;
				$contig->add_field('ASSEMBLEDREADS', $ar);
	
				# Initialize numerical properties
				$ar->add_field('RELSTART',$bsrelstart);
				$ar->add_field('RELEND', $bsrelend);
				$ar->add_field('RELLEFTEND', $bsrelleft);
				$ar->add_field('RELRIGHTEND', $bsrelright);
				$ar->add_field('RELDIRECTION', $bsreldir);
				$ar->add_field('RELLENGTH', $bsrellength);
				
				$ar->add_field('ABSSTART', $bsabsstart);
				$ar->add_field('ABSEND', $bsabsend);
				$ar->add_field('ABSLEFTEND', $bsabsleft);
				$ar->add_field('ABSRIGHTEND', $bsabsright);
				$ar->add_field('ABSDIRECTION', $bsabsdir);
				$ar->add_field('ABSLENGTH', $bsabslength);
			    }
			    else {
				$ar = $ASSEMBLEDREAD{$bsname};

				# Calculate new numerical properties in light of new Basesegment
				$arrelleft = &min($bsrelleft, $ar->get_field('RELLEFTEND'));
				$arrelright = &max($bsrelright, $ar->get_field('RELRIGHTEND'));
				$arreldir = $ar->get_field('RELDIRECTION');
				$arabsleft = &min($bsabsleft, $ar->get_field('ABSLEFTEND'));
				$arabsright = &max($bsabsright, $ar->get_field('ABSRIGHTEND'));
				$arabsdir = $ar->get_field('ABSDIRECTION');
				($arrelstart, $arrelend, $arrellength) =
				    &left_right_direction_to_start_end_length($arrelleft, $arrelright, $arreldir);
				($arabsstart, $arabsend, $arabslength) =
				    &left_right_direction_to_start_end_length($arabsleft, $arabsright, $arabsdir);
		    
				
				# Update numerical properties
				$ar->add_field('RELSTART',$arrelstart);
				$ar->add_field('RELEND', $arrelend);
				$ar->add_field('RELLEFTEND', $arrelleft);
				$ar->add_field('RELRIGHTEND', $arrelright);
				$ar->add_field('RELDIRECTION', $arreldir);
				$ar->add_field('RELLENGTH', $arrellength);
				
				$ar->add_field('ABSSTART', $arabsstart);
				$ar->add_field('ABSEND', $arabsend);
				$ar->add_field('ABSLEFTEND', $arabsleft);
				$ar->add_field('ABSRIGHTEND', $arabsright);
				$ar->add_field('ABSDIRECTION', $arabsdir);
				$ar->add_field('ABSLENGTH', $arabslength);
			    }
			    
			    # Add basesegment alignment to AssembledRead
			    $bs = [ [$bsrelstart, $bsrelend, $bsabsstart, $bsabsend] ];
			    $ar->add_field('BASESEGMENTS', $bs);

			    # Adjust left end and right end of contig
			    if ( $firstbasesegment ) {
				$contig->add_field('LEFTEND', $bsrelleft);
				$contig->add_field('RIGHTEND', $bsrelright);
				
				undef $firstbasesegment;
			    }
			    else {
				$contig->add_field('LEFTEND', &min ($contig->get_field('LEFTEND'), $bsrelleft) );
				$contig->add_field('RIGHTEND', &max ($contig->get_field('RIGHTEND'), $bsrelright) );
				$contig->add_field('LENGTH', $contig->get_field('RIGHTEND') - $contig->get_field('LEFTEND') + 1);
			    }
			    
			    undef $value;
			}
			# Special treatment :  Tag
			elsif ( $field eq 'Tag' ) {
			    
			    shift(@fields);
			    
			    # Create new Tag object
			    $obj = new Tag;
			    $obj->add_field('METHOD', $fields[0]);
			    $obj->add_field('LEFT', $fields[1]);
			    $obj->add_field('RIGHT', $fields[2]);
			    $obj->add_field('COMMENT', $fields[3]);
			    
			    $value = $obj;
			}
			else {
			    shift(@fields);			
			    $value = $fields[0];
			}

			# Register $value in $field of $read if rules are not defined or rules permit
			if ( defined $field && defined $value ) {
			    $contig->add_field($obj_caf{$field}, $value);
			}
		    }
		}
		$new_paragraph = 1;
	    }
	
	    # DNA class
	    elsif ($dna && ($class eq 'DNA') && ($name eq $contigname) ) {
		
		# Register presence of DNA paragraph
		$dna_ok = 1;
		
		# If needed, parse it
		if ( (! defined $rules) || $rules->{$obj_caf{'DNA'}} ) {

		    # Initialise DNA
		    undef $DNA;
		    
		    # Read and parse remaining lines of the paragraph
		    while ( ($line = <$file>) && ($line ne "\n") ) {
			chomp $line;
			
			# Discard white space
			$line =~ s/\s//g;
			
			# Concatenate to DNA
			$DNA .= $line;		
		    } 
		    
		    # Register DNA sequence
		    $contig->add_field($obj_caf{'DNA'}, $DNA);
		}
		# If not needed, just go to end of paragraph
		else {
		    while ( ($line = <$file>) && ($line ne "\n") ) {
		    } 
		}
		$new_paragraph = 1;
            } 
	    
	    # BaseQuality class
	    elsif ($basequality && ($class eq 'BaseQuality') && ($name eq $contigname) ) {
		
		# Register presence of BaseQuality  paragraph
		$basequality_ok = 1;
		
		# If needed, parse it
		if ( (! defined $rules) || $rules->{$obj_caf{'BaseQuality'}} ) {

		    # Initialise BaseQuality
		    $BASEQUALITY = [];

		    # Read and parse remaining lines of the paragraph
		    while ( ($line = <$file>) && ($line ne "\n") ) {
			chomp $line;	
	    
			$BASEQUALITY = [split(/\s+/, $line)];

			# Register BaseQuality sequence if rules are not defined or rules permit
			$contig->add_field($obj_caf{'BaseQuality'}, $BASEQUALITY);

			undef $BASEQUALITY;
		    } 
		}

		# If not needed, just go to end of paragraph
		else {
		    while ( ($line = <$file>) && ($line ne "\n") ) {
		    } 
		}
		$new_paragraph = 1;
	    } 
	    else {
		undef $new_paragraph;
	    }
	}

	# Check if contig is complete
	$complete = (! $sequence || $sequence_ok) && (! $dna || $dna_ok) && ( ! $basequality || $basequality_ok);
    }

    # Check if contig coordinates are defined. If not, make attempt to extract them
    if ( (! defined $contig->get_field('LEFTEND')) ||  (! defined $contig->get_field('RIGHTEND')) || (! defined $contig->get_field('LENGTH')) ) {
	if ( defined $contig->get_field('DNA') ) {
	    $contig->add_field('LEFTEND', '1');
	    $contig->add_field('RIGHTEND', length($contig->get_field('DNA')) );
	    $contig->add_field('LENGTH', $contig->get_field('RIGHTEND') - $contig->get_field('LEFTEND') + 1);
	}
    }

    # Return $contig if complete
    if ( $complete ) {
	$contig;
    }
    else {
	undef;
    }
}

############################################################
# Determine order of sequence entries in CAF file
############################################################

sub get_caf_order {

    # Input
    my ($file) = @_;

    # Internal
    my ($line, $new_paragraph);
    my ($class, $colon, $name);
    my (%PRESENT, %SUBCLASS);

    # Output
    my ($names, $reads, $contigs, $classes);

    # Initialize list
    $names = [];

    # Determine order of all relevant read classes 
    $new_paragraph = 1;
    while ( $line = <$file> ) {

	chomp $line;

	# Discard CAF file comments
        $line =~ s/\/\/.*//;

	# If line is empty, discard it, set new paragraph flag
	if ( $line eq "" ) {
	    $new_paragraph = 1;
	    next;
	}

	# If first line of paragraph is recognized as valid
	if ( $new_paragraph ) {
	    if ( ($class, $colon, $name) = $line =~ /^(\S+)(\s+:)?\s+(\S+)/ ) {

		###############################################################
		# Branch on class
		# Only recognise: Sequence, DNA, BaseQuality
		###############################################################
		# Sequence class: Register name and search for subclass specification
		if ( $class eq 'Sequence' ) {

		    # Register presence of Sequence
		    $classes->{$name}->{'Sequence'} = 1;
		
		    # Register presence of valid read
		    if ( ! $PRESENT{$name} ) {
			push(@{$names}, $name);
			$PRESENT{$name} = 1;
		    }

		    # Search for subclass specification
		    do {
			$line = <$file>;
			chomp $line;

			if ( $line =~ /^Is_contig/ ) {
			    $SUBCLASS{$name} = 'Contig';
			}
			elsif ( $line =~ /^Is_read/ ) {
			    $SUBCLASS{$name} = 'Read';
			}
		    } until ( $SUBCLASS{$name} || $line eq "" || eof($file) );
		    $new_paragraph = 1;
		}	

		# DNA class
		elsif ( $class eq 'DNA' ) {
		    
		    # Register presence of DNA paragraph
		    $classes->{$name}->{'DNA'} = 1;
		    
		    # Register presence of valid read
		    if ( ! $PRESENT{$name} ) {
			push(@{$names}, $name);
			$PRESENT{$name} = 1;
		    }	    

		    # Read rest of paragraph
		    do {
			$line = <$file>;
			chomp $line;
		    } until ( $line eq "" || eof($file) );
		    $new_paragraph = 1;
		} 
		# BaseQuality class
		elsif ( $class eq 'BaseQuality' ) {
		    
		    # Register presence of BaseQuality  paragraph
		    $classes->{$name}->{'BaseQuality'} = 1;
		    
		    # Register presence of valid read
		    if ( ! $PRESENT{$name} ) {
			push(@{$names}, $name);
			$PRESENT{$name} = 1;
		    }	    

		    # Read rest of paragraph
		    do {
			$line = <$file>;
			chomp $line;
		    } until ( $line eq "" || eof($file) );
		    $new_paragraph = 1;
		}
		# BasePosition class
		elsif ( $class eq 'BasePosition' ) {
		    
		    # Register presence of BaseQuality  paragraph
		    $classes->{$name}->{'BasePosition'} = 1;
		    
		    # Register presence of valid read
		    if ( ! $PRESENT{$name} ) {
			push(@{$names}, $name);
			$PRESENT{$name} = 1;
		    }	    

		    # Read rest of paragraph
		    do {
			$line = <$file>;
			chomp $line;
		    } until ( $line eq "" || eof($file) );
		    $new_paragraph = 1;
		}
	    }
	    else {
		undef $new_paragraph;
	    }
	}
    }

    # Separate out reads and contigs
    if ( $#{$names} >= 0 ) {
	foreach $name (@{$names}) {
	    if ( $SUBCLASS{$name} eq 'Read' ) {
		push(@{$reads}, $name);
	    }
	    elsif ( $SUBCLASS{$name} eq 'Contig' ) {
		push(@{$contigs}, $name);
	    }
	}
    }

    # Return output
    if ( $#{$names} >= 0 ) {
	($names, $reads, $contigs, $classes);
    }
    else {
	(undef, undef, undef, undef);
    }
}

############################################################
# Write perl raed object into CAF file
############################################################

sub write_read_into_caf {

    # Input
    my ($file, $line_length, $read) = @_;

    # Internal
    my (%caf_obj);
    my ($field, @fields, $obj, $obj);
    my ($i, @basequal, @basepos);

    # Translation hash from object fields to CAF fields
    %caf_obj = ( 
		'NAME' => 'Name',
		'STATE' => 'State',
		'DNA' => 'DNA',
		'BASEQUALITY' => 'BaseQuality',
		'BASEPOSITION' => 'BasePosition',
		'TRACEFILE' => 'SCF_File',
		'PHDVERSION' => 'Phd_version',
		'TEMPLATE' => 'Template',
		'TEMPLATETYPE' => 'Sequencing_vector',
		'PRIMER' => 'Primer',
		'DYE' => 'Dye',
		'STRAND' => 'Strand',
		'CLONE' => 'Clone',
		'SEQVECTOR' => 'Seq_vec',
		'CLONEVECTOR' => 'Clone_vec',
		'CLIPPING' => 'Clipping',
		'TRACEALIGNMENT' => 'Align_to_SCF',
		'TAGS' => 'Tag',
		);
    
    # These sequence objects are only valid if they have a name
    if ( defined $read->get_field('NAME') ) {
	print $file 'Sequence', " : ", $read->get_field('NAME'), "\n";
	print $file 'Is_read', "\n";
	if ( defined $read->get_field('STATE') ) {
	    print $file $read->get_field('STATE'), "\n";
	}

	# Get all used fields of $read
	@fields = $read->get_used_fieldnames;
		
	# Print all fields of sub sequence object
	foreach $field (@fields) {
	    # Special treatment: Tag
	    if ( $field eq 'TAGS' ) {
		if ( defined $read->get_field($field) ) {
		    foreach $obj (@{$read->get_field($field)}) {
			if ( defined $obj->get_field('METHOD')
			    && defined $obj->get_field('LEFT')
			    && defined $obj->get_field('RIGHT') ) {
			    print $file "$caf_obj{$field} ", 
			    $obj->get_field('METHOD'), " ",
			    $obj->get_field('LEFT'), " ",
			    $obj->get_field('RIGHT'), " ",
			    "\"", $obj->get_field('COMMENT'), "\"", "\n";
			}
		    }
		}
	    }
	    # Special treatment: Region objects (Sequencing vector, Cloning vector and Clipping)
	    elsif ( $field eq 'SEQVECTOR' || $field eq 'CLONEVECTOR' || $field eq 'CLIPPING' ) {
		foreach $obj (@{$read->get_field($field)}) {
		    if ( defined $obj->get_field('METHOD')
			&& defined $obj->get_field('LEFT')
			&& defined $obj->get_field('RIGHT') ) {
			print $file "$caf_obj{$field} ", 
			$obj->get_field('METHOD'), " ",
			$obj->get_field('LEFT'), " ",
			$obj->get_field('RIGHT'), " ",
			"\"", $obj->get_field('TYPE'), "\"", "\n";
		    }
		}
	    }
	    # Special treatment: Alignment to original base call
	    elsif ( $field eq 'TRACEALIGNMENT' ) {
		foreach $obj (@{$read->get_field($field)}) {
		    if ( defined $obj->get_field('RELSTART')
			&& defined $obj->get_field('RELEND')
			&& defined $obj->get_field('ABSSTART')
			&& defined $obj->get_field('ABSEND') ) {
			print $file "$caf_obj{$field} ", 
			$obj->get_field('RELSTART'), " ",
			$obj->get_field('RELEND'), " ",
			$obj->get_field('ABSSTART'), " ",
			$obj->get_field('ABSEND'), "\n";
		    }
		}
	    }
	    # The rest of the objects to be printed out
	    elsif ( defined $caf_obj{$field} && $field ne 'NAME' && $field ne 'STATE' 
		   && $field ne 'DNA' && $field ne 'BASEQUALITY' && $field ne 'BASEPOSITION' ) {
		print $file "$caf_obj{$field} ", $read->get_field($field), "\n";
	    }
	}
	print $file "\n";
	
	# Special treatment: DNA (needs to be in separate paragraph)
	if ( defined $read->get_field('DNA') ) {
	    print $file "DNA : ", $read->get_field('NAME'), "\n";
	    for ( $i = 0; $i < length($read->get_field('DNA')); $i += $line_length ) {
		print $file substr( $read->get_field('DNA'), $i, $line_length), "\n";
	    }
	    print $file "\n";
	}
	
	# Special treatment: BASEQUALITY (needs to be in separate paragraph)
	if ( defined $read->get_field('BASEQUALITY') ) {
	    print $file "BaseQuality : ", $read->get_field('NAME'), "\n";
	    @basequal = @{$read->get_field('BASEQUALITY')};
	    while ( $#basequal >= 0 ) {
		print $file join(' ', splice(@basequal, 0, $line_length)) , "\n";
	    }
	    print $file "\n";
	}
	
	# Special treatment: BASEPOSITION (needs to be in separate paragraph)
	if ( defined $read->get_field('BASEPOSITION') ) {
	    print $file "BasePosition : ", $read->get_field('NAME'), "\n";
	    @basepos = @{$read->get_field('BASEPOSITION')};
	    while ( $#basepos >= 0 ) {
		print $file join(' ', splice(@basepos, 0, $line_length)) , "\n";
	    }
	    print $file "\n";
	}
    }
}

############################################################
# Write perl contig object into CAF file
############################################################

sub write_contig_into_caf {

    # Input
    my ($file, $line_length, $contig) = @_;

    # Internal
    my (%caf_obj);
    my ($field, @fields, $obj, $bs, $tile);
    my ($i, @basequal);
    my ($name, $relstart, $relend);
    my ($absstart, $absend);

    # Translation hash from object fields to CAF fields
    %caf_obj = ( 
		'NAME' => 'Name',
		'STATE' => 'State',
		'DNA' => 'DNA',
		'BASEQUALITY' => 'BaseQuality',
		'TAGS' => 'Tag',
		'ASSEMBLEDREADS' => 'Assembled_from',
		'TILES' => 'Golden_path',
		);
    
    # These sequence objects are only valid if they have a name
    if ( defined $contig->get_field('NAME') ) {
	print $file 'Sequence', " : ", $contig->get_field('NAME'), "\n";
	print $file 'Is_contig', "\n";
	if ( defined $contig->get_field('STATE') ) {
	    print $file $contig->get_field('STATE'), "\n";
	}

	# Get all used fields of $read
	@fields = $contig->get_used_fieldnames;
		
	# Print all fields of sub sequence object
	foreach $field (@fields) {
	    my($value) = $contig->get_field($field);

	    # Special treatment: Tag
	    if ( $field eq 'TAGS' ) {
		if ( defined $contig->get_field($field) ) {
		    foreach $obj (@{$contig->get_field($field)}) {
			if ( defined $obj->get_field('METHOD')
			    && defined $obj->get_field('LEFT')
			    && defined $obj->get_field('RIGHT') ) {
			    print $file "$caf_obj{$field} ", 
			    $obj->get_field('METHOD'), " ",
			    $obj->get_field('LEFT'), " ",
			    $obj->get_field('RIGHT'), " ",
			    "\"", $obj->get_field('COMMENT'), "\"", "\n";
			}
		    }
		}
	    }
	    # Special treatment: Basesegments
	    elsif ( $field eq 'ASSEMBLEDREADS' ) {
		if ( defined $contig->get_field($field) ) {
		    foreach $obj (@{$contig->get_field($field)}) {
			$name = $obj->get_field('NAME');
			if ( defined $obj->get_field('BASESEGMENTS') ) {
			    foreach $bs (@{$obj->get_field('BASESEGMENTS')}) {
				print $file $caf_obj{$field}, " ", $name, " ", $bs->[0], " ", $bs->[1], " ", $bs->[2], " ",  $bs->[3], "\n";
			    }
			}
		    }
		}
	    }
	    # Special treatment: Golden path
	    elsif ( $field eq 'TILES' ) {
		if ( defined $contig->get_field($field) ) {
		    foreach $tile (@{$contig->get_field($field)}) {
			print $file $caf_obj{$field}, " ", $tile->[0], " ", $tile->[1], " ",  $tile->[2], "\n";
		    }
		}
	    }
	    # The rest of the objects to be printed out
	    elsif ( defined $caf_obj{$field} && $field ne 'NAME' && $field ne 'STATE' 
		   && $field ne 'DNA' && $field ne 'BASEQUALITY' ) {
		print $file "$caf_obj{$field} ", $contig->get_field($field), "\n";
	    }
	}
	print $file "\n";
	
	# Special treatment: DNA (needs to be in separate paragraph)
	if ( defined $contig->get_field('DNA') ) {
	    print $file "DNA : ", $contig->get_field('NAME'), "\n";
	    for ( $i = 0; $i < length($contig->get_field('DNA')); $i += $line_length ) {
		print $file substr( $contig->get_field('DNA'), $i, $line_length), "\n";
	    }
	    print $file "\n";
	}
	
	# Special treatment: BASEQUALITY (needs to be in separate paragraph)
	if ( defined $contig->get_field('BASEQUALITY') ) {
	    print $file "BaseQuality : ", $contig->get_field('NAME'), "\n";
	    @basequal = @{$contig->get_field('BASEQUALITY')};
	    while ( $#basequal >= 0 ) {
		print $file join(' ', splice(@basequal, 0, $line_length)) , "\n";
	    }
	    print $file "\n";
	}
    }
}

############################################################
# Parsing of cross_match alignment program output
############################################################

sub parse_crossmatch_output {

    # Parse arguments

    my ($file) = @_;

    # Internal
    my ($matchindex, $errorindex, $line, @FIELDS);

    # Output variables

    my (@matches);

    # Open cross_match output file

    open(FILE, $file) || die "couldnt open $file: $!\n";

    # Process cross_match output file

    $matchindex = -1;
    $errorindex = 0;

    while ($line = <FILE>) {
        chop($line);
        if ( ( $line =~ /\(/ ) && ( $line =~ /\)/ ) ) {

            $line =~ s/^\s+//g;         # Remove leading whitespace at beginning of line
            $line =~ s/[\(\)]/ /g;      # Remove parentheses
            $line =~ s/[\*]/ /g;        # Remove wild card characters

            @FIELDS = split(' ', $line);

            # If $line is likely to indicate an alignment

            if ( ( $#FIELDS >= 11 ) && ( $#FIELDS <= 13 ) ) {
            
                $matchindex++;
                $errorindex = 0;

                if ( $FIELDS[8] eq "C" ) {
                    $matches[$matchindex]->{'QUERY'}->{'NAME'} = $FIELDS[4];
                    $matches[$matchindex]->{'QUERY'}->{'START'} = $FIELDS[5];
                    $matches[$matchindex]->{'QUERY'}->{'END'} = $FIELDS[6];

                    $matches[$matchindex]->{'SUBJECT'}->{'NAME'} = $FIELDS[9];
                    $matches[$matchindex]->{'SUBJECT'}->{'START'} = $FIELDS[11];
                    $matches[$matchindex]->{'SUBJECT'}->{'END'} = $FIELDS[12];
                }
                else {
                    $matches[$matchindex]->{'QUERY'}->{'NAME'} = $FIELDS[4];
                    $matches[$matchindex]->{'QUERY'}->{'START'} = $FIELDS[5];
                    $matches[$matchindex]->{'QUERY'}->{'END'} = $FIELDS[6];

                    $matches[$matchindex]->{'SUBJECT'}->{'NAME'} = $FIELDS[8];
                    $matches[$matchindex]->{'SUBJECT'}->{'START'} = $FIELDS[9];
                    $matches[$matchindex]->{'SUBJECT'}->{'END'} = $FIELDS[10];
                }

                # Calculate LEFTEND, RIGHTEND and DIRECTION of QUERY and SUBJECT FIELDS

                if ( $matches[$matchindex]->{'QUERY'}->{'END'} >= $matches[$matchindex]->{'QUERY'}->{'START'} ) {
                    $matches[$matchindex]->{'QUERY'}->{'LEFTEND'} = $matches[$matchindex]->{'QUERY'}->{'START'};
                    $matches[$matchindex]->{'QUERY'}->{'RIGHTEND'} = $matches[$matchindex]->{'QUERY'}->{'END'};
                    $matches[$matchindex]->{'QUERY'}->{'DIRECTION'} = 1;
                }
                else {
                    $matches[$matchindex]->{'QUERY'}->{'LEFTEND'} = $matches[$matchindex]->{'QUERY'}->{'END'};
                    $matches[$matchindex]->{'QUERY'}->{'RIGHTEND'} = $matches[$matchindex]->{'QUERY'}->{'START'};
                    $matches[$matchindex]->{'QUERY'}->{'DIRECTION'} = -1;
                }

                if ( $matches[$matchindex]->{'SUBJECT'}->{'END'} >= $matches[$matchindex]->{'SUBJECT'}->{'START'} ) {
                    $matches[$matchindex]->{'SUBJECT'}->{'LEFTEND'} = $matches[$matchindex]->{'SUBJECT'}->{'START'};                    
		    $matches[$matchindex]->{'SUBJECT'}->{'RIGHTEND'} = $matches[$matchindex]->{'SUBJECT'}->{'END'};
                    $matches[$matchindex]->{'SUBJECT'}->{'DIRECTION'} = 1;
                }
                else {
                    $matches[$matchindex]->{'SUBJECT'}->{'LEFTEND'} = $matches[$matchindex]->{'SUBJECT'}->{'END'};
                    $matches[$matchindex]->{'SUBJECT'}->{'RIGHTEND'} = $matches[$matchindex]->{'SUBJECT'}->{'START'};
                    $matches[$matchindex]->{'SUBJECT'}->{'DIRECTION'} = -1;
                }

                # Calculate DIRECTION of $match

                $matches[$matchindex]->{'DIRECTION'} = $matches[$matchindex]->{'QUERY'}->{'DIRECTION'}
                * $matches[$matchindex]->{'SUBJECT'}->{'DIRECTION'};
            }

            # If $line is likely to indicate an error corresponding to an alignment

            if ( ($#FIELDS == 5) && ($FIELDS[1] =~ /^\d+$/) && ($FIELDS[4] =~ /^\d+$/) && (($FIELDS[0] =~ /^S/) || ($FIELDS[0] =~ /^D/) || ($FIELDS[0] =~ /^I/)) ) {

		# Check type & multiplicity of error
		if ( $FIELDS[0] =~ /(^\D)\-(\d+)/ ) {
		    $matches[$matchindex]->{'ERRORS'}->[$errorindex]->{'TYPE'} = $1;
		    $matches[$matchindex]->{'ERRORS'}->[$errorindex]->{'MULTIPLICITY'} = $2;
		}
		else {
		    $matches[$matchindex]->{'ERRORS'}->[$errorindex]->{'TYPE'} = $FIELDS[0];
		    $matches[$matchindex]->{'ERRORS'}->[$errorindex]->{'MULTIPLICITY'} = 1;
		}
		
                $matches[$matchindex]->{'ERRORS'}->[$errorindex]->{'QUERYPOSITION'} = $FIELDS[1];
                $matches[$matchindex]->{'ERRORS'}->[$errorindex]->{'NUCLEOTIDE'} = $FIELDS[2];
                $matches[$matchindex]->{'ERRORS'}->[$errorindex]->{'QUALITY'} = $FIELDS[3];
                $matches[$matchindex]->{'ERRORS'}->[$errorindex]->{'SUBJECTPOSITION'} = $FIELDS[4];
                $matches[$matchindex]->{'ERRORS'}->[$errorindex]->{'CONTEXT'} = $FIELDS[5];

                $errorindex++;
            }
        }
    }

    # Output
    @matches;
}
