#!/usr/local/bin/perl

# a program to detect de novo CNVs from birdseye output 

die "Usage: post_twinning_de_novos.pl <.out file> <output file>\n"    
        unless @ARGV == 2;

$output = "@ARGV[1]";

open OUT , ">$output" or die "Can't open $output: @ARGV[1].txt!\n";



###################################################################################################################################
######################################## open CNV calls and store info in hashes ##################################################

open (READ1, "@ARGV[0]") || die "Cannot open @ARGV[0]\n";

# The tab-delimited columns of the input file represent: 
# 0 - perID:			person ID - must be in format: 00001-01 -> family ID (5 nrs) + "-" + person ID (2 nrs)
# 1 - overlap_prop:		proportion overlap (Birdsuite & PennCNV)
# 2 - overlap_size:		size overlap (Birdsuite & PennCNV)
# 3 - start:			segment start (overlap Birdsuite & PennCNV)
# 4 - end:			segment end (overlap Birdsuite & PennCNV)
# 5 - chr_BS: 		chr (Birdsuite)
# 6 - start_BS:		segment start (Birdsuite)
# 7 - end_BS:			segment end (Birdsuite)
# 8 - size_BS:		size of CNV (Birdsuite)
# 9 - num_probes_BS:		number of loci (CN or SNP probes) in the segment (Birdsuite)
# 10 - cn_BS:			CN call (Birdsuite)
# 11 - chr_PC: 		chr (PennCNV)
# 12 - start_PC:		segment start (PennCNV)
# 13 - end_PC:		segment end (PennCNV)
# 14 - size_PC:		size of CNV (PennCNV)
# 15 - num_probes_PC:	number of loci (CN or SNP probes) in the segment (PennCNV)
# 16 - cn_PC:			CN call (PennCNV)
# 17 - lod_score_BS:		LOD score of the probability of the segment being the stated copy number versus the copy number of the flanking segments
# 18 - cn_agreement: 	has to be 1 for all CNVs (means Birdsuite & PennCNV call is the same) or -1 (means that there were no losses or gains -> come from Birdsuite)

$nrfam = -1;
$nrstart = -1;
$nrend = -1;

		while ($line=<READ1>) 
		{
				if ($line =~ "perID")				# skip the header
				{
				next;
				}

		@splitCNline = split (/\s+/  , $line);				# create array 

		if (@famIDraw[$nrfam] != substr(@splitCNline[0], 0, 5))
			{
			$nrfam++;
			}

		@famIDraw[$nrfam] = substr(@splitCNline[0], 0, 5);		# store family id (first 5 numbers of sample number)
		$perID = substr(@splitCNline[0], 6, 2);				# store person id (last 2 numbers of sample number)
		$chr = @splitCNline[5];						# store chromosome number (BS!)
		$start = @splitCNline[3];
		$end = @splitCNline[4];

		$CNV{$chr}{$start,$end} = "$chr\t" . "$start\t" . "$end" ;

		$CNVchr{$chr}{$start,$end} = $chr;
		$CNVstart{$chr}{$start,$end} = $start;
		$CNVend{$chr}{$start,$end} = $end;

		$size{$chr}{$start,$end} = @splitCNline[2];
		$numprobes_1{$chr}{$start,$end} = @splitCNline[9];
		$numprobes_2{$chr}{$start,$end} = @splitCNline[15];

		$cn{@famIDraw[$nrfam]}{$perID}{$CNV{$chr}{$start,$end}} = @splitCNline[10];

		$lod{@famIDraw[$nrfam]}{$perID}{$CNV{$chr}{$start,$end}} = @splitCNline[17];

		}

for($j=1;$j<=24;++$j)
	{
$totalcnvs += scalar(keys %{$CNV{$j}});
	}

$totalcnvs = 0;

##################################################################################################################################################


######## filter out duplicates from the family nr array #######

    	undef %saw;
    	@famID = grep(!$saw{$_}++, @famIDraw);  

###############################################################


################ print header of the output file ################

	print OUT "famID\t";

	print OUT "tw1_chr\t";
	print OUT "tw1_start\t";
	print OUT "tw1_end\t";
	print OUT "tw1_size\t";
	print OUT "tw1_num_probes_BS\t";
	print OUT "tw1_num_probes_PC\t";
	print OUT "tw1_cn\t";
	print OUT "tw1_lod_score\t";

	print OUT "tw2_chr\t";
	print OUT "tw2_start\t";
	print OUT "tw2_end\t";
	print OUT "tw2_size\t";
	print OUT "tw2_num_probes_BS\t";
	print OUT "tw2_num_probes_PC\t";
	print OUT "tw2_cn\t";
	print OUT "tw2_lod_score\t";

	print OUT "tw1tw2_cnv_overlap\t";
	print OUT "tw1tw2_cnv_overlap_start\t";
	print OUT "tw1tw2_cnv_overlap_end\t";
	print OUT "tw1tw2_cnv_overlap_prop\n";

