require "crispr_blast_modules.pm";
require "rna_struct.pm";
require "svm_classifier.pm";

my $h_constr = 0;
my $hc_start = 1;
my $hc_end = 3;
my $mfc_cutoff = 20.0;
my $p_cutoff = 0.5;
my $e_cutoff = 0.05;
my $cd_cutoff = 0;
my $around_cas9 = 0;
my $rr_unmask_start = 0;
my $rr_unmask_end = 0;
my $tt_skip = 25;
my $tail_cut = 204;
my $rpt_revcom = 1;
my $tail_pat_file = "pat_chh.txt";
my $rr_start_seq = "GYY";
my $junct_start_seq = "AR";
my $path = '';
my $outpath = '';
my $folder = "tracrRNA_search";
my $cas9_not_required = 0;
my $min_hybrid_mismatch = 3;
my $max_hybrid_mismatch = 15;
my $min_hybrid_helix_size = 9;
my $mask_annotated_genes = 0;

my $range = 10000;
my $ind = 0;
foreach(@ARGV) {
	if (@ARGV[$ind] eq '-in') {
		$path = @ARGV[$ind + 1];
		
		#if ($path !~ /.gbff/) {
		#	die "File must be a GBFF file";
		#}

		if (! (-e $path)) {
			die "cannot open file: " . $path;
		}
	}
	
	if (@ARGV[$ind] eq '-out') {
		$outpath = @ARGV[$ind + 1];
		if (! (-d $outpath)) {
			mkdir($outpath);
		}
	}
	
	if (@ARGV[$ind] eq '-helix_constr') {
		my $val = @ARGV[$ind + 1];
		if ($val =~ /\d+,\d+/) {
			($hc_start, $hc_end) = ($val =~ /(\d+),(\d+)/);
			if ($hc_start > $hc_end) {
				my $temp = $hc_start;
				$hc_start = $hc_end;
				$hc_end = $temp;
			}
			
			if ($hc_start < 1) {
				$hc_start = 1;
			}
			
			if ($hc_end < 1) {
				$hc_end = 1;
			}
			
			if ($hc_start < $hc_end) {
				$h_constr = 1;
			}
		} else {	
			if ($hc_start < $hc_end) {
				$h_constr = 1;
			}
		}
	}
	
	if (@ARGV[$ind] eq '-mfe_cutoff') {
		my $val = @ARGV[$ind + 1];
		
		if ($val =~ /-*\d+[\.]*[\d]*/) {
			$mfc_cutoff = $val;
			if ($mfc_cutoff < 0) {
				$mfc_cutoff = -1 * $mfc_cutoff;
			}
		}
	}
	
	if (@ARGV[$ind] eq '-p_cutoff') {
		my $val = @ARGV[$ind + 1];
		
		if ($val =~ /-*\d+[\.]*[\d]*/) {
			$p_cutoff = $val;
			if ($p_cutoff < 0) {
				$p_cutoff = 0;
			}
			
			if ($p > 1) {
				$p = 1;
			}
		}
	}
	
	if (@ARGV[$ind] eq '-e_cutoff') {
		my $val = @ARGV[$ind + 1];
		
		if ($val =~ /-*\d+[\.]*[\d]*/) {
			$e_cutoff = $val;
			if ($e_cutoff < 0) {
				$e_cutoff = 0;
			}
		}
	}
	
	if (@ARGV[$ind] eq '-cd_cutoff') {
		my $val = @ARGV[$ind + 1];
		
		if ($val =~ /-*\d+[\.]*[\d]*/) {
			$cd_cutoff = $val;
		}
	}
	
	if (@ARGV[$ind] eq '-tail_cut') {
		my $val = @ARGV[$ind + 1];
		
		if ($val =~ /\d+/) {
			$tail_cut = $val;
		}
	}
	
	if (@ARGV[$ind] eq '-tt_skip') {
		my $val = @ARGV[$ind + 1];
		
		if ($val =~ /\d+/) {
			$tt_skip = $val;
		}
	}
	
	if (@ARGV[$ind] eq '-rr_start_seq') {
		my $val = @ARGV[$ind + 1];
		
		if ($val =~ /[ACGTRYSWKMBDHVN]+/) {
			$rr_start_seq = $val;
		}
	}
	
	if (@ARGV[$ind] eq '-junct_start_seq') {
		my $val = @ARGV[$ind + 1];
		
		if ($val =~ /[ACGTRYSWKMBDHVN]+/) {
			$junct_start_seq = $val;
		}
	}
	
	if (@ARGV[$ind] eq '-rr_unmask_start') {
		my $val = @ARGV[$ind + 1];
		
		if ($val =~ /-*\d+/) {
			$rr_unmask_start = $val;
		}
	}
	
	if (@ARGV[$ind] eq '-rr_unmask_end') {
		my $val = @ARGV[$ind + 1];
		
		if ($val =~ /-*\d+/) {
			$rr_unmask_end = $val;
		}
	}
	
	if (@ARGV[$ind] eq '-around_cas9') {
		my $val = @ARGV[$ind + 1];
		
		if ($val =~ /\d+/) {
			$range = $val;
			if ($range < 100) {
				$range = 100;
			}
			$around_cas9 = 1;
		}
	}
	
	if (@ARGV[$ind] eq '-repeat_as_is') {
		$rpt_revcom = 0;
	}
	
	if (@ARGV[$ind] eq '-pat') {
		$tail_pat_file = @ARGV[$ind + 1];
	}
	
	if (@ARGV[$ind] eq '-cas9_not_required') {
		$cas9_not_required = 1;
	}
	
	if (@ARGV[$ind] eq '-min_hybrid_mismatch') {
		my $val = @ARGV[$ind + 1];
		
		if ($val =~ /-*\d+/) {
			$min_hybrid_mismatch = $val;
		}
	}
	
	if (@ARGV[$ind] eq '-max_hybrid_mismatch') {
		my $val = @ARGV[$ind + 1];
		
		if ($val =~ /-*\d+/) {
			$max_hybrid_mismatch = $val;
		}
	}
	
	if (@ARGV[$ind] eq '-min_hybrid_helix_size') {
		my $val = @ARGV[$ind + 1];
		
		if ($val =~ /-*\d+/) {
			$min_hybrid_helix_size = $val;
		}
	}
	
	if (@ARGV[$ind] eq '-mask_annotated_genes') {
		$mask_annotated_genes = 1;
	}
	
	$ind++;
}

