#!/usr/bin/perl
package GCG;

##########################################################################
##### Author :	Shuai Weng
##### Date   :  Nov. 2001
##### Description : This object contains methods for 
#####               retrieving gcg sequence for a given 
#####               sequence name or a given chrnum/start_coord/
#####               stop_coord.
#####
##### See documentation for usage details:
#####
##### http:///usr/local/dicty/www_dictybase/db/lib/staff/dictyBase/programmer/GCG.html
#####
#######################################################################
use strict;
use DBI;
use lib "/usr/local/dicty/www_dictybase/db/lib/common";
use TextUtil qw (DeleteUnwantedChar);
use lib "/usr/local/dicty/www_dictybase/db/lib/dictyBase/Objects";
use ConfigPathdictyBase;


#######################################################################
######################## global variables #############################
my $configPath;
my $GCGDataDir;
#######################################################################
sub new {      ############ constructor ###############################
#######################################################################
       
    my ($self, %args) = @_;

    $self = {};
    bless $self;

    $self->{'_map'} = $args{'map'} || 'amap';
    $self->{'_seqname'} = $args{'seqname'};
    $self->{'_chrnum'} = $args{'chrnum'};
    $self->{'_beg'} = $args{'beg'} || '1';
    $self->{'_end'} = $args{'end'} || '10000';
    $self->{'_rev'} = $args{'rev'};
    $self->{'_men'} = $args{'men'} || '-MEN=o'; ### o=open reading frames 
    $self->{'_sequence'} = $args{'sequence'};
   
    if (!$self->{'_seqname'} && !$self->{'_chrnum'} && 
	!$self->{'_sequence'}) {
	print "You must pass a sequence name, a chromosome number or DNA sequence to GCG object.";
	return;
    }
    return $self;
}

#######################################################################
sub getSequence {
#######################################################################
    my ($self) = @_;

    $configPath = ConfigPathdictyBase->new;
    my $tmpDir = $configPath->tmpDir;
    
    my $cmdtmp = $tmpDir."gcgcmd.tmp.$$";   # Command temp file    
    my $seqtmp = $tmpDir."gcgseq.tmp.$$";   # Sequence temp file 
    my $seqtmp2 = $tmpDir."gcgseq.tmp.2.$$";

    my $cmd = $self->getCmd($seqtmp, $seqtmp2);

    ### build the command file
    open (CMDTMP, ">$cmdtmp") ||
	die "GCG: Can't create cmd file '$cmdtmp' for writing:$!.\n";

    print CMDTMP "#!/bin/csh -f\n";
    print CMDTMP "cd /tmp\n";
    print CMDTMP "source /share/Dictyostelium/www-data/cgi-bin/gcg/gcgbatch.startup\n";
    print CMDTMP "$cmd\n";
    close(CMDTMP);

    system "rsh Dictyostelium /bin/csh $cmdtmp" || 
	print "GCG: problem running command: $!\n";

    unlink $cmdtmp;

    if ($self->{'_map'} =~ /^[rf]/i && $self->{'_seq'}) {

	return $seqtmp;

    }

    if (-e "$seqtmp2") {

	return $seqtmp2;

    }

    return $seqtmp;
    
}

