#!usr/bin/perl

#STEPS:
#	1. Get the sequences in 18bp sequence. Inner 18, then outer 18. They will overlap by one. You can get this from scaff_pos 
#		- 35 to - 18, then -17 to SNP, ref base. 
#	2. Put the sequences in a hash. 
#	3. Run the genome on 18bp sliding window. 1-18, 2-19, 3-20 etc.... It should go to the end of the scaffold. 

#################################################################################################################################################################
#1. Creating the database of probe kmers

	open (INFILE, $ARGV[0]);					# PROBES. NO IUPAC. REFERENCE SEQUENCE!!!

	while ($line = <INFILE>){
		$line =~ s/\s+$//;
		
		@split = split "\t", $line;
		@split_left = split "", $split[0];
	
		$split_left_1 = $split_left[0] . $split_left[1] .
			$split_left[2] . $split_left[3] .
			$split_left[4] . $split_left[5] .
			$split_left[6] . $split_left[7] .
			$split_left[8] . $split_left[9] .
			$split_left[10] . $split_left[11] .
			$split_left[12] . $split_left[13] .
			$split_left[14] . $split_left[15] .
			$split_left[16] . $split_left[17];

		$split_left_1 =~ tr/ATGCN/atgcn/;
	
		$split_left_2 = $split_left[18] . $split_left[19] .
                        $split_left[20] . $split_left[21] .
                        $split_left[22] . $split_left[23] .
                        $split_left[24] . $split_left[25] .
                        $split_left[26] . $split_left[27] .
                        $split_left[28] . $split_left[29] .
                        $split_left[30] . $split_left[31] .
                        $split_left[32] . $split_left[33] .
                        $split_left[34] . $split[1];

                $split_left_2 =~ tr/ATGCN/atgcn/;

		@split_right = split "", $split[2];

		 $split_right_1 = $split[1] .
			$split_right[0] . $split_right[1] .
                        $split_right[2] . $split_right[3] .
                        $split_right[4] . $split_right[5] .
                        $split_right[6] . $split_right[7] .
                        $split_right[8] . $split_right[9] .
                        $split_right[10] . $split_right[11] .
                        $split_right[12] . $split_right[13] .
                        $split_right[14] . $split_right[15] .
                        $split_right[16];

                $split_right_1 =~ tr/ATGCN/atgcn/;

		 $split_right_2 = $split_right[17] . 
			$split_right[18] . $split_right[19] .
                        $split_right[20] . $split_right[21] .
                        $split_right[22] . $split_right[23] .
                        $split_right[24] . $split_right[25] .
                        $split_right[26] . $split_right[27] .
                        $split_right[28] . $split_right[29] .
                        $split_right[30] . $split_right[31] .
                        $split_right[32] . $split_right[33] .
                        $split_right[34];

                $split_right_2 =~ tr/ATGCN/atgcn/;
		
		$Identifier_1 = $split[3] . "_" . $split[4] . "_" . "1";
		$Identifier_2 = $split[3] . "_" . $split[4] . "_" . "2";
		$Identifier_3 = $split[3] . "_" . $split[4] . "_" . "3";
		$Identifier_4 = $split[3] . "_" . $split[4] . "_" . "4";

		$kmer_hash{$split_left_1} = $Identifier_1;
		$kmer_hash{$split_left_2} = $Identifier_2;
		$kmer_hash{$split_right_1} = $Identifier_3;
		$kmer_hash{$split_right_2} = $Identifier_4;

		$INPUT = $INPUT + 4;

		$full_sequence = $split[0] . $split[1] . $split[2];
		$full_sequence =~ tr/ATGCN/atgcn/;

		$full_hash{$full_sequence} = $full_sequence;

			if ($full_sequence =~ /n/){
				#print $full_sequence, "\n";
			}
		$full_seq++;

	}

	$INPUT = $INPUT - 4;										# Because the first line is not adding info 

	print "Number of sequences added to has:\t", $INPUT, "\n";
	print "Scalar of input hash:\t", scalar keys %kmer_hash, "\n";
	print "Full sequences:\t", $full_seq, "\n";
	print "Scalar of all sequence:\t", scalar keys %full_hash, "\n";


	close INFILE;

#################################################################################################################################################################
#2. Querying the genome

open (FASTA_TABLE, "/UCHC/LABS/Wegrzyn/ConiferGenomes/Pita/Freebayes_test/genome/pita.v2.0.1.masked.3k.table")
        or die "CANNOT OPEN FASTA TABLE!!!";
        while ($line = <FASTA_TABLE>){
        $line =~ /\>(.+)/;
        @fasta_table = split "\t", $1;
        @fasta_sequence = split "", $fasta_table[1];
	
	for ($a = 1; $a < scalar @fasta_sequence - 18; $a++){
		$sequence = $fasta_sequence[$a] . 
				$fasta_sequence[$a+1] .
				$fasta_sequence[$a+2] .
				$fasta_sequence[$a+3] .
				$fasta_sequence[$a+4] .
				$fasta_sequence[$a+5] .
				$fasta_sequence[$a+6] .
				$fasta_sequence[$a+7] .
				$fasta_sequence[$a+8] .
				$fasta_sequence[$a+9] .
				$fasta_sequence[$a+10] .
				$fasta_sequence[$a+11] .
				$fasta_sequence[$a+12] .	
				$fasta_sequence[$a+13] .
				$fasta_sequence[$a+14] .
				$fasta_sequence[$a+15] .                              
                                $fasta_sequence[$a+16] .
                                $fasta_sequence[$a+17];
	$sequence =~ tr/ATCGN/atcgn/;						# Making everything lower case

		if ($kmer_hash{$sequence} =~ /./){
			if ($sequence =~ /nnnnnnnnnnnnnnnnnn/){
				next;
			}
			else{ 
				$occurance{$kmer_hash{$sequence}}++;
				#print $sequence, "\t", $occurance{$kmer_hash{$sequence}}, "\n";
			}
		}
	}


}
close FASTA_TABLE;
#################################################################################################################################################################
# 3. Evauluate

open (OUTPUT, ">overlap_table.txt");

	foreach $key (%kmer_hash){
		print OUTPUT $occurance{$kmer_hash{$key}}, "\t", $key, "\n";
	}

close OUTPUT;





