[2] | 1 | #!/usr/bin/env perl |
---|
| 2 | |
---|
| 3 | use strict; |
---|
| 4 | use warnings; |
---|
| 5 | use File::Basename; |
---|
| 6 | use File::Temp qw/ tempfile /; |
---|
| 7 | |
---|
| 8 | $ENV{'PATH'} .= ':' . dirname($0); |
---|
| 9 | |
---|
| 10 | #this is a wrapper for gpass that converts a linkage pedigree file to input |
---|
| 11 | #for this program |
---|
| 12 | |
---|
| 13 | my($map, $ped, $out, $fdr) = @ARGV; |
---|
| 14 | |
---|
| 15 | if (!$map or !$ped or !$out or !$fdr) { die "missing args\n"; } |
---|
| 16 | |
---|
| 17 | my($fh, $name) = tempfile(); |
---|
| 18 | #by default this file is removed when these variable go out of scope |
---|
| 19 | print $fh "map=$map ped=$ped\n"; |
---|
| 20 | close $fh; #converter will overwrite, just keep name |
---|
| 21 | |
---|
| 22 | #run converter |
---|
| 23 | system("lped_to_geno.pl $map $ped > $name") == 0 |
---|
| 24 | or die "system lped_to_geno.pl $map $ped > $name failed\n"; |
---|
| 25 | |
---|
| 26 | #system("cp $name tmp.middle"); |
---|
| 27 | |
---|
| 28 | #run GPASS |
---|
| 29 | system("gpass $name -o $out -fdr $fdr 1>/dev/null") == 0 |
---|
| 30 | or die "system gpass $name -o $out -fdr $fdr, failed\n"; |
---|
| 31 | |
---|
| 32 | #merge SNP data with results |
---|
| 33 | merge(); |
---|
| 34 | |
---|
| 35 | exit; |
---|
| 36 | |
---|
| 37 | ######################################## |
---|
| 38 | |
---|
| 39 | #merge the input and output files so have SNP data with result |
---|
| 40 | sub merge { |
---|
| 41 | open(FH, $out) or die "Couldn't open $out, $!\n"; |
---|
| 42 | my %res; |
---|
| 43 | my @ind; |
---|
| 44 | while (<FH>) { |
---|
| 45 | chomp; |
---|
| 46 | my $line = $_; |
---|
| 47 | if ($line =~ /^(\d+)/) { $res{$1} = $line; push(@ind, $1); } |
---|
| 48 | else { $res{'index'} = $line; } |
---|
| 49 | } |
---|
| 50 | close FH; |
---|
| 51 | if (!@ind) { return; } #no results, leave alone |
---|
| 52 | @ind = sort { $a <=> $b } @ind; |
---|
| 53 | $res{'index'} =~ s/Index/#ID\tchr\tposition/; |
---|
| 54 | #read input file to get SNP data |
---|
| 55 | open(FH, $name) or die "Couldn't open $name, $!\n"; |
---|
| 56 | my $i = 0; #index is 0 based not counting header line |
---|
| 57 | my $c = shift @ind; |
---|
| 58 | while (<FH>) { |
---|
| 59 | chomp; |
---|
| 60 | if (/^ID/) { next; } |
---|
| 61 | my @f = split(/\s+/); |
---|
| 62 | if ($i == $c) { |
---|
| 63 | $res{$i} =~ s/^$i/$f[0]\t$f[1]\t$f[2]/; |
---|
| 64 | if (!@ind) { last; } |
---|
| 65 | $c = shift @ind; |
---|
| 66 | } |
---|
| 67 | $i++; |
---|
| 68 | } |
---|
| 69 | close FH; |
---|
| 70 | #now reprint results with SNP data included |
---|
| 71 | open(FH, ">", $out) or die "Couldn't write to $out, $!\n"; |
---|
| 72 | print FH $res{'index'}, "\n"; |
---|
| 73 | delete $res{'index'}; |
---|
| 74 | foreach $i (keys %res) { |
---|
| 75 | print FH $res{$i}, "\n"; |
---|
| 76 | } |
---|
| 77 | close FH; |
---|
| 78 | } |
---|
| 79 | |
---|