clf 
clear

N=400;

w0=0.10;
randSeed=17;

startSite=401;
middleSite=601;
endSite=800;
startSwap=801;
endSwap=1000;


data= Sirens(randSeed,w0);

dna0= data(1:4,startSite:endSite);
dna_swap= data(1:4,startSwap:endSwap);

replaced_sites_1=0;
% Replacing 10 percent of the sites
dna1=dna0;
for i=middleSite:endSite
	if rand<0.1
		replaced_sites_1=replaced_sites_1+1;
		dna1(:,i-startSite+1)= dna_swap(:,i-middleSite+1);
	end
end
replaced_sites_1

replaced_sites_2=0;
% Replacing 20 percent of the sites
dna2=dna0;
for i=middleSite:endSite
	if rand<0.2
		replaced_sites_2=replaced_sites_2+1;
		dna2(:,i-startSite+1)= dna_swap(:,i-middleSite+1);
	end
end
replaced_sites_2


% -------------------------------------------------------------
%  	Write out alignment to FastA and Phylip file: dna0
% -------------------------------------------------------------
fid =fopen('dna0.fastA','w');
K=4;
for k=1:K
        strain=strcat('>strain',num2str(k));
        fprintf(fid,'%8s\n',strain);
        for n=1:N
                if dna0(k,n)==1
                        fprintf(fid,'%1s','A');
                elseif dna0(k,n)==2
                        fprintf(fid,'%1s','C');
                elseif dna0(k,n)==3
                        fprintf(fid,'%1s','G');
                elseif dna0(k,n)==4
                        fprintf(fid,'%1s','T');
                else
                        error('Unknown nucleotide');
                end
                if mod(n,50)==0 | n==N
                        fprintf(fid,'\n');
                end
        end
end  
fclose(fid);

!readseq -f=Phylip3.2 dna0.fastA -out=stercus -a 
!sed 's/ YF//g' stercus > dna0.phy 
!/bin/rm stercus


% -------------------------------------------------------------
%  	Write out alignment to FastA and Phylip file: dna1
% -------------------------------------------------------------
fid =fopen('dna1.fastA','w');
K=4;
for k=1:K
        strain=strcat('>strain',num2str(k));
        fprintf(fid,'%8s\n',strain);
        for n=1:N
                if dna1(k,n)==1
                        fprintf(fid,'%1s','A');
                elseif dna1(k,n)==2
                        fprintf(fid,'%1s','C');
                elseif dna1(k,n)==3
                        fprintf(fid,'%1s','G');
                elseif dna1(k,n)==4
                        fprintf(fid,'%1s','T');
                else
                        error('Unknown nucleotide');
                end
                if mod(n,50)==0 | n==N
                        fprintf(fid,'\n');
                end
        end
end  
fclose(fid);

!readseq -f=Phylip3.2 dna1.fastA -out=stercus -a 
!sed 's/ YF//g' stercus > dna1.phy 
!/bin/rm stercus


% -------------------------------------------------------------
%  	Write out alignment to FastA and Phylip file: dna2
% -------------------------------------------------------------
fid =fopen('dna2.fastA','w');
K=4;
for k=1:K
        strain=strcat('>strain',num2str(k));
        fprintf(fid,'%8s\n',strain);
        for n=1:N
                if dna2(k,n)==1
                        fprintf(fid,'%1s','A');
                elseif dna2(k,n)==2
                        fprintf(fid,'%1s','C');
                elseif dna2(k,n)==3
                        fprintf(fid,'%1s','G');
                elseif dna2(k,n)==4
                        fprintf(fid,'%1s','T');
                else
                        error('Unknown nucleotide');
                end
                if mod(n,50)==0 | n==N
                        fprintf(fid,'\n');
                end
        end
end  
fclose(fid);

!readseq -f=Phylip3.2 dna2.fastA -out=stercus -a 
!sed 's/ YF//g' stercus > dna2.phy 
!/bin/rm stercus