peaks2genes

#!/usr/bin/perl -w

# Copyright © 2015, Owen Marshall

# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or (at
# your option) any later version. 
# 
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# General Public License for more details. 
# 
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 
# USA

use strict;
use File::Basename;
use 5.010;
$|++;

my $version = "1.0.1";

print STDERR "\npeaks2genes v$version\nCopyright © 2015, Owen Marshall\n\n";

my %vars = (
	'pad' => 0,
	'genes_file' => '',
	'gene_pad' => 1000,
	'mouse' => 0,
	'promoters_only'=>0,
	'exclude_lincs' => 0,
);

my %vars_details = (
	'pad' => 'Add padding around peaks (bp)',
	'genes_file' => 'GFF file of annotated genes to use',
	'gene_pad' => 'Consider a peak associated with a gene if it lies within this distance of the gene start or end',
	'mouse' => 'Deals with Ensembl formatted GFFs for the mouse genome',
	'promoters_only' => 'Only check promoters for overlap',
	'exclude_lincs' => 'Ignore lncRNAs in output',
);

my $HOME = (getpwuid($<))[7];

my @in_files;
process_cli();

unless (@in_files) {
	help();
}

my @genes;
my @genes_f;
my @genes_r;
my @peaks;

unless ($vars{'genes_file'}) {
	die("Error: no genes file specified.  Please specify a valid genes file with the --genes_file option ...\n\n")
}

unless (-e $vars{'genes_file'}) {
	die("Error: Genes file does not exist.  Please specify a valid genes file with the --genes_file option ...\n\n")
}


### pipeline

foreach my $fn (@in_files) {
	my ($name,$dir,$ext) = fileparse($fn, qr/\.[^.]*/);
	load_file($fn, \@peaks);
	read_genes($vars{'genes_file'});
	report_genes($name,$dir);
	say '';
}

say "All done.\n";

exit 0;


