load results.mat%output file from hierBAPS analysis of the original alignment n_loci=length(c.snpPosition); sing_locus_freqs=zeros(1,n_loci); for i=1:n_loci     freq_table=tabulate(c.snpData(:,i));     here=freq_table(:,1)<16;     if length(find(freq_table(here,2)))==2 & length(find(freq_table(here,2)==1))==1         sing_locus_freqs(1,i)=1;     end end singleton_loci=c.snpPosition(find(sing_locus_freqs));      for i=1:620 [Header,Sequence] = fastaread('snp.aln','Blockread',i); Sequence(singleton_loci)='-'; %if length(find(strcmp(Header,outliers)))==0 fastawrite('snpfinal_nosingletonsnps.fas',Header,Sequence); %end

end