#!/usr/bin/env perl

# Autor:	Katrin Grothaus
# Berechnungen der durchschnittlichen Methylierung einer bestimmten Region
# Eingabe: 	Regions-Datei 
#			Methylierungs-Datei 
#			Variable letzte Chromosom

# Aufruf: 	perl Bestimmung_rel_Lage_CGI_retrogen.pl <Regions-Datei> <Methylierungs-Datei>

use warnings;
use strict;

# Variablen für die Eingabedateien
my $infile_region=$ARGV[0];
my $infile_meth=$ARGV[1];
my $last_meth_chr=$ARGV[2];

# Variablen für die Ausgabedateien
my $outfile="mean_".$infile_meth;

# benötigte Variablen für die unterschiedlichen Berechnungen und Parameter
my $counter_meth=0;
my $counter_all_cgs=0;
my $counter_all_cgs_w_reads=0;
my $counter_all_cgs_wo_reads=0;
my $counter_regions=0;

my $sum_meth=0;
my $mean_meth=0;
my $sum_reads=0;
my $mean_reads=0;

my $number100=0;
my $number0=0;

my @regions;
my @meths;

my @regions_chr;
my @regions_start;
my @regions_stop;

# Einlesen der Regions-Datei in ein Array
open (INCGI, $infile_region)or die "test Can't open file	 $!\n";
	@regions= <INCGI>;
close INCGI;

# Regionen in die einzelnen Parameter Chromosom, Start und Stopp aufsplitten
for(my $i=0; $i<@regions; $i++){
	my $region=$regions[$i];
	chomp($region);
	my @region_parts=split("\t", $region);
	push(@regions_chr,$region_parts[0]);
	push(@regions_start,$region_parts[1]);
	push(@regions_stop,$region_parts[2]);
}

# Öffnen der Ausgabedateien zum Schreiben der Ergebnisse
open(OUT,">".$outfile);
# Header der Ausgabedatei schreiben
print OUT "chr_region"."\t"."start_region"."\t"."stop_region"."\t"."#CGs"."\t"."#CGs w reads"."\t"."#CGs w/o reads"."\t"."mean meth"."\t"."mean coverage"."\t"."#CGs 100% meth"."\t"."#CGs 0% meth"."\n";

