<jupyter_start><jupyter_text>Single amplicon processing<jupyter_code># sample names per species
sps = dict()
for sp in SPP:
    sps[sp] = samples.loc[samples.species==sp, 'ox_code']
display([x.shape for x in sps.values()])
# read amplicon data
ampl = '60'
callset = zarr.open_group(AMPL_HAP_ZARR, 'r')
gt = allel.GenotypeArray(callset[ampl]['calldata/genotype'])
# extract samples - those differ between autosome and X
gs = pd.Series(callset[ampl]['samples'][:])
# nuber of variants
nvar = gt.shape[0]
# extract genotypes per species
gts = dict()
for sp in SPP:
    gts[sp] = gt[:, gs.isin(sps[sp]), :]
display([x.shape for x in gts.values()])
# convert to haplotypes
hts = dict()
for sp in SPP:
    hts[sp] = gts[sp].reshape(gts[sp].shape[0], 2 * gts[sp].shape[1]).transpose()
display([x.shape for x in hts.values()])
# pairwise distances
pds = dict()
# within sp
for sp in SPP:
    pds[sp] = pdist(hts[sp], metric=METRIC)
# between spp for pairwise combinations of species
pds['between_spp'] = np.array([])
for (sp1, sp2) in itertools.combinations(SPP, 2):
    pds['between_spp'] = np.append(pds['between_spp'],
                                   cdist(hts[sp1],hts[sp2], metric=METRIC))
# combine multiple
for (sp, dist) in pds.items():
    plt.hist(dist * nvar, bins=range(8), alpha=0.5, label=sp);
plt.legend();<jupyter_output><empty_output><jupyter_text>Batch processing<jupyter_code>callset = zarr.open_group(AMPL_HAP_ZARR, 'r')
# accumulate pairwise distances into dict
pds = dict()

for ampl in callset:
    sys.stdout.write('\r' + ampl)
    pds[ampl] = dict()
    # read genotypes
    gt = allel.GenotypeArray(callset[ampl]['calldata/genotype'])
    # extract samples - those differ between autosome and X amplicons
    gs = pd.Series(callset[ampl]['samples'][:])
    # nuber of variants
    pds[ampl]['nvar'] = gt.shape[0]
    # extract genotypes per species
    gts = dict()
    for sp in SPP:
        gts[sp] = gt[:, gs.isin(sps[sp]), :]
    # convert to haplotypes
    hts = dict()
    for sp in SPP:
        hts[sp] = gts[sp].reshape(gts[sp].shape[0], 2 * gts[sp].shape[1]).transpose()
    # pairwise distances
    # within sp
    for sp in SPP:
        pds[ampl][sp] = pdist(hts[sp], metric=METRIC)
    # between spp for pairwise combinations of species
    pds[ampl]['between_spp'] = np.array([])
    for (sp1, sp2) in itertools.combinations(SPP, 2):
        pds[ampl]['between_spp'] = np.append(pds[ampl]['between_spp'],
                                             cdist(hts[sp1],hts[sp2], metric=METRIC))
fig, axs = plt.subplots(8,8, figsize=(15,15))
for i, ampl in enumerate(AMPLS):
    ax = axs.flatten()[i]
    n = pds[ampl]['nvar']
    for sp in SPP + ['between_spp']:
        pd_var = pds[ampl][sp] * n
        bins = np.arange(0, int(pd_var.max()) + 1.5) - 0.5
        ax.hist(pd_var, 
                bins=bins, 
                alpha=0.4,
                label=sp);
        ax.set_xticks(bins + 0.5)
        if i == 0:
            ax.legend()
    ax.set_title(ampl)
#     ax.set_xscale('log')
    ax.set_yscale('log')
    ax.set_ylim(0.9)
    
fig.tight_layout()<jupyter_output><empty_output><jupyter_text>Within-species distance thresholdsAbove debug histogram shows that in most cases the distributions of pairwise distances between haplotypes are contiguous, thus taking maximum observed within-species distances will work as a threshold. The only doubtful amplicon is 33, where gambiae is far much more variable than both coluzzii and between-species.<jupyter_code>var_data = panel_mosquito.loc[:, ['start_insert','end_insert']].reset_index()
# sequence length in agam genome, both start and end positions taken into insert by subsetting code 
# (see previous notebook)
var_data['len_insert'] = var_data.end_insert - var_data.start_insert + 1
# find maximum within-species distances across both species
def max_sp_pd_var(ampl):
    max_sp_dist = [max(pds[ampl][sp]) for sp in SPP]
    return max(max_sp_dist) * pds[ampl]['nvar']
var_data['max_withinspecies_var'] = var_data.Primary_ID.apply(max_sp_pd_var)
# normalise by sequence length
var_data['max_withinspecies_dist'] = var_data.max_withinspecies_var / var_data.len_insert
var_data.head()
# write
var_data.to_csv(WSP_VAR_FILE, index=False)
! head {WSP_VAR_FILE}<jupyter_output>Primary_ID,start_insert,end_insert,len_insert,max_withinspecies_var,max_withinspecies_dist
0,10356439,10356587,149,5.0,0.03355704697986577
1,10357476,10357629,154,2.0,0.012987012987012988
2,12084517,12084667,151,6.0,0.039735099337748346
3,13862890,13863031,142,4.0,0.028169014084507043
4,18693503,18693693,191,4.0,0.020942408376963352
5,19065193,19065344,152,7.0,0.046052631578947366
6,20840408,20840599,192,5.0,0.026041666666666668
7,26760718,26760865,148,7.0,0.0472972972972973
8,27674665,27674814,150,6.0,0.04<jupyter_text>Compare to previous estimates<jupyter_code># estimate based on extended sample set
seq_wsp_file = '../../../data/phylo_ampl_dada2/comb1_5/3_thresholds.tsv'
seq_wsp = pd.read_csv(seq_wsp_file, sep='\t')
seq_wsp.head()
# convert to dict
seq_t = seq_wsp.set_index('target')['dist_threshold'].to_dict()
var_data['seq_threshold'] = var_data.Primary_ID.astype(int).replace(seq_t)
fig, ax = plt.subplots(1,1, figsize=(10,4))
ax.scatter(var_data.Primary_ID.astype(int), var_data.seq_threshold, c='b', alpha=0.5, label='sequencing')
ax.scatter(var_data.Primary_ID.astype(int), var_data.max_withinspecies_dist, c='r', alpha=0.5, label='Ag1000g')
ticks = range(62)
ax.set_xticks(ticks)
plt.xticks(rotation=45)
for t in ticks:
    ax.axvline(t - 0.5, c='grey')
ax.legend();
fig, ax = plt.subplots(1,1, figsize=(4,4))
ax.scatter(var_data.max_withinspecies_dist, 
           var_data.seq_threshold)
ax.set_xlabel('Ag1000g max dist')
ax.set_ylabel('Sequencing data threshold')
ax.plot([0, 0.1], [0, 0.1], c='grey', ls='--');<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14562275141339728, 'fraction_numerical': 0.07383930101079321, 'mean_word_length': 4.4868421052631575, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '146898', 'license': 'MIT', 'n_tokens_mistral': 2542, 'n_tokens_neox': 2243, 'n_words': 442, 'path': 'work/6_ag1k_extraction/1_ag1k_within_species_distances.ipynb', 'repo_name': 'malariagen/ano-spp', 'size': '5837'}