use strict;
use warnings;

#perl crispr.meyer.pl ref.fasta
#perl crispr.meyer.pl dmel-all-chromosome-r5.1.fasta
#  http://wfitch.bio.uci.edu/~tdlong/test.bed

my $ref = shift(@ARGV);
chomp $ref;

# read chromes into hash
my %genome = ();
my $id = '';
open F,$ref,or die $!;
while(<F>){
	chomp;
    if($_ =~ /^>(.+)/){
		$id = $1;
	}else{
        $genome{$id} .= $_;
    }
}
close F;

# forward

print("track name=MeyerCrispr description=\"Crispr Guides (Meyer)\" useScore=1\n");
foreach my $chr (keys (%genome)){
	if ($chr eq 'X' || $chr eq '2L' || $chr eq '2R' || $chr eq '3L' || $chr eq '3R'){
	my $seq = $genome{$chr};
	$seq =~ tr/acgt/ACGT/;
	while ($seq =~ /([ACGT]{18}GG)([ACGT]GG)/g){
		my $hit = $1;
		my $PAM = $2;
		my $pos = $-[0] + 1;
		print "chr$chr\t".($pos-1)."\t".($pos+20-1)."\t$hit\t750\t".'+'."\n";
		} # candidate hits

	# reverse
	while ($seq =~ /(CC[ACGT]{1})(CC[ACGT]{18})/g){
		my $hit = $2;
		my $rchit = reverse($hit);
		$rchit =~ tr/ACGT/TGCA/;
		my $PAM = $1;
		my $pos = $-[0] + 1;
		print "chr$chr\t".($pos-1+3)."\t".($pos+20-1+3)."\t$rchit\t250\t".'-'."\n";
		} # candidate hits
	}
	}
exit();

