Skip to content

Commit

Permalink
Add files via upload
Browse files Browse the repository at this point in the history
  • Loading branch information
@jorgenripa
jorgenripa committed Jan 23, 2020
1 parent da06c02 commit b9081de
Show file tree
Hide file tree
Showing 3 changed files with 319 additions and 81 deletions.
49 changes: 28 additions & 21 deletions TREES_Matlab/F_ST.m
View file
Expand Up @@ -26,38 +26,45 @@
a_tot = 0;
b_tot = 0;
c_tot = 0;

Group1_genes = [sample.G1(:,group1) sample.G2(:,group1)];
Group2_genes = [sample.G1(:,group2) sample.G2(:,group2)];
for li = 1:length(loci)
locus = loci(li);
% group genes:
G1 = [sample.G1(locus,group1); sample.G2(locus,group1)];
G2 = [sample.G1(locus,group2); sample.G2(locus,group2)];
G1 = reshape(Group1_genes(locus,:),2,n1); % [sample.G1(locus,group1); sample.G2(locus,group1)];
G2 = reshape(Group2_genes(locus,:),2,n2); % [sample.G1(locus,group2); sample.G2(locus,group2)];
a_locus = 0;
b_locus = 0;
c_locus = 0;

% Find all alleles:
all_alleles = unique([G1(:);G2(:)])';
if isa(G1,'int8')
all_alleles = [-1 1];
else
all_alleles = unique([G1(:);G2(:)])';
end
for allele = all_alleles
p1 = sum(G1(:)==allele)/2/n1;
p2 = sum(G2(:)==allele)/2/n2;
p_bar = (n1*p1+n2*p2)/N;
s2 = (n1*(p1-p_bar)^2 + n2*(p2-p_bar)^2)/(r-1)/n_bar;
% homozygotes:
homo1 = sum(G1(1,:)==allele & G1(2,:)==allele);
% proportion heterozygotes:
hetero1 = (p1*2*n1 - 2*homo1) / n1;
% homozygotes:
homo2 = sum(G2(1,:)==allele & G2(2,:)==allele);
% proportion heterozygotes:
hetero2 = (p2*2*n2 - 2*homo2) / n2;
h_bar = (n1*hetero1 + n2*hetero2)/N; % average heterozygosity
a = n_bar/n_c*( s2 - 1/(n_bar-1)*( p_bar*(1-p_bar) - (r-1)/r*s2 - h_bar/4));
a_locus = a_locus + a;
b = n_bar/(n_bar-1)*( p_bar*(1-p_bar) - (r-1)/r*s2 - (2*n_bar-1)/4/n_bar * h_bar );
b_locus = b_locus + b;
c = h_bar/2;
c_locus = c_locus + c;
if p1>0 || p2>0
p_bar = (n1*p1+n2*p2)/N;
s2 = (n1*(p1-p_bar)^2 + n2*(p2-p_bar)^2)/(r-1)/n_bar;
% homozygotes:
homo1 = sum(G1(1,:)==allele & G1(2,:)==allele);
% proportion heterozygotes:
hetero1 = (p1*2*n1 - 2*homo1) / n1;
% homozygotes:
homo2 = sum(G2(1,:)==allele & G2(2,:)==allele);
% proportion heterozygotes:
hetero2 = (p2*2*n2 - 2*homo2) / n2;
h_bar = (n1*hetero1 + n2*hetero2)/N; % average heterozygosity
a = n_bar/n_c*( s2 - 1/(n_bar-1)*( p_bar*(1-p_bar) - (r-1)/r*s2 - h_bar/4));
a_locus = a_locus + a;
b = n_bar/(n_bar-1)*( p_bar*(1-p_bar) - (r-1)/r*s2 - (2*n_bar-1)/4/n_bar * h_bar );
b_locus = b_locus + b;
c = h_bar/2;
c_locus = c_locus + c;
end
end
F_ST_i(li) = a_locus/(a_locus + b_locus + c_locus);
a_tot = a_tot + a_locus;
Expand Down
58 changes: 29 additions & 29 deletions TREES_Matlab/plot_sim.m
View file
Expand Up @@ -24,26 +24,26 @@ function plot_sim(sim)
bins = cell(1,nT);
if strcmpi(sim.Genetics.model,'diallelic')
maxx = zeros(1,nT);
for t = 1:nT
x = [sim.samples.(Trait_names{t})];
maxx(t) = max(abs(x(:)));
for tr = 1:nT
x = [sim.samples.(Trait_names{tr})];
maxx(tr) = max(abs(x(:)));
end

