#!usr/bin/perl
#################################################################################################################################
# PART ONE : Gather relevant information about strict and standard vcf file SNPs.

# How to run this script:
#	perl 70_mer_probes.pl [strict] [standard]
#		
# Outfiles created:
#	scoring_temporary_outfile.txt
#	scored_SNPs.txt
#	SNP_probes.txt	


open (STRICT, $ARGV[0])
        or die "Cannot open the strict file!!!";
while ($line = <STRICT>){								#Logic: Store all the SNPs that are of STRICT quality
	$line =~ s/\s+$//;
        if ($line =~ /^\#/){
                next;
        }
        else{
		@line = split "\t", $line;
                $unique_ID = "$line[0]_$line[1]";
		$SNP_ID_to_strict{$unique_ID} = "STRICT";
	}
}
close STRICT;


open (STANDARD, $ARGV[1])
	or die "Cannot open the standard file!!!";
while ($line = <STANDARD>){
	$line =~ s/\s+$//;
	if ($line =~ /^\#/){
		next;
	}
	else{
		@SNPeff_info = split /\|/, $line;					#Grab all the SNPeff information 
		@line = split "\t", $line;
		$unique_ID = "$line[0]_$line[1]";
		$scaffold_position_to_SNPeff{$unique_ID} = $SNPeff_info[1];		#Scaffold_position : SNPeff info
		
		push @unique_IDs, $unique_ID;						#Store all SNPs
		$scaffold_position_to_scaffold{$unique_ID} = $line[0];			#Scaffold_position : scaffold
		$scaffold_position_to_ref{$unique_ID} = $line[3];			#Scaffold_position : ref base 
		$scaffold_position_to_alt{$unique_ID} = $line[4];			#Scaffold_position : alt base
	}
}
close STANDARD;

open (FASTA_FAI_REFERENCE, "/UCHC/LABS/Wegrzyn/ConiferGenomes/Pita/Freebayes_test/genome/pita.v2.0.1.masked.3k.fasta.fai")
	or die "Cannot open the reference file!!!";
	while ($reference = <FASTA_FAI_REFERENCE>){
                if ($reference =~ /^\w+/){
			@reference_info = split "\t", $reference;
		        $scaffold_length{$reference_info[0]} = $reference_info[1];	#Scaffold : scaffold length
                }
        }
close FASTA_FAI_REFERENCE;

open (OUTFILE, ">scoring_temporary_outfile.txt")
        or die "Cannot create the strict temporary outfile";

for ($a=0; $a < scalar @unique_IDs; $a++){
	        $unique_IDs[$a] =~ /(.+)_(.+$)/;
                        $current_SNP_scaffold = $1;					#Current line = current SNP
                        $current_SNP_position = $2;
                $unique_IDs[$a+1] =~ /(.+)_(.+$)/;					#Next line = next SNP 
                        $next_SNP_scaffold = $1;
                        $next_SNP_position = $2;
                $unique_IDs[$a-1] =~ /(.+)_(.+$)/;					#Previous line = previous SNP
                        $previous_SNP_scaffold = $1;
                        $previous_SNP_position = $2;

	$separation_R = $next_SNP_position - $current_SNP_position;			#Positions (may be weird if first or last SNP on scaffold)
        $separation_L = $current_SNP_position - $previous_SNP_position;

print OUTFILE $scaffold_position_to_scaffold{$unique_IDs[$a]}, "\t",
                $current_SNP_position, "\t",
                $scaffold_length{$current_SNP_scaffold}, "\t",
                $previous_SNP_scaffold, "\t",
                $previous_SNP_position, "\t",
                $next_SNP_scaffold, "\t",
                $next_SNP_position, "\t",
                $separation_L, "\t",
                $separation_R, "\t",
                $scaffold_position_to_ref{$unique_IDs[$a]}, "\t",
                $scaffold_position_to_alt{$unique_IDs[$a]}, "\t",
                ".", "\t",
                $scaffold_position_to_SNPeff{$unique_IDs[$a]}, "\t",
		$SNP_ID_to_strict{$unique_IDs[$a]}, "\n";
}

close OUTFILE;

open (TEMP_INFILE, "scoring_temporary_outfile.txt")
        or die "Cannot open the strict temporary outfile!!!";
open (SCORES_OUTFILE, ">scored_SNPs.txt")
        or die "Cannot create scoring file!!!";

print SCORES_OUTFILE "Scaffold", "\t",
                "Position", "\t",
                "Length_of_Scaffold", "\t",
                "Previous_SNP_Scaffold", "\t",
                "Previous_SNP_Position", "\t",
                "Next_SNP_Scaffold", "\t",
                "Next_SNP_Position", "\t",
                "Bp_Left", "\t",
                "Bp_Right", "\t",
                "Ref", "\t",
                "Alt", "\t",
                "SNP", "\t",
		"IUPAC_Code", "\t";
                "Grade", "\t",
                "SNP_Location", "\t",
		"Score";

