#! /usr/bin/perl -w
# File: screen_YY1.pl
# This program is designed to identify clustered YY1 binding sites from genomes.
# This program requires one argument, the name of input file containing a chromosome sequence 
# within the same directory as this program.  More fine tuned searches can be possible by modifying 
# several parameters including searching strings or distances or density, marked by #md.
# If you need more assistance, feel free to contact Dr. Joomyeong Kim at 225-578-7692 or 
# jkim@lsu.edu.  This program is also available as part of Supplementary material for the paper of Genome 
# Research entitled "Identification of clustered YY1 binding sites within Imprinting Control Regions."

$target_seq='';

if (@ARGV) {

           $target_seq= $ARGV[0];

           }else{

           print "Input all the required ARGV.\n";

}

open(RF,"$target_seq") || die "cannot open file: $!";

$/ = '>'; 

$name='';

@sequence=();

$sequence_1=''; 

%data=();

$reg_exp_f="CGCCAT.TT.\{0,2000\}CGCCAT.TT";

$reg_exp_r="AA.ATGGCG.\{0,2000\}AA.ATGGCG";

while (<RF>) {

           chomp;

           ($name,@sequence) = split "\n"; 
 
           next unless $name =~ /^(\S+)/;
     
           $sequence_1= join("", @sequence);

           $sequence='';  

          for( $i=0; $i < (length($sequence_1)) ; $i += 100000000) {  

             $sequence= substr($sequence_1, $i, 100000000); 

              while ( $sequence =~ /$reg_exp_f/g ) {    #md
     
                   my $file_id = $name;

                   my $position =  pos($sequence) - length($&) + 1 +$i;

                   my $domain = "YY1";

                   my $orientation = "+";

                   my $nucleotide = substr($sequence, pos($sequence) - 1.5*length($&), 2*length($&)); #md

                   my $count = YY1_count ($nucleotide);

                   my $density = length($nucleotide)/$count; #md

                                      if ( ($count > 2) and ($density <300) ){

                                                 $data{$position} = {       "file_id" => $file_id,

                                                         "position" => $position,

                                                         "domain"   => $domain,

                                                       "orientation" => $orientation,

                                                        "nucleotide" => $nucleotide,

                                                         "count"     =>  $count,

                                                 };

                                        }

                    }

             while ( $sequence =~ /$reg_exp_r/g ) {      #md
     
                   my $file_id = $name;

                   my $position =  pos($sequence) - length($&) + 1 +$i;

                   my $domain = "YY1";

                   my $orientation = "-";

                   my $nucleotide = revcom (substr($sequence, pos($sequence) - 1.5*length($&), 2*length($&))); #md

                   my $count = YY1_count ($nucleotide);

                   my $density = length($nucleotide)/$count; #md

                                     if ( ($count > 2) and ($density <300) ) {                              

                                             $data{$position} = {       "file_id" => $file_id,

                                                         "position" => $position,

                                                         "domain"   => $domain,

                                                       "orientation" => $orientation,

                                                        "nucleotide" => $nucleotide,

                                                        "count"       => $count,

                                                 };

                                       }

                    }

      } 

}

close (RF);

@position = sort { $data{$a}->{position} <=> $data{$b}->{position} } keys %data;

$result_YY1 = "$name"."_YY1_result";

open (OUTNT, ">$result_YY1") || die "$result_YY1: $!";

foreach my $position (@position) {

      if ( $data{$position}->{domain} eq "YY1" ) {

        print OUTNT ">$data{$position}->{file_id}"."_$data{$position}->{position}"."_$data{$position}->{orientation}"."_$data{$position}->{count}\n";

        print OUTNT "$data{$position}->{nucleotide}\n\n";

       }

}

close (OUTNT);

sub revcom {

    my($dna) = @_;

    my $revcom = reverse $dna;

    $revcom =~ tr/ACGTacgt/TGCAtgca/;

    return $revcom;
}

sub YY1_count {

    my($dna) = @_;

    my $probe = "CCAT.TT"; #md

    $count = 0;

     while ($dna =~ /$probe/g) {  #md

     $count++;

     }

    return $count;

}