for i=1:length(sim.samples)
s = sim.samples(i);
for t = 1:nT
x = s.(Trait_names{t});
if maxx(t) > 0
x = round(x/maxx(t),7)*maxx(t);
for si=1:length(sim.samples)
s = sim.samples(si);
for tr = 1:nT
x = s.(Trait_names{tr});
if maxx(tr) > 0
x = round(x/maxx(tr),3)*maxx(tr);
end
bins{t} = unique([bins{t} x(:)']);
bins{tr} = unique([bins{tr} x(:)']);
end
end
for t=1:nT
if length(bins{t})>1
bins{t} = min(bins{t}):min(diff(bins{t})):max(bins{t});
for tr=1:nT
if length(bins{tr})>1
bins{tr} = min(bins{tr}):min(diff(bins{tr})):max(bins{tr});
else
bins{t} = bins{t}+[-1 0 1];
bins{tr} = bins{tr}+[-1 0 1];
end
end
else
Expand All @@ -52,7 +52,7 @@ function plot_sim(sim)
for i=1:length(sim.samples)
s = sim.samples(i);
for t = 1:length(Trait_names)
x = s.(Trait_names{t});
x = s.traits.(Trait_names{t});
m = max(x(:));
Smaxv(t) = max(Smaxv(t),m);
m = min(x(:));
Expand All @@ -69,17 +69,17 @@ function plot_sim(sim)
end
% plot all traits:
figure(fig), clf
tt = [sim.samples(:).gen];
nn = [sim.samples(:).sample_size];
tt = [sim.samples.gen];
nn = [sim.samples.size];
%sp = 1;
rows = max(Trait_dims);
for tr=1:nT
for d=1:Trait_dims(tr)
positions = zeros(length(bins{tr}),sim.sample_count);
for si = 1:sim.sample_count
s = sim.samples(si);
if s.sample_size > 0
xx = sim.samples(si).(Trait_names{tr});
if s.size > 0
xx = s.(Trait_names{tr});
xx = xx(d,:);
if 0 %isfield(s,'pos')
pp = s.pos;
Expand Down Expand Up @@ -111,16 +111,16 @@ function plot_sim(sim)

figtitle(sim.name)

figure(fig+1), clf
%[mean_Achoose,mean_Bchoose,std_colour,nnsp] = plot_RI_trends(sim,false,0.1);
subplot(2,1,1)
if tt(end) < sim.t_max % extinction?
nn = [nn 0];
tt = [tt tt(end)+sim.sample_interval];
end
plot(tt,nn,'b-')
xlabel('generations')
ylabel('total population size (blue)')
% figure(fig+1), clf
% %[mean_Achoose,mean_Bchoose,std_colour,nnsp] = plot_RI_trends(sim,false,0.1);
% subplot(2,1,1)
% if tt(end) < sim.t_max % extinction?
% nn = [nn 0];
% tt = [tt tt(end)+sim.sample_interval];
% end
% plot(tt,nn,'b-')
% xlabel('generations')
% ylabel('total population size (blue)')
% raxis
% plot(tt,nnsp,'r-');
% ylabel('#species (red)')
Expand Down

0 comments on commit b9081de

Please sign in to comment.