#!/usr/bin/perl -w
package GenomicsUtils;
require Exporter;

use strict;
use DBI;
use Gene;
our (@ISA, @EXPORT);
 
@ISA = qw(Exporter);	
@EXPORT = qw( connect disconnect dbh query queryList queryScalar queryTable queryHashList queryAttrs
	      getDebug setDebug warn
	      addToHashList );
{

sub translation($);
sub reverseComplement($);
sub readFastaData($);
# Reads data from FASTA file

sub readFASTA($);
# Read FASTA files from data directory

sub dbh();
sub getDebug;
sub setDebug;
sub query($);

sub queryList($);
sub queryScalar($);
sub queryAttrs($$$$$$);
sub queryTable($$);
sub queryHashList($);

sub connect;
sub disconnect();
sub getSequence($;$;$;$);
# Returns gene's sequence.
# Parameters:
#    first - sequence
#    second - strand
#    		

sub formatFASTA($;$);
###//sub get_seq($;$;$);

sub setDebug;
sub getDebug;
sub warn;


my $DATA_DIR = "/data/genomics/data/genomes/";
my $dbh_;
my ( $noDebug, $warnings, $verbose, $full ) = ( 0, 1, 2, 10 );
my $debug_ = $warnings;

#--------------------------------------------------------
sub data_dir()
{
    return $DATA_DIR;
}

#--------------------------------------------------------

sub reverseComplement($) 
{
    my $seq = shift;
    chomp $seq;

    die "Invalid sequence \"$seq\" in reverseComplement" if ( ($seq =~ 
tr/RYKMSWBDHVNATCGrykmswbdhvnatcg/YRMKWSVHDBNTAGCyrmkwsvhdbntagc/) != 
length($seq) );
    $seq = reverse $seq;
    return $seq;
}

#--------------------------------------------------------
sub readFastaData($)
{
    my $file = shift;
    open( F, "$file" ) or print STDERR "Cannot open file $file\n";
    my @lines = <F>;
    close F;
    my $dna = "";
    foreach my $l (@lines) {
	next if ( $l =~ /^>/ );
	chomp $l;  
	$dna .= $l;
    }

    return $dna;
}

sub readFASTA($)
{
    my $file = shift;
    $file = $DATA_DIR.$file;
    return readFastaData( $file );
}

#--------------------------------------------------------

sub dbh()
{
    return $dbh_;
}

sub getDebug{
    return $debug_;
}

sub setDebug{
    $debug_ = 1;
    $debug_ = $_[0] if defined $_[0];
}

sub connect
{

    my $host = (defined($_[0])&&($_[0] ne "localhost"))?$_[0].".lbl.gov":"localhost";
    my $db = "genomics_test";
    my $dsn = "DBI:mysql:$db:".$host;
    my $user_name="test";
    my $password="test";

    $dbh_ = DBI->connect( $dsn, $user_name, $password, {RaiseError => 1} );

    return $dbh_;
}

sub disconnect(){
    $dbh_->disconnect();
}

sub query($){
    my $temp_dbh = $dbh_;
    my $temp_statement = shift;
    print STDERR "GenomicsUtils(line 130): $temp_statement\n" if $debug_>=10;
    my $sth = $temp_dbh->prepare (qq{$temp_statement});
    $sth->execute();

    my $refArray = $sth->fetchall_arrayref();
    $sth->finish(); 
    &warn( "GenomicsUtils: query returned NULL result; $temp_statement" )
	unless $refArray->[0] || $debug_ <= 1;
    return $refArray;
}

sub queryList($){
    return map {$_->[0]} @{&query(@_)};
}

sub queryScalar($){
    my $qRef = &query(@_);
    return $qRef->[0][0];
}

# queryAttrs(table name, unique key name, unique key value, reference to list of attribute names,
#            hash to fill, suffix to put on names)
# returns: nothing
sub queryAttrs($$$$$$){
    my ($table, $keyName, $key, $refAttrs, $refHash, $suffix) = @_;
    my $result = &query("SELECT " . join(",", @$refAttrs) . " from $table where $keyName=$key");
    die "Missing or non-unique key $keyName=$key in $table" if scalar @$result != 1;
    $result = $result->[0];
    my @columns = @$refAttrs; # copy
    while(scalar @$result > 0) {
	$refHash->{(shift @columns) . $suffix} = shift @$result;
    }
}

# queryTable(attrNames, fromClause)
# returns reference to a list of hashes for each row, each of the form attr->value
sub queryTable($$) {
    my ($refAttrs, $from) = @_;
    my $result = &query("SELECT " . join(",", @$refAttrs) . " " . $from);
    my $list = [];
    foreach my $row (@$result) {
	my $hash = {};
	my @attrs = @$refAttrs;
	while (scalar @attrs > 0) {
	    $hash->{shift @attrs} = shift @$row;
	}
	push @$list, $hash;
    }
    return($list);
}

# queryHashList(query)
# returns a reference to a hash of key -> list of values
# query must return two-column results
sub queryHashList($) {
    my $query = shift;
    my $result = &query($query);
    my $hashRef = {};
    foreach my $row (@$result) {
	die "row size != 2: query $query" if scalar @$row != 2;
	addToHashList($hashRef, $row->[0], $row->[1]);
    }
    return($hashRef);
}

sub warn{
    return unless $debug_;
    my $warning = shift;
    my $level;
    $level = $warnings unless ($level = shift);
    print STDERR "$warning\n" if $level >= $debug_;
}

sub getSequence($;$;$;$) {
    my ($seq, $strand, $start, $end) = @_;
    my $sequence = substr(${$seq}, $start, $end - $start  + 1);

    if( $strand eq '-' ) {
	return reverseComplement( $sequence );
    }

    return $sequence;
}

#-------------------------------------------------------------------
sub formatFASTA($;$)
{
    my ($title, $sequence) = @_;
    my $s = qq/>$title\n/;
    my @whole = $sequence =~ /.{1,60}/g;
    my $line = join "\n", @whole;
    return $s."$line\n";
}

sub addToHashList($$$) {
    my ($hashRef, $key, $value) = @_;
    if (!defined $hashRef->{$key}) {
	$hashRef->{$key} = [];
    }
    push @{ $hashRef->{$key} }, $value;
    return 1;
}

## stop codon is included in subroutine translation as * symbol
sub translation($){
    my $transcript = shift;
    if (length($transcript)%3){&warn("Transcript is not multiple of 3");}
    my %codonTable = (Gene::codonUsageTable());
    my $protein = "";
    for(my $i = 0; $i < length($transcript); $i+=3){
	my $codon = substr($transcript, $i, 3);
	if ($codonTable{$codon}){
	    $protein .= $codonTable{$codon};
	}
	else{
	    $protein .= 'X';
	}

    } 
    
    return $protein;
    
}
#=======
}

return 1;
