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