%This code is used in "Inversions and genomic differentiation after %secondary contact: when drift contributes to maintenance, not loss, of %differentiation" by Rafajlovic et al. (2021) to analyse the simulated %data, specifically to compute the weighted time of differentiation (see %Rafajlovic et al, (2021) for further details). %Load the raw data. %Here, as an example: Model 2,half-half array, conserved-size version %of the model (vers=1), with inversion polymorphism (invYN=1) vers=1; invYN=1; %also state the specific parameters used s=0.1; mu=0; m=0.1; L=20; rInv=2*10^(-4); span=50; %Here: quoted are the values of N used in Rafajlovic et al. (2021) for %Model 2 NI1=[1000 5000 10000 20000 100000]; AVpTimeWeight=zeros(1,length(NI1)); up1=zeros(1,length(NI1)); down1=zeros(1,length(NI1)); for II=1:length(NI1) N=NI1(II); filename=sprintf('Model2_ArrHalf_vers%g_invYN%g_s%g_N%g_mu%g_m%g_rInv%g_L%g', vers, invYN, s, N, mu, m, rInv, L); filename=strrep(filename,'.','p'); load(filename); pTimeWeight1=zeros(size(frAl1,1),L); for j=1:size(frAl1,1) a1=reshape(frAl1(j,:,:),[],L); a2=reshape(frAl2(j,:,:),[],L); a3=a1-(1-a2); for i=1:L bb1=min(find(a3(:,i)==0)); if size(bb1,1)==0 bb1=size(a1,1); end pTimeWeight1(j,i)=span*sum(a3(1:bb1,i)); end end a=reshape(pTimeWeight1,1,[]); AVpTimeWeight(1,II)=mean(a); up1(1,II)=max(max(a)); down1(1,II)=min(min(a)); end %%%%%%%%%%%%%%%%%% %Plot the Average weighted time of differentiation for the model with %inversion polymorphism: log scale on x axis, linear scale on y-axis figure(1);clf errorbar(log10(NI1),AVpTimeWeight,AVpTimeWeight-down1,up1-AVpTimeWeight,'or','markersize',10,'markerfacecolor','w'); hold on %Here one can add the corresponding results obtained under the collinear model box off %Use desired settings for ticks etc. %%%%%%%%%%%%%%%%% %Plot the Average weighted time of differentiation for the model with %inversion polymorphism: log scale on both axes. figure(2);clf errorbar(log10(NI1),log10(max(1,AVpTimeWeight)),log10(max(1,AVpTimeWeight))-log10(max(1,down1)),log10(max(1,up1))-log10(max(1,AVpTimeWeight)),'or','markersize',10,'markerfacecolor','w'); hold on %Here one can add the corresponding results obtained under the collinear model. box off %Use desired settings for ticks etc.