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 $
17 genotypeizer - converts varoius genotype files into various formats
21 genotypeizer [options] [files to read]
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
30 The files to read should be either map/ped pairs or a single snpassoc
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.]
44 =item B<--output-format>
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
52 Debug verbosity. (Default 0)
56 Display brief usage information.
69 use POSIX qw(:fcntl_h);
74 my %options = (output_format => 'R',
85 'output_format|output-format=s',
87 'debug|d+','help|h|?','man|m');
89 pod2usage() if $options{help};
90 pod2usage({verbose=>2}) if $options{man};
92 $DEBUG = $options{debug};
96 if ($options{output_format} !~ /^(R|pedmap|mapped)$/i) {
97 push @USAGE_ERRORS,"--output-format must be one of R or pedmap.";
100 $options{format} = lc($options{format});
102 if ($options{format} eq 'mapped') {
103 $options{format} = 'pedmap';
106 pod2usage(join("\n",@USAGE_ERRORS)) if @USAGE_ERRORS;
109 # open up output filehandles
113 my %file_ext= (pedmap => {ped => '.ped',
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}};
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: $!";
129 $output_fh{$file_ext} = \*STDOUT;
133 # read it in; currently only support map/ped, but eventually do them
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: $!";
151 my ($chr,$rsid,$something,$position) = split /\t/;
152 push @map,{chr => $chr,
154 something => $something,
155 position => $position,
157 $snps{$rsid} = {chr => $chr,
159 something => $something,
160 position => $position,
163 push @snp_order,@map;
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,
172 family_id => $family_id,
173 individual_id => $individual_id,
174 maternal_id => $maternal_id,
175 paternal_id => $paternal_id,
177 for my $snp_info (@map) {
178 $samples{$individual_id}{$snp_info->{rsid}} = [shift @snps, shift @snps];
181 die "The map file and the number of snps don't appear to match".join(' ',@snps)."were left over";
188 @valid_snps{map {s/^rs//; "rs$_";} @{$options{snps}}} = (1) x @{$options{snps}};
191 my @snp_ordering = map {$_->{rsid}} @snp_order;
193 if (keys %valid_snps) {
194 @snp_ordering = grep {exists $valid_snps{$_}} @snp_ordering;
198 @valid_snps{@{$options{snps}}} = (1) x $options{snps} if defined $options{snps} and @{$options{snps}};
200 if ($options{output_format} eq 'pedmap') {
201 if (keys %valid_snps) {
202 %snps = map {($_,$snps{$_})} keys %valid_snps;
205 for my $snp (@snp_ordering) {
206 print {$output_fh{map}} join("\t",@{$snps{$snp}}{qw(chr rsid something position)})."\n";
209 for my $individual_id (sort keys %extra_info) {
210 print {$output_fh{ped}}
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}
219 print {$output_fh{main}}
221 qw(individual_id family_id paternal_id maternal_id sex phenotype),
224 print {$output_fh{main}} "\n";
225 for my $individual_id (keys %samples) {
226 print {$output_fh{main}}
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') {
233 join('',map{uc($_)} @$a);
237 print {$output_fh{main}} "\n";