|
4
|
1 #!/usr/bin/perl
|
|
|
2
|
|
|
3 use strict;
|
|
|
4 use Getopt::Long;
|
|
|
5 use File::Basename;
|
|
|
6
|
|
|
7 my $usage = qq~Usage:$0 <args> [<opts>]
|
|
|
8 where <args> are:
|
|
|
9 -i, --input <input basename (for PED, FAM, BIM)>
|
|
|
10 -o, --output <output>
|
|
|
11 -k, --kmin <K min. int>
|
|
|
12 -m, --maxK <K max. int>
|
|
|
13 -d, --directory <temporary directory>
|
|
5
|
14 -t, --threshold <threshold admixture proportion for group assignation>
|
|
4
|
15 ~;
|
|
|
16 $usage .= "\n";
|
|
|
17
|
|
5
|
18 my ($input,$output,$kmin,$kmax,$directory,$threshold);
|
|
4
|
19
|
|
|
20
|
|
|
21 GetOptions(
|
|
|
22 "input=s" => \$input,
|
|
|
23 "output=s" => \$output,
|
|
|
24 "kmin=s" => \$kmin,
|
|
|
25 "maxK=s" => \$kmax,
|
|
|
26 "directory=s" => \$directory,
|
|
5
|
27 "threshold" => \$threshold,
|
|
4
|
28 );
|
|
|
29
|
|
|
30
|
|
|
31 die $usage
|
|
|
32 if ( !$input || !$output || !$kmin || !$kmax || !$directory );
|
|
|
33
|
|
|
34 if ($kmin =~/^(\d+)\s*$/){
|
|
|
35 $kmin = $1;
|
|
|
36 }
|
|
|
37 else{
|
|
|
38 die "Error: kmin must be an integer\n";
|
|
|
39 }
|
|
|
40 if ($kmax =~/^(\d+)\s*$/){
|
|
|
41 $kmax = $1;
|
|
|
42 }
|
|
|
43 else{
|
|
|
44 die "Error: kmax must be an integer\n";
|
|
|
45 }
|
|
|
46
|
|
|
47
|
|
|
48
|
|
|
49 open(my $P2,"$input.fam");
|
|
|
50 my $n = 0;
|
|
|
51 my $ind_num = 0;
|
|
|
52 my @individus;
|
|
|
53 while(<$P2>)
|
|
|
54 {
|
|
|
55 my ($ind,$other) = split(/ /,$_);
|
|
|
56 push(@individus,$ind);
|
|
|
57 }
|
|
|
58 close($P2);
|
|
|
59
|
|
|
60 my $basename = basename($input);
|
|
|
61
|
|
|
62 ###################################
|
|
|
63 # launch admixture for different K
|
|
|
64 ###################################
|
|
|
65 my %errors;
|
|
|
66 for (my $k = $kmin; $k <= $kmax; $k++)
|
|
|
67 {
|
|
6
|
68 system("admixture --cv $input.bed $k >>$directory/log.$k 2>&1");
|
|
4
|
69 my $cv_error_line = `grep -h CV $directory/log.$k`;
|
|
|
70 if ($cv_error_line =~/: (\d+\.*\d*)$/)
|
|
|
71 {
|
|
|
72 $errors{$1} = $k;
|
|
|
73 }
|
|
|
74 system("cat $directory/log.$k >>$directory/logs");
|
|
|
75 system("echo '\n\n====================================\n\n' >>$directory/logs");
|
|
|
76
|
|
|
77 open(my $O2,">$basename.$k.final.Q");
|
|
5
|
78 open(my $O3,">$directory/groups.$k");
|
|
4
|
79 open(my $O,"$basename.$k.Q");
|
|
|
80 my %hash_groupes;
|
|
|
81 my %hash_indv;
|
|
|
82 my %group_of_ind;
|
|
|
83 my $i = 0;
|
|
|
84 while (<$O>){
|
|
|
85 $i++;
|
|
|
86 my $line = $_;
|
|
|
87 $line =~s/\n//g;
|
|
|
88 $line =~s/\r//g;
|
|
|
89 my @infos = split(/\s+/,$line);
|
|
|
90 my $group = "admix";
|
|
|
91 my $ind = $individus[$i];
|
|
|
92 for (my $j = 0; $j <$k; $j++){
|
|
|
93 my $val = $infos[$j];
|
|
5
|
94 if ($val > ($threshold/100)){$group = "Q$j";}
|
|
4
|
95 }
|
|
|
96 if ($ind){
|
|
|
97 $hash_indv{$ind} = join(" ",@infos);
|
|
|
98 $hash_groupes{$group}{"ind"} .= ",".$ind;
|
|
|
99 $group_of_ind{$ind} = $group;
|
|
|
100 }
|
|
|
101 }
|
|
|
102 close($O);
|
|
|
103
|
|
|
104 foreach my $group(sort keys(%hash_groupes)){
|
|
|
105 my @inds = split(",",$hash_groupes{$group}{"ind"});
|
|
|
106 foreach my $ind(@inds){
|
|
|
107 if ($ind =~/\w+/){
|
|
5
|
108 print $O3 "$ind;$group\n";
|
|
4
|
109 print $O2 $ind." ".$hash_indv{$ind}. "\n";
|
|
|
110 }
|
|
|
111 }
|
|
|
112 }
|
|
5
|
113 close($O3);
|
|
4
|
114 close($O2);
|
|
|
115
|
|
|
116 system("cat $basename.$k.final.Q >>$directory/outputs.Q");
|
|
|
117 system("echo '\n\n====================================\n\n' >>$directory/outputs.Q");
|
|
|
118 system("cat $basename.$k.P >>$directory/outputs.P");
|
|
|
119 system("echo '\n\n====================================\n\n' >>$directory/outputs.P");
|
|
|
120 }
|
|
|
121
|
|
|
122 my @sorted_errors = sort {$a<=>$b} keys(%errors);
|
|
|
123 my $best_K = $errors{@sorted_errors[0]};
|
|
|
124
|
|
|
125
|
|
|
126 open(BEST1,"$basename.$best_K.final.Q");
|
|
|
127 open(BEST2,">$directory/output");
|
|
|
128 print BEST2 "<Covariate>\n";
|
|
|
129 print BEST2 "<Trait>";
|
|
|
130 for (my $j=1;$j<=$best_K;$j++)
|
|
|
131 {
|
|
|
132 print BEST2 " Q" . $j;
|
|
|
133 }
|
|
|
134 print BEST2 "\n";
|
|
|
135 my $i = 0;
|
|
|
136 while(<BEST1>)
|
|
|
137 {
|
|
|
138 my $line = $_;
|
|
|
139 $line =~s/ /\t/g;
|
|
|
140 print BEST2 $line;
|
|
|
141 $i++;
|
|
|
142 }
|
|
|
143 close(BEST1);
|
|
|
144 close(BEST2);
|
|
|
145
|
|
|
146 system("cp -rf $directory/log.$best_K $directory/log");
|
|
5
|
147 system("cp -rf $directory/groups.$best_K $directory/groups");
|
|
4
|
148
|
|
|
149
|
|
|
150
|