print SCORES_OUTFILE "\n";

while ($annotation = <TEMP_INFILE>){							# THE IF STATEMENT HELL
	$annotation =~ s/\s+$//;
        @annotations = split "\t", $annotation;						# annotation = each line of scoring file
	if ($annotation =~ /^\w+\d+\t/){						# If line has something on it
		if ($annotations[5] ne $annotations[0]){				# If next SNP is on new scaffold...
			 $distance_to_end = $annotations[2] - $annotations[1];		# 	Right distance is distance to end of scaffold
			 $annotations[8] =~ s/.+/$distance_to_end/;			#	So replace what was a wrong number!
		}
		if ($annotations[3] ne $annotations[0]){				# Likewise, if SNP is first on scaffold...
			$position = $annotations[1];					# 	Left space is position of SNP on scaffold.
			$annotations[7] =~ s/.+/$position/;				# 	So replace what was a wrong number!
		}
		if ($annotations[7] >= 36 and $annotations[8] >= 36){			# Numberic scoring. L and R spacing > 35 to pass
                        $annotations[11] =~ s/./PASS/;					# Update with PASS
                }
		elsif ($annotations[7] <=35 or $annotations[8] <=35){
                        $annotations[11] =~ s/./FAIL/;					# Update with fail
                }
			

#IUPAC / SNP organizer section								#What is happening:
		if ($annotations[9] =~ /a/i and $annotations[10] =~ /g/i){			# If ref is a and alt is g
			$SNP = "A/G";								# SNP is [a/g]
			if ($annotations[13] =~ /STRICT/){$IUPAC = "R";}			# If this SNP is strict, IUPAC (R) is CAPITALIZED
			else {$IUPAC = "r";}							# If it is NOT strict, IUPAC (r) is lower case.
		}											# I do this for printing. Probe filtering favors
                if ($annotations[9] =~ /g/i and $annotations[10] =~ /a/i){				# SNPs in flanking regions that are low quality.
			$SNP = "A/G";									# This is how I tell them apart!
                        if ($annotations[13] =~ /STRICT/){$IUPAC = "R";}
                        else {$IUPAC = "r";}
		}
                if ($annotations[9] =~ /c/i and $annotations[10] =~ /t/i){
			$SNP = "C/T";
			if ($annotations[13] =~ /STRICT/){$IUPAC = "Y";}
                        else {$IUPAC = "y";}
		}
                if ($annotations[9] =~ /t/i and $annotations[10] =~ /c/i){
			$SNP = "C/T";
			if ($annotations[13] =~ /STRICT/){$IUPAC = "Y";}
                        else {$IUPAC = "y";}
		}
                if ($annotations[9] =~ /g/i and $annotations[10] =~ /c/i){
			$SNP = "G/C";
			if ($annotations[13] =~ /STRICT/){$IUPAC = "S";}
                        else {$IUPAC = "s";}
                }
                if ($annotations[9] =~ /c/i and $annotations[10] =~ /g/i){
			$SNP = "G/C";
			if ($annotations[13] =~ /STRICT/){$IUPAC = "S";}
                        else {$IUPAC = "s";}
		}
                if ($annotations[9] =~ /a/i and $annotations[10] =~ /t/i){
			$SNP = "A/T";
			if ($annotations[13] =~ /STRICT/){$IUPAC = "W";}
                        else {$IUPAC = "w";}
		}
                if ($annotations[9] =~ /t/i and $annotations[10] =~ /a/i){
			$SNP = "A/T";
			if ($annotations[13] =~ /STRICT/){$IUPAC = "W";}
                        else {$IUPAC = "w";}
		}
                if ($annotations[9] =~ /g/i and $annotations[10] =~ /t/i){
			$SNP = "G/T";
			if ($annotations[13] =~ /STRICT/){$IUPAC = "K";}
                        else {$IUPAC = "k";}
		}
                if ($annotations[9] =~ /t/i and $annotations[10] =~ /g/i){
			$SNP = "G/T";
			if ($annotations[13] =~ /STRICT/){$IUPAC = "K";}
                        else {$IUPAC = "k";}
		}
                if ($annotations[9] =~ /a/i and $annotations[10] =~ /c/i){
			$SNP = "A/C";
			if ($annotations[13] =~ /STRICT/){$IUPAC = "M";}
                        else {$IUPAC = "m";}
		}
                if ($annotations[9] =~ /c/i and $annotations[10] =~ /a/i){
			$SNP = "A/C";
			if ($annotations[13] =~ /STRICT/){$IUPAC = "M";}
                        else {$IUPAC = "m";}
		}


		
# Dealing with Indels
		@ref_array = split "", $annotations[9];					#Split up ref and alt bases. Indels have at least 2 bases:
		@alt_array = split "", $annotations[10];				# A and AT or AT and A
			if (scalar @ref_array > 1){					# If 2+ bases...
				if ($annotations[9] ne /\,/){				# But not a comma
					$SNP = "$annotations[9]\/$annotations[10]";		#SNP will be [AT/A];				
					if ($annotations[13] =~ /STRICT/){$IUPAC = "I";}	#IUPAC will be I for indel if strict
					else {$IUPAC = "i";}					# i if standard
				}
			}		
			if (scalar @alt_array > 1){                                     # Same thing but this would be deletions
                               if ($annotations[9] ne /\,/){                           
                                       $SNP = "$annotations[9]\/$annotations[10]";
                                       if ($annotations[13] =~ /STRICT/){$IUPAC = "I";}
                                        else {$IUPAC = "i";}
                                }
	                        }
		if ($annotations[9] =~ /\,/ or $annotations[10] =~ /\,/){               # A comma = multiallelic SNP
                        $SNP = "X";
                        $IUPAC = "X";
                        $annotations[13] =~ s/STRICT/STANDARD/;                         # Demotes any high quality mutlialleic indel to Standard.
                }  	
		

		if ($annotations[13] =~ /STRICT/){
			$score = "STRICT";
		}
		else{
			$score = "STANDARD";
		}
			
		$Unique_ID = "$annotations[0]_$annotations[1]";
		$Position_to_IUPAC{$Unique_ID} = $IUPAC;

		print SCORES_OUTFILE $annotations[0], "\t",
                $annotations[1], "\t",
                $annotations[2], "\t",
                $annotations[3], "\t",
                $annotations[4], "\t",
                $annotations[5], "\t",
                $annotations[6], "\t",
                $annotations[7], "\t",
                $annotations[8], "\t",
                $annotations[9], "\t",
                $annotations[10], "\t",
                $SNP, "\t",
		$IUPAC, "\t",
                $annotations[11], "\t",
                $annotations[12], "\t",
		$score, "\n";
        }
}