##################################################################

################################ print output sub ########################################

sub printCNVs {

my (@CNVout) = @_;
#return 

print OUT "$CNVout[0]\t";

print OUT "$CNV{$CNVout[2]}{$CNVout[3],$CNVout[4]}\t";
print OUT "$size{$CNVout[2]}{$CNVout[3],$CNVout[4]}\t";
print OUT "$numprobes_1{$CNVout[2]}{$CNVout[3],$CNVout[4]}\t";
print OUT "$numprobes_2{$CNVout[2]}{$CNVout[3],$CNVout[4]}\t";
print OUT "$cn{$CNVout[0]}{'01'}{$CNV{$CNVout[2]}{$CNVout[3],$CNVout[4]}}\t";
print OUT "$lod{$CNVout[0]}{'01'}{$CNV{$CNVout[2]}{$CNVout[3],$CNVout[4]}}\t";

print OUT "$CNV{$CNVout[5]}{$CNVout[6],$CNVout[7]}\t";
print OUT "$size{$CNVout[5]}{$CNVout[6],$CNVout[7]}\t";
print OUT "$numprobes_1{$CNVout[5]}{$CNVout[6],$CNVout[7]}\t";
print OUT "$numprobes_2{$CNVout[5]}{$CNVout[6],$CNVout[7]}\t";
print OUT "$cn{$CNVout[0]}{'02'}{$CNV{$CNVout[5]}{$CNVout[6],$CNVout[7]}}\t";
print OUT "$lod{$CNVout[0]}{'02'}{$CNV{$CNVout[5]}{$CNVout[6],$CNVout[7]}}\t";

if (exists $cnv_tw1tw2_overlap{$CNV{$CNVout[2]}{$CNVout[3],$CNVout[4]}}{$CNVout[0]})
{
print OUT "$cnv_tw1tw2_overlap{$CNV{$CNVout[2]}{$CNVout[3],$CNVout[4]}}{$CNVout[0]}\t";
print OUT "$cnv_tw1tw2_overlap_start{$CNV{$CNVout[2]}{$CNVout[3],$CNVout[4]}}{$CNVout[0]}\t";
print OUT "$cnv_tw1tw2_overlap_end{$CNV{$CNVout[2]}{$CNVout[3],$CNVout[4]}}{$CNVout[0]}\t";
print OUT "$cnv_tw1tw2_overlap_prop{$CNV{$CNVout[2]}{$CNVout[3],$CNVout[4]}}{$CNVout[0]}\n";
} else {
print OUT "-1\t";
print OUT "-1\t";
print OUT "-1\t";
print OUT "-1\n";
}

}

# print_CNVs(0$famID, 1$discrepancy, 2$tw1_chr, 3$tw1_start, 4$tw1_end, 5$tw2_chr, 6$tw2_start, 7$tw2_end);

##########################################################################################


############## get maximum & minimum value subs ##############

sub max {

my (@maxinput) = @_;

	if ($maxinput[0] > $maxinput[1])
	{
	return $maxinput[0];
	}

	if ($maxinput[0] <= $maxinput[1])
	{
	return $maxinput[1];
	}

}

sub min {

my (@mininput) = @_;

	if ($mininput[0] > $mininput[1])
	{
	return $mininput[1];
	}

	if ($mininput[0] <= $mininput[1])
	{
	return $mininput[0];
	}

}

#################################################################

############################################################################################################################################################
######################################################### print the rest of the output file ################################################################

