# 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 <args> [<opts>]
+where <args> are:
+    -i, --input         <input VCF>
+    -o, --output        <output>
+    -k, --kmin          <K min>
+    -m, --maxK          <K max>
+    -d, --directory     <temporary directory>
+    -t, --threshold     <threshold admixture proportion for group assignation>
+~;
+$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 @@
+<tool id="snmf" name="sNMF" version="1.2">
+	<description>a population structure from large SNP genotype datasets</description>
+	<requirements>
+		<requirement type="package" version="1.2">snmf</requirement>
+		<requirement type="package" version="1.90">plink</requirement>
+	</requirements>
+	<command interpreter="bash">./snmf.sh $vcf $outputs $logs $best_k_output $best_k_logfile $kmin $kmax $best_k_groups $threshold_group
+    </command>
+	<inputs>
+		<param format="vcf" name="vcf" type="data" label="VCF file" help="VCF file"/>
+		<param type="text" name="kmin" label="K min" value="2"/>
+		<param type="text" name="kmax" label="K max" value="5"/>
+		<param type="text" name="threshold_group" label="Minimum admixture proportion percentage for group assignation" value="50"/>
+	</inputs>
+	<outputs>
+		<data format="txt" name="best_k_output" label="Best K Output"/>
+		<data format="txt" name="best_k_groups" label="Best K Groups"/>
+		<data format="txt" name="best_k_logfile" label="Best K Logfile"/>
+		<data format="txt" name="outputs" label="All Outputs"/>
+		<data format="txt" name="logs" label="All Logs"/>
+	</outputs>
+<!--
+	<tests>
+		<test>
+			<param name="vcf" value="input.vcf" />
+			<param name="kmax" value="3" />
+			<param name="threshold_group" value="60" />
+			<output name="best_k_output" file="output" />
+			<output name="outputs" file="outputs.Q" />
+			<output name="best_k_groups" file="groups" />
+		</test>
+
+	</tests>
+-->
+	<help>
+	
+
+.. 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
+	</help>
+<citations>
+<citation type="doi" >10.1534/genetics.113.160572</citation>
+</citations>
+</tool>
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 @@
+<?xml version="1.0"?>
+<tool_dependency>
+    <package name="snmf" version="1.2">
+        <repository changeset_revision="fd556ed7ce46" name="package_snmf_1_2" owner="dereeper" toolshed="https://toolshed.g2.bx.psu.edu" />
+    </package>
+    <package name="plink" version="1.90">
+        <repository changeset_revision="bae25c2eb2f6" name="package_plink_1_90" owner="dereeper" toolshed="https://toolshed.g2.bx.psu.edu" />
+    </package>
+</tool_dependency>