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 | |
---|