#!/usr/bin/perl

use strict;

# kmer length 
my $K = 10;

# counters
my $c = 0; my $ct = 1; my $contig_number = 1;

my @sequences = ();

# read sequences from file given in first argument to command line

open(IN,$ARGV[0]);
# noting the original complete sequence (provided) 
my $full_seq = <IN>;
# cheating slightly by making the first read equal to the first 100 bases of the original
# sequence, just to make it easier to align the assembled and original sequences following 
# assembly 
$sequences[0] = substr($full_seq,0,100);

while(<IN>) {
$sequences[$ct] = <IN>;
$ct++;
}


print ">original\n$full_seq\n";

# assign verteces (i.e. all kmers that occur in the reads) to an array
my @verteces = &assign_verteces(\@sequences,$K);

# assign edges to a hash array 
my %out_edges = &assign_edges(\@sequences,$K);

# make contig starting at vertex 0 (could start at any vertex, but because of the
# cheat above this will ensure that the assembled sequence is alignable with original
&make_contig(0,\%out_edges,\@verteces); # construct a contig starting from vertex 0

for($ct = 0; $ct <= $#verteces; $ct++) {
	if($#{$out_edges{$verteces[$ct]}} > -1) { # out degree of ct^th vertex
	&make_contig($ct,\%out_edges,\@verteces);
	$ct = 0;
	}
} 



# subroutines

sub make_contig { 
my $vertex = $_[0]; 
my $ref = $_[1]; 
my %out_edges = %{$ref};
my $ref = $_[2]; 
my @verteces = @{$ref};
my  @cycle = ();
my  @new_cycle = ();

@cycle = &make_cycle($verteces[$vertex],\%out_edges);

my $path_growth = 1;
while($path_growth > 0) {
my $size0 = $#cycle;
	for($c = 0; $c <= $#cycle; $c++) {
		if($#{$out_edges{$cycle[$c]}} > -1) { # out degree of c^th vertex in cycle
		@new_cycle = &make_cycle($cycle[$c],\%out_edges);
		@cycle = (@cycle[0..$c],@new_cycle,@cycle[$c+1..$#cycle]);
		}
	}

$path_growth = $#cycle - $size0;
	if($path_growth > 0) {
	$c = 0;
	}
}

my $contig = &reconstruct_sequence(@cycle);
print ">contig $contig_number\n$contig\n";
$contig_number++;
}

sub reconstruct_sequence {
my @path = @_;
my $sequence = $path[0];
my $i = 0;
my $kmer_len = length($path[0]);
	for($i = 1; $i <= $#path; $i++) {
	$sequence = $sequence . substr($path[$i],$kmer_len-1,1);
	}
return($sequence);
}

sub make_cycle {
my $start_vertex = $_[0];
my $ref = $_[1];
my %out_edges = %{$ref};
my $ct = 0;
my $vertex = $start_vertex;
my @path = ($vertex);
	WHILE: while($vertex ne $start_vertex || $ct == 0) {
	$ct++;
		unless($vertex = pop(@{$out_edges{$vertex}})) {
		last WHILE;
		}
	push(@path,$vertex);
	}
return(@path);
}

sub assign_verteces {
my @sequences = @{$_[0]};
my $K = $_[1];;
my $sub = "";
my $seq = "";
my $i = 0;
my $ct = 0;
my @verteces = ();
my %seen = ();
	foreach $seq (@sequences) {
	chomp($seq);
		for($i = 0; $i <= length($seq) - $K; $i++) {
		$sub = substr($seq,$i,$K);
			if($seen{$sub} == 0) {
			$seen{$sub} = 1;
			push(@verteces,$sub);
			}
		}
	}
return(@verteces);
}

sub assign_edges {
my @sequences = @{$_[0]};
my $K = $_[1];;
my $seq = "";
my $sub1 = "";
my $sub2 = "";
my %edge_defined = ();
my $i = 0;
	foreach $seq (@sequences) {
		for($i = 0; $i < length($seq) - ($K+1); $i++) {
		$sub1 = substr($seq,$i,$K);
		$sub2 = substr($seq,$i+1,$K);
			if(!exists($edge_defined{$sub1}{$sub2})) {
			$edge_defined{$sub1}{$sub2} = 1;
			push(@{$out_edges{$sub1}},$sub2);
			}
		}
	}
return(%out_edges);
}