File Coverage

blib/lib/Shatterproof.pm

Criterion	Covered	Total	%
statement	797	1818	43.8
branch	174	490	35.5
condition	42	180	23.3
subroutine	34	51	66.6
pod	0	27	0.0
total	1047	2566	40.8

#!/usr/local/bin/perl
#The ShatterProof package is copyright (c) 2013 Ontario Institute for Cancer Research (OICR).
#
#This package and its accompanying libraries is free software; you can redistribute it and/or modify it under the terms of the GPL (either version 1, or at your option, any later version) or the Artistic License 2.0.  Refer to LICENSE for the full license text.

#OICR makes no representations whatsoever as to the SOFTWARE contained herein.  It is experimental in nature and is provided WITHOUT WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE OR ANY OTHER WARRANTY, EXPRESS OR IMPLIED.  CSHL MAKES NO REPRESENTATION OR WARRANTY THAT THE USE OF THIS SOFTWARE WILL NOT INFRINGE ANY PATENT OR OTHER PROPRIETARY RIGHT.

#By downloading this SOFTWARE, your Institution hereby indemnifies OICR against any loss, claim, damage or liability, of whatsoever kind or
#nature, which may arise from your Institution's respective use, handling or storage of the SOFTWARE.

#If publications result from research using this SOFTWARE, we ask that the Ontario Institute for Cancer Research be acknowledged and/or credit be given to OICR scientists, as scientifically appropriate.

### Shatterproof.pm ###############################################################################
# ShatterProof is a tool that can be used to analyze next generation sequencing data for signs
# of chromothripsis.
# See POD at end of file for more description
#

### INCLUDES ######################################################################################
use strict;
use warnings;
use Carp;

use vars qw($VERSIONS);

use feature 'switch';
#use Switch;
use File::Basename;
use List::Util qw[min max];
use Statistics::Distributions;
use POSIX;
###################################################################################################

### HISTORY #######################################################################################
# Version       Date            Coder   	Comments
# 0.001         2012/03/19      sgovind      	Versioning start point
# 0.002		2012/04/03	sgovind		moved input validation methods and run methods
#						from shatterproof.pl to here
# 0.003		2012/10/08	sgovind		Updated POD
# 0.04		2012/11/25	sgovind		Stable build before changing translocation scoring equation
# 0.05		2012/12/26	sgovind		See change log for details
# 0.06		2012/12/27	sgovind		Added additional documentation for new config variable
# 0.07		2012/12/27	sgovind		Added example guide for provided sample data
# 0.08		2013/05/22	sgovind		Added EXPORT code for test case, added minor error checking
# 0.09		2013/06/10	sgovind		Minor changes to accomodate testing
# 0.10		2013/06/19	sgovind		Minor changes to accomodate testing
# 0.11		2013/06/24	sgovind		Changed all sorts to stable sorts. Reduced number of posix
#						calculations.
# 0.12		2013/06/28	sgovind		Changed sort order of interchromosomal translocation output
# 0.13		2013/0716	sgovind		Corrected logical error in calculate_loh_score

our $VERSION = '0.14';

package Shatterproof;
use Exporter;
our @ISA = 'Exporter';
our @EXPORT = qw(
			$bin_size
			$genome_cnv_data_hash_ref
			$chromosome_copy_number_count_hash_ref
			$chromosome_cnv_breakpoints_hash_ref
			$tp53_mutation_found
			$genome_trans_data_hash_ref
			$chromosome_translocation_count_hash_ref
			$genome_trans_breakpoints_hash_ref
			$genome_mutation_density_hash_ref
			$suspect_regions_array_ref
			$likely_regions_array_ref
			$genome_cnv_data_windows_hash_ref
			$genome_trans_data_windows_hash_ref
			$genome_mutation_data_windows_hash_ref
			$localization_window_size
		);

### Global Variables ##############################################################################
my $pos = 0;	#used to parse command line variables
my $ARGC;	#stores the number of command line arguments provided

my %chromosome_length = (			#stores the sequence length of each chromosome
			 X => 154913754,
			 Y => 57741652,
			 1 => 247199719,
			 2 => 242751149,
			 3 => 199446827,
			 4 => 191263063,
			 5 => 180837866,
			 6 => 170896993,
			 7 => 158821424,
			 8 => 146274826,
			 9 => 140442298,
			 10 => 135374737,
			 11 => 134452384,
			 12 => 132289534,
			 13 => 114127980,
			 14 => 106360585,
			 15 => 100338915,
			 16 => 88822254,
			 17 => 78654742,
			 18 => 76117153,
			 19 => 63806651,
			 20 => 62435965,
			 21 => 46944323,
			 22 => 49528953
			);

my $TP53_start 	= 1000000*7.57;	#start base pair of the TP53 gene
my $TP53_end 	= 1000000*7.59;	#end base pair of the TP53 gene

my $insertion_data_present = 0;
my $LOH_data_present = 0;


#The values for the following 13 variables are defined in the config.pl file
our $bin_size;			#number of bases pairs that will be compressed into 1 region when analyzing the genome
						#this value defines how many base pairs are included in one array element in the data_hash_ref varaibles

our $localization_window_size;			#number of regions to sum together when performing sliding window analysis of the genome

our $expected_mutation_density;			#the expected mutation density of translocations in a highly mutated region
						#used to calculate spread factor of translocations
our $low_mutation_density_threshold;		#the mutation density that will be used to call likely regions

our $collapse_regions;				#flag variable
						#value 1:	merge overlapping CNV regions that have the same copy number
						#value 0:	do not merge overlapping CNV regions that have the same copy number. If such
						#		regions are found an error is thrown

our $outlier_deviation;				#the number of standard deviations away from the mean a value has to be in order to be considered non-significant

our $translocation_cut_off_count;		#the max number of translocation chromosomes that will be tolerated before translocation score is set to 0

our $chromosome_localization_weight;		#the scoring formula weight given to the localization of mutations to a specific region on the chromosome
our $genome_localization_weight; 		#the scoring formula weight given to the localization of mutations to a specific chromosome
our $cnv_weight;				#the scoring formula weight given to the aberrant CNV hallmark
our $translocation_weight; 			#the scoring formula weight given to the localization of translocations
our $insertion_breakpoint_weight; 		#the scoring formula weight given to the number of insertions found at translocation breakpoints
our $loh_weight;				#the scoring formula weight given to the amount of heterozygosity that is retained in a mutated region
our $tp53_mutated_weight; 			#the scoring formula weight given to the presents or absence of a TP53 mutation

### SUB METHODS ###################################################################################

#=head2 Sub-Method: run

### run ###########################################################################################
# Description:
#	Main method called by shatterproof.pl
#	Calls primary sub methods
#
# Input variables:
#	$argv_ref:	reference to @ARGV

#=cut
sub run {

	my $argv_ref = shift;		#parse parameters
	my @argv = @{$argv_ref};	#dereference array reference

	my $cnv_directory;		#stores the path to the directory where the CNV input files are found
	my $trans_directory;		#stores the path to the directory where the translocation input files are found
	my $insertion_directory;		#stores the path to the directory where the insertion input files are found
	my $loh_directory;		#stores the path to the directory where the loss of hetrozygosity input files are found
	my $output_directory;		#stores the path to the directory where output files will be placed

	my $config_file_path;		#stores the path to the configuration file

	my $tp53_mutated = 0;		#flag variable to indicate if the TP53 gene should be considered to be mutated
	my $tp53_mutation_found = 0;	#flag variable to indicate if a mutation was found in the TP53 region. This does not affect scoring

	my @cnv_files;			#list of CNV input files
	my @trans_files;		#list of translocation input files
	my @insertion_files;			#list of insertion input files
	my @loh_files;			#list of LOH input files


	my $chromosome_copy_number_count_hash_ref;	#hash
							#key1: 	chromosome eg. 1,2,X,Y
							#key2:	copy-number state eg. 0,1,3,20
							#value:	number of regions with copy number key2

	my $chromosome_cnv_breakpoints_hash_ref;	#hash
							#key: 	chromosome eg. 1,2,X,Y
							#value:	an array storing the start and end points of all cnv regions on key

	my $chromosome_translocation_count_hash_ref;	#hash
							#key1: 	chromosome eg. 1,2,X,Y
							#key2: 	chromosome eg. 1,2,X,Y
							#value:	number of translocations between key1 and key2

	my $chromosome_insertion_count_hash_ref;		#hash
							#key: 	chromosome eg. 1,2,X,Y
							#value:	number of insertions on key

	my $chromosome_loh_breakpoints_hash_ref;	#hash
							#key: 	chromosome eg. 1,2,X,Y
							#value:	an array storing the start and end points of all loh regions on key


	my $genome_trans_breakpoints_hash_ref;		#hash
							#key: 	chromosome eg. 1,2,X,Y
							#value:	an array storing all the translocation breakpoints on key

	my $genome_trans_insertion_breakpoints_hash_ref;	#hash
							#key: 	chromosome eg. 1,2,X,Y
							#value:	an array storing only the translocation breakpoints on key that have a insertion with 10 base pairs

	my $genome_mutation_density_hash_ref;		#hash
							#key: 	chromosome eg. 1,2,X,Y
							#value:	the total number of mutation on key divided by the sequence length of key


	my $genome_cnv_data_hash_ref = initialize_genome_hash();	#hash {key1}[index]{key2}
									#key1: 	chromosome eg. 1,2,X,Y
									#value:	an array storing references to hashes which contain information about the
									#	CNVs in each region of key1. The index of the array indicated the region
									#key2:	'BPcount' -> gives the number of CNV breakpoints in the region.
									#	a number, eg. '1' -> gives the number of subregions within the region that
									#	have a copy number of 1

	my $genome_trans_data_hash_ref = initialize_genome_hash();	#hash {key1}[index]{key2}{key3}
									#key1: 	chromosome eg. 1,2,X,Y
									#value1:an array storing references to hashes which contain information about the
									#	translocations in each region of key1. The index of the array indicated the region.
									#key2:	'BPcount' -> gives the number of translocation breakpoints in the region
									#	'in' -> gives a reference to a hash that contains information about translocations
									#	into the region
									#	'out' -> gives a reference to a hash that contains information about translocations
									#	out of the region
									#key3:	chromosome eg. 1,2,X,Y
									#value2:the number of subregions in the region that were translocated to key1 from key3 if key2 = 'in'
									#	or
									#	the number of subregions in the region that were translocated from key1 to key3 if key2 = 'out'

	my $genome_insertion_data_hash_ref = initialize_genome_hash();	#hash {key1}[index]
									#key1: 	chromosome eg. 1,2,X,Y
									#value:	an array storing the number of insertions in each region of key1
									#	The index of the array indicated the region


	my $genome_cnv_data_windows_hash_ref;		#hash {key1}[index]
							#key1: 	chromosome eg. 1,2,X,Y
							#value:	an array storing the number of CNV breakpoints in each window of the genome.
							#	Each window begins at the region indicated by the array index

	my $genome_trans_data_windows_hash_ref;		#hash {key1}[index]
							#key1: 	chromosome eg. 1,2,X,Y
							#value:	an array storing the number of translocation breakpoints in each window of the genome.
							#	Each window begins at the region indicated by the array index


	my $genome_mutation_data_windows_hash_ref;	#hash {key1}[index]
							#key1: 	chromosome eg. 1,2,X,Y
							#value:	an array storing the total number of mutation breakpoints in each window of the genome.
							#	Each window begins at the region indicated by the array index

	my $suspect_regions_array_ref;		#reference to array that stores regions where chromothripsis most likely occured. Format: chr start end
	my $likely_regions_array_ref;		#reference to array that stores regions where chromothripsis may have occured. Format: chr start end


	#Validate input arguements and parse them to the correct variables
	validate_input(\@argv, \$cnv_directory, \$trans_directory, \$insertion_directory, \$loh_directory, \$tp53_mutated, \$output_directory, \$config_file_path);

	#Load the values from the config file
	if(load_config_file($config_file_path)!=1){
		die("ERROR - could not load config file\n");
		}

	print "CNV dir:\t$cnv_directory\n";
	print "Trans dir:\t$trans_directory\n";

	if(defined($insertion_directory)){
		print "insertion dir:\t$insertion_directory\n";
		}

	if(defined($loh_directory)){
		print "LOH dir:\t$loh_directory\n";
		}

	print "Output dir:\t$output_directory\n";

	print "Force TP53 Mutation:\t$tp53_mutated\n\n";

	#Get a list of files for each of the provided input directories
	@cnv_files = glob ("$cnv_directory"."*.spc");
	@trans_files = glob ("$trans_directory"."*.spt");

	if(scalar(@cnv_files)==0 || scalar(@trans_files)==0){
		die "ERROR: no CNV or translocation input files found\n";
		}

	if(defined($insertion_directory)){
		@insertion_files = glob ("$insertion_directory"."*.vcf");
		}

	if(defined($loh_directory)){
		@loh_files = glob ("$loh_directory"."*.spl");
		}

	#Echo a list of all the input files
	$" = "\n\t\t";
	print "CNV files:\t@cnv_files\n\n";
	print "Trans files:\t@trans_files\n\n";
	$" = "\n\t\t";
	if(scalar(@insertion_files)==0){
		print "Indel files:\t-none\n\n";
		}
	else{
		print "Indel files:\t@insertion_files\n\n";
		}
	$" = "\n\t";
	if(scalar(@loh_files)==0){
		print "LOH files:\t-none\n\n";
		}
	else{
		print "LOH files:\t@loh_files\n\n";
		}
	$" = " ";

	#Create the output directory if it does not exist
	mkdir ("$output_directory",0770) unless (-d "$output_directory");

	#Check that the output directory exists
	if(!(-e $output_directory)){
		die "ERROR: could not create directory: $output_directory\n";
		}

	print "\n--analyzing CNV data\n";
	($genome_cnv_data_hash_ref, $chromosome_copy_number_count_hash_ref, $chromosome_cnv_breakpoints_hash_ref) = analyze_cnv_data($output_directory, \@cnv_files, $bin_size, \$tp53_mutation_found);
	print "---done analyzing CNV data\n\n";

	print "--analyzing translocation data\n";
	($genome_trans_data_hash_ref, $chromosome_translocation_count_hash_ref, $genome_trans_breakpoints_hash_ref) = analyze_trans_data($output_directory, \@trans_files, $bin_size, \$tp53_mutation_found);
	print "---done analyzing translocation data\n\n";

	#If insertion data was provided then analyze it
	if(defined($insertion_directory)){
		print "--analyzing insertion data\n";
		($genome_insertion_data_hash_ref, $chromosome_insertion_count_hash_ref, $genome_trans_insertion_breakpoints_hash_ref) = analyze_insertion_data($output_directory, \@insertion_files, $bin_size, $genome_trans_breakpoints_hash_ref, \$tp53_mutation_found);
		print "---done analyzing insertion data\n\n";
		}

	#Delete intermediate storage
	%$genome_trans_breakpoints_hash_ref = ();
	undef $genome_trans_breakpoints_hash_ref;

	#If LOH data was provided then analyze it
	if(defined($loh_directory)){
		print "--analyzing LOH data\n";
		($chromosome_loh_breakpoints_hash_ref) = analyze_loh_data($output_directory, \@loh_files, \$tp53_mutation_found);
		print "---done analyzing LOH data\n\n";

		#Check that the correct format of the LOH hash has been preserved
		my %loh_hash = %{$chromosome_loh_breakpoints_hash_ref};
		for my $key1 (keys %loh_hash){
			my @a = @{$loh_hash{$key1}};
			my $size = @a;

			if($size % 2 != 0){
				die "ERROR: odd number of loh breakpoints recorded for chromosome $key1\n";
				}
			}

		}

	print "--calculating chromosome mutation densities\n";
	$genome_mutation_density_hash_ref = calculate_genome_localization($output_directory, $chromosome_copy_number_count_hash_ref, $chromosome_translocation_count_hash_ref);
	print "---done calculating chromosome mutation densities\n\n";

	print "--calculating chromosome region mutation densities\n";
	($suspect_regions_array_ref, $likely_regions_array_ref, $genome_cnv_data_windows_hash_ref, $genome_trans_data_windows_hash_ref, $genome_mutation_data_windows_hash_ref) = calculate_chromosome_localization($output_directory, $genome_cnv_data_hash_ref, $genome_trans_data_hash_ref, $bin_size, $localization_window_size);
	print "---done calculating chromosome region mutation densities\n\n";

	print "--analyzing suspect regions\n";
	analyze_suspect_regions($output_directory, $suspect_regions_array_ref, $genome_mutation_density_hash_ref, $genome_cnv_data_hash_ref, $genome_trans_data_hash_ref, $genome_trans_insertion_breakpoints_hash_ref, $bin_size, $localization_window_size, $tp53_mutated, $tp53_mutation_found, $chromosome_cnv_breakpoints_hash_ref, $chromosome_loh_breakpoints_hash_ref);
	print "---done analyzing suspect regions\n\n";

	print "--analyzing likely regions\n";
	analyze_likely_regions($output_directory, $likely_regions_array_ref, $genome_mutation_density_hash_ref, $genome_cnv_data_hash_ref, $genome_trans_data_hash_ref, $bin_size, $localization_window_size);
	print "---done analyzing likely regions\n\n";

	print "--calculating copy number count\n";
	check_copy_number_count($output_directory, $chromosome_copy_number_count_hash_ref);
	print "---done calculating copy number count\n\n";

	print "--calculating switch count\n";
	check_copy_number_switches($output_directory, $chromosome_copy_number_count_hash_ref);
	print "---done calculating switch count\n\n";

	print "--calculating interchromosomal translocation rate\n";
	calculate_interchromosomal_translocation_rate($output_directory, $chromosome_translocation_count_hash_ref);
	print "---done calculating interchromosomal translocation rate\n";

	}#sub run

#=head2 Sub-Method: validate_input

### validate_input ################################################################################
# Description:
#	Validates command line arguments. Prints error messages if some input if invalid.
#
# Input variables:
#	$argv_ref:		reference to @ARGV
#	$cnv_directory_ref:	reference to variable storing the CNV input directory
#	$trans_directory_ref:	reference to variable storing the translocation input directory
#	$insertion_directory_ref:	reference to variable storing the insertion input directory
#	$loh_directory_ref:	reference to variable storing the LOH input directory
#	$tp53_mutated_ref:	reference to variable storing the tp53 mutated flag
#	$output_directory_ref:	reference to variable storing the output directory
#	$config_file_path_ref:	reference to variable storing the path to the config file

