#!/usr/bin/perl -w

# Take a (zipped) plot file and return number of unique insertion sites and total read count

use strict;

print ("Please enter a plot file:\n");
my $infile = <STDIN>;
chomp $infile;

# list of variables
my ($plus, $minus);
my $i = 0; my $j = 0; my $insT = 0;
my $readsP = 0; my $readsM = 0; my $total = 0;

# If it's a zipped plot file
if ($infile =~ m/\.gz$/){
open (ZIP, "gunzip -c $infile |") or die ("Couldn't open pipe to $infile:$!");

	while (<ZIP>){
   	# look for the two numbers relating to reads on the forward/reverse strands
		if (m/(\d+)\s+(\d+)/){
      	$plus = $1; $minus = $2; #print ("$plus\t$minus\n"); <STDIN>;
      }
      
      # count insertion sites on the forward strand
      if ($plus > 0) {
      	$i++; #print ("$i"); <STDIN>; 
         #print ("$plus"); <STDIN>
      }
      
      # count insertion sites on the reverse strand
      if ($minus > 0){
      	$j++; #print ("$j"); <STDIN>;
      }
      
      # count total unique insertion sites
      if (($plus > 0) or ($minus > 0)){
      	$insT++; #print ("$insT"); <STDIN>;
      }
      
      # add up total number of reads per strand
      $readsP += $plus; #print ("Forward strand: $readsP"); <STDIN>;
      $readsM += $minus; #print ("Reverse strand: $readsM"); <STDIN>;
	}
close ZIP;

}

# if it's not zipped
else {
open (IN, "<$infile")  or die ("Couldn't open $infile:$!");

	while (<IN>){
		# look for the two numbers relating to reads on the forward/reverse strands
      if (m/(\d+)\s+(\d+)/){
      	$plus = $1; $minus = $2; #print ("$plus\t$minus\n"); <STDIN>;
      }
      
      # count insertion sites on the forward strand
      if ($plus > 0) {
      	$i++; #print ("$i"); <STDIN>; 
         #print ("$plus"); <STDIN>
      }
      
      # count insertion sites on the reverse strand
      if ($minus > 0){
      	$j++; #print ("$j"); <STDIN>;
      }
      
      # count total unique insertion sites
      if (($plus > 0) or ($minus > 0)){
      	$insT++; #print ("$insT"); <STDIN>;
      }
      
      # add up total number of reads per strand
      $readsP += $plus; #print ("Forward strand: $readsP"); <STDIN>;
      $readsM += $minus; #print ("Reverse strand: $readsM"); <STDIN>;
	}
close IN;	

}

# calculate total read count
$total = ($readsP + $readsM);

print ("There are $i insertion sites on the forward strand, $j on the reverse ($insT total unique sites)\n");
print ("There are $readsP reads on the forward strand, $readsM on the reverse ($total total reads)\n");
