genotypeizer

   1 #! /usr/bin/perl
   2 # genotypeizer converts various genotype files into various formats,
   3 # and is released under the terms of the GPL version 2, or any later
   4 # version, at your option. See the file README and COPYING for more
   5 # information. Copyright 2009 by Don Armstrong <don@donarmstrong.com>.
   6 # $Id: perl_script 1432 2009-04-21 02:42:41Z don $
   7
   8
   9 use warnings;
  10 use strict;
  11
  12 use Getopt::Long;
  13 use Pod::Usage;
  14
  15 =head1 NAME
  16
  17 genotypeizer - converts varoius genotype files into various formats
  18
  19 =head1 SYNOPSIS
  20
  21 genotypeizer [options] [files to read]
  22
  23  Options:
  24   --output file to output to; defaults to stdout
  25   --output-format format to output to; defaults to R
  26   --debug, -d debugging level (Default 0)
  27   --help, -h display this help
  28   --man, -m display manual
  29
  30 The files to read should be either map/ped pairs or a single snpassoc
  31 R file.
  32
  33 =head1 OPTIONS
  34
  35 =over
  36
  37 =item B<--output>
  38
  39 File to output to. If the file format that you choose requires
  40 multiple files, the appropriate extension(s) will be appended to the
  41 name provided. [If --output isn't provided, STDOUT will be used, which
  42 will produce suboptimal results in the case of pedmap format.]
  43
  44 =item B<--output-format>
  45
  46 Output format to output. Currently valid values are R or pedmap. (R is
  47 the format used by the R snpassoc module, and pedmap is the format
  48 used by phylip.
  49
  50 =item B<--debug, -d>
  51
  52 Debug verbosity. (Default 0)
  53
  54 =item B<--help, -h>
  55
  56 Display brief usage information.
  57
  58 =item B<--man, -m>
  59
  60 Display this manual.
  61
  62 =back
  63
  64 =head1 EXAMPLES
  65
  66
  67 =cut
  68
  69 use POSIX qw(:fcntl_h);
  70 use IO::File;
  71
  72 use vars qw($DEBUG);
  73
  74 my %options = (output_format   => 'R',
  75                output          => undef,
  76                snps            => undef,
  77                debug           => 0,
  78                help            => 0,
  79                man             => 0,
  80                force           => 0,
  81                );
  82
  83 GetOptions(\%options,
  84            'output=s',
  85            'output_format|output-format=s',
  86            'snps=s@',
  87            'debug|d+','help|h|?','man|m');
  88
  89 pod2usage() if $options{help};
  90 pod2usage({verbose=>2}) if $options{man};
  91
  92 $DEBUG = $options{debug};
  93
  94 my @USAGE_ERRORS;
  95
  96 if ($options{output_format} !~ /^(R|pedmap|mapped)$/i) {
  97      push @USAGE_ERRORS,"--output-format must be one of R or pedmap.";
  98 }
  99
 100 $options{format} = lc($options{format});
 101
 102 if ($options{format} eq 'mapped') {
 103     $options{format} = 'pedmap';
 104 }
 105
 106 pod2usage(join("\n",@USAGE_ERRORS)) if @USAGE_ERRORS;
 107
 108
 109 # open up output filehandles
 110
 111 my %output_fh;
 112
 113 my %file_ext= (pedmap => {ped => '.ped',
 114                           map => '.map',
 115                          },
 116                R      => {main => '',
 117                          },
 118               );
 119 die "format $options{output_format} doesn't appear to actually be supported; this is a bug."
 120     unless exists $file_ext{$options{output_format}};
 121
 122 for my $file_ext (keys %{$file_ext{$options{output_format}}}) {
 123     if (defined $options{output}) {
 124         $output_fh{$file_ext} = IO::File->new($options{output}.$file_ext{$options{output_format}}{$file_ext},
 125                                               O_CREAT|O_EXCL|O_WRONLY)
 126             or die "Unable to open $options{output}"."$file_ext{$options{output_format}}{$file_ext} for writing: $!";
 127     }
 128     else {
 129         $output_fh{$file_ext} = \*STDOUT;
 130     }
 131 }
 132
 133 # read it in; currently only support map/ped, but eventually do them
 134 # all.
 135
 136 my %samples;
 137 my %extra_info;
 138 my %snps;
 139
 140 my @map_ped = @ARGV;
 141
 142 my @snp_order;
 143
 144 my ($map_file,$ped_file);
 145 while (($map_file,$ped_file) = splice @map_ped,0,2) {
 146      my $map_fh = IO::File->new($map_file,'r') or die "Unable to open $map_file for reading: $!";
 147      my $ped_fh = IO::File->new($ped_file,'r') or die "Unable to open $ped_file for reading: $!";
 148      my @map = ();
 149      while (<$map_fh>) {
 150           chomp;
 151           my ($chr,$rsid,$something,$position) = split /\t/;
 152           push @map,{chr  => $chr,
 153                      rsid => $rsid,
 154                      something => $something,
 155                      position => $position,
 156                     };
 157           $snps{$rsid} = {chr  => $chr,
 158                           rsid => $rsid,
 159                           something => $something,
 160                           position => $position,
 161                          };
 162      }
 163      push @snp_order,@map;
 164
 165      while (<$ped_fh>) {
 166           chomp;
 167           next unless /^\d+/;
 168           my @row = split /\s+/;
 169           my ($family_id,$individual_id,$paternal_id,$maternal_id,$sex,$phenotype,@snps) = @row;
 170           $extra_info{$individual_id} = {phenotype     => $phenotype,
 171                                          sex           => $sex,
 172                                          family_id     => $family_id,
 173                                          individual_id => $individual_id,
 174                                          maternal_id   => $maternal_id,
 175                                          paternal_id   => $paternal_id,
 176                                         };
 177           for my $snp_info (@map) {
 178               $samples{$individual_id}{$snp_info->{rsid}} = [shift @snps, shift @snps];
 179           }
 180           if (@snps) {
 181                die "The map file and the number of snps don't appear to match".join(' ',@snps)."were left over";
 182           }
 183      }
 184 }
 185
 186
 187 my %valid_snps;
 188 @valid_snps{map {s/^rs//; "rs$_";} @{$options{snps}}} = (1) x @{$options{snps}};
 189
 190
 191 my @snp_ordering = map {$_->{rsid}} @snp_order;
 192
 193 if (keys %valid_snps) {
 194     @snp_ordering = grep {exists $valid_snps{$_}} @snp_ordering;
 195 }
 196
 197
 198 @valid_snps{@{$options{snps}}} = (1) x $options{snps} if defined $options{snps} and @{$options{snps}};
 199
 200 if ($options{output_format} eq 'pedmap') {
 201     if (keys %valid_snps) {
 202         %snps = map {($_,$snps{$_})} keys %valid_snps;
 203     }
 204     # output the map
 205     for my $snp (@snp_ordering) {
 206         print {$output_fh{map}} join("\t",@{$snps{$snp}}{qw(chr rsid something position)})."\n";
 207     }
 208     # output the ped
 209     for my $individual_id (sort keys %extra_info) {
 210         print {$output_fh{ped}}
 211             join("\t",
 212                  @{$extra_info{$individual_id}}{qw(family_id individual_id),
 213                                                     qw(paternal_id maternal_id sex phenotype)},
 214                  map {@{$_}} @{$samples{$individual_id}}{@snp_ordering}
 215                 )."\n";
 216     }
 217 }
 218 else {
 219     print {$output_fh{main}}
 220         join("\t",
 221              qw(individual_id family_id paternal_id maternal_id sex phenotype),
 222              sort keys %snps
 223             );
 224     print {$output_fh{main}} "\n";
 225     for my $individual_id (keys %samples) {
 226         print {$output_fh{main}}
 227             join("\t",
 228                  @{$extra_info{$individual_id}}{qw(individual_id family_id paternal_id maternal_id sex phenotype)},
 229                  map {my $a = $samples{$individual_id}{$_};
 230                       if (not defined $a or not @$a or $a->[0] eq '0' or $a->[1] eq '0') {
 231                           'NA'
 232                       } else {
 233                           join('',map{uc($_)} @$a);
 234                       }
 235                   } sort keys %snps
 236                 );
 237         print {$output_fh{main}} "\n";
 238     }
 239 }
 240
 241
 242 __END__