#! /usr/bin/perl
# genotypeizer converts various genotype files into various formats,
# and is released under the terms of the GPL version 2, or any later
# version, at your option. See the file README and COPYING for more
# information. Copyright 2009 by Don Armstrong <don@donarmstrong.com>.
# $Id: perl_script 1432 2009-04-21 02:42:41Z don $


use warnings;
use strict;

use Getopt::Long;
use Pod::Usage;

=head1 NAME

genotypeizer - converts varoius genotype files into various formats

=head1 SYNOPSIS

genotypeizer [options] [files to read]

 Options:
  --output file to output to; defaults to stdout
  --output-format format to output to; defaults to R
  --debug, -d debugging level (Default 0)
  --help, -h display this help
  --man, -m display manual

The files to read should be either map/ped pairs or a single snpassoc
R file.

=head1 OPTIONS

=over

=item B<--output>

File to output to. If the file format that you choose requires
multiple files, the appropriate extension(s) will be appended to the
name provided. [If --output isn't provided, STDOUT will be used, which
will produce suboptimal results in the case of pedmap format.]

=item B<--output-format>

Output format to output. Currently valid values are R or pedmap. (R is
the format used by the R snpassoc module, and pedmap is the format
used by phylip.

=item B<--debug, -d>

Debug verbosity. (Default 0)

=item B<--help, -h>

Display brief usage information.

=item B<--man, -m>

Display this manual.

=back

=head1 EXAMPLES


=cut

use POSIX qw(:fcntl_h);
use IO::File;

use vars qw($DEBUG);

my %options = (output_format   => 'R',
	       output          => undef,
	       snps            => undef,
	       debug           => 0,
	       help            => 0,
	       man             => 0,
	       force           => 0,
	       );

GetOptions(\%options,
	   'output=s',
	   'output_format|output-format=s',
	   'snps=s@',
	   'debug|d+','help|h|?','man|m');

pod2usage() if $options{help};
pod2usage({verbose=>2}) if $options{man};

$DEBUG = $options{debug};

my @USAGE_ERRORS;

if ($options{output_format} !~ /^(R|pedmap|mapped)$/i) {
     push @USAGE_ERRORS,"--output-format must be one of R or pedmap.";
}

$options{format} = lc($options{format});

if ($options{format} eq 'mapped') {
    $options{format} = 'pedmap';
}

pod2usage(join("\n",@USAGE_ERRORS)) if @USAGE_ERRORS;


# open up output filehandles

my %output_fh;

my %file_ext= (pedmap => {ped => '.ped',
			  map => '.map',
			 },
	       R      => {main => '',
			 },
	      );
die "format $options{output_format} doesn't appear to actually be supported; this is a bug."
    unless exists $file_ext{$options{output_format}};

for my $file_ext (keys %{$file_ext{$options{output_format}}}) {
    if (defined $options{output}) {
	$output_fh{$file_ext} = IO::File->new($options{output}.$file_ext{$options{output_format}}{$file_ext},
					      O_CREAT|O_EXCL|O_WRONLY)
	    or die "Unable to open $options{output}"."$file_ext{$options{output_format}}{$file_ext} for writing: $!";
    }
    else {
	$output_fh{$file_ext} = \*STDOUT;
    }
}

# read it in; currently only support map/ped, but eventually do them
# all.

my %samples;
my %extra_info;
my %snps;

my @map_ped = @ARGV;

my @snp_order;

my ($map_file,$ped_file);
while (($map_file,$ped_file) = splice @map_ped,0,2) {
     my $map_fh = IO::File->new($map_file,'r') or die "Unable to open $map_file for reading: $!";
     my $ped_fh = IO::File->new($ped_file,'r') or die "Unable to open $ped_file for reading: $!";
     my @map = ();
     while (<$map_fh>) {
	  chomp;
	  my ($chr,$rsid,$something,$position) = split /\t/;
	  push @map,{chr  => $chr,
		     rsid => $rsid,
		     something => $something,
		     position => $position,
		    };
	  $snps{$rsid} = {chr  => $chr,
			  rsid => $rsid,
			  something => $something,
			  position => $position,
			 };
     }
     push @snp_order,@map;

     while (<$ped_fh>) {
	  chomp;
	  next unless /^\d+/;
	  my @row = split /\s+/;
	  my ($family_id,$individual_id,$paternal_id,$maternal_id,$sex,$phenotype,@snps) = @row;
	  $extra_info{$individual_id} = {phenotype     => $phenotype,
					 sex           => $sex,
					 family_id     => $family_id,
					 individual_id => $individual_id,
					 maternal_id   => $maternal_id,
					 paternal_id   => $paternal_id,
					};
	  for my $snp_info (@map) {
	      $samples{$individual_id}{$snp_info->{rsid}} = [shift @snps, shift @snps];
	  }
	  if (@snps) {
	       die "The map file and the number of snps don't appear to match".join(' ',@snps)."were left over";
	  }
     }
}


my %valid_snps;
@valid_snps{map {s/^rs//; "rs$_";} @{$options{snps}}} = (1) x @{$options{snps}};


my @snp_ordering = map {$_->{rsid}} @snp_order;

if (keys %valid_snps) {
    @snp_ordering = grep {exists $valid_snps{$_}} @snp_ordering;
}


@valid_snps{@{$options{snps}}} = (1) x $options{snps} if defined $options{snps} and @{$options{snps}};

if ($options{output_format} eq 'pedmap') {
    if (keys %valid_snps) {
	%snps = map {($_,$snps{$_})} keys %valid_snps;
    }
    # output the map
    for my $snp (@snp_ordering) {
	print {$output_fh{map}} join("\t",@{$snps{$snp}}{qw(chr rsid something position)})."\n";
    }
    # output the ped
    for my $individual_id (sort keys %extra_info) {
	print {$output_fh{ped}}
	    join("\t",
		 @{$extra_info{$individual_id}}{qw(family_id individual_id),
						    qw(paternal_id maternal_id sex phenotype)},
		 map {@{$_}} @{$samples{$individual_id}}{@snp_ordering}
		)."\n";
    }
}
else {
    print {$output_fh{main}}
	join("\t",
	     qw(individual_id family_id paternal_id maternal_id sex phenotype),
	     sort keys %snps
	    );
    print {$output_fh{main}} "\n";
    for my $individual_id (keys %samples) {
	print {$output_fh{main}}
	    join("\t",
		 @{$extra_info{$individual_id}}{qw(individual_id family_id paternal_id maternal_id sex phenotype)},
		 map {my $a = $samples{$individual_id}{$_};
		      if (not defined $a or not @$a or $a->[0] eq '0' or $a->[1] eq '0') {
			  'NA'
		      } else {
			  join('',map{uc($_)} @$a);
		      }
		  } sort keys %snps
		);
	print {$output_fh{main}} "\n";
    }
}


__END__