#=cut
sub validate_input {

	#Parse parameters
	my $argv_ref = shift;
	my @argv = @{$argv_ref};

	my $cnv_directory_ref = shift;
	my $trans_directory_ref = shift;
	my $insertion_directory_ref = shift;
	my $loh_directory_ref = shift;
	my $tp53_mutated_ref = shift;
	my $output_directory_ref = shift;
	my $config_file_path_ref = shift;

	#Determine number of command line arguements
	$ARGC = @argv;

	#Parse the command line arguements
	given ($ARGC) {
		when (/^0$/) { usage(0); }	#Print error message if no arguements were entered

		when (/^1$/) {		#Check for help option
			if($argv[0] eq "--help"){
				man_text();
				}
			else{
				usage(1);
				}
		       }#case 1

		default   {

			if($argv[$pos] eq "--cnv"){	#Check for the cnv input directory option, this field is mandatory
				next_arg(2);
				$$cnv_directory_ref = $argv[$pos];
				if(!(substr($$cnv_directory_ref,-1,1) eq '/')){
					$$cnv_directory_ref = $$cnv_directory_ref.'/';
					}
				next_arg(3);
				}
			else {
				usage(4);
				}


			if($argv[$pos] eq "--trans"){	#Check for the translocation input directory option, this field is mandatory
				next_arg(5);
				$$trans_directory_ref = $argv[$pos];
				if(!(substr($$trans_directory_ref,-1,1) eq '/')){
					$$trans_directory_ref = $$trans_directory_ref.'/';
					}
				next_arg(6);
				}
			else {
				usage(7);
				}

			if($argv[$pos] eq "--insrt"){	#Check for the insertion input directory option
				next_arg(8);
				$$insertion_directory_ref = $argv[$pos];
				if(!(substr($$insertion_directory_ref,-1,1) eq '/')){
					$$insertion_directory_ref = $$insertion_directory_ref.'/';
					}
				$insertion_data_present = 1;
				next_arg(9);
				}

			if($argv[$pos] eq "--loh"){	#Check for the LOH input directory option
				next_arg(10);
				$$loh_directory_ref = $argv[$pos];
				if(!(substr($$loh_directory_ref,-1,1) eq '/')){
					$$loh_directory_ref = $$loh_directory_ref.'/';
					}
				$LOH_data_present = 1;
				next_arg(11);
				}

			if($argv[$pos] eq "--tp53"){	#Check for the TP53 gene mutation check option
				$$tp53_mutated_ref = 1;
				next_arg(12);
				}

			if($argv[$pos] eq "--config"){	#Check for the config file option, this field is mandatory
				next_arg(13);
				$$config_file_path_ref = $argv[$pos];
				next_arg(14);
				}
			else{
				usage(15);
				}

			if($argv[$pos] eq "--output"){	#Check for the output directory option, this field is mandatory
				next_arg(16);
				$$output_directory_ref = $argv[$pos];
				if(!(substr($$output_directory_ref,-1,1) eq '/')){
					$$output_directory_ref = $$output_directory_ref.'/';
					}
				}
			else {
				usage(17);
				}

			#Check that there are no other command line arguments
			if($pos != $ARGC-1){
				usage(18);
				}
		       }#default case
		}#given ($ARGC)

	}#sub validate_input

#=head2 Sub-Method: analyze_cnv_data

### analyze_cnv_data ##############################################################################
# Description:
#	Reads data from files located in the CNV input directory and populates:
#		$genome_cnv_data_hash_ref
#		$chromosome_copy_number_count_hash_ref
#		$chromosome_cnv_breakpoints_hash_ref
#
# Input variables:
#	$output_directory:	stores the path to the output directory
#	$cnv_files_array_ref:	reference to array containing all the CNV input files
#	$bin_size:		stores the size of the bins which the chromosome will be divided into
#	$tp53_mutation_found_ref:	reference to the tp53 mutation found flag