my $flag = 0;
my $count = 0;
if (! (-d $folder)) {
	$flag = system("mkdir $folder");
}

my $fname_ext = '';
if ($path =~ m|[/\\]|) {
	my ($dir, $fname_ext_0) = ($path =~ m|^(.*[/\\])([^/\\]+?)$|);
	$fname_ext = $fname_ext_0;
} else {
	$fname_ext = $path;
}

my $dir_name = '';
my $session = get_unique_id();

$flag = system("perl sfm_run.pl $tail_pat_file test_seq_in.fa test_seq_out$session.txt >sfmerr$session.txt 2>&1");

open(SFMERR, "sfmerr$session.txt");
my @errs = <SFMERR>;
close(SFMERR);

if (scalar(@errs) > 0) {
	$tail_pat_file = "pat_chh.txt";
}

if (-e "test_seq_out$session.txt") {
	unlink("test_seq_out$session.txt");
}

if (-e "sfmerr$session.txt") {
	unlink("sfmerr$session.txt");
}

$count = process_gbff_psiblast($path, $session);
my $n_contig = $count;

my $hmmflag = run_hmmsearch($session);
	
if (($count == 0) && ($hmmflag == 0) && ($cas9_not_required == 0)) {
	if (-d "psiblast_files$session") {
		$flag = system("rm -rf psiblast_files$session");
	}

	if (-d "work_01$session") {
		$flag = system("rm -rf work_01$session");
	}
	exit;
}

$fname_ext =~ s/\./_/g;
$fname_ext =~ s/_gbff//g;

$dir_name = "$folder/" . $fname_ext;