close TEMP_INFILE;
close SCORES_OUTFILE;

#die "Remove me to move onto probe creation";

########################################################################################################################################
# Part 2 : Using the scoring information to build probes

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{$fasta_table[0]} = $fasta_table[1];			# scaffold : scaffold sequence
	
}
close FASTA_TABLE;

open (SCORED_SNPS, "scored_SNPs.txt") 
	or die "CANNOT OPEN SCORED SNPS FILE!!!";
open (PROBES_OUTFILE, ">SNP_probes.txt")
        or die "CANNOT CREATE PROBES OUTPUT!!!";

print PROBES_OUTFILE "Sequence left", "\t", 	
	"SNP", "\t", 
	"Sequence right", "\t", 
	"Scaffold", "\t", 
	"Position", "\n";

while ($scores = <SCORED_SNPS>){
	$scores =~ s/\s+$//;
	if ($scores =~ /STRICT/){
        	@pass = split "\t", $scores;				
        	@sequence = split "", $fasta{$pass[0]};			# Break up the sequence of the matching scaffold.

                        $more = $pass[1] + 34;
                        $start = $pass[1];
                        $less = $pass[1] - 36;
                        $stop = $pass[1]- 2;
                        $current = $pass[1] -1;

			# TAKING INDELS INTO CONSIDERATION

			@split_for_indel = split /\//, $pass[11];					# Split up the XXXXX/XX array
			@split_ref = split "", $split_for_indel[0];					# How big is the first part (>1 = indel)

			if (scalar @split_ref > 1){
				$spaces_to_add = scalar @split_ref - 1;					#1234567/12 means I must add 6 spaces to back end start
				$start = $start + $spaces_to_add;
				$more = $more + $spaces_to_add;
			}

                        for ($a = $less; $a <= $stop; $a++){						# Left base probe printing
				$position = $a + 1;
				$Unique_ID = "$pass[0]_$position";
				
				if ($Position_to_IUPAC{$Unique_ID}){					
					print PROBES_OUTFILE $Position_to_IUPAC{$Unique_ID};		# If the position matches a IUPAC made before, print it!
				}
				else{
                                	print PROBES_OUTFILE $sequence[$a];
                        	}
			}

                        print PROBES_OUTFILE "\t", $pass[11], "\t";					# This is the actual SNP.

                        for ($a = $start; $a <= $more; $a++){						# Same thing, but now right end of probe.
                                $position = $a + 1;
                                $Unique_ID = "$pass[0]_$position";

                        	if ($Position_to_IUPAC{$Unique_ID}){
                                        print PROBES_OUTFILE $Position_to_IUPAC{$Unique_ID};
                                }
                                else{
                                        print PROBES_OUTFILE $sequence[$a];
                                }

			}

                        print PROBES_OUTFILE "\t", $pass[0], 						# extra information
				"\t", $pass[1],
				"\t", $pass[12], 
				"\t", $pass[13], 
				"\t", $pass[14], 
				"\t", $pass[15],
				"\n";

                        $more = 0;
                        $less = 0;
                        $start = 0;
                        $stop = 0;

	}
}


close PROBES_OUTFILE;
close SCORED_SNPS;

