Mercurial > repos > dereeper > snmf
changeset 0:9a1729b89405 draft
Uploaded
| author | dereeper |
|---|---|
| date | Tue, 17 May 2016 10:51:34 -0400 |
| parents | |
| children | 395c7ab476a6 |
| files | Snmf.pl snmf.sh snmf.xml tool_dependencies.xml |
| diffstat | 4 files changed, 226 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /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"); + +
--- /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 + +
--- /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>
--- /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>