# HG changeset patch # User dereeper # Date 1463496694 14400 # Node ID 9a1729b8940568310aeba25c860975bd678013d0 Uploaded diff -r 000000000000 -r 9a1729b89405 Snmf.pl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Snmf.pl Tue May 17 10:51:34 2016 -0400 @@ -0,0 +1,129 @@ +#!/usr/bin/perl + +use strict; +use Switch; +use Getopt::Long; +use Bio::SeqIO; + +my $usage = qq~Usage:$0 [] +where are: + -i, --input + -o, --output + -k, --kmin + -m, --maxK + -d, --directory + -t, --threshold +~; +$usage .= "\n"; + +my ($input,$output,$kmin,$kmax,$directory,$threshold); + + +GetOptions( + "input=s" => \$input, + "output=s" => \$output, + "kmin=s" => \$kmin, + "maxK=s" => \$kmax, + "directory=s" => \$directory, + "threshold=s" => \$threshold +); + + +die $usage + if ( !$input || !$output || !$kmin || !$kmax || !$directory || !$threshold); + + +my $PLINK_EXE = "plink"; + +system("$PLINK_EXE --vcf $input --allow-extra-chr --recode-vcf --out $directory/input >>$directory/plink.log 2>&1"); + +system("vcf2geno $directory/input.vcf $directory/polymorphisms.geno >>$directory/vcf2geno.log 2>&1"); + + +my $ind_cmd = `grep '#CHROM' $input`; +chomp($ind_cmd); +my @individuals = split(/\t/,$ind_cmd);shift @individuals;shift @individuals;shift @individuals;shift @individuals;shift @individuals;shift @individuals;shift @individuals;shift @individuals; + +################################### +# launch admixture for different K +################################### +my %errors; +for (my $k = $kmin; $k <= $kmax; $k++) +{ + system("sNMF -x $directory/polymorphisms.geno -K $k -c >>$directory/log.$k 2>&1"); + + open(my $O3,">$directory/out.$k.group"); + open(my $O2,">$directory/out.$k.final.Q"); + + my $ent; + open(my $LOG,"$directory/log.$k"); + while(<$LOG>){ + if (/Cross-Entropy \(masked data\).*(\d+\.\d+)$/){ + $ent = $1; + $errors{$ent} = $k; + } + } + close($LOG); + + open(E,">>$directory/entropy"); + print E "K=$k $ent\n"; + close(E); + + print $O2 "Indiv"; + print $O3 "Indiv;Group\n"; + for (my $j = 0; $j <$k; $j++){ + print $O2 " Q$j"; + } + print $O2 "\n"; + + open(my $O,"$directory/polymorphisms.$k.Q"); + my %hash_groupes; + my %hash_indv; + my %group_of_ind; + my $i = 0; + while (<$O>){ + $i++; + my $line = $_; + $line =~s/\n//g; + $line =~s/\r//g; + my @infos = split(/\s+/,$line); + my $group = "admix"; + my $ind = $individuals[$i]; + for (my $j = 0; $j <$k; $j++){ + my $val = $infos[$j]; + if ($val > 0.5){$group = "Q$j";} + } + if ($ind){ + $hash_indv{$ind} = join(" ",@infos); + $hash_groupes{$group}{"ind"} .= ",".$ind; + $group_of_ind{$ind} = $group; + } + } + close($O); + + foreach my $group(sort keys(%hash_groupes)){ + my @inds = split(",",$hash_groupes{$group}{"ind"}); + foreach my $ind(@inds){ + if ($ind =~/\w+/){ + print $O3 "$ind;$group\n"; + print $O2 $ind." ".$hash_indv{$ind}. "\n"; + } + } + } + + system("cat $directory/log.$k >>$directory/logs"); + system("echo '\n\n====================================\n\n' >>$directory/logs"); + system("cat $directory/out.$k.final.Q >>$directory/outputs.Q"); + system("echo '\n\n====================================\n\n' >>$directory/outputs.Q"); +} + +my @sorted_errors = sort {$a<=>$b} keys(%errors); +my $best_K = $errors{@sorted_errors[0]}; + + +system("cp -rf $directory/out.$best_K.final.Q $directory/output"); + +system("cp -rf $directory/log.$best_K $directory/log"); +system("cp -rf $directory/out.$best_K.group $directory/groups"); + + diff -r 000000000000 -r 9a1729b89405 snmf.sh --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/snmf.sh Tue May 17 10:51:34 2016 -0400 @@ -0,0 +1,24 @@ +#!/bin/bash +vcf=$1 +outputs=$2 +logs=$3 +best_k_output=$4 +best_k_logfile=$5 +kmin=$6 +kmax=$7 +groups=$8 +threshold_group=$9 + +directory=`dirname $0` +mkdir tmpdir$$ + + +perl $directory/Snmf.pl -i $vcf -o $outputs -k $kmin -m $kmax -d tmpdir$$ -t $threshold_group + +mv tmpdir$$/output $best_k_output +mv tmpdir$$/log $best_k_logfile +mv tmpdir$$/outputs.Q $outputs +mv tmpdir$$/logs $logs +mv tmpdir$$/groups $groups + + diff -r 000000000000 -r 9a1729b89405 snmf.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/snmf.xml Tue May 17 10:51:34 2016 -0400 @@ -0,0 +1,64 @@ + + a population structure from large SNP genotype datasets + + snmf + plink + + ./snmf.sh $vcf $outputs $logs $best_k_output $best_k_logfile $kmin $kmax $best_k_groups $threshold_group + + + + + + + + + + + + + + + + + + +.. class:: infomark + +**Program encapsulated in Galaxy by Southgreen** + +.. class:: infomark + +**sNMF version 1.2** + +----- + +=========== + Overview: +=========== + +Fast and efficient program for estimating individual admixture coefficients based on sparse non-negative matrix factorization and population genetics. + +----- + +For further informations, please visite the sNMF_ website. + + +.. _sNMF: http://membres-timc.imag.fr/Olivier.Francois/snmf/index.htm + + +10.1534/genetics.113.160572 + + diff -r 000000000000 -r 9a1729b89405 tool_dependencies.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tool_dependencies.xml Tue May 17 10:51:34 2016 -0400 @@ -0,0 +1,9 @@ + + + + + + + + +