sub load_file {
	# Load peaks file (gff)
	my ($fn, $aref) = @_;
	say "Loading $fn ...";
	open (IN, "<$fn") || die "Unable to read input file $fn: $!\n";
	
	while (<IN>) {
		chomp;
		next if m/^#/;
		
		my @tmp = split("\t");
		
		if ($#tmp > 1) {
			push @$aref, [@tmp[0,3,4,5]];
		} else {
			# allow a list of chr:start-end ...
			my ($chr, $start, $end) = $tmp[0] =~ m/(.*?)\:(\d+)-(\d+)/;
			push @$aref, [$chr,$start,$end,1];
		}
	}
	close (IN);
	
	map { @{$_}[0] =~ s/chr//g } @$aref;
}

sub read_genes {
	my $fn = shift;
	
	if ($fn =~ m/\.gz$/) {
		# gzipped file
		open (GENES, '-|', "gunzip -c $fn") || die "Error: cannot open genes file $fn: $!\n\n";
	} else {
		open(GENES, "<$fn")|| die "Error: cannot open genes file $fn: $!\n\n";
	}
	
	while (<GENES>) {
		my ($chr, $a, $b, $start, $end, $c, $pol, $d, $info) = split(/\t/);
		next unless $info;
		
		my $name;
		if ($vars{'mouse'}) {
			($name) = $info =~ m/gene_name \"(.*?)\"/;
		} else {
			($name) = $info =~ m/Name=(.*?);/;
		}
		
		if ($vars{'exclude_lincs'}) {
			next if $name =~ m/^CR\d+/;
		}
		
		$chr =~ s/chr//g;
		
		push @genes, [$chr, $start, $end, $pol, $name];
		if ($pol eq '+') {
			push @genes_f, [$chr, $start, $end, $pol, $name];
		} else {
			push @genes_r, [$chr, $start, $end, $pol, $name];
		}
		
	}
	close GENES;
}

sub report_genes {
	my ($name,$dir) = @_;
	$dir ||= '.';
	
	say "  Searching for peaks overlapping genes ...";
	my ($p1inc, @genes_overlap) = find_genes(\@peaks);	
	
	say "  Writing genes list ...";
	write_gene_lists(\@genes_overlap, "$dir/$name.genes.in.peak.csv");
}

sub write_gene_lists {
	my ($ar, $name) = @_;
	open (GENEOUT, '>', $name) || die "Cannot open gene lists for writing: $!\n";
	say GENEOUT join("\t",@$_) foreach @$ar;
	close GENEOUT;
}

sub find_genes {
	my $ar = shift;
	my @matched;
	my $overlap=0;
	foreach (@$ar) {
		my ($c,$s,$e,$sc) = @{$_};
		if ($vars{'promoters_only'}) {
			my @mgf = grep {
				@{$_}[0] eq $c
				&& $e >= @{$_}[1] - $vars{'gene_pad'}
				&& $s <= @{$_}[1] + $vars{'gene_pad'}
			} @genes_f;
			
			if (@mgf) {
				$overlap++;
				push @matched, [@{$_}[4], $sc] foreach @mgf;
			}
			
			my @mgr = grep {
				@{$_}[0] eq $c
				&& $e >= @{$_}[2] - $vars{'gene_pad'}
				&& $s <= @{$_}[2] + $vars{'gene_pad'}
			} @genes_r;
			
			if (@mgr) {
				$overlap++;
				push @matched, [@{$_}[4], $sc] foreach @mgr;
			}
			
		} else {
			my @mg = grep {
				@{$_}[0] eq $c
				&& $e >= @{$_}[1] - $vars{'gene_pad'}
				&& $s <= @{$_}[2] + $vars{'gene_pad'}
			} @genes;
			
			if (@mg) {
				$overlap++;
				push @matched, [@{$_}[4], $sc] foreach @mg;
			}
		}
	}
	return ($overlap, sort {@$b[1] <=> @$a[1]} uniq(@matched));
}

sub printout {
	my ($hr, $fh, $ar) = @_;
	foreach (keys %$hr) {
		my ($c,$s,$e,$sc) = split(/-/);
		push @$ar, [$c,$s,$e,$sc] if $ar;
		print $fh join("\t",$c,'.','.',$s,$e,$sc,'.','.','.'),"\n";
	}
}

sub excl {
	my ($phr, $ar, $fh, $aer) = @_;
	foreach (@$ar) {
		my ($c,$a,$b,$s,$e,$sc) = @{$_};
		next if $phr -> {"$c-$s-$e-$sc"};
		push @$aer, [$c,$s,$e,$sc] if $aer;
		print $fh join("\t",$c,'.','.',$s,$e,$sc,'.','.','.'),"\n";
	}
}


sub peak_score {
	my ($href, $aref) = @_;
	my ($peaksum, @peakscore);
	foreach (keys %$href) {
		my ($chr, $start, $end, $score) = split("-");
		$peaksum += $end-$start;
		push @peakscore, $score;
	}
	
	my $avg_peak_score = mean(@peakscore);
	return ($peaksum, $avg_peak_score);
}

sub max {
    my ($max, @vars) = @_;
	my $index=0;
	$max||=0;
    for my $i (0..$#vars) {
        ($max, $index) = ($vars[$i], $i+1) if $vars[$i] > $max;
    }
    return ($index, $max);
}

sub min {
    my ($min, @vars) = @_;
	my $index=0;
	$min||=0;
    for my $i (0..$#vars) {
        ($min, $index) = ($vars[$i],$i+1) if $vars[$i] < $min;
    }
    return ($index, $min);
}

sub uniq {
	return unless @_;
	my %seen;
	return grep { !$seen{$_}++ } @_;
}


sub mean {
	my @a = @_;
	my $n = @a;
	return $a[0] if $n<2;
	
	my $total;
	map {$total += $_} @a;
	my $mean = $total/$n;
	
	return $mean;
}

sub median {
	my @a = sort {$a <=> $b} @_;
	my $index = int($#a/2);
	return $a[$index];
}

sub stdev {
	my @a = @_;
	my $n = @a;
	return (0, 0) if $n<2;

	my $mean = mean(@a);
	
	my $sum = 0;
	map {$sum += abs($_-$mean)} @a;
	
	my $sd= $sum/($n-1);
	my $se = $sd/sqrt($n);
	
	return $sd;
}



sub process_cli {
	foreach (@ARGV) {
		if (/--(.*)=(.*)/) {
			unless (defined($vars{$1})) {
				print STDERR "Did not understand $_ ...\n";
				help();
			}
			my ($v, $opt) = ($1,$2);
			$opt =~ s/~/$HOME/;
			$vars{$v} = $opt;
			next;
		} elsif (/--h[elp]*/) {
			help();
		} elsif (/--(.*)/) {
			if (defined($vars{$1})) {
				$vars{$1} = 1;
			} else {
				print STDERR "Did not understand $_ ...\n";
				help();
			}
			next;
		}
		push @in_files, $_;
	}
}


sub help {
	print STDOUT <<EOT;
peaks2genes.pl -- calls genes in close proximity to peaks
Usage: peaks2genes [options] [peak files in GFF format]
	
Options:
EOT
	
	my $opt_len = 0;
	foreach (keys %vars) {
		my $l = length($_);
		$opt_len = $l if $l > $opt_len;
	}
	
	$opt_len+=2;
	
	my $cols= `tput cols` || 80;
	
	my ($v, $val, $def, $def_format);
	my $help_format = "format STDOUT =\n"
		.' '.'^'.'<'x$opt_len . ' '. '^' . '<'x($cols-$opt_len-4) . "\n"
		.'$v, $def_format'."\n"
		.' '.'^'.'<'x$opt_len . '   '. '^' . '<'x($cols-$opt_len-6) . "~~\n"
		.'$v, $def_format'."\n"
		.".\n";
		
	eval $help_format;
	die $@ if $@;
	
	foreach my $k (sort (keys %vars)) {
		($v, $val, $def) = ($k, $vars{$k}, $vars_details{$k});
		$def||="";
		$def_format = $val ? "$def\n\r[Current value: $val]" : $def;
		$v = "--$v";
#		format =
# ^<<<<<<<<<<<<<<<<<<<< ^<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
#$v, $def_format
# ^<<<<<<<<<<<<<<<<<<<<   ^<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ~~
#$v, $def_format
#.

		write();
		
	}
	print STDOUT "\n";
	exit 1;
}