# Einlesen der Methylierungs-Datei zeilenweise mit Aufsplitten jeder Zeile in die unterschiedlichen Variablen
open (INMETH, $infile_meth)or die "test Can't open file	 $!\n";
while(<INMETH>){
	my $meth= $_;
	chomp($meth);
	my @meth_parts=split("\t", $meth);
	my $meth_chr=$meth_parts[0];
	my $meth_pos=$meth_parts[1];
	my $meth_meth=$meth_parts[2];
	my $meth_reads=$meth_parts[3];
	# letztes Chromosom gleich aktuelles Chromosom?
	if($last_meth_chr eq $meth_chr){
		# überprüfen ob die Position innerhalb der Region liegt
		if ($meth_pos<=$regions_stop[$counter_regions]){
			$counter_all_cgs++;
			# sind Reads an dieser Position vorhanden? -> Berechnung
			if ($meth_reads > 0){
				$counter_all_cgs_w_reads++;
				$sum_meth=$sum_meth + $meth_meth;
				$sum_reads=$sum_reads + $meth_reads;
				# überprüfen ob CGs vollständig methyliert oder vollständig unmethyliert sind
				if ($meth_meth==1.0){
					$number100++;
				}
				if ($meth_meth==0.0){
					$number0++;
				}
			}
			# keine Reads sind vorhanden
			else{
				$counter_all_cgs_wo_reads++;
			}
		}
		# Position liegt nicht innerhalb der Region 
		else{
			
			if($counter_all_cgs_w_reads>0){
				$mean_meth= $sum_meth/$counter_all_cgs_w_reads;
				$mean_reads=$sum_reads/$counter_all_cgs_w_reads;
			}
			# Ausgabe
			print OUT $regions_chr[$counter_regions]."\t".$regions_start[$counter_regions]."\t".$regions_stop[$counter_regions]."\t".$counter_all_cgs."\t".$counter_all_cgs_w_reads."\t".$counter_all_cgs_wo_reads."\t".$mean_meth."\t".$mean_reads."\t".$number100."\t".$number0."\n";
			$counter_regions++;
			$mean_meth=0;
			$mean_reads=0;
			$counter_all_cgs=0;
			$counter_all_cgs_w_reads=0;
			$counter_all_cgs_wo_reads=0;
			$number100=0;
			$number0=0;
			$sum_meth=0;
			$sum_reads=0;
			# überprüfen ob die Position innerhalb der Region liegt
			if ($meth_pos<=$regions_stop[$counter_regions]){
				$counter_all_cgs++;
				if ($meth_reads > 0){
					$counter_all_cgs_w_reads++;
					$sum_meth=$sum_meth + $meth_meth;
					$sum_reads=$sum_reads + $meth_reads;
					# überprüfen ob CGs vollständig methyliert oder vollständig unmethyliert sind
					if ($meth_meth==1.0){
						$number100++;
					}
					if ($meth_meth==0.0){
						$number0++;
					}
				}
				else{
					$counter_all_cgs_wo_reads++;
				}
			}
		}
	}
	# letztes Chromosom ungleich aktuelles Chromosom?
	if($last_meth_chr ne $meth_chr){
		if($counter_all_cgs_w_reads>0){
				$mean_meth= $sum_meth/$counter_all_cgs_w_reads;
				$mean_reads=$sum_reads/$counter_all_cgs_w_reads;
		}
		$last_meth_chr=$meth_chr;
		# Ausgabe
		print OUT $regions_chr[$counter_regions]."\t".$regions_start[$counter_regions]."\t".$regions_stop[$counter_regions]."\t".$counter_all_cgs."\t".$counter_all_cgs_w_reads."\t".$counter_all_cgs_wo_reads."\t".$mean_meth."\t".$mean_reads."\t".$number100."\t".$number0."\n";
		$counter_regions++;
		$mean_meth=0;
		$mean_reads=0;
		$counter_all_cgs=0;
		$counter_all_cgs_w_reads=0;
		$counter_all_cgs_wo_reads=0;
		$number100=0;
		$number0=0;
		$sum_meth=0;
		$sum_reads=0;
		# überprüfen ob die Position innerhalb der Region liegt
		if ($meth_pos<=$regions_stop[$counter_regions]){
			$counter_all_cgs++;
			# sind Reads an dieser Position vorhanden? -> Berechnung
			if ($meth_reads > 0){
				$counter_all_cgs_w_reads++;
				$sum_meth=$sum_meth + $meth_meth;
				$sum_reads=$sum_reads + $meth_reads;
				# überprüfen ob CGs vollständig methyliert oder vollständig unmethyliert sind
				if ($meth_meth==1.0){
					$number100++;
				}
				if ($meth_meth==0.0){
					$number0++;
				}
			}
			# keine Reads sind vorhanden
			else{
				$counter_all_cgs_wo_reads++;
			}
		}
	}
}
# Wenn es Reads für eine Region gibt werden nun die durchschnittliche Anzahl an Reads und die durchschnittliche Methylierung berechnet
if($counter_all_cgs_w_reads>0){
		$mean_meth= $sum_meth/$counter_all_cgs_w_reads;
		$mean_reads=$sum_reads/$counter_all_cgs_w_reads;
}
# Ausgabe
print OUT $regions_chr[$counter_regions]."\t".$regions_start[$counter_regions]."\t".$regions_stop[$counter_regions]."\t".$counter_all_cgs."\t".$counter_all_cgs_w_reads."\t".$counter_all_cgs_wo_reads."\t".$mean_meth."\t".$mean_reads."\t".$number100."\t".$number0."\n";
close(OUT);

