root/galaxy-central/tools/human_genome_variation/gpass.pl

リビジョン 2, 2.0 KB (コミッタ: hatakeyama, 14 年 前)

import galaxy-central

  • 属性 svn:executable の設定値 *
行番号 
1#!/usr/bin/env perl
2
3use strict;
4use warnings;
5use File::Basename;
6use 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
13my($map, $ped, $out, $fdr) = @ARGV;
14
15if (!$map or !$ped or !$out or !$fdr) { die "missing args\n"; }
16
17my($fh, $name) = tempfile();
18#by default this file is removed when these variable go out of scope
19print $fh "map=$map ped=$ped\n";
20close $fh;  #converter will overwrite, just keep name
21
22#run converter
23system("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
29system("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
33merge();
34
35exit;
36
37########################################
38
39#merge the input and output files so have SNP data with result
40sub 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
Note: リポジトリブラウザについてのヘルプは TracBrowser を参照してください。