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