Mercurial > repos > dereeper > sniplay

comparison hapmap2mlmm/HapmapToMLMMFiles.pl @ 1:420b57c3c185 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author dereeper
date Fri, 10 Jul 2015 04:39:30 -0400
parents
children
comparison
equal deleted inserted replaced
0:3e19d0dfcf3e 1:420b57c3c185
1 #!/usr/bin/perl
2
3 use strict;
4 use Switch;
5 use Getopt::Long;
6
7 my $usage = qq~Usage:$0 <args> [<opts>]
8 where <args> are:
9 -h, --hapmap <Hapmap input file>
10 -m, --map <Map output file>
11 -g, --geno <Genotype output file>
12 -p, --path <Path for transpose executable>
13 ~;
14 $usage .= "\n";
15
16 my ($hapmap,$map,$geno,$path);
17
18
19 GetOptions(
20 "geno=s" => \$geno,
21 "map=s" => \$map,
22 "hapmap=s" => \$hapmap,
23 "path=s" => \$path,
24 );
25
26
27 die $usage
28 if ( !$geno || !$map || !$hapmap || !$path);
29
30 my $TRANSPOSE_EXE = "$path/transpose.awk";
31
32 my @snps;
33 my %chrom_pos;
34 my $num_line = 0;
35 open(my $O,">geno_transposed");
36 open(my $H,$hapmap);
37 while(<$H>)
38 {
39 $num_line++;
40 my $line = $_;
41 chomp($line);
42 $line =~s/\r//g;
43 $line =~s/\n//g;
44 my @infos = split(/\t/,$line);
45 if ($num_line == 1)
46 {
47 print $O "Ind_id";
48 for (my $i = 11; $i <= $#infos; $i++)
49 {
50 my $individual = $infos[$i];
51 print $O " " . $individual;
52 }
53 print $O "\n";
54 }
55 elsif ($num_line > 1)
56 {
57 my $snp = $infos[0];
58 my $variation = $infos[1];
59 my %scores;
60 if ($variation =~/(\w)\/(\w)/)
61 {
62 my $allele1 = $1;
63 my $allele2 = $2;
64 $scores{$allele1} = 0;
65 $scores{$allele2} = 1;
66 }
67 my $chrom = $infos[2];
68 my $pos = $infos[3];
69 $chrom_pos{$snp}{"chrom"} = $chrom;
70 $chrom_pos{$snp}{"pos"} = $pos;
71 push(@snps,$snp);
72 print $O "$snp";
73 for (my $i = 11; $i <= $#infos; $i++)
74 {
75 my $genotype = $infos[$i];
76 my @alleles = split("",$genotype);
77 if ($genotype ne "NN")
78 {
79 my $score = $scores{$alleles[0]} + $scores{$alleles[1]};
80 print $O " $score";
81 }
82 else
83 {
84 print $O " NA";
85 }
86 }
87 print $O "\n";
88 }
89 }
90 close($H);
91 close($O);
92
93 open(my $M,">$map");
94 print $M "SNP Chr Pos\n";
95 foreach my $snp(@snps)
96 {
97 print $M "$snp " . $chrom_pos{$snp}{"chrom"} . " ". $chrom_pos{$snp}{"pos"} . "\n";
98 }
99 close($M);
100
101 system("$TRANSPOSE_EXE geno_transposed >geno_transposed2");
102
103 open(my $F,">$geno");
104 open(my $G,"geno_transposed2");
105 while(<$G>)
106 {
107 my $line = $_;
108 $line =~s/ /\t/g;
109 print $F $line;
110 }
111 close($G);
112 close($F);
113
114 unlink("geno_transposed");
115 unlink("geno_transposed2");
116
117