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