%%%%%%%%%%%%%%% %Ravinet et al. (2017) "Interpreting the genomic landscape of speciation: %a road map for finding barriers to gene flow", JEB. %This code simulates the secondary contact after a hard sweep in isolation, %explained in detail in Ravinet et al. (2017). %This code is for the process after the initial "neutral-evolution phase" %Note: the results obtained during the first 500 generations in this code are not %shown in Fig. 3C. RandStream.setDefaultStream (RandStream('mt19937ar','seed',sum(100*clock))); %the input data are obtained from the simulation called: %"Ravinet_et_al_Neutral_Evolution_Phase.m" load('InConds.mat') %choose the initial conditions from the ones generated in a separate %simulation. the following is to go over all initial conditions generated, %but this can be varied condi=1:length(pop1Neutr(:,1,1)); %effective population size in either population N=500; %selection coefficient at the target of selection si=0.2; %mutation model: no mutations at the selected %locus and symmetric two-allele model for neutral loci %mutation rate at the target of selection mu=0; %mutation rate at neutral loci mu_neut=4*10^(-5); %recombination rate between a pair of nearby loci: denote by r1 and r2; but %also assume that the first locus is fully linked to the selected locus %(r=0) r1=10^(-3); r0=r1; %for some loci allow for a larger recombination rate r2=10*r1 r2=10*r0; %number of selected loci L=1; %number of neutral loci (denoted by Ln) linked to the selected one at %increasing recombination distances Ln=61; %of these Ln=61 loci, the recombination rate between a pair of adjacent loci %numbered 51 to 61 is r2 dLn=10; %add also LnU neutral loci unlinked to the selected one, to see the %background pattern LnU=10; %total number of loci in the genetic region analysed K=L+Ln+LnU; %vector of recombination rates between pairs of nearby loci (denote by ri); %here assumed that the total number of loci (K) is larger than 1 ri=r0*ones(1,K-1); %the first locus is tightly linked to the target of selection ri(1,1)=0; %remember that the recombination rate is r2 for a subset of dLn loci ri(1,Ln-(dLn-1):Ln)=r2; %the last LnU loci are fully unlinked ri(1,K-1-(LnU-1):end)=1/2; %secondary contact: no migration during the first N generations after the %initial neutral phase, but now without migration TimeNoMig=N; %migration rate per generation, individual, population mi=[0 0.004]; n=50; neutr_time=0; %the duration of the "selection phase": sel_time sel_time=10^4; %total duration: repeat repeat=neutr_time+sel_time;%*10^4; %eq_time=0; eq_time=0; stat_time=sel_time-eq_time; add_mutALL=zeros(length(condi),1); trialALL=zeros(length(condi),1); timesweep=zeros(length(condi),1); %measure relevant variables every 5th generation (denote this as "span") span=5; CensusPoints=span:span:sel_time; FstALL=zeros(length(condi)*length(CensusPoints),K); ksi1ALL=zeros(length(condi)*length(CensusPoints),K); ksi2ALL=zeros(length(condi)*length(CensusPoints),K); alphaALL=zeros(length(condi)*length(CensusPoints),K); %statistics for samples taken %sample size: n (the same number taken from each population) FstSampleALL=zeros(length(condi)*length(CensusPoints),K); ksi1SampleALL=zeros(length(condi)*length(CensusPoints),K); ksi2SampleALL=zeros(length(condi)*length(CensusPoints),K); alphaSampleALL=zeros(length(condi)*length(CensusPoints),K); dxyALL=zeros(length(condi)*length(CensusPoints),K); dxySampleALL=zeros(length(condi)*length(CensusPoints),K); fr1ALL=zeros(length(condi)*length(CensusPoints),1); fr2ALL=zeros(length(condi)*length(CensusPoints),1); for J=1:length(condi) cond=condi(J); %initialisation at the selected locus. note: no varioation. selection %initially has nothing to act on, so the evolution is neutral pop1=zeros(2*L,N); pop2=zeros(2*L,N); %genetics at neutral loci: taken from the initial conditions generated by a %separate simulations (called "Neutral_evolution_phase.m") pop_neutr1=reshape(pop1Neutr(cond,:,:),2*(K-L),[]); pop_neutr2=reshape(pop2Neutr(cond,:,:),2*(K-L),[]); %statistics for the whole population Fst=zeros(length(CensusPoints),K); ksi1=zeros(length(CensusPoints),K); ksi2=zeros(length(CensusPoints),K); alpha=zeros(length(CensusPoints),K); %statistics for samples taken %sample size: n (the same number taken from each population) n=50; FstSample=zeros(length(CensusPoints),K); ksi1Sample=zeros(length(CensusPoints),K); ksi2Sample=zeros(length(CensusPoints),K); alphaSample=zeros(length(CensusPoints),K); count=0; countNoMig=0; dxy=zeros(stat_time/span,K); dxySample=zeros(stat_time/span,K); fr1=zeros(stat_time/span,1); fr2=zeros(stat_time/span,1); success=0; if neutr_time==0 s=si; end I=0; trial=1; add_mut=0; success_sweep=0; while I1 b=[0 r0 1]; a=rand(K-LnU-dLn-1,2*N); [~,c]=histc(a,b); c(1,:)=2;%first neutral locus completely linked to the selected one b=[0 r2 1]; a=rand(dLn,2*N); [~,c1]=histc(a,b); c=[c;c1]; end b1=[0 1/2 1];%for unlinked loci a=rand(LnU,2*N); [~,c2]=histc(a,b1); a1=pp1([K-LnU+1:K K+K-LnU+1:2*K],:); for i=1:LnU x=i:length(a1(:,1)):2*N*length(a1(:,1)); x=x+(c2(i,:)-1)*LnU; a1(i,:)=a1(x); end pp1([K-LnU+1:K K+K-LnU+1:2*K],:)=a1; %whenever c==1, recombination occurs if Ln>1 a=find(c==1); for i=1:length(a) %find matrix position: individual number is ind=ceil(a(i)/length(c(:,1))) %locus number is a(i)-(ind-1)*(K-1)+1 ind=ceil(a(i)/length(c(:,1))); locN=a(i)-(ind-1)*length(c(:,1))+1; da=pp1(locN:K-LnU,ind); pp1(locN:K-LnU,ind)=pp1(K+locN:end-LnU,ind); pp1(K+locN:end-LnU,ind)=da; end end prog_all_new1=pp1(1:K,:); prog1_new=pp1(1,:); %second population a1=cumsum([0 fit2/sum(fit2)]); a=rand(1,2*N); [~,b]=histc(a,a1); prog_all2=pop_all_new2(:,b);%these are the parents %decide for the first locus whether it comes from the first or the %second chromosome pp2=zeros(2*K,2*N); a=rand(1,2*N); [~,c1]=histc(a,[0 0.5 1]); cc=find(c1==1); pp2(:,cc)=prog_all2(:,cc); cc=find(c1==2); pp2(:,cc)=prog_all2([K+1:end 1:K],cc); if Ln>1 b=[0 r0 1]; a=rand(K-LnU-dLn-1,2*N); [~,c]=histc(a,b); c(1,:)=2;%first neutral locus completely linked to the selected one b=[0 r2 1]; a=rand(dLn,2*N); [~,c1]=histc(a,b); c=[c;c1]; end b1=[0 1/2 1];%for unlinked loci a=rand(LnU,2*N); [~,c2]=histc(a,b1); a1=pp2([K-LnU+1:K K+K-LnU+1:2*K],:); for i=1:LnU x=i:length(a1(:,1)):2*N*length(a1(:,1)); x=x+(c2(i,:)-1)*LnU; a1(i,:)=a1(x); end pp2([K-LnU+1:K K+K-LnU+1:2*K],:)=a1; %whenever c==1, recombination occurs if Ln>1 a=find(c==1); for i=1:length(a) %find matrix position: individual number is ind=ceil(a(i)/length(c(:,1))) %locus number is a(i)-(ind-1)*(K-1)+1 ind=ceil(a(i)/length(c(:,1))); locN=a(i)-(ind-1)*length(c(:,1))+1; da=pp2(locN:K-LnU,ind); pp2(locN:K-LnU,ind)=pp2(K+locN:end-LnU,ind); pp2(K+locN:end-LnU,ind)=da; end end prog_all_new2=pp2(1:K,:); prog2_new=pp2(1,:); prog_neutr1=prog_all_new1(L+1:end,:); prog_neutr2=prog_all_new2(L+1:end,:); %introduce mutations; first selected loci if mu>0 a=rand(1,2*L*N); b=[0 mu 1+mu]; [~,c]=histc(a,b); c=find(c==1); prog1_new(c)=1-prog1_new(c); a=rand(1,2*L*N); b=[0 mu 1+mu]; [~,c]=histc(a,b); c=find(c==1); prog2_new(c)=1-prog2_new(c); end %mutations at neutral loci if mu_neut>0 a=rand(1,2*N*(K-L)); b=[0 mu_neut 1+mu_neut]; [~,c]=histc(a,b); c=find(c==1); prog_neutr1(c)=1-prog_neutr1(c); a=rand(1,2*N*(K-L)); b=[0 mu_neut 1+mu_neut]; [~,c]=histc(a,b); c=find(c==1); prog_neutr2(c)=1-prog_neutr2(c); end %lastly, the new population; selected part pop1=[prog1_new(:,1:N);prog1_new(:,N+1:2*N)]; pop2=[prog2_new(:,1:N);prog2_new(:,N+1:2*N)]; %neutral part pop_neutr1=[prog_neutr1(:,1:N);prog_neutr1(:,N+1:2*N)]; pop_neutr2=[prog_neutr2(:,1:N);prog_neutr2(:,N+1:2*N)]; pop1_all=[prog1_new;prog_neutr1]; pop2_all=[prog2_new;prog_neutr2]; %summary statistics if mod(I,span)==0 count=count+1; %for a sample b=randperm(N); popS1=[pop1_all(:,b(1:n)) pop1_all(:,N+b(1:n))]; b=randperm(N); popS2=[pop2_all(:,b(1:n)) pop2_all(:,N+b(1:n))]; a=popS1'; a1=mean(a); ksi1Sample(count,:)=a1.^2+(1-a1).^2; a2=popS2'; a1=mean(a2); ksi2Sample(count,:)=a1.^2+(1-a1).^2; a=[a;a2]; a1=mean(a); alphaSample(count,:)=a1.^2+(1-a1).^2; a=pop1_all'; a1=mean(a); fr1(count,1)=1-a1(1); ksi1(count,:)=a1.^2+(1-a1).^2; a2=pop2_all'; a1=mean(a2); fr2(count,1)=a1(1); ksi2(count,:)=a1.^2+(1-a1).^2; a=[a;a2]; a1=mean(a); alpha(count,:)=a1.^2+(1-a1).^2; Fst(count,:)=((ksi1(count,:)+ksi2(count,:))/2-alpha(count,:))./(1-alpha(count,:)); FstSample(count,:)=((ksi1Sample(count,:)+ksi2Sample(count,:))/2-alphaSample(count,:))./(1-alphaSample(count,:)); dxy(count,:)=1-2*(alpha(count,:)-ksi1(count,:)/4-ksi2(count,:)/4); dxySample(count,:)=1-2*(alphaSample(count,:)-ksi1Sample(count,:)/4-ksi2Sample(count,:)/4); end %this is where variation is introduced at the selected locus! before %this point, selection is "turned on" but there is no variation for %selection to act on, and so the evolution before "TimeNoMig" is fully %neutral. %with the following part of the code, introduce only a single "mutant" %allele (denoted "1"); all others are "0". Allele "1" is introduced in %the second population where it is beneficial. if I==TimeNoMig success_sweep=1; a=1+floor(rand(1,1).*2*L*N); pop2(a)=1; end if success_sweep==1 if I>TimeNoMig a=length(find(prog2_new==1)); if a==0 %this means that the mutant is lost; then start everythign %again I=0; count=0; pop1=zeros(2*L,N); pop2=zeros(2*L,N); pop_neutr1=reshape(pop1Neutr(cond,:,:),2*(K-L),[]); pop_neutr2=reshape(pop2Neutr(cond,:,:),2*(K-L),[]); trial=trial+1; success=0; success_sweep=0; else if length(find(prog2_new==1))==2*N if length(find(prog1_new==0))==2*N %this means that the locally beneficial allele is %fixated in its corresponding population, and the %sweep is successful. after this, migration is %reintroduced between the two populations success_sweep=2; timesweep(J,1)=I; success=1; end end end end end end fr1ALL((J-1)*length(Fst(:,1))+1:J*length(Fst(:,1)),:)=fr1; fr2ALL((J-1)*length(Fst(:,1))+1:J*length(Fst(:,1)),:)=fr2; add_mutALL(J)=add_mut; trialALL(J)=trial; FstALL((J-1)*length(Fst(:,1))+1:J*length(Fst(:,1)),:)=Fst; FstSampleALL((J-1)*length(Fst(:,1))+1:J*length(Fst(:,1)),:)=FstSample; ksi1ALL((J-1)*length(Fst(:,1))+1:J*length(Fst(:,1)),:)=ksi1; ksi2ALL((J-1)*length(Fst(:,1))+1:J*length(Fst(:,1)),:)=ksi2; alphaALL((J-1)*length(Fst(:,1))+1:J*length(Fst(:,1)),:)=alpha; alphaSampleALL((J-1)*length(Fst(:,1))+1:J*length(Fst(:,1)),:)=alphaSample; ksi1SampleALL((J-1)*length(Fst(:,1))+1:J*length(Fst(:,1)),:)=ksi1Sample; ksi2SampleALL((J-1)*length(Fst(:,1))+1:J*length(Fst(:,1)),:)=ksi2Sample; dxyALL((J-1)*length(Fst(:,1))+1:J*length(Fst(:,1)),:)=dxy; dxySampleALL((J-1)*length(Fst(:,1))+1:J*length(Fst(:,1)),:)=dxySample; end %finally, the data are to be saved in a convenient form and with a %convenient name, for example, as follows: save('ALL_NoMutSel_SecHardSweepNew_N500_s0p2_m0p004_Ln61_dLn10_r0p001_r0p01_LnU10_mu4p10(-5)_00','timesweep','success_sweep','trialALL','add_mutALL','condi','fr1ALL','fr2ALL','eq_time','stat_time','dxyALL','si','dxySampleALL','mu_neut','neutr_time','sel_time','FstALL','FstSampleALL','n','alphaALL','alphaSampleALL','ksi1ALL','ksi1SampleALL','ksi2ALL','ksi2SampleALL','r1','r0','Ln','L','LnU','repeat','ri','m','N','mu','span','r2','cond') %%%%%%%%%%%%%%%%%%%%%%%% %note: to obtain good statistics, it is advisable to repeat this code %several times (with a random seed for the random-number generator). %In Ravinet et al. (2017), this code (or its alternative with a subset of %initial conditions) was repeated so that the total number of individual %realisations of the process is 5000.