########################################################################
sub getCmd {
########################################################################
    my ($self, $seqtmp, $seqtmp2) = @_;

    my $gcgbin = "/tools/gcg/10.3/gcgbin/execute/";


    my $cmd;

    if ($self->{'_sequence'}) { 
 
	$self->{'_seq'} = $self->{'_sequence'};

        ### translate DNA sequence to protein sequence

	$self->{'_sequence'} =~ s/[^a-zA-Z]//g;

	my $tmpfile = $seqtmp.".tmp";
	my $seqtmp3 = $seqtmp2;
	$seqtmp3 =~ s/\.2\./\.3\./;
	my $translateFile = $seqtmp;
	$translateFile =~ s/gcgseq\./gcgseq\.translated\./; 
	

	open(TMP, ">$tmpfile") || 
	    die "GCG: Can't open '$tmpfile' for writing:$!\n";
	print TMP "tmp sequence file ..\n";

	if (length($self->{'_sequence'}) > 200) {
	    my $subseq = substr($self->{'_sequence'}, 0, 200);
	    $self->{'_sequence'} = substr($self->{'_sequence'}, 200);
	    while($subseq) {
		print TMP "$subseq\n";
		$subseq = substr($self->{'_sequence'}, 0, 200);
		$self->{'_sequence'} = substr($self->{'_sequence'}, 200);
	    }
	    print TMP "\n";
	}
	else {
	    print TMP $self->{'_sequence'}, "\n";
	}
	close (TMP);	
	
	$cmd = $gcgbin."reformat -D -in=$tmpfile -out=$seqtmp2\n";
	if ($self->{'_map'} =~ /^rmap/i) { ## six frame
    
	    $cmd .= $gcgbin."map -D -in=$seqtmp2 -out=$seqtmp ".$self->{'_men'};
 
	}
	elsif ($self->{'_map'} =~ /^fmap/i) {## restriction fragments

	    $cmd .= $gcgbin."mapsort -D -in=$seqtmp2 -out=$seqtmp";

	}
	else { ## translate DNA to protein

	    $cmd .= $gcgbin."translate -D -in=$seqtmp2 -out=$translateFile";

	}

	return $cmd;

    }
 
    my $infile;
    $self->{'_chrnum'} =~ s/^17$/Mt/;
    if ($self->{'_chrnum'}) {
	$infile = $configPath->DictyosteliumFtpPubDir."data_download/sequence/genomic_sequence/chromosomes/gcg/chr".$self->{'_chrnum'}.".gcg";
    }
    else {
	$infile = "gb_sc:";
    }

    if ($self->{'_map'} =~ /^p/i) { ### protein sequence
	if ($self->{'_chrnum'}) {  ### chromosome region/seqment
	    if ($self->{'_seqname'} =~ /^Y[A-Z](L|R)[0-9]/i || 
		$self->{'_seqname'} =~ /^Y[A-Z](L|R)X[0-9]/i || 
		$self->{'_seqname'} =~ /^Q[0-9][0-9][0-9][0-9]/i) { ### orf name
		$infile = "orfp:";
		$cmd = $gcgbin."fetch -D -in=$infile".$self->{'_seqname'}." -out=$seqtmp\n"; 
	    }
	    else {
		$cmd = $gcgbin."translate -D -in=$infile -out=$seqtmp -beg=",$self->{'_beg'}." -end=".$self->{'_end'}." ".$self->{'_rev'}."\n";
	    }
	}
	else {
	    $cmd = $gcgbin."translate -D -in=$infile".$self->{'_seqname'}." -out=$seqtmp -beg=",$self->{'_beg'}." -end=".$self->{'_end'}." ".$self->{'_rev'}."\n";
	}
    }
    elsif ($self->{'_map'} =~ /^n/i) { 
	$infile = "orfn:";
	$cmd = $gcgbin."assemble -D -in=$infile".$self->{'_seqname'}." -out=$seqtmp ".$self->{'_rev'}."\n"; 
    }
    else {
	if ($self->{'_chrnum'}) {  ### chromosome region/seqment
	    $cmd = $gcgbin."assemble -D -in=$infile -out=$seqtmp -beg=".$self->{'_beg'}." -end=".$self->{'_end'}." ".$self->{'_rev'}."\n";
	}
	else { ### GenBank sequence
	    $cmd = $gcgbin."assemble -D -in=$infile".$self->{'_seqname'}." -out=$seqtmp -beg=".$self->{'_beg'}." -end=".$self->{'_end'}." ".$self->{'_rev'}."\n";
	}
    }
    
    #################################################
    #################################################

    ######## translate dna sequence into protein
    if ($self->{'_map'} eq 'pmap') {
	return $cmd;
    }

    ######## translate dna sequence into protein (p2map=No header)
    if ($self->{'_map'} eq 'p2map') {
	$cmd .= $gcgbin."tostaden -D -in=$seqtmp -out=$seqtmp2";
	return $cmd;
    }

    ######## translate dna sequence into protein (p3map=fasta format)
    if ($self->{'_map'} eq 'p3map') {
	$cmd .= $gcgbin."tofasta -D -in=$seqtmp -out=$seqtmp2 >& /dev/null\n";
	return $cmd;
    }

    ######## dna sequence without introns for splice orfs
    if ($self->{'_map'} eq 'nmap') {
	return $cmd;
    }
    
    ######## dna sequence without introns for splice orfs (n2map=No header)
    if ($self->{'_map'} eq 'n2map') {
	$cmd .= $gcgbin."tostaden -D -in=$seqtmp -out=$seqtmp2";
	return $cmd;
    }
    ######## dna sequence without introns for splice orfs (n3map=fasta format)

    if ($self->{'_map'} eq 'n3map') {
	$cmd .= $gcgbin."tofasta -D -in=$seqtmp -out=$seqtmp2 >& /dev/null\n";
	return $cmd;
    }
    
    ######## restriction map
    if ($self->{'_map'} eq 'rmap') {  
	if ($self->{'_chrnum'} eq "Mt") {
	    $cmd .= $gcgbin."map -D -in=$seqtmp -out=$seqtmp2 ".$self->{'_men'}." -TRANS=/tools/gcg/10.3/gcgcore/data/moredata/transl_table_03.txt";
	}
	else {
	    $cmd .= $gcgbin."map -D -in=$seqtmp -out=$seqtmp2 ".$self->{'_men'};
	}
	return $cmd;
    }
    
    ######## restriction fragments
    if ($self->{'_map'} eq "fmap") {
	$cmd .= $gcgbin."mapsort -D -in=$seqtmp -out=$seqtmp2";
	return $cmd;
    }

    ######## assemble single strand sequence (amap=GCG) 
    if ($self->{'_map'} eq "amap") {
	return $cmd;
    }
    
    ######## assemble single strand sequence (a2map=No header)
    if ($self->{'_map'} eq "a2map") {
	$cmd .= $gcgbin."tostaden -D -in=$seqtmp -out=$seqtmp2";
	return $cmd;
    }

    ######## assemble sequence to fasta format (a3map=fasta)
    if ($self->{'_map'} eq "a3map") {
	$cmd .= $gcgbin."tofasta -D -in=$seqtmp -out=$seqtmp2 >& /dev/null\n";
	return $cmd;
    }
}