if (-d $dir_name) {
	$flag = system("rm -rf " . $dir_name);
}

$flag = system('mkdir ' . $dir_name);

my @dir_w = glob("psiblast_files$session/*.*");
	
for (my $i = 0; $i < scalar(@dir_w); $i++) {
	my $path_w = @dir[0];
	$flag = system("mv " . @dir_w[$i] . " " . $dir_name);
}
	
$flag = system("rm -rf psiblast_files$session");

my @gbff_file_paths = glob($dir_name . '/*.gbff');
my @gbff_names = ();
for (my $i = 0; $i < scalar(@gbff_file_paths); $i++) {
	my $line = @gbff_file_paths[$i];
	my ($dir, $fname_ext) = ($line =~ m|^(.*[/\\])([^/\\]+?)$|);
	my @arr = split(/\./, $fname_ext);
	my $fname = @arr[0];
	push(@gbff_names, $fname);
	
	print "run CRISPRDetect on contig " . $fname . "\n";
	$flag = crispr_call($dir_name, $fname, $cd_cutoff, $session);
	
	print "process CRISPRDetect output for contig " . $fname . "\n";
	$flag = handle_crispr_repeats_data_txt($dir_name, $fname, $rpt_revcom);
}

open(CASPSI, ">>$dir_name/cas9_genome_psi_summary.csv");

for (my $k = 0; $k < scalar(@gbff_names); $k++) {
	my @coords = get_cas9_range_psiblast($dir_name, @gbff_names[$k]);
	
	my $fap = "$dir_name/" . @gbff_names[$k] . ".fa";
	open (FAGBFF, $fap);
	my @contents = <FAGBFF>;
	close(FAGBFF);
	
	my $seq = @contents[1];
	$seq =~ s/\n//g;
	my $seqlen = length($seq);
	
	print CASPSI @gbff_names[$k] . "," . $n_contig . "," . $seqlen . ",";
	if (scalar(@coords) > 0) {
		if (scalar(@{$coords[0]} > 0)) {
			for (my $j = 0; $j < 6; $j++) {
				if ($j == 5) {
					print CASPSI $coords[$j][0] . "\n";
				} else {
					print CASPSI $coords[$j][0] . ",";
				}
			}
		} else {
			print CASPSI "NA,NA,NA,NA,NA,NA\n";
		}
	} else {
		print CASPSI "NA,NA,NA,NA,NA,NA\n";
	}
}
close(CASPSI);

open(CASHMM, ">>$dir_name/cas9_genome_hmm_summary.csv");

for (my $k = 0; $k < scalar(@gbff_names); $k++) {
	my @coords = get_cas9_range_hmm($dir_name, @gbff_names[$k]);
	
	my $fap = "$dir_name/" . @gbff_names[$k] . ".fa";
	open (FAGBFF, $fap);
	my @contents = <FAGBFF>;
	close(FAGBFF);
	
	my $seq = @contents[1];
	$seq =~ s/\n//g;
	my $seqlen = length($seq);
	
	print CASHMM @gbff_names[$k] . "," . $n_contig . "," . $seqlen . ",";
	if (scalar(@coords) > 0) {
		if (scalar(@{$coords[0]} > 0)) {
			for (my $j = 0; $j < 6; $j++) {
				if ($j == 5) {
					print CASHMM $coords[$j][0] . "\n";
				} else {
					print CASHMM $coords[$j][0] . ",";
				}
			}
		} else {
			print CASHMM "NA,NA,NA,NA,NA,NA\n";
		}
	} else {
		print CASHMM "NA,NA,NA,NA,NA,NA\n";
	}
}
close(CASHMM);

open(CASPSI_R, "$dir_name/cas9_genome_psi_summary.csv");
my @caspsi_contents = <CASPSI_R>;
close(CASPSI_R);

