Mercurial > repos > xuebing > sharplabtool

comparison tools/human_genome_variation/gpass.pl @ 0:9071e359b9a3

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author xuebing
date Fri, 09 Mar 2012 19:37:19 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:9071e359b9a3
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