Mercurial > repos > dereeper > admixture

comparison Admixture.pl @ 5:97c9c8daa3c3 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload
author dereeper
date Wed, 13 Apr 2016 07:51:13 -0400
parents 58df6910f1c3
children 3f66f32dc5d9
comparison
equal deleted inserted replaced
4:58df6910f1c3 5:97c9c8daa3c3
1 #!/usr/bin/perl 1 #!/usr/bin/perl
2 2
3 use strict; 3 use strict;
4 use Switch;
4 use Getopt::Long; 5 use Getopt::Long;
6 use Bio::SeqIO;
5 use File::Basename; 7 use File::Basename;
6 8
7 my $usage = qq~Usage:$0 <args> [<opts>] 9 my $usage = qq~Usage:$0 <args> [<opts>]
8 where <args> are: 10 where <args> are:
9 -i, --input <input basename (for PED, FAM, BIM)> 11 -i, --input <input basename (for PED, FAM, BIM)>
10 -o, --output <output> 12 -o, --output <output>
11 -k, --kmin <K min. int> 13 -k, --kmin <K min. int>
12 -m, --maxK <K max. int> 14 -m, --maxK <K max. int>
13 -d, --directory <temporary directory> 15 -d, --directory <temporary directory>
16 -t, --threshold <threshold admixture proportion for group assignation>
14 ~; 17 ~;
15 $usage .= "\n"; 18 $usage .= "\n";
16 19
17 my ($input,$output,$kmin,$kmax,$directory); 20 my ($input,$output,$kmin,$kmax,$directory,$threshold);
18 21
19 22
20 GetOptions( 23 GetOptions(
21 "input=s" => \$input, 24 "input=s" => \$input,
22 "output=s" => \$output, 25 "output=s" => \$output,
23 "kmin=s" => \$kmin, 26 "kmin=s" => \$kmin,
24 "maxK=s" => \$kmax, 27 "maxK=s" => \$kmax,
25 "directory=s" => \$directory, 28 "directory=s" => \$directory,
29 "threshold" => \$threshold,
26 ); 30 );
27 31
28 32
29 die $usage 33 die $usage
30 if ( !$input || !$output || !$kmin || !$kmax || !$directory ); 34 if ( !$input || !$output || !$kmin || !$kmax || !$directory );
61 # launch admixture for different K 65 # launch admixture for different K
62 ################################### 66 ###################################
63 my %errors; 67 my %errors;
64 for (my $k = $kmin; $k <= $kmax; $k++) 68 for (my $k = $kmin; $k <= $kmax; $k++)
65 { 69 {
66 system("admixture --cv $input.bed $k >>$directory/log.$k 2>&1"); 70 system("/apps/www/sniplay.cirad.fr/tools/admixture/admixture_linux-1.23/admixture --cv $input.bed $k >>$directory/log.$k 2>&1");
67 my $cv_error_line = `grep -h CV $directory/log.$k`; 71 my $cv_error_line = `grep -h CV $directory/log.$k`;
68 if ($cv_error_line =~/: (\d+\.*\d*)$/) 72 if ($cv_error_line =~/: (\d+\.*\d*)$/)
69 { 73 {
70 $errors{$1} = $k; 74 $errors{$1} = $k;
71 } 75 }
72 system("cat $directory/log.$k >>$directory/logs"); 76 system("cat $directory/log.$k >>$directory/logs");
73 system("echo '\n\n====================================\n\n' >>$directory/logs"); 77 system("echo '\n\n====================================\n\n' >>$directory/logs");
74 78
75 open(my $O2,">$basename.$k.final.Q"); 79 open(my $O2,">$basename.$k.final.Q");
80 open(my $O3,">$directory/groups.$k");
76 open(my $O,"$basename.$k.Q"); 81 open(my $O,"$basename.$k.Q");
77 my %hash_groupes; 82 my %hash_groupes;
78 my %hash_indv; 83 my %hash_indv;
79 my %group_of_ind; 84 my %group_of_ind;
80 my $i = 0; 85 my $i = 0;
86 my @infos = split(/\s+/,$line); 91 my @infos = split(/\s+/,$line);
87 my $group = "admix"; 92 my $group = "admix";
88 my $ind = $individus[$i]; 93 my $ind = $individus[$i];
89 for (my $j = 0; $j <$k; $j++){ 94 for (my $j = 0; $j <$k; $j++){
90 my $val = $infos[$j]; 95 my $val = $infos[$j];
91 if ($val > 0.5){$group = "Q$j";} 96 if ($val > ($threshold/100)){$group = "Q$j";}
92 } 97 }
93 if ($ind){ 98 if ($ind){
94 $hash_indv{$ind} = join(" ",@infos); 99 $hash_indv{$ind} = join(" ",@infos);
95 $hash_groupes{$group}{"ind"} .= ",".$ind; 100 $hash_groupes{$group}{"ind"} .= ",".$ind;
96 $group_of_ind{$ind} = $group; 101 $group_of_ind{$ind} = $group;
100 105
101 foreach my $group(sort keys(%hash_groupes)){ 106 foreach my $group(sort keys(%hash_groupes)){
102 my @inds = split(",",$hash_groupes{$group}{"ind"}); 107 my @inds = split(",",$hash_groupes{$group}{"ind"});
103 foreach my $ind(@inds){ 108 foreach my $ind(@inds){
104 if ($ind =~/\w+/){ 109 if ($ind =~/\w+/){
105 #print $O3 "$ind;$group\n"; 110 print $O3 "$ind;$group\n";
106 print $O2 $ind." ".$hash_indv{$ind}. "\n"; 111 print $O2 $ind." ".$hash_indv{$ind}. "\n";
107 } 112 }
108 } 113 }
109 } 114 }
110 #close($O3); 115 close($O3);
111 close($O2); 116 close($O2);
112 117
113 system("cat $basename.$k.final.Q >>$directory/outputs.Q"); 118 system("cat $basename.$k.final.Q >>$directory/outputs.Q");
114 system("echo '\n\n====================================\n\n' >>$directory/outputs.Q"); 119 system("echo '\n\n====================================\n\n' >>$directory/outputs.Q");
115 system("cat $basename.$k.P >>$directory/outputs.P"); 120 system("cat $basename.$k.P >>$directory/outputs.P");
139 } 144 }
140 close(BEST1); 145 close(BEST1);
141 close(BEST2); 146 close(BEST2);
142 147
143 system("cp -rf $directory/log.$best_K $directory/log"); 148 system("cp -rf $directory/log.$best_K $directory/log");
149 system("cp -rf $directory/groups.$best_K $directory/groups");
144 150
145 151
146 152
147