Mercurial > repos > xuebing > sharplabtool
diff tools/human_genome_variation/gpass.pl @ 0:9071e359b9a3
Uploaded
author | xuebing |
---|---|
date | Fri, 09 Mar 2012 19:37:19 -0500 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tools/human_genome_variation/gpass.pl Fri Mar 09 19:37:19 2012 -0500 @@ -0,0 +1,79 @@ +#!/usr/bin/env perl + +use strict; +use warnings; +use File::Basename; +use File::Temp qw/ tempfile /; + +$ENV{'PATH'} .= ':' . dirname($0); + +#this is a wrapper for gpass that converts a linkage pedigree file to input +#for this program + +my($map, $ped, $out, $fdr) = @ARGV; + +if (!$map or !$ped or !$out or !$fdr) { die "missing args\n"; } + +my($fh, $name) = tempfile(); +#by default this file is removed when these variable go out of scope +print $fh "map=$map ped=$ped\n"; +close $fh; #converter will overwrite, just keep name + +#run converter +system("lped_to_geno.pl $map $ped > $name") == 0 + or die "system lped_to_geno.pl $map $ped > $name failed\n"; + +#system("cp $name tmp.middle"); + +#run GPASS +system("gpass $name -o $out -fdr $fdr 1>/dev/null") == 0 + or die "system gpass $name -o $out -fdr $fdr, failed\n"; + +#merge SNP data with results +merge(); + +exit; + +######################################## + +#merge the input and output files so have SNP data with result +sub merge { + open(FH, $out) or die "Couldn't open $out, $!\n"; + my %res; + my @ind; + while (<FH>) { + chomp; + my $line = $_; + if ($line =~ /^(\d+)/) { $res{$1} = $line; push(@ind, $1); } + else { $res{'index'} = $line; } + } + close FH; + if (!@ind) { return; } #no results, leave alone + @ind = sort { $a <=> $b } @ind; + $res{'index'} =~ s/Index/#ID\tchr\tposition/; + #read input file to get SNP data + open(FH, $name) or die "Couldn't open $name, $!\n"; + my $i = 0; #index is 0 based not counting header line + my $c = shift @ind; + while (<FH>) { + chomp; + if (/^ID/) { next; } + my @f = split(/\s+/); + if ($i == $c) { + $res{$i} =~ s/^$i/$f[0]\t$f[1]\t$f[2]/; + if (!@ind) { last; } + $c = shift @ind; + } + $i++; + } + close FH; + #now reprint results with SNP data included + open(FH, ">", $out) or die "Couldn't write to $out, $!\n"; + print FH $res{'index'}, "\n"; + delete $res{'index'}; + foreach $i (keys %res) { + print FH $res{$i}, "\n"; + } + close FH; +} +