open(CAS, ">>$dir_name/cas9_genome_summary.csv");
for (my $k = 0; $k < scalar(@caspsi_contents); $k++) {
	my $line = @caspsi_contents[$k];
	$line =~ s/\n//g;
	
	my @toks = split(/\,/, $line);
	my $start = @toks[3];
	my $contig = @toks[0];
	
	if ($start eq "NA") {
		my $cmd = "grep -w $contig $dir_name/cas9_genome_hmm_summary.csv";
		my @rs = `$cmd`;
		
		if (scalar(@rs) > 0) {
			my $rline = @rs[0];
			$rline =~ s/\n//g;
			
			my @rtoks = split(/\,/, $rline);
			my $rstart = @rtoks[3];
			
			if ($rstart eq "NA") {
				print CAS "$line\n";
			} else {
				print CAS "$rline\n";
			}
		} else {
			print CAS "$line\n";
		}
	} else {
		print CAS "$line\n";
	}
}
close(CAS);

open(LOOKUP, ">>$dir_name/blast_lookup.csv");
for (my $i = 0; $i < scalar(@gbff_names); $i++) {
	my $gname = @gbff_names[$i];
	my @repeat_summary_paths = glob($dir_name . '/*repeat_summary.csv');
	for (my $j = 0; $j < scalar(@repeat_summary_paths); $j++) {
		my $rfile = @repeat_summary_paths[$j];
		open(RPT, $rfile);
		my @rpt_contents = <RPT>;
		close(RPT);
		
		my ($dir, $fname_ext) = ($rfile =~ m|^(.*[/\\])([^/\\]+?)$|);
		my ($rcname) = ($fname_ext =~ /(.*)_repeat_summary\.csv/);
		
		my $org = get_organism($dir_name, $gname);
		
		for (my $k = 0; $k < scalar(@rpt_contents); $k++) {
			my $line = @rpt_contents[$k];
			$line =~ s/\n//g;
			my @toks = split(/\,/, $line);
			my $rptid = @toks[0];
			
			my $bfname = $gname . "_" . $rptid . "_BLAST"; 
			
			print LOOKUP "$org,$gname,$rcname,$line,$dir_name/$gname.fa,$dir_name/$bfname.txt\n";
		}
	}
}
close(LOOKUP);

$session = get_unique_id();
$flag = blast_call($dir_name, $e_cutoff, $around_cas9, $range, $rr_unmask_start, $rr_unmask_end, $mask_annotated_genes);
$flag = run_struct_calc_dir($dir_name, $mfc_cutoff, $h_constr, $hc_start, $hc_end, $tt_skip, $tail_cut, $session);
$flag = add_Cas9_data($dir_name);
$flag = add_distance_data($dir_name);
$flag = add_additional_data($dir_name, $rr_start_seq, $junct_start_seq);
$flag = parse_tail_fold($dir_name);
$flag = filterCSV($dir_name, $rr_start_seq, $junct_start_seq, $min_hybrid_mismatch, $max_hybrid_mismatch, $min_hybrid_helix_size);
$flag = svm_predict($dir_name);
$flag = search_tail_pattern($dir_name, $tail_pat_file);
$flag = add_gene_overlapping_flag($dir_name);
$flag = orderCSV($dir_name);
$flag = add_cas_classification($dir_name);
$flag = generate_tracr_txt($dir_name, $p_cutoff, $tail_pat_file, $session);
$flag = generate_tracr_gff($dir_name, $p_cutoff);

my @files = glob($dir_name . "/*.*");
		
for (my $j = 0; $j < scalar(@files); $j++) {
	my $file = @files[$j];
			
	if (index($file, "genome_tracr_") >= 0) {
		next;
	}
			
	if (index($file, "_repeat.gff") >= 0) {
		next;
	}
			
	if (index($file, "_repeat.txt") >= 0) {
		next;
	}
	
	if (index($file, "cas9_genome") >= 0) {
		next;
	}
	
	if (index($file, "repeat_summary") >= 0) {
		next;
	}
			
	unlink($file);
}

if ((length($outpath) > 0) && (-d $outpath)) {
	my @files = glob($dir_name . "/*.*");
	for (my $j = 0; $j < scalar(@files); $j++) {
		my $file = @files[$j];
		system("cp $file $outpath/");
	}
	system("rm -rf $dir_name");
}