#!/usr/bin/perl

# This script reads the result of gensimreads_dRNAseq.py (BED-format) counts the
# number of reads starting at and covering a certain position and prints this
# information to stdout. The columns are tab separated and contain the following
# information:
# STARTS_plus_strand COVERAGE_plus-strand STARTS_minus_strand COVERAGE_minus_strand
#
# IMPORTANT: In its current version the script can only be used for BED-files
#            containing information for a single genome/chromosome/plasmid/...!
#            The script does not take care about the ID in the first column!

use strict;

my @starts;
my @cov;

my $max_pos = 0;
my $c = 0;
while(<>) {
	$c++;	
	if ($c%1000 == 0) { 
		print STDERR "Processed $c lines of the input!\r";
	}

	my @a = split(/\s+/);

	if ($a[5] eq "-") {
		$starts[1][$a[2]]--;
		for (my $i = $a[1]; $i<=$a[2]; $i++) { $cov[1][$i]--; if ($i > $max_pos) {$max_pos = $i;}}
 	} else {
		$starts[0][$a[1]]++;
		for (my $i = $a[1]; $i<=$a[2]; $i++) { $cov[0][$i]++; if ($i > $max_pos) {$max_pos = $i;}}
	}

}

for (my $j=1; $j<=$max_pos; $j++) {

	if (not defined $starts[0][$j]) {$starts[0][$j] = 0;}
	if (not defined $starts[1][$j]) {$starts[1][$j] = 0;}
	if (not defined $cov[0][$j]) {$cov[0][$j] = 0;}
	if (not defined $cov[1][$j]) {$cov[1][$j] = 0;}
	print $starts[0][$j]."\t".$cov[0][$j]."\t".$starts[1][$j]."\t".$cov[1][$j]."\n";

}
exit;
