Mercurial > repos > dereeper > admixture
diff Admixture.pl @ 5:97c9c8daa3c3 draft
planemo upload
author | dereeper |
---|---|
date | Wed, 13 Apr 2016 07:51:13 -0400 |
parents | 58df6910f1c3 |
children | 3f66f32dc5d9 |
line wrap: on
line diff
--- a/Admixture.pl Tue Apr 12 09:31:56 2016 -0400 +++ b/Admixture.pl Wed Apr 13 07:51:13 2016 -0400 @@ -1,7 +1,9 @@ #!/usr/bin/perl use strict; +use Switch; use Getopt::Long; +use Bio::SeqIO; use File::Basename; my $usage = qq~Usage:$0 <args> [<opts>] @@ -11,10 +13,11 @@ -k, --kmin <K min. int> -m, --maxK <K max. int> -d, --directory <temporary directory> + -t, --threshold <threshold admixture proportion for group assignation> ~; $usage .= "\n"; -my ($input,$output,$kmin,$kmax,$directory); +my ($input,$output,$kmin,$kmax,$directory,$threshold); GetOptions( @@ -23,6 +26,7 @@ "kmin=s" => \$kmin, "maxK=s" => \$kmax, "directory=s" => \$directory, + "threshold" => \$threshold, ); @@ -63,7 +67,7 @@ my %errors; for (my $k = $kmin; $k <= $kmax; $k++) { - system("admixture --cv $input.bed $k >>$directory/log.$k 2>&1"); + system("/apps/www/sniplay.cirad.fr/tools/admixture/admixture_linux-1.23/admixture --cv $input.bed $k >>$directory/log.$k 2>&1"); my $cv_error_line = `grep -h CV $directory/log.$k`; if ($cv_error_line =~/: (\d+\.*\d*)$/) { @@ -73,6 +77,7 @@ system("echo '\n\n====================================\n\n' >>$directory/logs"); open(my $O2,">$basename.$k.final.Q"); + open(my $O3,">$directory/groups.$k"); open(my $O,"$basename.$k.Q"); my %hash_groupes; my %hash_indv; @@ -88,7 +93,7 @@ my $ind = $individus[$i]; for (my $j = 0; $j <$k; $j++){ my $val = $infos[$j]; - if ($val > 0.5){$group = "Q$j";} + if ($val > ($threshold/100)){$group = "Q$j";} } if ($ind){ $hash_indv{$ind} = join(" ",@infos); @@ -102,12 +107,12 @@ my @inds = split(",",$hash_groupes{$group}{"ind"}); foreach my $ind(@inds){ if ($ind =~/\w+/){ - #print $O3 "$ind;$group\n"; + print $O3 "$ind;$group\n"; print $O2 $ind." ".$hash_indv{$ind}. "\n"; } } } - #close($O3); + close($O3); close($O2); system("cat $basename.$k.final.Q >>$directory/outputs.Q"); @@ -141,7 +146,7 @@ close(BEST2); system("cp -rf $directory/log.$best_K $directory/log"); +system("cp -rf $directory/groups.$best_K $directory/groups"); -