function BFAC = bayenv2_wrapper(haplo,matrixfile,envfile, popids) % Provide Bayenv values for haplotypes given in P, a matrix file, and environmental files. % % Inputs: % haplo : haplotype matrix % matrixfile : covariance matrix file (Bayenv). % envfile : environmental matrix file (Bayenv)- % % Outputs: % BFAC: Bayes factors. % bayenv2_simuwrapper.m % Written by Andreas Wollstein, Ludwig Maximilians University of Munich. % Copyright (c) 2015 % % This file is part of a set of Matlab scripts published with: % Bozicevic et al. 2015 Molecular Ecology % % This program is free software: you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by % the Free Software Foundation, either version 3 of the License, or % (at your option) any later version. % % This program is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. % % You should have received a copy of the GNU General Public License % along with this program. If not, see . BFAC = []; upop = unique(popids); for i=1:size(haplo,2); fprintf(' bf %d / %d\n',i,size(haplo,2)); g = haplo(:,i); %% count ancestral and derived alleles nanc = []; nder = []; for j=1:length(upop) ix = popids==upop(j); nanc(j) = sum(g(ix)==0); nder(j) = sum(g(ix)==1); end if sum(nder)>0 && sum(nanc)>0 %% write out line 1 % if sum(nder) > sum(nanc) % t = nder; nder = nanc; nanc = t; % end rndn=['snp_r' num2str(rand*1e10,'%010.10f') 'r']; rndn = strrep(rndn,'.',''); f = fopen(['/tmp/' rndn '.txt'],'w'); for j=1:length(upop) fprintf(f,'%d\t',nder(j)); end fprintf(f,'\n'); for j=1:length(upop) fprintf(f,'%d\t',nanc(j)); end fprintf(f,'\n'); fclose(f); %% write out %% apply 10 replicates BF = []; for j=1:10 cmd = ['bayenv2 -i /tmp/' rndn '.txt -m ' matrixfile ' -e ' envfile ' -p 4 -k 10000 -n 6 -t -o /tmp/' rndn '.out']; [s w] = unix(['rm /tmp/' rndn '.out.bf']); [s w] = unix(cmd); try f = fopen(['/tmp/' rndn '.out.bf'],'r'); l = fgetl(f); fclose(f); X = textscan(l,'%s'); bf = str2double(X{1}(2:end)); BF(j,:) = bf'; catch end end BFAC(i,:) = median(BF); [s w] = unix(['rm /tmp/' rndn '.*']); end end