Results for Appendix 3

This notebook explores the expected time to first coalescence of an ancient lineage into a modern panel. We explore two particular approximations for this to make claims on the scales for haplotype copying.

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm
import os
import sys
sys.path.append('../../src/')
from first_coal_anc_samples import *
from plot_utils import *

%matplotlib inline
In [2]:
plt.rcParams['font.sans-serif'] = "Arial"
plt.rcParams['figure.facecolor'] = "w"
plt.rcParams['figure.autolayout'] = True
plt.rcParams['pdf.fonttype'] = 3

# Making the relevant figure directories that we want
main_figdir = '../../plots/supp_figs/appendix3/'
os.makedirs(main_figdir, exist_ok=True)
In [3]:
nt = 20
tas = 10**np.linspace(-5,-0.5, nt)

# Reference Panel Size and colors for points...
Ks = [1000, 5000, 10000]
colors = ['blue','orange', 'green']

f, ax = plt.subplots(1,1,figsize=(4,3))
for i in range(len(Ks)):
    e_t = np.zeros(nt)
    for j in tqdm(range(nt)):
        t = tas[j]
        # calculating the time to the first coalescent using griffiths approximation
        e_t[j] = time_first_coal_griffiths(Ks[i],t)

    ax.loglog(tas, e_t, lw=2, color=colors[i], label=r'$K = %d$' % Ks[i])

ax.legend(fontsize=10)
ax.set_xlabel(r'$t_a$', fontsize=14);
ax.set_ylabel(r'$E[T^*]$', fontsize=14);
debox(ax);

# Save both the png and pdf versions of the plot 
plt.savefig(main_figdir + 'exp_time_to_coal_anc.pdf', bbox_inches='tight', dpi=300)
plt.savefig(main_figdir + 'exp_time_to_coal_anc.png', bbox_inches='tight', dpi=400)

100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 20/20 [00:00<00:00, 22.92it/s]
100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 20/20 [00:11<00:00,  1.68it/s]
100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 20/20 [00:48<00:00,  2.42s/it]
In [4]:
print(time_first_coal_griffiths(1000,1e-2)/time_first_coal_griffiths(10000,1e-2))
print(time_first_coal_griffiths(1000,1e-4)/time_first_coal_griffiths(10000,1e-4))

1.1726186634543034
6.99580269588812
In [5]:
# Testing asymptotic calculations against exact values from Griffiths 1984 (Table 1)
nt = 100
Ne = 1e4
gens = np.arange(1,10000)
tas = gens/Ne
K = 50
e_nt_asymptotic = np.array([appx_num_lineages_mean_var(K,t)[0] for t in tas])
var_nt_asymptotic = np.array([appx_num_lineages_mean_var(K,t)[1] for t in tas])

fig, ax = plt.subplots(2,1,figsize=(3, 2*2), sharex=True)
ax[0].plot(tas, e_nt_asymptotic, lw=2)
# ax[0].set_xlabel(r't', fontsize=14)
ax[0].scatter(0.01, 40.161, color='red', marker='x')
ax[0].scatter(0.1,14.61, color='red', marker='x')
ax[0].scatter(1.0,2.289, color='red', marker='x', label=r'Exact')
ax[0].set_ylabel(r'$\mathbb{E}[A_n(t)]$',fontsize=14)
ax[0].legend()

ax[1].plot(tas, var_nt_asymptotic, lw=2)
ax[1].scatter(0.01,6.443, color='red', marker='x',label=r'Exact')
ax[1].scatter(0.1, 4.668, color='red', marker='x')
ax[1].scatter(1.0, 0.643, color='red', marker='x')
ax[1].set_ylabel(r'$Var[A_n(t)]$',fontsize=14)
ax[1].set_xlabel(r't',fontsize=14)
ax[1].legend()

debox(ax[0]); debox(ax[1]);
fig.tight_layout();
plt.savefig(main_figdir + 'verify_An.pdf', bbox_inches='tight', dpi=300)