########################################################################
1;
########################################################################

=pod

=head1 Name

GCG.pm    

=head1 Description

This perl object (GCG.pm) contains method for retrieving gcg DNA/protein sequences or for translating DNA sequence to protein sequence. 

=head1 Instantiating a New GCG Object

To instantiate a new GCG object, you may use one of following syntaxes:

my $obj = GCG->new(map=>$map,
		   seqname=>$seqname,
		   chrnum=>$chrnum,
		   beg=>$beg,
		   end=>$end,
		   rev=>$rev);


OR


my $obj = GCG->new(sequence=>$sequence);



###### $map = amap, a2map, a3map, nmap, n2map, n3map, pmap, p2map, p3map, rmap, or fmap


***** amap --- assemble single strand sequence (amap=GCG) 

***** a2map --- assemble single strand sequence (a2map=No header) 

***** a3map --- assemble single strand sequence (a3map=fasta format) 

***** nmap --- dna sequence without introns for orf (nmap=GCG) 

***** n2map --- dna sequence without introns for orf (n2map=No header)

***** n3map --- dna sequence without introns for orf (n3map=fasta format)

***** pmap --- translate dna sequence into protein (pmap=GCG) 

***** p2map --- translate dna sequence into protein (p2map=No header)

***** p3map --- translate dna sequence into protein (p3map=fasta format)

***** rmap --- restriction map

***** fmap --- restriction fragments

###### $seqname can be orf name or genBank sequence name.

###### $chrnum = chromosome number (1->17)

###### $beg = start_coord

###### $end = stop_coord

###### $rev = '-REV' or ''

###### $beg always less than $end

###### $sequence = DNA sequence -- for protein translation


=head1 Accessor Methods


=head2 getSequence

Usage:

my $seqfile = $obj->getSequence;

This method returns the temp protein/dna sequence file name for a given sequence name or for a given chromosome number plus start_coord and stop_coord or for the translated protein sequence.

=head1 Author

shuai@genome.stanford.edu

=cut