for($j=1;$j<=24;++$j)
	{
		foreach $key1 (sort keys %{$CNV{$j}})
		{

			foreach $key2 (sort keys %{$CNV{$j}})
			{

				if (($CNVstart{$j}{$key2} < $CNVstart{$j}{$key1}) && ($CNVend{$j}{$key2} < $CNVstart{$j}{$key1}))
				{
				next;
				}


				if (($CNVstart{$j}{$key2} > $CNVend{$j}{$key1}) && ($CNVend{$j}{$key2} > $CNVend{$j}{$key1}))
				{
				next;
				}
			

				for($i=0;$i<@famID;++$i)
				{

				if 	(	($CNVchr{$j}{$key1} == $CNVchr{$j}{$key2}) && 
						(	(($CNVstart{$j}{$key1} >= $CNVstart{$j}{$key2}) && ($CNVstart{$j}{$key1} <= $CNVend{$j}{$key2})) ||
							(($CNVend{$j}{$key1} >= $CNVstart{$j}{$key2}) && ($CNVend{$j}{$key1} <= $CNVend{$j}{$key2})) ||
							(($CNVstart{$j}{$key2} >= $CNVstart{$j}{$key1}) && ($CNVstart{$j}{$key2} <= $CNVend{$j}{$key1})) ||
							(($CNVend{$j}{$key2} >= $CNVstart{$j}{$key1}) && ($CNVend{$j}{$key2} <= $CNVend{$j}{$key1}))	
						)
					) 			
				{
				
					######## set missing copy nrs to -1 ########	$key1
	
					if ($cn{$famID[$i]}{'01'}{$CNV{$j}{$key1}} eq "")
						{
						$cn{$famID[$i]}{'01'}{$CNV{$j}{$key1}} = -1;
						}	

					if ($cn{$famID[$i]}{'02'}{$CNV{$j}{$key1}} eq "")
						{
						$cn{$famID[$i]}{'02'}{$CNV{$j}{$key1}} = -1;
						}

					#############################################
					######## set missing copy nrs to -1 ########	$key2
	
					if ($cn{$famID[$i]}{'01'}{$CNV{$j}{$key2}} eq "")
						{
						$cn{$famID[$i]}{'01'}{$CNV{$j}{$key2}} = -1;
						}	

					if ($cn{$famID[$i]}{'02'}{$CNV{$j}{$key2}} eq "")
						{
						$cn{$famID[$i]}{'02'}{$CNV{$j}{$key2}} = -1;
						}

					#############################################


					####### set missing lod scores to -1 ########	$key1

					if ($lod{$famID[$i]}{'01'}{$CNV{$j}{$key1}} eq "")
						{
						$lod{$famID[$i]}{'01'}{$CNV{$j}{$key1}} = -1;
						}

					if ($lod{$famID[$i]}{'02'}{$CNV{$j}{$key1}} eq "")
						{
						$lod{$famID[$i]}{'02'}{$CNV{$j}{$key1}} = -1;
						}

					#############################################
					####### set missing lod scores to -1 ########	$key2

					if ($lod{$famID[$i]}{'01'}{$CNV{$j}{$key2}} eq "")
						{
						$lod{$famID[$i]}{'01'}{$CNV{$j}{$key2}} = -1;
						}

					if ($lod{$famID[$i]}{'02'}{$CNV{$j}{$key2}} eq "")
						{
						$lod{$famID[$i]}{'02'}{$CNV{$j}{$key2}} = -1;
						}

					#############################################


					if (($cn{$famID[$i]}{'01'}{$CNV{$j}{$key1}} != $cn{$famID[$i]}{'02'}{$CNV{$j}{$key2}}) && ($cn{$famID[$i]}{'01'}{$CNV{$j}{$key1}} != -1) && ($cn{$famID[$i]}{'02'}{$CNV{$j}{$key2}} != -1))
					{

					# store CNV info for CNV 2 that has a discrepancy with CNV 1
					$cnv_tw1tw2_chr{$CNV{$j}{$key1}}{$famID[$i]} = $CNVchr{$j}{$key2};
					$cnv_tw1tw2_start{$CNV{$j}{$key1}}{$famID[$i]} = $CNVstart{$j}{$key2};
					$cnv_tw1tw2_end{$CNV{$j}{$key1}}{$famID[$i]} = $CNVend{$j}{$key2};

					$max = max($CNVstart{$j}{$key1},$cnv_tw1tw2_start{$CNV{$j}{$key1}}{$famID[$i]});
					$min = min($CNVend{$j}{$key1},$cnv_tw1tw2_end{$CNV{$j}{$key1}}{$famID[$i]});
					$size = $min - $max;

					$cnv_tw1tw2_overlap{$CNV{$j}{$key1}}{$famID[$i]} = $size;
					$cnv_tw1tw2_overlap_start{$CNV{$j}{$key1}}{$famID[$i]} = $max;
					$cnv_tw1tw2_overlap_end{$CNV{$j}{$key1}}{$famID[$i]} = $min;
					if (($size{$j}{$key1} + $size{$j}{$key2}) == $size)
						{
						$cnv_tw1tw2_overlap_prop{$CNV{$j}{$key1}}{$famID[$i]} = 1;
						} else {
						$cnv_tw1tw2_overlap_prop{$CNV{$j}{$key1}}{$famID[$i]} = $size / ($size{$j}{$key1} + $size{$j}{$key2} - $size);
						}

					# print the discrepancy to the output
					printCNVs($famID[$i],'twin1-twin2',
					$CNVchr{$j}{$key1},$CNVstart{$j}{$key1},$CNVend{$j}{$key1},
					$cnv_tw1tw2_chr{$CNV{$j}{$key1}}{$famID[$i]},$cnv_tw1tw2_start{$CNV{$j}{$key1}}{$famID[$i]},$cnv_tw1tw2_end{$CNV{$j}{$key1}}{$famID[$i]});

					}


				}
				########################################################################################################
				}
			}
	
		}
$totalcnvs += scalar(keys %{$CNV{$j}});
}


############################################################################################################################################################

close READ1;