#=cut
sub analyze_cnv_data {

	#Parse the parameters
	my $output_directory = shift;

	my $cnv_files_array_ref = shift;
	my @cnv_files = @$cnv_files_array_ref;

	my $bin_size = shift;
	my $tp53_mutation_found_ref = shift;

	my %genome_cnv_data = ();	#hash
					#key:	chromosome eg. 1,2,X,Y
					#value:	a reference to an array where each element corresponds to a bin along the
					#	chromosome

	my @file_data;		#an array storing all the entries from every input file

	my $CURRENT_FILE;	#file handle to the current file that is open
	my $TP53_FILE;		#file handle to the TP53 CNV mutation output file

	my $line;		#stores raw line read in from file
	my @line_data;		#stores tokenized line read in from file

	my %chromosome_copy_number_count = (); 		#hash  {chr}{copy number}{count}
							#key1: chromosome eg. 1,2,X,Y
							#key2: a copy-number state eg 0,1,3,15
							#value: the number of region on key1 that have a copy number of key2

	my %chromsome_cnv_breakpoints = (		#hash 	{chr}[start and end pairs]
					 X => [],	#key: 	chromosome eg. 1,2,X,Y
					 Y => [],	#value: an array that stored an ordered list of CNV breakpoints on key
					 1 => [],
					 2 => [],
					 3 => [],
					 4 => [],
					 5 => [],
					 6 => [],
					 7 => [],
					 8 => [],
					 9 => [],
					 10 => [],
					 11 => [],
					 12 => [],
					 13 => [],
					 14 => [],
					 15 => [],
					 16 => [],
					 17 => [],
					 18 => [],
					 19 => [],
					 20 => [],
					 21 => [],
					 22 => []
					);

	#Read the contents of the cnv files into memory
	if($#cnv_files==-1){
		die "ERROR: no cnv files found in analyze_cnv_data\n";
		}

	foreach my $file (@cnv_files){
		#Open the file
		open ($CURRENT_FILE, "<", $file) or die "ERROR: could not open file at path $file\n";

		#Check that the file is not empty
		if(eof($CURRENT_FILE)){
			close ($CURRENT_FILE);
			die "ERROR: $file is empty\n";
			}

		#Read header line and validate
		$line = <$CURRENT_FILE>;
		chomp($line);

		#Check that the format of the header line is correct
		if(!($line =~ m/^#chr\tstart\tend\tnumber\tquality$/)){
			close ($CURRENT_FILE);
			die "ERROR: header of cnv file $file is invalid\n";
			}

		#Read all the data lines in the file
		while( !(eof($CURRENT_FILE)) ){
			#read data line
			$line = <$CURRENT_FILE>;
			chomp($line);

			#Validate the data line
			if(!($line =~ m/^(chr)?(1[0-9]|2[0-2]|X|Y|[1-9])\t[0-9]+\t[0-9]+\t[0-9]+\t([0-9]+|\.)$/)){
				die "ERROR: invalid line found ($line) in file $file\n";
				}

			#Split the data line and add it to the file_data array
			@line_data = (split (/\t/,$line));

			push(@file_data,[@line_data]);
			}

		close ($CURRENT_FILE);

		}#foreach my $file (@cnv_files)

	#Check that there are no overlapping CNV regions with different copy-numbers
	@file_data = check_for_overlaps("cnv", \@file_data);


	#Create TP53 directory and output folder
	mkdir ("$output_directory"."TP53",0770) unless (-d "$output_directory"."TP53");
	if(!(-e "$output_directory"."TP53")){
		die "ERROR: could not create folder $output_directory"."TP53\n";
		}

	#Create the TP53 CNV mutation file
	open ($TP53_FILE, ">", "$output_directory"."TP53/TP53.spc") or die "ERROR: could not create file: $output_directory"."TP53/TP53.spc";

	#Print the header of the file (same format as a .spc file)
	print $TP53_FILE "#chr\tstart\tend\tnumber\tquality";

	#For every data line that was read in from an input file,
	#record the CNV mutation in the genome_cnv_data_hash,
	#record the exact breakpoints of the CNV,
	#update the chromosome_copy_number_count hash and
	#check if the CNV is in the TP53 region
	for (my $n = 0; $n < scalar(@file_data); $n++){
		my $hash = {};

		#Ensure that the chromosome value is valid
		if(!($file_data[$n][0] =~ m/^(chr)?(1[0-9]|2[0-2]|X|Y|[1-9])$/)){
			die "ERROR: invalid line found in CNV input file: @{$file_data[$n]}\n";
			}

		#Parse the chromosome
		my $chr = $2;

		#Increment the copy-number count hash based on the line data
		$chromosome_copy_number_count {$chr}{$file_data[$n][3]}++;

		#Record the exact breakpoints of the CNV
		push (@{$chromsome_cnv_breakpoints{$chr}}, ($file_data[$n][1],$file_data[$n][2]));

		#Calculate the bin for the start and end breakpoint
		my $start_index = int($file_data[$n][1]/$bin_size);
		my $end_index = int($file_data[$n][2]/$bin_size);

		my $update_bin = sub {
			my $source_chr 	= shift;
			my $index 	= shift;
			my $copy_num	= shift;

			my $genome_hash_ref	= shift;
			my %genome_hash		= %{$genome_hash_ref};

			my $hash = ();

			if(!defined(@{$genome_hash{$source_chr}}[$index])){
				$hash->{'BPcount'} = 1;
				$hash->{$copy_num} = 0.5;
				@{$genome_hash{$source_chr}}[$index] = $hash;
				}
			else{	#if a bin does exist then increment the counts
				${@{$genome_hash{$source_chr}}[$index]}{'BPcount'}++;
				${@{$genome_hash{$source_chr}}[$index]}{$copy_num}+=0.5;
				}
			};

		#Check if a bin exists at the start index
		#if not, create one
		if(!defined(@{$genome_cnv_data{$chr}}[$start_index])){
			$hash->{'BPcount'} = 1;
			$hash->{$file_data[$n][3]} = 0.5;
			@{$genome_cnv_data{$chr}}[$start_index] = $hash;
			}
		else{ 	#If one does exist increment the counts
			${@{$genome_cnv_data{$chr}}[$start_index]}{'BPcount'}++;
			${@{$genome_cnv_data{$chr}}[$start_index]}{$file_data[$n][3]}+=0.5;
			}

		$hash = {};

		#Check if a bin exists at the end index
		#if not, create one
		if(!defined(@{$genome_cnv_data{$chr}}[$end_index])){
			$hash->{'BPcount'} = 1;
			$hash->{$file_data[$n][3]} = 0.5;
			@{$genome_cnv_data{$chr}}[$end_index] = $hash;
			}
		else{ 	#If one does exist increment the counts
			${@{$genome_cnv_data{$chr}}[$end_index]}{'BPcount'}++;
			${@{$genome_cnv_data{$chr}}[$end_index]}{$file_data[$n][3]}+=0.5;
			}


		#Check if the variation was in the TP53 gene
		if(
		   ( $chr ne 'X' && $chr ne 'Y' ) &&
		   ( $chr==17 )
		  ){
			if(
			   ( $file_data[$n][3] != 2 ) &&
			   ( ( $file_data[$n][1] >= $TP53_start && $file_data[$n][1] <= $TP53_end ) || ( $file_data[$n][2] >= $TP53_start && $file_data[$n][2] <= $TP53_end ) )
			  ){
				#If a CNV was found in the TP53 region
				$$tp53_mutation_found_ref = 1;

				#Record the mutation in the TP53 CNV file
				print $TP53_FILE "\n";
				for (my $i = 0; $i < scalar(@{$file_data[$n]}); $i++){
					print $TP53_FILE "$file_data[$n][$i]";
					if($i != scalar(@{$file_data[$n]})-1){
						print $TP53_FILE "\t";
						}#if
					}#for (my $i = 0; $i < scalar(@{$file_data[$n]}); $i++)
				}#if
			}#if
		}#for (my $n = 0; $n < scalar(@file_data); $n++)

	close($TP53_FILE);

	#return hash
	return(\%genome_cnv_data, \%chromosome_copy_number_count, \%chromsome_cnv_breakpoints);

	}#sub analyze_cnv_data

#=head2 Sub-Method: check_for_overlaps

### check_for_overlaps ############################################################################
# Description:
#	Checks if there were overlapping CNV regions with different copy-numbers in the input files.
#	Also checks if there are any overlapping translocation destinations or overlapping LOH
#	regions.
#
# Input variables:
#	$type:		flag variable indicating which type of overlap to check for "cnv", "trans",
#			or "loh"
#	$file_data_ref:	reference to an array storing all the data lines read in from the specific
#			type of input file

#=cut
sub check_for_overlaps {

	my $type = shift;
	my $file_data_ref = shift;
	my @file_data = @$file_data_ref;

	my $start_overlap = 0;	#Flag variable indicating if the start position of one region is overlapping with
				#another region

	my $end_overlap = 0;	#Flag variable indicating if the end position of one region is overlapping with
				#another region

	#Check for overlapping regions
	#Compares each entry in the array to every following entry
	for (my $n = 0; $n < scalar(@file_data); $n++){
		for (my $k = $n+1; $k < scalar(@file_data); $k++){

			#Check if the 2 regions in question are on the same chromosome
			if($file_data[$n][0] eq $file_data[$k][0]) {

				#Check if the end point of region 1 is within region 2
				if($file_data[$n][2]>=$file_data[$k][1] && $file_data[$n][2]<=$file_data[$k][2]){
					$end_overlap = 1;
					}

				#Check if the start point of region 1 is within region 2
				if($file_data[$n][1]>=$file_data[$k][1] && $file_data[$n][1]<=$file_data[$k][2]){
					$start_overlap = 1;
					}

				#If an overlap was detected
				if($start_overlap==1 || $end_overlap==1) {

					#if it was a translocation overlap then throw an error
					if($type eq "trans"){
						die "ERROR: found overlapping translocation source regions:\n\t@{$file_data[$n]}\n\t@{$file_data[$k]}\n";
						}

					#if it was a LOH overlap then throw an error
					if($type eq "loh"){
						die "ERROR: found overlapping LOH regions:\n\t@{$file_data[$n]}\n\t@{$file_data[$k]}\n";
						}

					#if it was a CNV overlap check if the copy numbers are the same
					#If they are different throw an error
					if(
					   ( $type eq "cnv") &&
					   ( $file_data[$n][3] != $file_data[$k][3] )
					  ){
						die "ERROR: found overlapping regions with different copy number values:\n\t@{$file_data[$n]}\n\t@{$file_data[$k]}\n";
						}

					#if they are the same
				 	#and the user does not wish to collapse overlapping regions with the same copy number then throw an error
					elsif ($collapse_regions==0) {
						die "ERROR: found overlapping copy number regions:\n\t@{$file_data[$n]}\n\t@{$file_data[$k]}\n";
						}

					#if the user wishes to collapses overlapping regions with the same copy number then do so
					elsif ($collapse_regions==1) {
						#Region 2 completely encompasses region 1
						#So replace region 1 with region 2
						if($start_overlap==1 && $end_overlap==1){
							$file_data[$n][1] = $file_data[$k][1];
							$file_data[$n][2] = $file_data[$k][2];
							}
						#The start point of region 1 is within region 2
						#So replace the start point of region 1 with the start point of region 2
						elsif($start_overlap==1){
							$file_data[$n][1] = $file_data[$k][1];
							}
						#The end point of region 1 is within region 2
						#So replace the end point of region 1 with the end point of region 2
						elsif($end_overlap==1){
							$file_data[$n][2] = $file_data[$k][2];
							}
						#If region 1 was modified then remove region 2 and re-check for overlaps
						if($start_overlap==1 || $end_overlap==1){
							@file_data = (@file_data[0..($k-1),($k+1)..(scalar(@file_data)-1)]);
							$start_overlap = 0;
							$end_overlap = 0;
							$k = $n+1;
							redo;
							}

						}#elsif ($collapse_regions==1)
					}#if($start_overlap==1 || $end_overlap==1)

				#Check if region 1 completely encompasses region 2
				elsif($file_data[$n][1]<=$file_data[$k][1] && $file_data[$n][2]>=$file_data[$k][2]){

					if($type eq "trans"){
						die "ERROR: found overlapping translocation source regions:\n\t@{$file_data[$n]}\n\t@{$file_data[$k]}\n";
						}

					if($type eq "loh"){
						die "ERROR: found overlapping LOH regions:\n\t@{$file_data[$n]}\n\t@{$file_data[$k]}\n";
						}

					#If the copy numbers are different throw an error
					if(
					   ( $type eq "cnv") &&
					   ( $file_data[$n][3] != $file_data[$k][3] )
					  ){
						die "ERROR: found overlapping regions with different copy number values:\n\t@{$file_data[$n]}\n\t@{$file_data[$k]}\n";
						}

					#If the copy numbers are the same but the user does not want to collapse then throw an error
					elsif ($collapse_regions==0) {
						die "ERROR: found overlapping copy number regions:\n\t@{$file_data[$n]}\n\t@{$file_data[$k]}\n";
						}

					#If the user does wish to collapse then remove the second region
					elsif ($collapse_regions==1) {
						@file_data = (@file_data[0..($k-1),($k+1)..(scalar(@file_data)-1)]);
						}
					}#elsif($file_data[$n][1]<=$file_data[$k][1] && $file_data[$n][2]>=$file_data[$k][2])
				}#if($file_data[$n][0] eq $file_data[$k][0])
			}#for (my $n = 0; $n < scalar(@file_data); $n++)
		}#for (my $n = 0; $n < scalar(@file_data); $n++)

	#Return the updated entries
	return (@file_data);

	}#sub check_for_overlaps

#=head2 Sub-Method: analyze_trans_data

### analyze_trans_data ############################################################################
# Description:
#	Reads data from files located in the trans input directory and popultates:
#		$genome_trans_data_hash_ref
#		$chromosome_translocation_count_hash_ref
#		$genome_trans_breakpoints_hash_ref
#
# Input variables:
#	$output_directory:	stores the path to the output directory
#	$trans_files_array_ref:	reference to array containing all the translocation
#				input files
#	$bin_size:		stores the size of the bins which the chromosome will be divided into
#	$tp53_mutation_found_ref:      reference to the tp53 mutation found flag

#=cut
sub analyze_trans_data {

	#Parse the parameters
	my $output_directory = shift;
	my $trans_files_array_ref = shift;
	my @trans_files = @$trans_files_array_ref;

	my $bin_size = shift;

	my $tp53_mutation_found_ref = shift;

	my %genome_trans_data = ();	#hash
					#key:   chromosome eg. 1,2,X,Y
					#value: a reference to an array where each element corresponds to a bin along the
					#	chromosome

	my @file_data;		#an array storing all the entries from every input file

	my $CURRENT_FILE;	#file handle to the current file that is open
	my $TP53_FILE;		#file handle to the TP53 translocation mutation output file

	my $line;		#stores raw line read in from file
	my @line_data;		#stores tokenized line read in from file

	my $chr1;
	my $chr2;

	my %chromosome_trans_count = (); 	#hash 	{chr}{chr}{count}
						#key1:	chromosome eg. 1,2,X,Y
						#key2: 	chromosome eg. 1,2,X,Y
						#value:	the number of translocations between key1 and key2

	my %genome_trans_breakpoints = 	(		#hash {chr}[array]
					 X => [],	#key:	chromosome eg. 1,2,X,Y
					 Y => [],	#value:	an array storing all the translocation breakpoints
					 1 => [],	#	on key
					 2 => [],
					 3 => [],
					 4 => [],
					 5 => [],
					 6 => [],
					 7 => [],
					 8 => [],
					 9 => [],
					 10 => [],
					 11 => [],
					 12 => [],
					 13 => [],
					 14 => [],
					 15 => [],
					 16 => [],
					 17 => [],
					 18 => [],
					 19 => [],
					 20 => [],
					 21 => [],
					 22 => []
					);


	#Read the contents of the cnv files into memory
	if($#trans_files==-1){
		die "ERROR: no trans files found in analyze_trans_data\n";
		}

	foreach my $file (@trans_files){
		#open the file
		open ($CURRENT_FILE, "<", $file) or die "ERROR: could not open file at path $file\n";

		#check that the file is not empty
		if(eof($CURRENT_FILE)){
			close ($CURRENT_FILE);
			die "ERROR: $file is empty\n";
			}

		#read header line and validate
		$line = <$CURRENT_FILE>;
		chomp($line);

		#validate the header line
		if(!($line =~ m/^#chr1\tstart\tend\tchr2\tstart\tend\tquality$/)){
			close ($CURRENT_FILE);
			die "ERROR: header of translocation file $file is invalid\n";
			}

		#read in every data line
		while( !(eof($CURRENT_FILE)) ){
			#read data line
			$line = <$CURRENT_FILE>;
			chomp($line);

			#validate the format of the data line
			if(!($line =~ m/^(chr)?(1[0-9]|2[0-2]|X|Y|[1-9])\t[0-9]+\t[0-9]+\t(chr)?(1[0-9]|2[0-2]|X|Y|[1-9])\t[0-9]+\t[0-9]+\t([0-9]+|\.)$/)){
				die "ERROR: invalid line found in translocation input file: ($line) in file $file\n";
				next;
				}

			#Split the data line and add it to the array
			@line_data = (split (/\t/,$line));

			#if the start position is  greater than the end position for either the source or destination throw an error
			if(
			   ( $line_data[1] >= $line_data[2] ) ||
			   ( $line_data[4] >= $line_data[5] )
			  ){
				#warn "ERROR: invalid line found ($line) in file $file. Start or end values invalid\n";
				#next;
				}

			#Add the data line to the file_data array
			push(@file_data,[@line_data]);
			}

		close ($CURRENT_FILE);

		}#foreach my $file (@trans_files)

	#ignoring overlapping translocations for now
	#@file_data = check_for_overlaps("trans", \@file_data);

	#Create TP53 directory and output folder
	mkdir ("$output_directory"."TP53",0770) unless (-d "$output_directory"."TP53");
	if(!(-e "$output_directory"."TP53")){
		die "ERROR: could not create folder $output_directory"."TP53\n";
		}

	#create TP53 translocation mutation file
	open ($TP53_FILE, ">", "$output_directory"."TP53/TP53.spt") or die "ERROR: could not create file: $output_directory"."TP53/TP53.spt";

	#Print the header of the file (same format as a .spt file)
	print $TP53_FILE "#chr1\tstart\tend\tchr2\tstart\tend\tquality";

	#For every data line that was read in from an input file,
	#record the translocation mutation in the genome_trans_data_hash,
	#record the exact breakpoints of the translocation,
	#update the chromosome_trans_count hash and
	#check if the translocation is in the TP53 region
	for (my $n = 0; $n < scalar(@file_data); $n++){

		my $hash = {};

		#verify that the chromosome 1 is valid
		if(!($file_data[$n][0] =~ m/^(chr)?(1[0-9]|2[0-2]|X|Y|[1-9])$/)){
			die "ERROR: invalid chromosome field detected in translocation file\n";
			}

		#parse out chromosome 1
		my $chr1 = $2;

		#verify that the chromosome 2 is valid
		if(!($file_data[$n][3] =~ m/^(chr)?(1[0-9]|2[0-2]|X|Y|[1-9])$/)){
			die "ERROR: invalid chromosome field detected in translocation file\n";
			}

		#parse out chromosome 2
		my $chr2 = $2;

		#calculate the bin where each breakpoint will be placed
		my $start_index1 = int($file_data[$n][1]/$bin_size);
		my $end_index1 = int($file_data[$n][2]/$bin_size);

		my $start_index2 = int($file_data[$n][4]/$bin_size);
		my $end_index2 = int($file_data[$n][5]/$bin_size);

		my $update_bin = sub {
			my $source_chr 	= shift;
			my $dest_chr	= shift;
			my $index 	= shift;
			my $data	= shift;
			my $type	= shift;

			my $genome_hash_ref	= shift;
			my %genome_hash		= %{$genome_hash_ref};

			my $hash = ();

			if(!defined(@{$genome_hash{$source_chr}}[$index])){
				$hash->{'BPcount'} = 1;
				push (@{$hash->{$type}{$dest_chr}}, $data);
				@{$genome_hash{$source_chr}}[$index] = $hash;
				}
			else{	#if a bin does exist then increment the counts
				${@{$genome_hash{$source_chr}}[$index]}{'BPcount'}++;
				push (@{${@{$genome_hash{$source_chr}}[$index]}{$type}{$dest_chr}}, $data);
				}

			};

		#check if a bin exists at $start_index1
		#if not, create one
		if(!defined(@{$genome_trans_data{$chr1}}[$start_index1])){
			$hash->{'BPcount'} = 1;
			push (@{$hash->{'out'}{$chr2}}, $file_data[$n][4]);
			@{$genome_trans_data{$chr1}}[$start_index1] = $hash;
			}
		else{	#if a bin does exist then increment the counts
			${@{$genome_trans_data{$chr1}}[$start_index1]}{'BPcount'}++;
			push (@{${@{$genome_trans_data{$chr1}}[$start_index1]}{'out'}{$chr2}}, $file_data[$n][4]);
			}

		$hash = {};
		if(!defined(@{$genome_trans_data{$chr1}}[$end_index1])){
			$hash->{'BPcount'} = 1;
			push (@{$hash->{'out'}{$chr2}}, $file_data[$n][5]);
			@{$genome_trans_data{$chr1}}[$end_index1] = $hash;
			}
		else{
			${@{$genome_trans_data{$chr1}}[$end_index1]}{'BPcount'}++;
			push (@{${@{$genome_trans_data{$chr1}}[$end_index1]}{'out'}{$chr2}}, $file_data[$n][5]);
			}

		$hash = {};
		if(!defined(@{$genome_trans_data{$chr2}}[$start_index2])){
			$hash->{'BPcount'} = 1;
			push (@{$hash->{'in'}{$chr1}}, $file_data[$n][1]);
			@{$genome_trans_data{$chr2}}[$start_index2] = $hash;
			}
		else{
			${@{$genome_trans_data{$chr2}}[$start_index2]}{'BPcount'}++;
			push (@{${@{$genome_trans_data{$chr2}}[$start_index2]}{'in'}{$chr1}}, $file_data[$n][1]);
			}

		$hash = {};
		if(!defined(@{$genome_trans_data{$chr2}}[$end_index2])){
			$hash->{'BPcount'} = 1;
			push (@{$hash->{'in'}{$chr1}}, $file_data[$n][2]);
			@{$genome_trans_data{$chr2}}[$end_index2] = $hash;
			}
		else{
			${@{$genome_trans_data{$chr2}}[$end_index2]}{'BPcount'}++;
			push (@{${@{$genome_trans_data{$chr2}}[$end_index2]}{'in'}{$chr1}}, $file_data[$n][2]);
			}

		#Increment hash translocation counts
		$chromosome_trans_count{$chr1}{$chr2}++;
		#if the translocation is intra-chromosomal then don't count it twice
		if($chr1 ne $chr2){
			$chromosome_trans_count{$chr2}{$chr1}++;
			}

		#store the breakpoints in their bins
		push (@{$genome_trans_breakpoints{$chr1}}, $file_data[$n][1]);
		push (@{$genome_trans_breakpoints{$chr1}}, $file_data[$n][2]);

		push (@{$genome_trans_breakpoints{$chr2}}, $file_data[$n][4]);
		push (@{$genome_trans_breakpoints{$chr2}}, $file_data[$n][5]);

		#Check if the translocation origin was in the TP53 gene
		if(
		   ( $chr1 ne 'X' && $chr1 ne 'Y' ) &&
		   ( $chr1==17 )
		  ){
			if(
			   ( ( $file_data[$n][1] >= $TP53_start && $file_data[$n][1] <= $TP53_end ) || ( $file_data[$n][2] >= $TP53_start && $file_data[$n][2] <= $TP53_end ) )
			  ){
				#if a mutation was found, set the TP53 mutated flag
				$$tp53_mutation_found_ref = 1;
				print $TP53_FILE "\n";
				#print the translocation data line to the TP53 translocation mutation output file
				for (my $i = 0; $i < scalar(@{$file_data[$n]}); $i++){
					print $TP53_FILE "$file_data[$n][$i]";
					if($i != scalar(@{$file_data[$n]})-1){
						print $TP53_FILE "\t";
						}
					}#for (my $i = 0; $i < scalar(@{$file_data[$n]}); $i++)
				}#if
			}#if

		#Check if the translocation destination was in the TP53 gene
		if(
		   ( $chr2 ne 'X' && $chr2 ne 'Y' ) &&
		   ( $chr2==17 )
		  ){
			if(
			   ( ( $file_data[$n][4] >= $TP53_start && $file_data[$n][4] <= $TP53_end ) || ( $file_data[$n][5] >= $TP53_start && $file_data[$n][5] <= $TP53_end ) )
			  ){
				$$tp53_mutation_found_ref = 1;
				print $TP53_FILE "\n";
				for (my $i = 0; $i < scalar(@{$file_data[$n]}); $i++){
					print $TP53_FILE "$file_data[$n][$i]";
					if($i != scalar(@{$file_data[$n]})-1){
						print $TP53_FILE "\t";
						}
					}#for (my $i = 0; $i < scalar(@{$file_data[$n]}); $i++)
				}#if
			}#if
		}#for (my $n = 0; $n < scalar(@file_data); $n++)

	close($TP53_FILE);

	#return hash
	return(\%genome_trans_data, \%chromosome_trans_count, \%genome_trans_breakpoints);

	}#sub analyze_trans_data


#=head2 Sub-Method: analyze_insertion_data

### analyze_insertion_data ##############################################################################
# Description:
#	Reads data from files located in the insertion input directory and populates:
#		$genome_insertion_data_hash_ref
#		$chromosome_insertion_count_hash_ref
#		$genome_trans_insertion_breakpoints_hash_ref
#
# Input variables:
#	$output_directory:			stores the path to the output directory
#	$insertion_files_array_ref:		reference to array containing all the insertion input files
#	$bin_size:				stores the size of the bins which the chromosome will be divided into
#	$genome_trans_breakpoints_hash_ref:	store reference to hash that contains the translocation breakpoints on
#						each chromosome
#	$tp53_mutation_found_ref:      		reference to the tp53 mutation found flag
#

#=cut
sub analyze_insertion_data {

	#Parse Parameters
	my $output_directory = shift;
	my $insertion_files_array_ref = shift;
	my @insertion_files = @$insertion_files_array_ref;

	my $bin_size = shift;

	my $genome_trans_breakpoints_hash_ref = shift;
	my %genome_trans_breakpoints = %$genome_trans_breakpoints_hash_ref;

	my $tp53_mutation_found_ref = shift;

	my %genome_insertion_data = ();	#hash
					#key:   chromosome eg. 1,2,X,Y
					#value: a reference to an array where each element corresponds to a bin along the
					#	chromosome

	my $CURRENT_FILE;	#file handle to the current file that is open
	my $TP53_FILE;		#file handle to the TP53 insertion mutation output file

	my $line;		#stores raw line read in from file
	my @line_data;		#stores tokenized line read in from file

	my $chr;

	my $file_name;
	my $path;
	my $suffix;

	my $rm_insertion_file_result;

	my $insertion_found = 0;

	my %chromosome_insertion_count = (); 	#hash {chr}{count}
						#key: chromosome eg. 1,2,X,Y
						#value: the number of insertions found on key

	my %genome_trans_insertion_breakpoints = (	#hash
					 X => [],	#key:	chromosome eg. 1,2,X,Y
					 Y => [],	#value: an array storing all the insertion start positions on key
					 1 => [],
					 2 => [],
					 3 => [],
					 4 => [],
					 5 => [],
					 6 => [],
					 7 => [],
					 8 => [],
					 9 => [],
					 10 => [],
					 11 => [],
					 12 => [],
					 13 => [],
					 14 => [],
					 15 => [],
					 16 => [],
					 17 => [],
					 18 => [],
					 19 => [],
					 20 => [],
					 21 => [],
					 22 => []
					);

	#Create TP53 directory and output folder
	mkdir ("$output_directory"."TP53",0770) unless (-d "$output_directory"."TP53");
	if(!(-e "$output_directory"."TP53")){
		die "ERROR: could not create folder $output_directory"."TP53\n";
		}

	#for each file in the insertion input file array
	foreach my $file (@insertion_files){

		#Parse the file name, path and file type
		( $file_name, $path, $suffix ) = File::Basename::fileparse( $file, "\.[^.]*");

		#open the file
		open ($CURRENT_FILE, "<", $file) or die "ERROR: could not open file at path $file\n";

		#ensure that the file is not empty
		if(eof($CURRENT_FILE)){
			die "ERROR: $file is empty\n";
			}

		#create the TP53 insertion mutation output file
		open ($TP53_FILE, ">", "$output_directory"."TP53/$file_name"."$suffix") or die "ERROR: could not create file: $output_directory"."TP53/$file_name"."$suffix";

		#read header lines
		$line = <$CURRENT_FILE>;
		chomp($line);

		#print the VCF header lines to the TP53 insertion mutation output file
		while ($line =~ m/^#(.*?)/){
			print $TP53_FILE "$line\n";
			$line = <$CURRENT_FILE>;
			chomp($line);
			}

		#read all the data lines in the file
		while(1){
			@line_data = (split (/\t/,$line));

			#verify that the chromosome is valid and that the mutation is an insertion type
			if(
			   ( !($line_data[1] =~ m/^[0-9]+$/) 				 ) ||
			   ( !($line_data[0] =~ m/^(chr)?(1[0-9]|2[0-2]|X|Y|x|y|[1-9])/) ) ||
			   ( length($line_data[4]) <= length($line_data[3])		 )
			  ){
				warn "ERROR: invalid chromosome or non-insertion VCF data line found and skipped:\t$line\n";
				$line = <$CURRENT_FILE>;
				unless($line){last;}
				chomp($line);
				next;
				}

			#parse the chromosome
			$chr = $2;

			#change to uppercase if 'x' or 'y' is found
			if($chr eq 'x'){
				$chr = 'X';
				}
			if($chr eq 'y'){
				$chr = 'Y';
				}

			#increment the insertion count of the chromosome
			$chromosome_insertion_count{$chr}++;

			#check if a bin exists at the insertion start position
			#if one does not, then create one
			if(!defined(@{$genome_insertion_data{$chr}}[int($line_data[1]/$bin_size)])){
				${$genome_insertion_data{$chr}}[int($line_data[1]/$bin_size)] = 1;
				}
			else{	#if one does, then increment the count
				${$genome_insertion_data{$chr}}[int($line_data[1]/$bin_size)]++;
				}

			#Search through the list of translocation breakpoints on the same chromosome
			foreach my $bp (@{$genome_trans_breakpoints{$chr}}){
				#if the insertion is within 10bps of the breakpoints and the insertion to the list stored in
				#the genome_trans_insertion_breakpoints hash
				if( $line_data[1] < $bp+10 && $line_data[1] > $bp-10){
					push (@{$genome_trans_insertion_breakpoints{$chr}}, $bp);
					}
				}

			#Check if the insertion was in the TP53 gene
			if(
			   ( $chr ne 'X' && $chr ne 'Y' ) &&
			   ( $chr==17 )
			  ){
				if($line_data[1] >= $TP53_start && $line_data[1] <= $TP53_end){
					$$tp53_mutation_found_ref = 1; #if a mutation was found in the region set the TP53 mutated flag
					$insertion_found = 1;
					print $TP53_FILE "$line\n"; 	#print the culprit data line to the TP53 insertion
									#mutation output file
					}
				}

			#read the next line
			$line = <$CURRENT_FILE>;
			#check that the end of file has not been reached
			unless($line){last;}
			chomp($line);
			}#while(1)

		#close file
		close($CURRENT_FILE);
		close ($TP53_FILE);

		#if an insertion was not found in the current file then delete the created TP53 insertion mutation output file
		if($insertion_found!=1){
			my $dir = "$output_directory"."TP53/$file_name"."$suffix";
			$rm_insertion_file_result = `rm $dir`;
			}

		$insertion_found = 0;
		}#foreach my $file (@insertion_files)

	#return hash
	return(\%genome_insertion_data, \%chromosome_insertion_count, \%genome_trans_insertion_breakpoints);

	}#sub analyze_insertion_data


#=head2 Sub-Method: analyze_loh_data

### analyze_loh_data ##############################################################################
# Description:
#	Reads data from files located in the LOH input directory and populates:
#		$chromosome_loh_breakpoints_hash_ref
#
# Input variables:
#	$output_directory:	stores the path to the output directory
#	$loh_files_array_ref:	reference to array containing all the LOH input files
#	$tp53_mutation_found_ref: reference to tp53 mutation found flag
#

#=cut

sub analyze_loh_data {

	#parse the parameters
	my $output_directory = shift;
	my $loh_files_array_ref = shift;
	my @loh_files = @$loh_files_array_ref;

	my $tp53_mutation_found_ref = shift;

	my @file_data;	#an array storing all the entries from every input file

	my $CURRENT_FILE;	#file handle to the current file that is open
	my $TP53_FILE;		#file handle to the TP53 translocation mutation output file

	my $line;		#stores raw line read in from file
	my @line_data;		#stores tokenized line read in from files

	my %chromsome_loh_breakpoints = (		#hash {chr}[start and end pairs]
					 X => [],	#key: 	chromosome eg. 1,2,X,Y
					 Y => [],	#value:	an array that stores all the LOH breakpoints on key
					 1 => [],
					 2 => [],
					 3 => [],
					 4 => [],
					 5 => [],
					 6 => [],
					 7 => [],
					 8 => [],
					 9 => [],
					 10 => [],
					 11 => [],
					 12 => [],
					 13 => [],
					 14 => [],
					 15 => [],
					 16 => [],
					 17 => [],
					 18 => [],
					 19 => [],
					 20 => [],
					 21 => [],
					 22 => []
					);

	#Read the contents of the cnv files into memory
	foreach my $file (@loh_files){
		#open the file
		open ($CURRENT_FILE, "<", $file) or die "ERROR: could not open file at path $file\n";

		#Ensure that the file is not empty
		if(eof($CURRENT_FILE)){
			close ($CURRENT_FILE);
			die "ERROR: $file is empty\n";
			}

		#read header line and validate
		$line = <$CURRENT_FILE>;
		chomp($line);

		#Validate the header line
		if(!($line =~ m/^#chr\tstart\tend\tquality$/)){
			close ($CURRENT_FILE);
			die "ERROR: header of loh file $file is invalid\n";
			}

		#Read all the data lines
		while( !(eof($CURRENT_FILE)) ){

			#read data line
			$line = <$CURRENT_FILE>;
			chomp($line);

			#validate the data line
			if(!($line =~ m/^(chr)?(1[0-9]|2[0-2]|X|Y|[1-9])\t[0-9]+\t[0-9]+\t([0-9]+|\.)$/)){
				die "ERROR: invalid line found ($line) in file $file\n";
				}

			#Split the data line and add it to the array
			@line_data = (split (/\t/,$line));
			push(@file_data,[@line_data]);
			}

		close ($CURRENT_FILE);

		}#foreach my $file (@cnv_files)

	#Ensure that there are no overlapping LOH regions, or join them if the user indicated
	@file_data = check_for_overlaps("loh", \@file_data);

	#Create TP53 directory and output folder
	mkdir ("$output_directory"."TP53",0770) unless (-d "$output_directory"."TP53");
	if(!(-e "$output_directory"."TP53")){
		die "ERROR: could not create folder $output_directory"."TP53\n";
		}

	#Create the TP53 LOH mutation output data file
	open ($TP53_FILE, ">", "$output_directory"."TP53/TP53.spl") or die "ERROR: could not create file: $output_directory"."TP53/TP53.spl";
	#Print the header for the output file (same format as a .spl file)
	print $TP53_FILE "#chr\tstart\tend\tquality";

	#For every data line that was read in
	for (my $n = 0; $n < scalar(@file_data); $n++){

		#Validate the chromosome field
		if(!($file_data[$n][0] =~ m/^(chr)?(1[0-9]|2[0-2]|X|Y|[1-9])$/)){
			die "ERROR: invalid chromosome field detected\n";
			}

		#Parse the chromosome
		my $chr = $2;

		#Add the breakpoints to the array for the chromosome
		push (@{$chromsome_loh_breakpoints{$chr}}, ($file_data[$n][1],$file_data[$n][2]));


		#Check if the variation was in the TP53 gene
		if(
		   ( $chr ne 'X' && $chr ne 'Y' ) &&
		   ( $chr==17 )
		  ){
			if(
			   ( ( $file_data[$n][1] >= $TP53_start && $file_data[$n][1] <= $TP53_end ) || ( $file_data[$n][2] >= $TP53_start && $file_data[$n][2] <= $TP53_end ) )
			  ){
				#If a mutation was found in the TP53 region then set the TP53 mutated flag
				$$tp53_mutation_found_ref = 1;

				#Print the LOH data line to the TP53 LOH mutation output file
				print $TP53_FILE "\n";
				for (my $i = 0; $i < scalar(@{$file_data[$n]}); $i++){
					print $TP53_FILE "$file_data[$n][$i]";
					if($i != scalar(@{$file_data[$n]})-1){
						print $TP53_FILE "\t";
						}#if
					}#for (my $i = 0; $i < scalar(@{$file_data[$n]}); $i++)
				}#if
			}#if

		}#for (my $n = 0; $n < scalar(@file_data); $n++)

	close($TP53_FILE);

	#return hash
	return(\%chromsome_loh_breakpoints);

	}#sub analyze_loh_data



#=head2 Sub-Method: calculate_genome_localization

### calculate_genome_localization #################################################################
# Description:
#	Caculates the mutation density for each chromosome
#
# Input variables:
#	$output_directory:			  stores the path to the output directory
#	$chromosome_copy_number_count_hash_ref:	  stores a reference to the hash storing the number
#						  of CNV events on each chromosome
#	#chromosome_translocation_count_hash_ref: stores a reference to the hash storing the number
#						  of translocation events on each chromosome
#

#=cut

sub calculate_genome_localization {

	#parse the parameters
	my $output_directory = shift;
	my $chromosome_copy_number_count_hash_ref = shift;
	my $chromosome_translocation_count_hash_ref = shift;

	my %chromosome_mutation_count;	#hash
					#key:	chromosome eg. 1,2,X,Y
					#value: the density of translocation and CNV events on key

	my $density;	#store the mutation density for a chromosome

	my $OUTPUT_FILE;	#file handle to the output file

	#initialize all the counts to 0
	for (my $i=1; $i<23; $i++){
		$chromosome_mutation_count{$i} = 0;
		}
	$chromosome_mutation_count{'X'} = 0;
	$chromosome_mutation_count{'Y'} = 0;

	#add the number of CNV events on each chromosome
	for my $cnv_key1 ( keys %$chromosome_copy_number_count_hash_ref){
		for my $cnv_key2 (keys %{$chromosome_copy_number_count_hash_ref->{$cnv_key1}}){
			$chromosome_mutation_count{$cnv_key1} += $chromosome_copy_number_count_hash_ref->{$cnv_key1}->{$cnv_key2};
			}
		}

	#add the number of translocation events on each chromosome
	for my $trans_key1 ( keys %$chromosome_translocation_count_hash_ref){
		for my $trans_key2 (keys %{$chromosome_translocation_count_hash_ref->{$trans_key1}}){
			$chromosome_mutation_count{$trans_key1} += $chromosome_translocation_count_hash_ref->{$trans_key1}->{$trans_key2};
			}
		}

	#Create the output file
	open ($OUTPUT_FILE, ">", "$output_directory/genome_localization.log") or die "ERROR: could not create file $output_directory/genome_localization.log\n";
	#Print the header
	print $OUTPUT_FILE "#chr\tcount\tdensity\n";

	#for each chromosome print the count and overall density
	{use sort 'stable';
	for my $chr ( sort keys %chromosome_mutation_count){
			$density = $chromosome_mutation_count{$chr}/$chromosome_length{$chr};

			print $OUTPUT_FILE "$chr";
			print $OUTPUT_FILE "\t";
			print $OUTPUT_FILE $chromosome_mutation_count{$chr};
			print $OUTPUT_FILE "\t";
			print $OUTPUT_FILE "$density";
			print $OUTPUT_FILE "\n";

			#Replace the count with the density
			$chromosome_mutation_count{$chr} = $density;
			}
	}#use sort 'stable'

	close ($OUTPUT_FILE);

	#return the hash containing the densities
	return(\%chromosome_mutation_count);

	}#sub calculate_genome_localization


#=head2 Sub-Method: calculate_chromosome_localization

### calculate_chromosome_localization #############################################################
# Description:
#	Performs a sliding window analysis on the CNV and translocation data. Identifies regions
#	that have a density of mutation much greater than the average rate of mutation of the
#	genome.
#
# Input variables:
#	$output_directory:		stores the directory where output files are created
#	$genome_cnv_data_hash_ref:	reference to hash that stores position of all CNV breakpoints in
#					the genome
#	$genome_trans_data_hash_ref:	reference to hash that stores position of all the
#					translocation breakpoints in the genome
#	$bin_size:			size of the bins that divide up the genome
#	$window_size:			number of bins to evaluate in each window
#

#=cut
sub calculate_chromosome_localization {

	#parse parameters
	my $output_directory = shift;
	my $genome_cnv_data_hash_ref = shift;
	my $genome_trans_data_hash_ref = shift;
	my $bin_size = shift;
	my $window_size = shift;

	my @suspect_regions;	#array storing the start position, end position and chromosome
				#of very highly mutated regions
	my @likely_regions;	#array storing the start position, end position and chromosome
				#of somewhat highly mutated regions

	my $in_suspect_region = 0;	#flag variables used in identifying highly mutated regions
	my $in_likely_region = 0;

	my $suspect_chr = -1;
	my $suspect_start = -1;
	my $suspect_end = -1;

	my $likely_chr = -1;
	my $likely_start = -1;
	my $likely_end = -1;

	my %genome_cnv_data_windows = 	(		#hash
					 X => [],	#key: chromosome eg. 1,2,X,Y
					 Y => [],	#value: an array storing the count of CNVs
					 1 => [],	#	in each window along the chromosome
					 2 => [],
					 3 => [],
					 4 => [],
					 5 => [],
					 6 => [],
					 7 => [],
					 8 => [],
					 9 => [],
					 10 => [],
					 11 => [],
					 12 => [],
					 13 => [],
					 14 => [],
					 15 => [],
					 16 => [],
					 17 => [],
					 18 => [],
					 19 => [],
					 20 => [],
					 21 => [],
					 22 => []
					);

	my %genome_trans_data_windows = (		#hash
					 X => [],	#key:	chromosome eg. 1,2,X,Y
					 Y => [],	#value:	an array storing a count of translocation
					 1 => [],	#	in each window along the chromosome
					 2 => [],
					 3 => [],
					 4 => [],
					 5 => [],
					 6 => [],
					 7 => [],
					 8 => [],
					 9 => [],
					 10 => [],
					 11 => [],
					 12 => [],
					 13 => [],
					 14 => [],
					 15 => [],
					 16 => [],
					 17 => [],
					 18 => [],
					 19 => [],
					 20 => [],
					 21 => [],
					 22 => []
					);

	my %genome_mutation_data_windows = (		#hash
					 X => [],	#key:	chromosome eg. 1,2,X,Y
					 Y => [],	#value:	an array storing a count of all mutations
					 1 => [],	#	in each window along the chromosome
					 2 => [],
					 3 => [],
					 4 => [],
					 5 => [],
					 6 => [],
					 7 => [],
					 8 => [],
					 9 => [],
					 10 => [],
					 11 => [],
					 12 => [],
					 13 => [],
					 14 => [],
					 15 => [],
					 16 => [],
					 17 => [],
					 18 => [],
					 19 => [],
					 20 => [],
					 21 => [],
					 22 => []
					);


	my $current_chr;		#current chromosome being analyzed
	my @current_chr_data = ();	#array storing the bins for the current chromosome

	my $genome_mean_mutation_density = 0;	#average density of all the windows across the genome
	my $total_genome_windows = 0;		#total number of windows across the genome
	my $genome_mutation_density_standard_deviation = 0;	#standard deviation of the mutation densities for
								#all the windows

	my $OUTPUT_FILE;	#file handle to output file

	$output_directory = $output_directory."mutation_clustering";

	#create output directories
	mkdir ("$output_directory",0770) unless (-d "$output_directory");
	if(!(-e "$output_directory")){
		die "ERROR: could not create folder $output_directory\n";
		}

	mkdir ("$output_directory/cnv",0770) unless (-d "$output_directory/cnv");
	if(!(-e "$output_directory")){
		die "ERROR: could not create folder $output_directory/cnv\n";
		}

	mkdir ("$output_directory/translocations",0770) unless (-d "$output_directory/translocations");
	if(!(-e "$output_directory")){
		die "ERROR: could not create folder $output_directory/translocations\n";
		}

	mkdir ("$output_directory/all_types",0770) unless (-d "$output_directory/all_types");
	if(!(-e "$output_directory")){
		die "ERROR: could not create folder $output_directory/all_types\n";
		}

	#compute the density of CNV mutations in each window
	for my $cnv_key ( keys %$genome_cnv_data_hash_ref){
		#get the array storing the CNV bins for the current chromosome
		@current_chr_data = @{$genome_cnv_data_hash_ref->{$cnv_key}};

		#check that the array is not empty
		if(scalar(@current_chr_data) > 0){

			#create an output file for this chromosome
			open ($OUTPUT_FILE, ">", "$output_directory/cnv/chr$cnv_key"."_cnv_localization.log") or die "ERROR: could not create file $output_directory/cnv/chr$cnv_key"."_cnv_localization.log";
			#print the header for the output file
			print $OUTPUT_FILE "#chr\tstart\tend\tdensity";

			@{$genome_cnv_data_windows{$cnv_key}}[0] = 0;	#initialize the count in the first window

			#Calculate the mutation count in the first window
			for(my $chr_pos = 0; $chr_pos < $window_size; $chr_pos++){
				my %region_hash;
				if(!defined($current_chr_data[$chr_pos])){
					next;
					}
				%region_hash = %{$current_chr_data[$chr_pos]};

				@{$genome_cnv_data_windows{$cnv_key}}[0] += $region_hash{'BPcount'};
				}

			#print the values from the first window to the output file
			print $OUTPUT_FILE "\n";
			print $OUTPUT_FILE "$cnv_key";
			print $OUTPUT_FILE "\t";
			print $OUTPUT_FILE "0";
			print $OUTPUT_FILE "\t";
			print $OUTPUT_FILE ($window_size)*$bin_size;
			print $OUTPUT_FILE "\t";

			if(!defined(@{$genome_cnv_data_windows{$cnv_key}}[0])){
				print $OUTPUT_FILE "0";
				}
			else{
				my $rounded = POSIX::ceil((@{$genome_cnv_data_windows{$cnv_key}}[0])/2);
				print $OUTPUT_FILE ($rounded)/($window_size*$bin_size);
				#add the cnv count to the total mutation count for the region
				@{$genome_mutation_data_windows{$cnv_key}}[0] += $rounded;
				}

			#perform the sliding window analysis for the rest of the chromosome
			for(my $chr_pos = 1; $chr_pos < scalar(@current_chr_data); $chr_pos++){

				#check that the window will not overshoot the length of the chromosome
				if( (($chr_pos+($window_size-1))*$bin_size) > $chromosome_length{$cnv_key}  ){
					last;
					}

				@{$genome_cnv_data_windows{$cnv_key}}[$chr_pos] = 0;	#initialize the count for the current window

				my %past_region_hash;
				my %next_region_hash;

				my $prev_value = 0;
				my $next_value = 0;

				#get the count of the from the first bin from the previous window
				if(defined($current_chr_data[$chr_pos-1])){
					%past_region_hash = %{$current_chr_data[$chr_pos-1]};
					$prev_value = $past_region_hash{'BPcount'};
					}

				#get the count from the bin following the last bin in the previous window
				if(defined($current_chr_data[$chr_pos+($window_size-1)])){
					%next_region_hash = %{$current_chr_data[$chr_pos+($window_size-1)]};
					$next_value = $next_region_hash{'BPcount'};
					}

				#the count for the current window = the count from the previous window - the first bin of the previous window + the next bin along the chromosome
				@{$genome_cnv_data_windows{$cnv_key}}[$chr_pos] += (@{$genome_cnv_data_windows{$cnv_key}}[$chr_pos-1]) - ($prev_value) + ($next_value);

				#print the values for this window
				print $OUTPUT_FILE "\n";
				print $OUTPUT_FILE "$cnv_key";
				print $OUTPUT_FILE "\t";
				print $OUTPUT_FILE $chr_pos*$bin_size;
				print $OUTPUT_FILE "\t";
				print $OUTPUT_FILE ($chr_pos+$window_size)*$bin_size;
				print $OUTPUT_FILE "\t";

				if(!defined(@{$genome_cnv_data_windows{$cnv_key}}[$chr_pos])){
					print $OUTPUT_FILE "0";
					}
				else{
					my $rounded = POSIX::ceil((@{$genome_cnv_data_windows{$cnv_key}}[$chr_pos])/2);
					print $OUTPUT_FILE ($rounded)/($window_size*$bin_size);
					#add the cnv count to the total mutation count for the region
					@{$genome_mutation_data_windows{$cnv_key}}[$chr_pos] += $rounded;
					}
				}#for(my $chr_pos = 1; $chr_pos < scalar(@current_chr_data); $chr_pos++)

			close ($OUTPUT_FILE);
			}#if(scalar(@current_chr_data) > 0)

		@current_chr_data = ();
		}#for my $cnv_key ( keys %$genome_cnv_data_hash_ref)

	#perform the sliding window analysis on the translocation mutation data
	for my $trans_key ( keys %$genome_trans_data_hash_ref){

		#get the array storing the translocation bins for the current chromosome
		@current_chr_data = @{$genome_trans_data_hash_ref->{$trans_key}};

		#check that the array is not empty
		if(scalar(@current_chr_data) > 0){

			#create the output file
			open ($OUTPUT_FILE, ">", "$output_directory/translocations/chr$trans_key"."_translocation_localization.log") or die "ERROR: could not create file $output_directory/translocations/chr$trans_key"."_translocation_localization.log";
			#print the header for the output file
			print $OUTPUT_FILE "#chr\tstart\tend\tdensity";

			#initialize the translocation count for the first window
			@{$genome_trans_data_windows{$trans_key}}[0] = 0;

			#calculate the translocation mutation count in the first window
			for(my $chr_pos = 0; $chr_pos < $window_size; $chr_pos++){
				my %region_hash;
				if(!defined($current_chr_data[$chr_pos])){
					next;
					}
				%region_hash = %{$current_chr_data[$chr_pos]};

				my %trans_hash_in;
				my %trans_hash_out;
				my $size = 0;

				#calculate the number of inbound translocation breakpoints
				if(defined($region_hash{'in'})){
					%trans_hash_in = %{$region_hash{'in'}};

					for my $key (keys %trans_hash_in){
						$size = @{$trans_hash_in{$key}};
						$size = $size/2;
						@{$genome_trans_data_windows{$trans_key}}[0] += $size;
						}
					}

				#calculate the number of outbound translocation breakpoints
				if(defined($region_hash{'out'})){
					%trans_hash_out = %{$region_hash{'out'}};

					for my $key (keys %trans_hash_out){
						if($key eq $trans_key){
							next;
							}
						$size = @{$trans_hash_out{$key}};
						$size = $size/2;
						@{$genome_trans_data_windows{$trans_key}}[0] += $size;
						}
					}

				}#for(my $chr_pos = 0; $chr_pos < $window_size; $chr_pos++)

			#print the values from the first window to the output file
			print $OUTPUT_FILE "\n";
			print $OUTPUT_FILE "$trans_key";
			print $OUTPUT_FILE "\t";
			print $OUTPUT_FILE "0";
			print $OUTPUT_FILE "\t";
			print $OUTPUT_FILE ($window_size)*$bin_size;
			print $OUTPUT_FILE "\t";

			if(!defined(@{$genome_trans_data_windows{$trans_key}}[0])){
				print $OUTPUT_FILE "0";
				}
			else{
				my $rounded = POSIX::ceil(@{$genome_trans_data_windows{$trans_key}}[0]);
				print $OUTPUT_FILE ($rounded)/($window_size*$bin_size);
				#add the translocation mutation count to the total mutation count for the region
				@{$genome_mutation_data_windows{$trans_key}}[0] += $rounded;
				}

			#perform the sliding window analysis for the rest of the chromosome
			for(my $chr_pos = 1; $chr_pos < scalar(@current_chr_data); $chr_pos++){

				if( (($chr_pos+($window_size-1))*$bin_size) > $chromosome_length{$trans_key} ){
					last;
					}

				@{$genome_trans_data_windows{$trans_key}}[$chr_pos] = 0;
				my %prev_region_hash;
				my %next_region_hash;

				my $prev_value = 0;
				my $next_value = 0;

				#Caculate the number of mutations in the first bin of the previous window
				if(defined($current_chr_data[$chr_pos-1])){
					%prev_region_hash = %{$current_chr_data[$chr_pos-1]};

					my $size = 0;
					my %prev_trans_hash_in;
					my %prev_trans_hash_out;

					if(defined($prev_region_hash{'in'})){
						%prev_trans_hash_in = %{$prev_region_hash{'in'}};

						for my $key (keys %prev_trans_hash_in){
							$size = @{$prev_trans_hash_in{$key}};
							$size = $size/2;
							$prev_value += $size;
							}
					}

					if(defined($prev_region_hash{'out'})){
						%prev_trans_hash_out = %{$prev_region_hash{'out'}};

						for my $key (keys %prev_trans_hash_out){
							if($key eq $trans_key){
								next;
								}
							$size = @{$prev_trans_hash_out{$key}};
							$size = $size/2;
							$prev_value += $size;
							}
						}

					}

				#Caculate the number of mutations in the last bin of the current window
				if(defined($current_chr_data[$chr_pos+($window_size-1)])){
					%next_region_hash = %{$current_chr_data[$chr_pos+($window_size-1)]};

					my $size = 0;
					my %next_trans_hash_in;
					my %next_trans_hash_out;

					if(defined($next_region_hash{'in'})){
						%next_trans_hash_in = %{$next_region_hash{'in'}};

						for my $key (keys %next_trans_hash_in){
							$size = @{$next_trans_hash_in{$key}};
							$size = $size/2;
							$next_value += $size;
							}
					}

					if(defined($next_region_hash{'out'})){
						%next_trans_hash_out = %{$next_region_hash{'out'}};

						for my $key (keys %next_trans_hash_out){
							if($key eq $trans_key){
								next;
								}
							$size = @{$next_trans_hash_out{$key}};
							$size = $size/2;
							$next_value += $size;
							}
						}


					}

				#total number of translocation mutations in the current window = number of mutations in previous window - the first bin of the previous window + next bin along the chromosome
				@{$genome_trans_data_windows{$trans_key}}[$chr_pos] += (@{$genome_trans_data_windows{$trans_key}}[$chr_pos-1]) - ($prev_value) + ($next_value);

				#print values from this window
				print $OUTPUT_FILE "\n";
				print $OUTPUT_FILE "$trans_key";
				print $OUTPUT_FILE "\t";
				print $OUTPUT_FILE $chr_pos*$bin_size;
				print $OUTPUT_FILE "\t";
				print $OUTPUT_FILE ($chr_pos+$window_size)*$bin_size;
				print $OUTPUT_FILE "\t";

				if(!defined(@{$genome_trans_data_windows{$trans_key}}[$chr_pos])){
					print $OUTPUT_FILE "0";
					}
				else{
					my $rounded = POSIX::ceil(@{$genome_trans_data_windows{$trans_key}}[$chr_pos]);
					print $OUTPUT_FILE ($rounded)/($window_size*$bin_size);
					@{$genome_mutation_data_windows{$trans_key}}[$chr_pos] += $rounded;
					}
				}#for(my $chr_pos = 1; $chr_pos < scalar(@current_chr_data); $chr_pos++)


			close ($OUTPUT_FILE);
			}
		@current_chr_data = ();
		}

	#calculate the density of both types of mutations in each window in the genome
	for my $mutation_key ( keys %genome_mutation_data_windows){

		#check that some data exisits for the current chromosome
		if(scalar(@{$genome_mutation_data_windows{$mutation_key}}) > 0){

			#create the output file
			open ($OUTPUT_FILE, ">", "$output_directory/all_types/chr$mutation_key"."_mutation_localization.log") or die "ERROR: could not create file $output_directory/all_types/chr$current_chr"."_mutation_localization.log";
			#print the header for the output file
			print $OUTPUT_FILE "#chr\tstart\tend\tdensity";

			my $density;

			#for every bin along the chromosome
			for(my $chr_pos = 0; $chr_pos < scalar(@{$genome_mutation_data_windows{$mutation_key}}); $chr_pos++){
				$total_genome_windows++;	#increment the total number of windows in the genome

				#calculate the density of mutations in the window
				if(!defined(@{$genome_mutation_data_windows{$mutation_key}}[$chr_pos])){
					$density = 0;
					}
				else{
					$density = (@{$genome_mutation_data_windows{$mutation_key}}[$chr_pos])/($window_size*$bin_size);
					}

				#sum the density values to calculate the mean
				$genome_mean_mutation_density += $density;

				#print the values for the window
				print $OUTPUT_FILE "\n";
				print $OUTPUT_FILE "$mutation_key";
				print $OUTPUT_FILE "\t";
				print $OUTPUT_FILE $chr_pos*$bin_size;
				print $OUTPUT_FILE "\t";
				print $OUTPUT_FILE ($chr_pos+$window_size)*$bin_size;
				print $OUTPUT_FILE "\t";
				print $OUTPUT_FILE $density;

				}#for(my $chr_pos = 0; $chr_pos < scalar(@{$genome_mutation_data_windows{$mutation_key}}); $chr_pos++)

			close ($OUTPUT_FILE);
			}
		}#for my $mutation_key ( keys %genome_mutation_data_windows)

	#calculate the mean mutation density for the windows in the genome
	$genome_mean_mutation_density = $genome_mean_mutation_density/$total_genome_windows;

	#find the sum of squared difference between the density of mutation of each window and the mean density of mutation
	for my $mutation_key ( keys %genome_mutation_data_windows){
		if(scalar(@{$genome_mutation_data_windows{$mutation_key}}) > 0){
			my $density;

			for(my $chr_pos = 0; $chr_pos < scalar(@{$genome_mutation_data_windows{$mutation_key}}); $chr_pos++){
				if(!defined(@{$genome_mutation_data_windows{$mutation_key}}[$chr_pos])){
					$density = 0;
					}
				else{
					$density = (@{$genome_mutation_data_windows{$mutation_key}}[$chr_pos])/($window_size*$bin_size);
					}

				#sum the squared differences
				$genome_mutation_density_standard_deviation += ($density-$genome_mean_mutation_density)**2;
				}
			}
		}#for my $mutation_key ( keys %genome_mutation_data_windows)

	#divided the sum of the squared differnces by the total number of windows and take the square root
	$genome_mutation_density_standard_deviation = ($genome_mutation_density_standard_deviation/$total_genome_windows)**0.5;

	#calculate z scores for each window and check if the window is greater than 2 SDs away from genome mean
	#use this value to identify highly mutated regions
	for my $mutation_key ( keys %genome_mutation_data_windows){
		if(scalar(@{$genome_mutation_data_windows{$mutation_key}}) > 0){
			my $density;
			my $region_z_score = 0;

			for(my $chr_pos = 0; $chr_pos < scalar(@{$genome_mutation_data_windows{$mutation_key}}); $chr_pos++){

				if(!defined(@{$genome_mutation_data_windows{$mutation_key}}[$chr_pos])){
					$density = 0;
					}
				else{
					$density = (@{$genome_mutation_data_windows{$mutation_key}}[$chr_pos])/($window_size*$bin_size);
					}

				#calculate z score for the window
				$region_z_score = ($density-$genome_mean_mutation_density)/$genome_mutation_density_standard_deviation;

				#check if the z score is above the threshold
				if( $region_z_score >= $outlier_deviation ) {
					if($in_suspect_region!=1){
						$suspect_start = $chr_pos*$bin_size;
						$in_suspect_region = 1;
						}
					$suspect_chr = $mutation_key;
					$suspect_end = ($chr_pos+$window_size)*$bin_size;
					}
				elsif ($in_suspect_region==1){
					#once a region has been called push the chromosome, start and end positions into the suspect region array
					push (@suspect_regions, ($suspect_chr,$suspect_start,$suspect_end));
					$suspect_chr = -1;
					$suspect_start = -1;
					$suspect_end = -1;
					$in_suspect_region = 0;
					}

				#check if the z score is below the threshold but still suspicously high
				if(
				   ( $region_z_score < $outlier_deviation ) &&
				   ( $region_z_score >= ($outlier_deviation-1) )
				   ){
					if($in_likely_region!=1){
						$likely_start = $chr_pos*$bin_size;
						$in_likely_region = 1;
						}
					$likely_chr = $mutation_key;
					$likely_end = ($chr_pos+$window_size)*$bin_size;
					}
				elsif ($in_likely_region==1){
					push (@likely_regions, ($likely_chr,$likely_start,$likely_end));
					$likely_chr = -1;
					$likely_start = -1;
					$likely_end = -1;
					$in_likely_region = 0;
					}

				}#for(my $chr_pos = 0; $chr_pos < scalar(@{$genome_mutation_data_windows{$mutation_key}}); $chr_pos++)

				#check if the last region analyzed was suspicious and push its data into the appropriate array
				if ($in_suspect_region==1){
					push (@suspect_regions, ($suspect_chr,$suspect_start,$suspect_end));
					$suspect_chr = -1;
					$suspect_start = -1;
					$suspect_end = -1;
					$in_suspect_region = 0;
					}

				if ($in_likely_region==1){
					push (@likely_regions, ($likely_chr,$likely_start,$likely_end));
					$likely_chr = -1;
					$likely_start = -1;
					$likely_end = -1;
					$in_likely_region = 0;
					}
			}
		}#for my $mutation_key ( keys %genome_mutation_data_windows)


	return (\@suspect_regions, \@likely_regions, \%genome_cnv_data_windows, \%genome_trans_data_windows, \%genome_mutation_data_windows);

	}#sub calculate_chromosome_localization


#=head2 Sub-Method: check_copy_number_count

### check_copy_number_count #######################################################################
# Description:
#	Produces an output file that records the number of regions of copy-number variation that
#	are present in each chromosome.
#
# Input variables:
#	$output_directory:			stores the path to the output directory
#	$chromosome_copy_number_count_hash_ref:	reference to hash that stores the count of regions
#						of copy-number variation on each chromosome
#

#=cut

sub check_copy_number_count {

	#parse parameters
	my $output_directory = shift;
	my $chromosome_copy_number_count_hash_ref = shift;

	my $OUTPUT_FILE;	#file handle to output file

	#open the output file
	open ($OUTPUT_FILE, ">", "$output_directory/copy_number_count.log") or die "ERROR: could not create file $output_directory/copy_number_count.log\n";

	#print the header
	print $OUTPUT_FILE "#chr\tcopy_number\tnumber_of_regions";

	#for each chromosome
	#print out the number of regions with the given copy-number
	{use sort 'stable';
	for my $chr (sort keys %$chromosome_copy_number_count_hash_ref){
		my %intermediate_hash = %{$chromosome_copy_number_count_hash_ref->{$chr}};
		for my $CN (sort {$a <=> $b} keys %intermediate_hash){
			print $OUTPUT_FILE "\n";
			print $OUTPUT_FILE "$chr";	#chromosome
			print $OUTPUT_FILE "\t";
			print $OUTPUT_FILE "$CN";	#copy-number
			print $OUTPUT_FILE "\t";
			print $OUTPUT_FILE $chromosome_copy_number_count_hash_ref->{$chr}->{$CN};	#number of regions with copy-number $CN
			}
		}
	}#use sort 'stable'
	close ($OUTPUT_FILE);

	}#sub check_copy_number_count

#=head2 Sub-Method: check_copy_number_switches

### check_copy_number_switches ####################################################################
# Description:
#	Creates an output file that records the number of breakpoints between CNV regions on each
#	chromosome
#
# Input variables:
#	$output_directory:			stores path to output directory
#	$chromosome_copy_number_count_hash_ref:	reference to hash that stores the count of regions
#						of copy-number variation on each chromosome
#

#=cut

sub check_copy_number_switches {

	#parse parameters
	my $output_directory = shift;
	my $chromosome_copy_number_count_hash_ref = shift;

	my $switch_count = 0;

	my $OUTPUT_FILE;	#file handle to output file

	#open output file
	open ($OUTPUT_FILE, ">", "$output_directory/copy_number_switches.log") or die "ERROR: could not create file $output_directory/copy_number_switches.log\n";

	#print header
	print $OUTPUT_FILE "#chr\tswitch_count";

	#for each chromosome
	#sum the total number of CNV events
	#multiply the sum by 2 to get the number of CNV breakpoints
	{use sort 'stable';
	for my $chr (sort keys %$chromosome_copy_number_count_hash_ref){
		my %intermediate_hash = %{$chromosome_copy_number_count_hash_ref->{$chr}};
		for my $CN (sort {$intermediate_hash{$b} <=> $intermediate_hash{$a} } keys %intermediate_hash){
			#only count regions that have abberant copy numbers
			if($CN != 2){
				$switch_count += ($chromosome_copy_number_count_hash_ref->{$chr}->{$CN} * 2);
				}
			}

		#print values to output file
		print $OUTPUT_FILE "\n";
		print $OUTPUT_FILE "$chr";
		print $OUTPUT_FILE "\t";
		print $OUTPUT_FILE $switch_count;

		$switch_count = 0;
		}
	}#use sort 'stable'

	close ($OUTPUT_FILE);

	}#sub check_copy_number_switches

#=head2 Sub-Method: calculate_interchromosomal_translocation_rate

### calculate_interchromosomal_translocation_rate #################################################
# Description:
#	Create an output file that records the number of translocations between each and every
#	chromosome
#
# Input variables:
#	$output_directory:				stores path to output directory
#	$chromosome_translocation_count_hash_ref:	reference to hash that stores the count of
#							translocations between each chromosome
#

#=cut

sub calculate_interchromosomal_translocation_rate {

	#parse parameters
	my $output_directory = shift;
	my $chromosome_translocation_count_hash_ref = shift;

	my $OUTPUT_FILE;	#file handle to output file

	#open output file
	open ($OUTPUT_FILE, ">", "$output_directory/interchromosomal_translocation_rate.log") or die "ERROR: could not create file $output_directory/interchromosomal_translocation_rate.log\n";

	#print header
	print $OUTPUT_FILE "#chr1\tchr2\tcount";

	#for each chromosome
	#print the number of translocations between every other chromosome
	{use sort 'stable';
	for my $chr1 (sort keys %$chromosome_translocation_count_hash_ref){
		my %intermediate_hash = %{$chromosome_translocation_count_hash_ref->{$chr1}};

		for my $chr2 (sort keys %intermediate_hash){
			print $OUTPUT_FILE "\n";
			print $OUTPUT_FILE $chr1;
			print $OUTPUT_FILE "\t";
			print $OUTPUT_FILE $chr2;
			print $OUTPUT_FILE "\t";
			print $OUTPUT_FILE $chromosome_translocation_count_hash_ref->{$chr1}->{$chr2};
			}
		}
	}#use sort 'stable'
	close ($OUTPUT_FILE);

	}#sub calculate_interchromosomal_translocation_rate


#=head2 Sub-Method: analyze_suspect_regions

### analyze_suspect_regions #######################################################################
# Description:
#	Produces the final report output file, that includes the chromothriptic scores for each of
#	the highly mutated regions
#
# Input variables:
#	$output_directory:				stores path to the output directory
#	$suspect_regions_array_ref:			reference to array storing the chromosome,
#							start, and end position of highly mutated
#							regions
#	$genome_mutation_density_hash_ref:		stores the average mutation density of each
#							chromosome
#	$genome_cnv_data_hash_ref:			stores the position of CNV mutations on
#							each chromosome
#	$genome_trans_data_hash_ref:			stores the position of translocation events
#							on each chromosome
#	$genome_trans_insertion_breakpoints_hash_ref:	stores the position of insertions on each
#							chromosome
#	$bin_size:					stores the size of a single bin
#	$localization_window_size:			stores the number of bins to include in a
#							window
#	$tp53_mutated:					stores whether the TP53 gene is mutatated
#							or not
#	$tp53_mutation_found:				stores whether or not a mutation was found
#							in the TP53 loci
#	$chromosome_cnv_breakpoints_hash_ref:		stores the breakpoints of CNV mutations on
#							each chromosome
#	$chromosome_loh_breakpoints_hash_ref:		stores the breakpoints of LOH regions on
#							each chromosome
#

#=cut

sub analyze_suspect_regions {

	#parse parameters
	my $output_directory 				= shift;

	my $suspect_regions_array_ref 			= shift;
	my @suspect_regions 				= @$suspect_regions_array_ref;

	my $genome_mutation_density_hash_ref 		= shift;
	my %genome_mutation_density_hash		= %{$genome_mutation_density_hash_ref};

	my $genome_cnv_data_hash_ref 				= shift;
	my $genome_trans_data_hash_ref				= shift;
	my $genome_trans_insertion_breakpoints_hash_ref 	= shift;

	my $bin_size 					= shift;
	my $localization_window_size			= shift;

	my $tp53_mutated				= shift;
	my $tp53_mutation_found				= shift;

	my $chromosome_cnv_breakpoints_hash_ref		= shift;
	my $chromosome_loh_breakpoints_hash_ref		= shift;


	my $suspect_regions_size = @suspect_regions;

	my $OUTPUT_FILE;

	my @suspect_region_data = ();

	my $header_string;

	#check that the suspect region data array is not malformed, should contain sets of 3 elements
	if($suspect_regions_size % 3 != 0){
		die "ERROR: suspect_regions_array has $suspect_regions_size entries. Value must be divisible by 3.\n";
		}

	#create output directory for report files
	mkdir ("$output_directory"."suspect_regions",0770) unless (-d "$output_directory"."suspect_regions");
	if(!(-e "$output_directory"."suspect_regions")){
		die "ERROR: could not create folder $output_directory"."suspect_regions";
		}

	#open final report output file
	open ($OUTPUT_FILE, ">", "$output_directory"."suspect_regions/suspect_regions.yml") or die "ERROR: could not create file: $output_directory"."suspect_regions/suspect_regions.yml\n";

	#construct and print header
	$header_string = "file: Suspect Chromothriptic Regions\n";
	$header_string .= "bin_size:\t\t\t$bin_size\n";
	$header_string .= "localization_window_size:\t$localization_window_size\n";
	$header_string .= "\n";
	$header_string .= "genome_localization_score_weight:\t$genome_localization_weight\n";
	$header_string .= "chromosome_localization_score_weight:\t$chromosome_localization_weight\n";
	$header_string .= "cnv_score_weight:\t\t\t$cnv_weight\n";
	$header_string .= "translocation_score_weight:\t\t$translocation_weight\n";
	$header_string .= "insertion_breakpoint_score_weight:\t$insertion_breakpoint_weight\n";
	$header_string .= "loh_score_weight:\t\t\t$loh_weight\n";
	$header_string .= "tp53_mutation_score_weight:\t\t$tp53_mutated_weight\n";
	$header_string .= "\n";
	$header_string .= "min_mutation_density_z_score:\t$outlier_deviation\n";
	$header_string .= "---\n";
	$header_string .= "\n";
	print $OUTPUT_FILE $header_string;

	#calculate a chromothripsis score for each region that was present in the suspect region array
	#store the results of the score calculation in a 2d array where elements in the first dimension correspond to each suspect region
	#and the elements in the second dimension are the results of the score calculation
	for (my $i = 0; $i < $suspect_regions_size; $i+=3){
		my @region_data = ();
		$region_data[0]   = $suspect_regions[$i];	#chr
		$region_data[1] = $suspect_regions[$i+1];	#start
		$region_data[2]   = $suspect_regions[$i+2]; 	#end

		($region_data[3], $region_data[4], $region_data[5], $region_data[6], $region_data[7], $region_data[8], $region_data[9], $region_data[10], $region_data[11], $region_data[12], $region_data[13]) = calculate_score($region_data[0], $region_data[1], $region_data[2], $genome_cnv_data_hash_ref, $genome_trans_data_hash_ref, $genome_mutation_density_hash_ref, $genome_trans_insertion_breakpoints_hash_ref, $tp53_mutated, $chromosome_cnv_breakpoints_hash_ref, $chromosome_loh_breakpoints_hash_ref, $bin_size);

		#add the results of the score calculation for this region to the array storing all the results
		push @suspect_region_data, [@region_data];
		}

	#sort the results so that the region with the highest chromothriptic score will be printed
	#to the final report output file first
	{use sort 'stable';
	@suspect_region_data = sort {$b->[3] <=> $a->[3] } @suspect_region_data;
	}#use sort 'stable'A

	#for each score that is generated print the score and the related statistics for that region
	foreach my $score_data (@suspect_region_data){
		my $chr   = $score_data->[0];	#chr
		my $start = $score_data->[1];	#start
		my $end   = $score_data->[2]; 	#end

		my $score			= sprintf("%.5f",$score_data->[3]);
		my $chr_z_score			= $score_data->[4];
		my $region_density		= sprintf("%e",$score_data->[5]);

		my $cnv_number_hash_ref		= $score_data->[6];
		my %cnv_number_hash;
		my $num_copy_num;
		if(defined($cnv_number_hash_ref)){
			%cnv_number_hash = %{$cnv_number_hash_ref};
			$num_copy_num = keys %cnv_number_hash;
			}
		else{
			$num_copy_num = 0;
			}

		my $cnv_density			= sprintf("%e",$score_data->[7]);

		my $intertranslocation_hash_ref	= $score_data->[8];
		my $translocation_density	= $score_data->[9];
		my %intertranslocation_hash;
		my $num_trans_chr;
		if(defined($intertranslocation_hash_ref)){
			%intertranslocation_hash = %{$intertranslocation_hash_ref};
			$num_trans_chr = keys %intertranslocation_hash;
			}
		else{
			$num_trans_chr = 0;
			}

		my $breakpoint_insertions_array_ref	= $score_data->[10];
		my @breakpoint_insertions_array;
		my $breakpoint_percentage;
		if(defined($breakpoint_insertions_array_ref)){
			@breakpoint_insertions_array = @$breakpoint_insertions_array_ref;
			$breakpoint_percentage	= sprintf("%.2f",($breakpoint_insertions_array[0]/$breakpoint_insertions_array[1])*100);
			}

		my $loh_size = $score_data->[11];
		my $hz_size = $score_data->[12];
		my $percent_hz_lost;
		if(defined($loh_size) && defined($hz_size)){
			$percent_hz_lost = sprintf("%.2f",($loh_size/$hz_size)*100);
			}

		my @score_array = @{$score_data->[13]};

		my $chr_density = $genome_mutation_density_hash{$chr};

		my $print_string;

		$print_string = "chromosome:\t$chr\n";
		$print_string .= "start:\t\t$start\n";
		$print_string .= "end:\t\t$end\n";
		$print_string .= "\n";

		$print_string .= "final_score:\t\t\t$score\n";
		$print_string .= "genome_localization_score:\t".$score_array[2]*$genome_localization_weight."\t(".$score_array[2].")"."\n";
		$print_string .= "chromosome_localization_score:\t".$score_array[1]*$chromosome_localization_weight."\t(".$score_array[1].")"."\n";
		$print_string .= "cnv_score:\t\t\t".$score_array[0]*$cnv_weight."\t(".$score_array[0].")"."\n";
		$print_string .= "translocation_score:\t\t".$score_array[3]*$translocation_weight."\t(".$score_array[3].")"."\n";
		$print_string .= "insertion_breakpoint_score:\t".$score_array[4]*$insertion_breakpoint_weight."\t(".$score_array[4].")"."\n";
		$print_string .= "loh_score:\t\t\t".$score_array[5]*$loh_weight."\t(".$score_array[5].")"."\n";
		$print_string .= "tp53_score:\t\t\t".$score_array[6]*$tp53_mutated_weight."\t(".$score_array[6].")\n";

		$print_string .= "\n";

		$print_string .= "mutation_density_of_region:\t$region_density\n";
		$print_string .= "mutation_density_of_chromosome:\t$chr_density\n";
		$print_string .= "standard_deviations_from_mean_of_chromosome_mutation_density:\t$chr_z_score\n";
		$print_string .= "\n";

		$print_string .= "density_of_copy_number_switches: $cnv_density\n";
		$print_string .= "number_of_aberrant_copy_number_states:\t$num_copy_num\n";
		if($num_copy_num>0){
			$print_string .= "aberrant_copy_number_states:\n";
			{use sort 'stable';
			foreach my $key (sort {$cnv_number_hash{$b} <=> $cnv_number_hash{$a} } keys %cnv_number_hash){
				$print_string .= "\t$key:\t$cnv_number_hash{$key}\n";
				}
			}#use sort 'stable'
			}

		$print_string .= "\n";

		$print_string .= "density_of_translocation_breakpoints: $translocation_density\n";
		$print_string .= "number_of_translocation_chromosomes:\t$num_trans_chr\n";

		if($num_trans_chr>0){
			$print_string .= "translocation_chromosomes:\n";
			{use sort 'stable';
			foreach my $key (sort {$intertranslocation_hash{$b} <=> $intertranslocation_hash{$a} } keys %intertranslocation_hash){
				$print_string .= "\t$key:\t$intertranslocation_hash{$key}\n";
				}
			}#use sort 'stable'
			}
		$print_string .= "\n";

		if(defined($breakpoint_insertions_array_ref)){
			$print_string .= "insertion_data:\n";
			$print_string .= "\tinsertions_found_at_translocation_breakpoints:\t$breakpoint_insertions_array[0]\n";
			$print_string .= "\ttotal_translocation_breakpoints:\t$breakpoint_insertions_array[1]\n";
			$print_string .= "\tpercentage:\t$breakpoint_percentage"."%\n";
			$print_string .= "\n";
			}

		if($loh_size!=-1){
			$print_string .= "loh_data:\n";
			$print_string .= "\ttotal_size_of_loh:\t$loh_size\n";
			$print_string .= "\ttotal_size_of_original_heterozygosity:\t$hz_size\n";
			$print_string .= "\tpercent_heterozygosity_lost:\t$percent_hz_lost"."%\n";
			$print_string .= "\n";
			}

		if($tp53_mutated && $tp53_mutation_found){
			$print_string .= "tp53_mutation_present:\t1 (forced and mutations found)\n";
			}
		elsif($tp53_mutated){
			$print_string .= "tp53_mutation_present:\t1 (forced)\n";
			}
		elsif($tp53_mutation_found){
			$print_string .= "tp53_mutation_present:\t1 (mutations found)\n";
			}
		else{
			$print_string .= "tp53_mutation_present:\t0\n";
			}


		$print_string .= "---\n";
		$print_string .= "\n";
		print $OUTPUT_FILE $print_string;
		$print_string = "";
		}

	print $OUTPUT_FILE "...";
	close($OUTPUT_FILE);
	}#sub analyze_suspect_regions


#=head2 Sub-Method: analyze_likely_regions

### analyze_likely_regions ########################################################################
# Description:
#	Generates an output file that lists the regions that have a mutation density that is less
#	than the outlier cut off but greater than 1 - the outlier cut off
#
# Input variables:
#	$output_directory:			stores path to the output directory
#	$likely_regions_array_ref:		reference to array storing the chromosome, start,
#						and end position of highly mutated regions
#	$genome_mutation_density_hash_ref:	stores the average mutation density of each
#						chromosome
#	$genome_cnv_data_hash_ref:		stores the position of CNV mutations on each
#						chromosome
#	$genome_trans_data_hash_ref:		stores the position of translocation events on each
#						chromosome
#	$bin_size:				stores the size of a single bin
#

#=cut

sub analyze_likely_regions {

	#parse parameters
	my $output_directory 			= shift;

	my $likely_regions_array_ref 		= shift;
	my @likely_regions 			= @$likely_regions_array_ref;

	my $genome_mutation_density_hash_ref 	= shift;
	my $genome_cnv_data_hash_ref 		= shift;
	my $genome_trans_data_hash_ref		= shift;
	my $bin_size 				= shift;

	my $likely_regions_size = @likely_regions;

	my $OUTPUT_FILE;

	my @return_vals;

	my @likely_region_data = ();	#stores start, end, chromsome and mutation density for each region

	#check that the likely region array is not malformed, should contain sets of 3 elements
	if($likely_regions_size % 3 != 0){
		die "ERROR: suspect_regions_array has $likely_regions_size entries. Value must be divisible by 3.\n";
		}


	#create output directory
	mkdir ("$output_directory"."suspect_regions",0770) unless (-d "$output_directory"."suspect_regions");
	if(!(-e "$output_directory"."suspect_regions")){
		die "ERROR: could not create folder $output_directory"."suspect_regions";
		}

	#create output file
	open ($OUTPUT_FILE, ">", "$output_directory"."suspect_regions/likely_regions.log") or die "ERROR: could not create file: $output_directory"."suspect_regions/likely_regions.log\n";

	#print file header
	print $OUTPUT_FILE "Likely Chromothriptic Regions\n";
	print $OUTPUT_FILE "High Mutation Density Z-Score:\t$outlier_deviation\n";
	print $OUTPUT_FILE "Min Mutation Density Z-Score:\t$outlier_deviation-1\n";
	print $OUTPUT_FILE "---------------------------------------\n";
	print $OUTPUT_FILE "#chr\tstart\tend\tmutation_density";

	#for each likely region calculate the mutation density for the region and store it in the likely_region_data array
	for (my $i = 0; $i < $likely_regions_size; $i+=3){
		my @region_data = ();
		$region_data[0]   = $likely_regions[$i];	#chr
		$region_data[1]   = $likely_regions[$i+1];	#start
		$region_data[2]   = $likely_regions[$i+2]; 	#end

		($region_data[3],$region_data[4]) = calculate_region_mutation_density_score($region_data[0], $region_data[1], $region_data[2], $genome_cnv_data_hash_ref, $genome_trans_data_hash_ref, $genome_mutation_density_hash_ref, $bin_size);

		push @likely_region_data, [@region_data];
		}

	#sort the regions by density, largest to smallest
	{use sort 'stable';
	@likely_region_data = sort {$b->[3] <=> $a->[3] } @likely_region_data;
	}#use sort 'stable'

	#print the density for each region to the output file
	foreach my $i (@likely_region_data){
		my $chr   = $i->[0];	#chr
		my $start = $i->[1];	#start
		my $end   = $i->[2]; 	#end

		my $region_density	= $i->[3];	#mutation density

		print $OUTPUT_FILE "\n";
		print $OUTPUT_FILE "$chr\t$start\t$end\t$region_density";
		push (@return_vals,$chr,$start,$end,$region_density);
		}

	close($OUTPUT_FILE);
	return(\@return_vals);
	}#sub analyze_likely_regions


#=head2 Sub-Method: calculate_score

### calculate_score ###############################################################################
# Description:
#	Calculates the chromothripic score for the given region. Calls sub methods to generate the
#	score for each hallmark
#
# Input variables:
#	$chr:						stores the chromosome on which the region
#							is found
#	$start:						stores the start base pair of the region
#	$end:						stores the end base pair of the region
#	$genome_cnv_data_hash_ref:			stores the position of CNV mutations on
#							each chromosome
#	$genome_trans_data_hash_ref:			stores the position of translocation events
#							on each chromosome
#	$genome_mutation_density_hash_ref:		stores the average mutation density of each
#							chromosome
#	$genome_trans_insertion_breakpoints_hash_ref:	stores the position of insertions on each
#							chromosome
#	$tp53_mutated:					stores whether the TP53 gene is mutatated
#							or not
#	$chromosome_cnv_breakpoints_hash_ref:		stores the breakpoints of CNV mutations on
#							each chromosome
#	$chromosome_loh_breakpoints_hash_ref:		stores the breakpoints of LOH regions on
#							each chromosome
#	$bin_size:					stores the size of a single bin
#

#=cut

sub calculate_score{

	#parse parameters
	my $chr 	= shift;
	my $start 	= shift;
	my $end		= shift;

	my $genome_cnv_data_hash_ref			= shift;
	my $genome_trans_data_hash_ref			= shift;
	my $genome_mutation_density_hash_ref		= shift;
	my $genome_trans_insertion_breakpoints_hash_ref	= shift;
	my $tp53_mutated				= shift;
	my $chromosome_cnv_breakpoints_hash_ref		= shift;
	my $chromosome_loh_breakpoints_hash_ref		= shift;
	my $bin_size					= shift;

	#initialize variable to store scores for each hallmark
	my $cnv_score				= 0;
	my $mutation_density_score		= 0;
	my $genome_localization_score		= 0;
	my $translocation_score			= 0;
	my $insertion_breakpoint_score		= 0;
	my $loh_score				= 0;
	my $final_score				= 0;

	my @score_array;	#array in which hallmark scores will be returned

	my $chr_mutation_density;	#stores the average mutation density of the chromosome where the region is found
	my $chr_z_score;		#stores the z_score of the mutation density of the chromosome where the region is found vs
					#all the other chromosomes
	my $cnv_number_hash_ref;	#stores a hash that contains the number of regions of each abberant copy-number
	my $cnv_density;		#stores the density of cnv mutations in the region
	my $translocation_density;	#stores the density of translocation mutations in the region
	my $mutation_density;		#stores the density of all mutations in the region
	my $intertranslocation_hash_ref;	#stores the number of translocations between all other chromosomes and the region
	my $breakpoint_insertions_array_ref;	#stores the total number of translocation breakpoints, and the number that have insertions nearby
	my $loh_size = -1;			#stores the amount of heterozygosity that was lost in the region
	my $heterozygous_size;			#stores the original amount of heterozygosity in the region

	($cnv_score, $cnv_number_hash_ref, $cnv_density) = calculate_copy_number_scores($chr, $start, $end, $genome_cnv_data_hash_ref, $bin_size);
	($genome_localization_score, $chr_z_score, $chr_mutation_density) = calculate_genome_localization_score($chr, $genome_mutation_density_hash_ref);
	($translocation_score, $intertranslocation_hash_ref, $translocation_density) = calculate_translocation_score($chr, $start, $end, $genome_trans_data_hash_ref, $bin_size);

	if(defined($genome_trans_insertion_breakpoints_hash_ref)){
		($insertion_breakpoint_score, $breakpoint_insertions_array_ref)	= calculate_insertion_breakpoint_score($chr, $start, $end, $genome_trans_data_hash_ref, $genome_trans_insertion_breakpoints_hash_ref, $bin_size);
		}

	($mutation_density, $mutation_density_score) = calculate_region_mutation_density_score($chr, $start, $end, $genome_cnv_data_hash_ref, $genome_trans_data_hash_ref, $genome_mutation_density_hash_ref, $bin_size);

	if(defined($chromosome_loh_breakpoints_hash_ref)){
		($loh_score, $loh_size, $heterozygous_size) = calculate_loh_score($chr, $start, $end, $chromosome_cnv_breakpoints_hash_ref, $chromosome_loh_breakpoints_hash_ref);
		}

	#calculate overall score for region based on hallmark weights and scores
	$final_score = ($cnv_score*$cnv_weight) + ($mutation_density_score*$chromosome_localization_weight) + ($genome_localization_score*$genome_localization_weight) + ($translocation_score*$translocation_weight) + ($insertion_breakpoint_score*$insertion_breakpoint_weight) + ($tp53_mutated*$tp53_mutated_weight) + ($loh_score*$loh_weight);

	#push the hallmark scores into the score array
	push (@score_array, ($cnv_score, $mutation_density_score, $genome_localization_score, $translocation_score, $insertion_breakpoint_score, $loh_score, $tp53_mutated));

	#return the scores and other region statistics
	return ($final_score, $chr_z_score, $mutation_density, $cnv_number_hash_ref, $cnv_density, $intertranslocation_hash_ref , $translocation_density, $breakpoint_insertions_array_ref, $loh_size, $heterozygous_size, \@score_array);
	}#sub calculate_score



#=head2 Sub-Method: calculate_copy_number_score

### calculate_copy_number_score  ##################################################################
# Description:
#	Calculates the score for the copy-number variation hallmark
#
# Input variables:
#	$chr:				stores the chromsome where the region is located
#	$start:				stores the starting location of the region
#	$end:				stores the end location of the region
#	$genome_cnv_data_hash_ref:	stores the position of CNV mutations on each chromosome
#	$bin_size:			stores the size of single bin
#

#=cut

sub calculate_copy_number_scores {

	#parse parameters
	my $chr = shift;
	my $start = shift;
	my $end = shift;

	my $genome_cnv_data_hash_ref = shift;
	my %genome_cnv_data_hash = %$genome_cnv_data_hash_ref;

	my $bin_size = shift;

	#calculate array index where the data for the first and last bins of the region are located
	my $start_index = $start / ($bin_size);
	my $end_index = $end / ($bin_size);

	my $cnv_score = 0;	#stores final score to return

	my %cnv_number_hash;	#hash
				#key: copy number eg 0,1,3,4
				#value: the number of regions with the given copy number

	my $cnv_switch_count = 0;	#number of switches between different copy numbers
	my $cnv_switch_density = 0;	#density of cnv events in the region

	my @chr_data;	#stores all the bins for the chromosome where the region is located
	my %cnv_hash;	#stores the cnv hash from each bin

	my $mean = 0;	#stores the average number of regions of aberrant copy-number
	my $SD = 0;	#stores the standard deviation of the number of regions of aberrant copy-number

	my %cnv_significant;		#hash
					#key: copy number (but only significant ones are stored)
					#value: the number of regions with the given copy number

	my $significant_count = 0;	#stores the number of unique significant copy-numbers

	#check if there is cnv data for the chromosome
	if(defined($genome_cnv_data_hash{$chr})){

		#extract the bin data for the chromosome
		@chr_data = @{$genome_cnv_data_hash{$chr}};

		#collect the data from the bins that contain the region
		for (my $i = $start_index; $i < $end_index+1; $i++){
			if(!defined($chr_data[$i])){
				next;
				}
			%cnv_hash = %{$chr_data[$i]};

			for my $key (keys %cnv_hash){
				if($key eq 'BPcount'){
					$cnv_switch_count += $cnv_hash{$key};
					}
				else{
					$cnv_number_hash{$key}+= $cnv_hash{$key};
					}
				}
			}

		#calculate the breakpoint density of cnv mutations for the region
		$cnv_switch_density = $cnv_switch_count/ ($end-$start);

		#calculate the number of cnv events in the region (half the number of breakpoints)
		for my $key (keys %cnv_number_hash){
			$cnv_number_hash{$key} = POSIX::ceil($cnv_number_hash{$key});
			$mean += $cnv_number_hash{$key};
			}

		#if no cnv mutations were found return a score of 0
		if(scalar(keys %cnv_number_hash)==0){
			$cnv_score = 0;
			return ($cnv_score, \%cnv_number_hash, $cnv_switch_density);
			}

		#calculate the mean of the number regions of each copy-number
		$mean = $mean/(scalar(keys %cnv_number_hash));

		#calculate the standard deviation of the number of regions of each copy-number
		for my $key (keys %cnv_number_hash){
			$SD += ($cnv_number_hash{$key}-$mean)**2;
			}
		$SD = $SD/(scalar(keys %cnv_number_hash));
		$SD = $SD**0.5;

		#determine which copy-numbers are significant (ie are not low out liers)
		for my $key (keys %cnv_number_hash){
			if(
			   ( $SD==0 )||
			   ( (($cnv_number_hash{$key}-$mean)/$SD) >= -1*$outlier_deviation )
			  ){
				$cnv_significant{$key} = $cnv_number_hash{$key};
				$cnv_score = $cnv_score + $cnv_significant{$key}**2;
				}
			}

		#score calculation
		$cnv_score = $cnv_score / (scalar(keys %cnv_significant));
		$cnv_score = log($cnv_score)/log(2);
		$cnv_score += 1;
		$cnv_score = 1 - (1/$cnv_score);
		$cnv_score = $cnv_score/(scalar(keys %cnv_significant));
		#$cnv_score = (1/(scalar(keys %cnv_significant)))*(0.25) + (1-(1/(1+log($cnv_score/(scalar(keys %cnv_significant)))/log(2))))*(0.75);
		}

	return ($cnv_score, \%cnv_number_hash, $cnv_switch_density);

	}#sub calculate_copy_number_scores


#=head2 Sub-Method: calculate_genome_localization_score

### calculate_genome_localization_score  ##########################################################
# Description:
#	Calculates the genome localization hallmark score
#
# Input variables:
#	$chr:					store the chromosome where the region is located
#	$genome_mutation_density_hash_ref:	stores the average mutation density of each
#						chromosome
#

#=cut

sub calculate_genome_localization_score {

	#parse parameters
	my $chr = shift;

	my $genome_mutation_density_hash_ref = shift;
	my %genome_mutation_density_hash = %{$genome_mutation_density_hash_ref};

	#read mutation density for the chromosome
	my $chr_mutation_density = $genome_mutation_density_hash{$chr};

	my $mean_density = 0;	#stores the average density of mutations for the chromosomes

	my $standard_deviation = 0;	#stores the standard deviation of the mutation densities of the chromosomes

	my $z_score = 0;	#stores the z-score for the suspect chromosome
	my $p_val = 0;		#stores the p-value calculated from the z-score for the suspect chromosome

	#sum mutation densities of all the chromosomes
	for my $key (keys %genome_mutation_density_hash){
		$mean_density += $genome_mutation_density_hash{$key}/$chromosome_length{$chr};
		}

	#calculate the mean
	$mean_density = $mean_density / 24;

	#calculate the standard deviation
	for my $key (keys %genome_mutation_density_hash){
		$standard_deviation += ((($genome_mutation_density_hash{$key}) - ($mean_density) )**2);
		}
	$standard_deviation = $standard_deviation / 24;
	$standard_deviation = ($standard_deviation)**0.5;

	#check for case where the standard deviation or mean comes back as 0
	if($mean_density == 0 || $standard_deviation == 0){
		return (0, 0, $chr_mutation_density);
		}

	#calculate the z-score for suspect chromosome
	$z_score = ($chr_mutation_density - $mean_density) / $standard_deviation;

	#calculate the p-value from the z-score
	$p_val = Statistics::Distributions::uprob($z_score);
	#only consider top tail
	$p_val = 0.5-$p_val;
	$p_val = $p_val/0.5;

	#check for case where z-score comes back as 0
	if($z_score < 0){
		$p_val = 0;
		}

	#p_val is the score for this hallmark
	return ($p_val, $z_score, $chr_mutation_density);

	}#sub calculate_genome_localization_score


#=head2 Sub-Method: calculate_region_mutation_density_score

### calculate_region_mutation_density_score  ######################################################
# Description:
#	Calculates the chromosome localization hallmark score
#
# Input variables:
#	$chr:					chromosome where the region is located
#	$start:					starting location of the region
#	$end:					end location of the region
#	$genome_cnv_data_hash_ref:		stores the position of CNV mutations on each
#						chromosome
#	$genome_trans_data_hash_ref:		stores the position of translocation events
#						on each chromosome
#	$genome_mutation_density_hash_ref:	stores the average mutation density of each
#						chromosome
#	$bin_size:				stores the size of single bin
#

#=cut

sub calculate_region_mutation_density_score {

	#parse parameters
	my $chr = shift;
	my $start = shift;
	my $end = shift;

	my $genome_cnv_data_hash_ref = shift;
	my %genome_cnv_data_hash = %{$genome_cnv_data_hash_ref};

	my $genome_trans_data_hash_ref = shift;
	my %genome_trans_data_hash = %{$genome_trans_data_hash_ref};

	my $genome_mutation_density_hash_ref = shift;
	my %genome_mutation_density_hash = %{$genome_mutation_density_hash_ref};

	my $bin_size = shift;

	#calculate array index where the data for the first and last bins of the region are located
	my $start_index = $start / $bin_size;
	my $end_index = $end / $bin_size;

	#get the mean mutation density for the suspect chromosome
	my $mean_chr_mutation_density = $genome_mutation_density_hash{$chr};

	my $chis_stat;	#stores the chi squared statistic value

	my $mutation_count = 0;		#stores the total mutation count in the region
	my $cnv_count = 0;		#stores the number of cnv breakpoints in the region
	my $trans_count = 0;		#stores the number of translocation breakpoints in the region
	my $mutation_density = 0;	#store the mutation density of the region
	my $mutation_density_score;	#stores the final score for the hallmark

	my @chr_data;	#stores the bin data for the chromosome
	my %cnv_hash;	#stores the cnv hash for each bin
	my %trans_hash;	#stores the translocation hash for each bin

	#check that there is cnv data for the chromosome
	if(defined($genome_cnv_data_hash{$chr})){
		#get the cnv data for the chromosome
		@chr_data = @{$genome_cnv_data_hash{$chr}};

		#sum the number of cnv breakpoints in the suspect region
		for (my $i = $start_index; $i < $end_index+1; $i++){
			if(!defined($chr_data[$i])){
				next;
				}

			%cnv_hash = %{$chr_data[$i]};
			$cnv_count += $cnv_hash{'BPcount'};
			}
		}

	@chr_data = ();

	#check that there is translocation data for the chromosome
	if(defined($genome_trans_data_hash{$chr})){
		#get the translocation data for the chromosome
		@chr_data = @{$genome_trans_data_hash{$chr}};

		#sum the number of translocation breakpoints in the suspect region
		for (my $i = $start_index; $i < $end_index+1; $i++){
			if(!defined($chr_data[$i])){
				next;
				}
			%trans_hash = %{$chr_data[$i]};
			$trans_count += $trans_hash{'BPcount'};
			}
		}

	#divide the count by 2 to get the number of translocation events
	$trans_count = POSIX::ceil($trans_count/2);

	#calculate the total number events in the region
	$mutation_count = $cnv_count + $trans_count;

	#calculate the mutation density for the region
	$mutation_density = $mutation_count / ($end-$start);

	#calculate the chi squared statistic
	$chis_stat = abs(((log($mutation_density)-log($mean_chr_mutation_density))**2)/(log($mean_chr_mutation_density)));

	#generate a p-value using the chi squared test
	$mutation_density_score = 1-(Statistics::Distributions::chisqrprob(1,$chis_stat));

	return ($mutation_density, $mutation_density_score);
	}#calculate_region_mutation_density_score


#=head2 Sub-Method: calculate_translocation_score

### calculate_translocation_score #################################################################
# Description:
#	Calculates the translocation hallmark score
#
# Input variables:
#	$chr:				chromosome where the region is located
#	$start:				starting location of the region
#	$end:				end location of the region
#	$genome_trans_data_hash_ref:	stores the position of translocation events on each
#					chromosome
#	$bin_size:			stores the size of single bin
#

#=cut

sub calculate_translocation_score {

	#parse parameters
	my $chr = shift;
	my $start = shift;
	my $end = shift;

	my $genome_trans_data_hash_ref = shift;
	my %genome_trans_data_hash = %{$genome_trans_data_hash_ref};

	my $bin_size = shift;

	#calculate array index where the data for the first and last bins of the region are located
	my $start_index = $start / ($bin_size);
	my $end_index = $end / ($bin_size);

	my $translocation_density = 0;	#stores the density of translocation events in the region

	my @chr_data;		#stores the bin data for the chromosome
	my %trans_breakpoints;	#hash
				#key: chromosome eg 1,2,X,Y
				#value: an array storing the position of the translocation breakpoints

	my %trans_breakpoint_spreads;	#stores the average distance between translocation breakpoints
	my @significant_chrs;		#stores a list of chromosomes that have a significant number of
					#translocation to or from the region
	my @diffs;			#stores the distance between adjacent translocation breakpoints on
					#one chromosome
	my $diff_sum;			#stores the sum of the distances
	my $diff_count;			#store the number regions between translocation breakpoints

	my %trans_number_hash;		#hash
					#key: chromosome eg 1,2,X,Y
					#value: the number of events between the chromosome and the region

	my $mean = 0;			#stores the average number of translocations between each chromosome
					#and the region
	my $SD = 0;			#stores the standard deviation of the above value

	my $weighted_sum = 0.00;	#component of the score calculation

	my $size = 0.00;		#intermediate variable used to collect sums
	my $count = 0;			#intermeidate variable

	my $translocation_score = 0;	#final hallmark score

	my $spread_factor = 0;

	my $translocation_count = 0;	#stores the total number of translocations from significant chromosomes

	#check that there is translocation data for the chromosome
	if(defined($genome_trans_data_hash{$chr})){
		#get the translocation data for the chromosome
		@chr_data = @{$genome_trans_data_hash{$chr}};

		#for each bin sum the number of translocation events
		#and record the position of breakpoints in the trans_breakpoints hash
		for (my $i = $start_index; $i < $end_index+1; $i++){
			if(!defined($chr_data[$i])){
				next;
				}
			my %trans_hash = %{$chr_data[$i]};
			my %trans_hash_in;
			my %trans_hash_out;

			#analyze the breakpoints from translocation into the region
			if(defined($trans_hash{'in'})){
				%trans_hash_in = %{$trans_hash{'in'}};

				for my $key (keys %trans_hash_in){
					$size = @{$trans_hash_in{$key}};

					#calculate the number of translocation events
					$size = $size/2;

					#add the count to the appropriate hash
					$trans_number_hash{$key} += $size;

					#add the breakpoints to the trans_breakpoints hash
					push(@{$trans_breakpoints{$key}},(@{$trans_hash_in{$key}}));
					}
				}

			#analyze the breakpoints from translocation out of the region
			if(defined($trans_hash{'out'})){
				%trans_hash_out = %{$trans_hash{'out'}};

				for my $key (keys %trans_hash_out){
					$size = @{$trans_hash_out{$key}};
					$count = $size;
					if($key eq $chr){
						foreach my $val (@{$trans_hash_out{$key}}){
							if($val > $start && $val < $end){
								$count--;
								}
							}
						}
					#calculate the number of translocation events
					$count = $count/2;

					#add the count to the appropriate hash
					$trans_number_hash{$key} += $count;

					#add the breakpoints to the trans_breakpoints hash
					push(@{$trans_breakpoints{$key}},(@{$trans_hash_out{$key}}));
					}
				}
			}

		$count = 0;

		#check that some translocation events were found else return a score of 0
		if (keys(%trans_number_hash) == 0){
			$translocation_score = 0;
			$translocation_density = 0;
			return ($translocation_score, \%trans_number_hash, $translocation_density);
			}

		for my $key (keys %trans_number_hash){
			#round the event count up since if only one breakpoint was present at the end
			#of the region we still count that as a whole translocation event
			$trans_number_hash{$key} = POSIX::ceil($trans_number_hash{$key});

			#sum the number of translocation events in the region
			$count += $trans_number_hash{$key};
			}

		#calculate the translocation density for the region
		$translocation_density = $count / ($end-$start);

		#calculate the mean and standard deviation of the number of translocation between the region
		#and each chromosome
		($SD, $mean) = standard_deviation_and_mean(\%trans_number_hash,0);

		#identify chromosomes that have a high number of translocations to or from the region and
		#add them to the significant chromosome list
		for my $key (keys %trans_number_hash){
			if(
			   ( $SD==0 )||
			   #( (($trans_number_hash{$key}-$mean)/$SD)>-2*$outlier_deviation)
			   ( (($trans_number_hash{$key}-$mean)/$SD)>-1*$outlier_deviation)
			  ){
				push (@significant_chrs, $key);
				}
			}

		#calculate the total number of translocations from significant chromosomes
		foreach my $key (@significant_chrs){
			$translocation_count += $trans_number_hash{$key};
			}

		#for each significant chromosome calculate the spread between translocation events
		foreach my $key (@significant_chrs){
			#sort the breakpoints
			{use sort 'stable';
			@{$trans_breakpoints{$key}} = sort {$a <=> $b} @{$trans_breakpoints{$key}};
			}#use sort 'stable'

			$size = @{$trans_breakpoints{$key}};
			@diffs = ();
			$diff_sum = 0;
			$diff_count = 0;

			#calculate and store the distance between adjacent breakpoints
			for (my $i = 1; $i<$size; $i++){
				push (@diffs, @{$trans_breakpoints{$key}}[$i]-@{$trans_breakpoints{$key}}[$i-1]);
				}

			#check that more that one distance was calculated
			if($size==1){
				$trans_breakpoint_spreads{$key} = 0;
				$diff_count = 1;
				}
			else{	#calculate the standard deviation and mean for the distance between breakpoints
				($SD,$mean) = standard_deviation_and_mean(\@diffs,1);

				#sum the distances that are not high outliers, indicating distance between 2 translocation
				#events and not distance between breakpoints of the same event
				foreach my $val (@diffs){
					if(
					   ( $SD==0 )||
					   ( (($val-$mean)/$SD)<$outlier_deviation )
					  ){
						$diff_sum += $val;
						$diff_count++;
						}
					}
				}

			#calculate the average spread of translocation breakpoints
			$trans_breakpoint_spreads{$key} = $diff_sum / $diff_count;

			#calculate the spread factor for the chromosome
			#my $spread_factor = (log($trans_breakpoint_spreads{$key}+1)/log(10))/((log($expected_mutation_density)/log(10))*-1);
			$spread_factor = (log($trans_breakpoint_spreads{$key}+1)/log(10))/((log($expected_mutation_density)/log(10))*-1);
			if($spread_factor==0){
			$spread_factor = 1;
				}


			#increase the weighted sum based on the number of translocation events and their spread multiplied by the proportion of translocations
			#from this specific chromosome relative to the total number of translocation events
			$weighted_sum += ($trans_number_hash{$key}/$spread_factor)*($trans_number_hash{$key}/$translocation_count);
			}#foreach my $key (@significant_chrs)

		my $t2 = (1-(1/(log(1+$weighted_sum)/log(2))));
		#final hallmark score calculation
		$size = @significant_chrs;

		#calculate second term of score
		$translocation_score = (1-(1/(log(1+$weighted_sum)/log(2))));

		#calculate first term of score
		if($size<$translocation_cut_off_count && $size>2){
			$translocation_score = (1-(0.10*($size-2)))*$translocation_score;
			}
		if($size>=$translocation_cut_off_count){
			$translocation_score = 0;
			}

		if($translocation_score>=1){
			print "score: ".$translocation_score."\n";
			print "ws: ".$weighted_sum."\n";
			print "sf: ".$spread_factor."\n";
			print "term 1: ";
			print (1/(1+(log($size)/log(4))));
			print "\n";
			print "term 2: ";
			print (1-(1/(log($weighted_sum)/log(2))));
			print "\n";
			}

		}#if(defined($genome_trans_data_hash{$chr}))

	return ($translocation_score, \%trans_number_hash, $translocation_density);

	}#sub calculate_translocation_score


#=head2 Sub-Method: calculate_insertion_breakpoint_score

### calculate_insertion_breakpoint_score ##########################################################
# Description:
#	Calculates the insertions at translocation breakpoints hallmark score
#
# Input variables:
#	$chr:						chromosome where the region is located
#	$start:						starting location of the region
#	$end:						end location of the region
#	$genome_trans_data_hash_ref:			stores the position of translocation events
#							on each chromosome
#	$genome_trans_insertion_breakpoints_hash_ref:	stores the position of breakpoints with
#							insertions nearby
#	$bin_size:					stores the size of single bin
#

#=cut

sub calculate_insertion_breakpoint_score {

	#parse parameters
	my $chr = shift;
	my $start = shift;
	my $end = shift;

	my $genome_trans_data_hash_ref = shift;
	my %genome_trans_data_hash = %{$genome_trans_data_hash_ref};

	my $genome_trans_insertion_breakpoints_hash_ref = shift;
	my %genome_trans_insertion_breakpoints_hash = %{$genome_trans_insertion_breakpoints_hash_ref};

	my $bin_size = shift;

	#calculate array index where the data for the first and last bins of the region are located
	my $start_index = $start / ($bin_size);
	my $end_index = $end / ($bin_size);

	my $total_breakpoints = 0;	#total number of breakpoints in the region
	my $inserted_breakpoints = 0;	#total number of breakpoints with nearby insertions in the region

	my $insertion_breakpoint_score = 0;

	my @chr_data;		#stores the bin data for the chromosome
	my %trans_hash;		#stores translocation hash for each bin

	my @inserted_breakpoint_list;	#stores the breakpoints that have insertions nearby
	my @breakpoint_data;		#stores $total_breakpoints and $inserted_breakpoints for return

	#check if there is translocation data for the region
	if(defined($genome_trans_data_hash{$chr})){

		#get the translocation data
		@chr_data = @{$genome_trans_data_hash{$chr}};

		#for each bin in the region sum the number of breakpoints
		for (my $i = $start_index; $i < $end_index+1; $i++){
			if(!defined($chr_data[$i])){
				next;
				}
			%trans_hash = %{$chr_data[$i]};
			$total_breakpoints += $trans_hash{'BPcount'};
			}

		#get the list of breakpoints with insertions nearby on the chromosome
		@inserted_breakpoint_list = @{$genome_trans_insertion_breakpoints_hash{$chr}};

		#sort the above list
		{use sort 'stable';
		@inserted_breakpoint_list = sort {$a <=> $b} @inserted_breakpoint_list;
		}#use sort 'stable'

		#calculate how many of the breakpoints with insertions are in the region
		foreach my $breakpoint (@inserted_breakpoint_list){
			if($breakpoint > $end){
				last;
				}
			if($breakpoint > $start){
				$inserted_breakpoints++;
				}
			}

		#calculate the hallmark score
		if ($total_breakpoints > 0) {
			$insertion_breakpoint_score = $inserted_breakpoints/$total_breakpoints;
			}

		if($insertion_breakpoint_score > 1){
			die "ERROR: found a insertion_breakpoint_score greater than 1\n";
			}

		}

	push (@breakpoint_data, $inserted_breakpoints);
	push (@breakpoint_data, $total_breakpoints);

	return ($insertion_breakpoint_score, \@breakpoint_data);
	}


#=head2 Sub-Method: calculate_loh_score

### calculate_loh_score ###########################################################################
# Description:
#	Calculates the loss of heterozgozity hallmark score
#
# Input variables:
#	$chr:					chromosome where the region is located
#	$start:					starting location of the region
#	$end:					end location of the region
#	$chromosome_cnv_breakpoints_hash_ref:	stores the breakpoints of CNV mutations on each
#						chromosome
#	$chromosome_loh_breakpoints_hash_ref:	stores the breakpoints of LOH regions on each
#						chromosome
#

#=cut

sub calculate_loh_score {

	#parse parameters
	my $chr		= shift;
	my $start	= shift;
	my $end		= shift;

	my $chromosome_cnv_breakpoints_hash_ref	= shift;
	my %chromosome_cnv_breakpoints_hash	= %{$chromosome_cnv_breakpoints_hash_ref};

	my $chromosome_loh_breakpoints_hash_ref	= shift;
	my %chromosome_loh_breakpoints_hash	= %{$chromosome_loh_breakpoints_hash_ref};

	my @cnv_breakpoints;		#stores cnv breakpoints in the region
	my $cnv_breakpoints_size;	#stores the number of cnv breakpoints in the region

	my @loh_breakpoints;		#stores the LOH breakpoints in the region
	my $loh_breakpoints_size;	#stores the number of LOH breakpoints in the region

	#calculate maximum potential amount of heterozygosity
	my $original_heterozygous_size = $end - $start;

	#calculate maximum pontentail amount of heterozygosity that can remain
	my $remaining_heterozygous_size 	= $end - $start;

	my $loh_size = 0;	#stores the size of all LOH regions in the region

	my $loh_score = 0;	#final hallmark score

	#check if there is any cnv data for the chromosome
	if(
	   ( !defined($chromosome_cnv_breakpoints_hash{$chr}) ) ||
	   ( !defined($chromosome_loh_breakpoints_hash{$chr}) )
	  ){
		$loh_score = 0;
		$loh_size = -1;
		$remaining_heterozygous_size = -1;
		return($loh_score, $loh_size, $remaining_heterozygous_size);
		}

	#get the cnv breakpoints for the chromosome
	@cnv_breakpoints = @{$chromosome_cnv_breakpoints_hash{$chr}};

	#sort the list
	{use sort 'stable';
	@cnv_breakpoints = sort {$a <=> $b} @cnv_breakpoints;
	}#use sort 'stable'

	#get the number of breakpoints
	$cnv_breakpoints_size = @cnv_breakpoints;

	#find all the cnv events that occur in the region and subtract the size of these regions
	#from the $original_heterozygous_size value
	for (my $i = 0; $i< $cnv_breakpoints_size; $i+=2){
		my $cnv_start = $i;
		my $cnv_end = $i+1;

		my $end_overlap = 0;
		my $start_overlap = 0;

		if($cnv_breakpoints[$cnv_start] > $end){
			last;
			}

		#Check if the end point of the cnv region is within the suspect region
		if($cnv_breakpoints[$cnv_end] >= $start && $cnv_breakpoints[$cnv_end] <= $end){
			$end_overlap = 1;
			}

		#Check if the start point of the region cnv is within the suspect region
		if($cnv_breakpoints[$cnv_start] >= $start && $cnv_breakpoints[$cnv_start] <= $end){
			$start_overlap = 1;
			}

		#If an overlap was detected
		if($start_overlap==1 && $end_overlap==1) {
			$remaining_heterozygous_size -= ($cnv_breakpoints[$cnv_end] - $cnv_breakpoints[$cnv_start]);
			}
		elsif($start_overlap==1){
			$remaining_heterozygous_size -= ($end-$cnv_breakpoints[$cnv_start]);
			}
		elsif($end_overlap==1){
			$remaining_heterozygous_size -= ($cnv_breakpoints[$cnv_end]-$start);
			}
		elsif($cnv_breakpoints[$cnv_start] < $start && $cnv_breakpoints[$cnv_end] > $end){
			$remaining_heterozygous_size = 0;
			}
		}

	#check if there is no potential heterozygous regions in the suspect region or if there are no cnv events
	#in the suspect region, if either is the case return a score of 0
	if($remaining_heterozygous_size == 0 || $remaining_heterozygous_size == $original_heterozygous_size){
		$loh_score = 0;
		$loh_size = -1;
		$remaining_heterozygous_size = -1;
		return($loh_score, $loh_size, $remaining_heterozygous_size);
		}

	#get a list of the LOH breakpoints on the chromosome
	@loh_breakpoints = @{$chromosome_loh_breakpoints_hash{$chr}};

	#sort the list
	{use sort 'stable';
	@loh_breakpoints = sort {$a <=> $b} @loh_breakpoints;
	}#use sort 'stable'

	#get the number of breakpoints
	$loh_breakpoints_size = @loh_breakpoints;

	#determine which LOH regions are in the suspect region
	for (my $i = 0; $i< $loh_breakpoints_size; $i+=2){
		my $start_overlap_region_loh = 0;
		my $end_overlap_region_loh = 0;

		my $loh_start = $i;
		my $loh_end = $i+1;

		my $loh_start_breakpoint = $loh_breakpoints[$loh_start];
		my $loh_end_breakpoint = $loh_breakpoints[$loh_end];
		my $loh_region_size;

		if($loh_breakpoints[$loh_start] > $end){
			last;
			}

		#Check if the end point of the cnv region is within the loh region
		if($loh_breakpoints[$loh_end] >= $start && $loh_breakpoints[$loh_end] <= $end){
			$end_overlap_region_loh = 1;
			}

		#Check if the start point of region 1 is within region 2
		if($loh_breakpoints[$loh_start] >= $start && $loh_breakpoints[$loh_start] <= $end){
			$start_overlap_region_loh = 1;
			}

		#If an overlap was detected
		if($start_overlap_region_loh==1 && $end_overlap_region_loh!=1){
			$loh_end_breakpoint = $end;
			}
		elsif($end_overlap_region_loh==1 && $start_overlap_region_loh!=1){
			$loh_start_breakpoint = $start;
			}
		elsif($loh_breakpoints[$loh_start] < $start && $loh_breakpoints[$loh_end] > $end){
			$loh_start_breakpoint = $start;
			$loh_end_breakpoint = $end;
			$start_overlap_region_loh=1;
			$end_overlap_region_loh=1;
			}
		if($start_overlap_region_loh != 1 && $end_overlap_region_loh != 1){	#if the loh region is not in the suspect region go to the next loh region
			next;
			}

		#if the loh region is in the suspect region then reduce the size of the loh by subtracting the size of cnv regions that over lap with it
		#this will tell us how much original heterozygosity remains
		$loh_region_size = $loh_end_breakpoint - $loh_start_breakpoint;

		#check for overlaps between the LOH region and cnv regions
		for (my $k = 0; $k< $cnv_breakpoints_size; $k+=2){
			my $start_overlap_loh_cnv = 0;
			my $end_overlap_loh_cnv = 0;

			my $cnv_start = $k;
			my $cnv_end   = $k+1;

			if($cnv_breakpoints[$cnv_start] > $loh_end_breakpoint){
				last;
				}

			#Check if the end point of the cnv region is within the loh region
			if($cnv_breakpoints[$cnv_end] >= $loh_start_breakpoint && $cnv_breakpoints[$cnv_end] <= $loh_end_breakpoint){
				$end_overlap_loh_cnv = 1;
				}

			#Check if the start point of region 1 is within region 2
			if($cnv_breakpoints[$cnv_start] >= $loh_start_breakpoint && $cnv_breakpoints[$cnv_start] <= $loh_end_breakpoint){
				$start_overlap_loh_cnv = 1;
				}

			#If an overlap was detected
			if($start_overlap_loh_cnv==1 && $end_overlap_loh_cnv==1) {
				$loh_region_size -= ($cnv_breakpoints[$cnv_end] - $cnv_breakpoints[$cnv_start]);
				}
			elsif($start_overlap_loh_cnv==1){
				$loh_region_size -= ($loh_end_breakpoint-$cnv_breakpoints[$cnv_start]);
				}
			elsif($end_overlap_loh_cnv==1){
				$loh_region_size -= ($cnv_breakpoints[$cnv_end]-$loh_start_breakpoint);
				}
			elsif($cnv_breakpoints[$cnv_start] < $loh_start_breakpoint && $cnv_breakpoints[$cnv_end] > $loh_end_breakpoint){
				$loh_region_size = 0;
				}
			}
		$loh_size += $loh_region_size;
		}


	#calculate the LOH score
	$loh_score = 1 - ($loh_size/$remaining_heterozygous_size);

	if($loh_size> $remaining_heterozygous_size){
		die "ERROR: invalid LOH size value found\n";
		}

	return($loh_score, $loh_size, $remaining_heterozygous_size);

	}#sub calculate_loh_score


#=head2 Sub-Method: standard_deviation_and_mean

### standard_deviation_and_mean ###################################################################
# Description:
#	Calculates the standard deviation and mean for a given set of values
#
# Input variables:
#	$data_ref:	reference to either a hash or an array
#	$type:		0 indicates a hash, 1 indicates an array
#

#=cut

sub standard_deviation_and_mean{

	#parse parameters
	my $data_ref = shift;
	my $type = shift;

	my %hash;
	my @array;
	my $size;

	my $mean = 0;
	my $SD = 0;

	if($type==0){
		%hash = %{$data_ref};

		if((scalar(keys %hash))==0){
			die"Found sample size of 0 when calculating SD-hash\n";
			}

		#calculate mean
		for my $key (keys %hash){
			$mean += $hash{$key};
			}

		$mean = $mean/(scalar(keys %hash));

		#calculate sum of squared differences
		for my $key (keys %hash){
			$SD += ($hash{$key}-$mean)**2;
			}

		#calculate final standard deviation value
		$SD = $SD/(scalar(keys %hash));
		$SD = $SD**0.5;
	}
	elsif($type==1){
		@array = @{$data_ref};
		$size = @array;

		if($size==0){
			die"Found sample size of 0 when calculating SD-array\n";
			}

		#calculate mean
		foreach my $val (@array){
			$mean += $val;
			}

		$mean = $mean/$size;

		#calculate sum of squared differences
		foreach my $val (@array){
			$SD += ($val-$mean)**2;
			}

		#calculate final standard deviation value
		$SD = $SD/$size;
		$SD = $SD**0.5;

		}
	else{
		die"ERROR: invalid SD/mean type found\n";
		}
	return ($SD, $mean);

	}#sub standard_deviation_and_mean


### next_arg ######################################################################################
# Parse the next arguement from the command line
#
sub next_arg {
	my $code = shift;
	$pos++;
	if($pos == $ARGC){
		usage($code);
		}
	}#sub next_arg

### man_text ######################################################################################
# Print the manual help text
#
sub man_text {
	print "Main Usage:\n";
	print "\tperl -w shatterproof.pl --cnv  --trans  [--insrt ] [--loh ] [--tp53] --config  --output \n";
	print "\n";
	print "\tArguments:\n";
	print "\t\t--cnv\t\tDefine the path to the directory containing the CNV input files\n";
	print "\t\t--trans\t\tDefine the path to the directory containing the Translocation input files\n";
	print "\t\t--insrt\t\tDefine the path to the directory containing the insertion VCF input files\n";
	print "\t\t--loh\t\tDefine the path to the directory containing the LOH input files\n";
	print "\t\t--tp53\t\tIndicate that TP53 should be considered mutated, regardless of data\n";
	print "\t\t--config\tDefine the path to the ShatterProof config file\n";
	print "\t\t--output\tDefine the path to the directory where output should be placed\n";
	print "\t\tdir\t\tPath to a directory\n";
	print "\t\tpath\t\tPath to a file\n";
	print "\n";
	print "Help Usage:\n";
	print "\tperl -w shatterproof.pl --help\t\tThis help message.\n";
	print "\n";
	exit 0;
	}#sub man_text

### usage #########################################################################################
# Prints an error message when invalid command line arguements are found
#
sub usage {
	my $usage_msg = shift;

	print "u $usage_msg \n";

	given($usage_msg){
		when (/^0/) 	{ print "ERROR: missing arguments\n"; 			}
		when (/^1/) 	{ print "ERROR: 2nd argument missing\n"; 		}
		when (/^2/) 	{ print "ERROR: CNV directory missing\n"; 		}
		when (/^3/) 	{ print "ERROR: --trans option missing\n"		}
		when (/^4/) 	{ print "ERROR: --cnv option missing\n"			}
		when (/^5/) 	{ print "ERROR: Translocation directory missing\n"	}
		when (/^6/) 	{ print "ERROR: --config option missing\n"		}
		when (/^7/) 	{ print "ERROR: --trans option missing\n"		}
		when (/^8/) 	{ print "ERROR: insertion directory missing\n"		}
		when (/^9/) 	{ print "ERROR: --config option missing\n"		}
		when (/^10/) 	{ print "ERROR: LOH directory missing \n"		}
		when (/^11/) 	{ print "ERROR: --config option missing\n"		}
		when (/^12/) 	{ print "ERROR: --config option missing\n"		}
		when (/^13/) 	{ print "ERROR: Path to config file missing\n"		}
		when (/^14/) 	{ print "ERROR: --output option missing\n"		}
		when (/^15/) 	{ print "ERROR: --config option missing\n"		}
		when (/^16/) 	{ print "ERROR: Output directory missing\n"		}
		when (/^17/) 	{ print "ERROR: --output option missing\n"		}
		when (/^18/) 	{ print "ERROR: too many arguments\n"			}
		}
	print "Try perl -w shatteproof.pl --help\n";
	exit 0;
	}#sub usage

### initialize_genome_hash ########################################################################
# Description:
#	Initializes a hash to store an array for each chromosome
#
sub initialize_genome_hash {

	my %genome_region_data = (  	#{chr}[region_num]->%region_data
			 X => [],
			 Y => [],
			 1 => [],
			 2 => [],
			 3 => [],
			 4 => [],
			 5 => [],
			 6 => [],
			 7 => [],
			 8 => [],
			 9 => [],
			 10 => [],
			 11 => [],
			 12 => [],
			 13 => [],
			 14 => [],
			 15 => [],
			 16 => [],
			 17 => [],
			 18 => [],
			 19 => [],
			 20 => [],
			 21 => [],
			 22 => []
			 );

	return (\%genome_region_data);
	}#sub initialize_genome_hash

### load_config_file #########################################################################################
# Description:
#	Opens the config file and reads the parameter values from it
#
# Input variables:
#	$path:	path to config file
#
sub load_config_file {

        my $path = shift;
        print "\nLoading configuration file";

        #Load the configuration file config.pl
        my $CONFIG;
        open($CONFIG, "<","$path") or die "COULD NOT OPEN CONFIG FILE at path: $path \n";
        eval (<$CONFIG>) while (!eof($CONFIG));
	close($CONFIG);

        print " - Done\n";
	1;
        }#sub load_config_file

=head1 NAME

ShatterProof - a script for analyzing next-generation sequencing data

=head1 SYNOPSIS

use Shatterproof

See "shatterproof.pl" in the scripts directory for a simple perl script which calls the ShatterProof module

Call ShatterProof via:

	ShatterProof::run(\@ARGV);

=head1 DESCRIPTION

ShatterProof is a tool that can be used to analyze next generation sequencing data for signs of chromothripsis. ShatterProof is implemented as a Perl module that processes input files and produces output files in both tab-delimited and YAML format. Perl version 5.0 or greater is required to run ShatterProof. Link to publication will be posted soon.

=head1 README

=head2 Installing ShatterProof

To install this module type the following:

   perl Makefile.PL
   make
   make test
   make install

Make sure that you have admin permission rights when running the previous commands.

=head2 Input File Types

ShatterProof bases its analysis of genomic data on calls of translocations, copy number variations (CNV), loss of heterozygosity (LOH) and insertions.
ShatterProof can takes as input 4 different types of input files.  See the scripts/conversion_scripts directory for some Perl scripts which will convert some common tools' output to the required input formats.

=head3 Translocation Input Files (.spt)

Tab delimited columns
First line is header line:
#chr1   start   end     chr2    start   end     quality

Example data entry line:

	1       1000    2000    4       4000    5000    78

If no value is available for quality, use a "." eg.:

	1       1000    2000    4       4000    5000    .


=head3 Copy-Number Input Files (.spc)

Tab delimited columns
First line is header line:
#chr    start   end     number  quality


Example data entry line:
	12      2000    3000    2       63

If no value is available for quality, use a "." eg.:

	12      2000    3000    2       .

=head3 Loss of Heterozygozity Input Files (.spl)

Tab delimited columns
First line is header line:
#chr    start   end	quality


Example data entry line:

	12      2000    3000	63

If no value is available for quality, use a "." eg.:

	12      2000    3000	.

=head3 Insertion Input Files (.vcf)

Additionally, ShatterProof accepts insertion calls in VCF files as input. See http://www.1000genomes.org/node/101 for details on the VCF file format.
ShatterProof analyzes the CHROM and POS fields of these files.


=head2 Configuring ShatterProof

See the config.pl file in the scripts directory for a sample ShatterProof configuration file.

$bin_size: number (integer) of base pairs to include in each bin of the sliding window analysis

$localization_window_size: number (integer) of bins to include in each window of the sliding window analysis

$expected_mutation_density: a reference value (double) used in determining if the concentration of translocation events on a particular chromosome is higher than expected.

$collapse_regions:

	flag variable

	value 1: merge overlapping CNV regions that have the same copy number

	value 0: do not merge overlapping CNV regions that have the same copy number. 	If such regions are found an error is thrown

$outlier_deviations: the number of standard deviations away from the mean a value has to be in order to be considered non-significant. Used to identify highly mutated regions.

$translocation_cut_off_count: the maximum number of translocation chromosomes to tolerate before the translocation score for a region is set to 0.

$genome_localization_weight: weight given to the localization of mutations to one chromosome hallmark

$chromosome_localization_weight: weight given to the localization of mutations to one area of a particular chromosome hallmark

$cnv_weight: weight given to the concentrated CNV hallmark

$translocation_weight: weight give to the concentrated translocations hallmark

$insertion_breakpoint_weight: weight given the the short breakpoint insertions hallmark

$loh_weight: weight given to the loss/retention of heterozygosity hallmark

$tp53_mutated_weight: weight given to the TP53 mutation hallmark


=head2 Running ShatterProof

From the scripts directory run execute the shatterproof.pl file using Perl.

Main Usage:

perl -w shatterproof.pl --cnv  --trans  [--insrt ] [--loh ] [--tp53] --config  --output

 Arguments:

--cnv		Define the path to the directory containing the CNV input files

--trans		Define the path to the directory containing the Translocation input files

--insrt		Define the path to the directory containing the insertion VCF input files

--loh		Define the path to the directory containing the LOH input files

--tp53		Indicate that TP53 should be considered mutated, regardless of data

--config	Define the path to the ShatterProof config file

--output	Define the path to the directory where output should be placed

dir		Path to a directory

path		Path to a file

=head1 PREREQUISITES

strict;
warnings;
Carp;
Switch;
File::Basename;
List::Util qw[min max];
Statistics::Distributions;
POSIX

=pod OSNAMES

any

=pod SCRIPT CATEGORIES

CPAN

=cut

1;
__END__