<jupyter_start><jupyter_text>build zarr from merged vcf<jupyter_code># def build_sample_zarr(input, output):
#     allel.vcf_to_zarr(
#         input=input,
#         output=output,
#         compressor=zarr.Blosc(cname='zstd', clevel=5, shuffle=0),
#         overwrite=True,
#         fields=['*'],
#         alt_number=3,
#         chunk_length=2**18,
#         log=sys.stdout,
#     )
input_vcf = '/home/chris/Universal_ref/variation/merged.vcf'
output_path = '../data/merged.zarr'
# build_sample_zarr(input_vcf, output_path)<jupyter_output><empty_output><jupyter_text>check zarr<jupyter_code>callset = zarr.open_group('merged.zarr/', mode='r')
callset
list(callset['calldata'].keys())
callset['calldata/GT'][:].shape
gt = allel.GenotypeArray(callset['calldata/GT'][:])
gt
chrom = callset['variants/CHROM']
chrom[:]<jupyter_output><empty_output><jupyter_text>get meta data from csv<jupyter_code>df_meta = pd.read_csv('1175-VO-KH-STLAURENT-S01-samples.csv', sep=',')
df_meta.columns, len(df_meta)
#keep only sequenced
df_meta_seq = df_meta[df_meta.Status == 'sequenced']
len(df_meta_seq)
all_128_samples = np.asarray(df_meta_seq.ROMA_ID)<jupyter_output><empty_output><jupyter_text>setup data for distance analysis 1. remove the Pf gene positions 2. look as missingness per variant and per sample, we only really want to keep informative positions/samples3. what happens to missing data if we leave it there then calculate distance?<jupyter_code>gt.shape
#1. remove Pf gene
pfbool = chrom[:]!='Pf_M76611'
gt_nopf = gt.compress(pfbool)
gt_nopf.shape
#2. look at per sample missingness
ismiss = gt_nopf.is_missing()
persam_ismiss = np.sum(ismiss, axis=0)
persam_ismiss.shape
#plot
persam_rate = persam_ismiss/260581*100
plt.hist(persam_rate);
#eyeballing plot - keep samples with <40% missingness - loose 25 samples
persam_bool = persam_rate <= 40
np.count_nonzero(persam_bool)
#keep track of the samples we have left
gud_samples = np.asarray(callset['samples']).compress(persam_bool)
len(gud_samples)
#compress genotype array to keep just gud samples
gt_gudsam = gt_nopf.compress(persam_bool, axis=1)
gt_gudsam.shape
#compress meta data to just keep the info for the gud samples
gudsam_bool = np.in1d(all_128_samples, gud_samples)
df_meta_gud = df_meta_seq[gudsam_bool]
df_meta_gud.shape
#3. look at per var missingness
ismissvar = gt_gudsam.is_missing()
pervar_ismiss = np.sum(ismissvar, axis=1)
pervar_ismiss.shape
pervar_rate = pervar_ismiss/96*100
plt.hist(pervar_rate);
#how many var do we have with 0 missingness?
full = pervar_rate == 0
np.count_nonzero(full)<jupyter_output><empty_output><jupyter_text>let's look at the dendrogram using just the 0% missing variants<jupyter_code>#if we just use the 0 missingness sites what happens?
gt_nomiss = gt_gudsam.compress(full)
gt_nomiss
#measure distance and cluster samples
d = allel.stats.pairwise_distance(gt_nomiss.to_n_alt(), metric='cityblock')
z = scipy.cluster.hierarchy.linkage(d, method='average')
#palette
palette = sns.color_palette('hls', n_colors=12, desat=0.8)
sns.palplot(palette);
# #use names and palette to define plotting colours
# colors = {
#     'An. annularis group': palette[0],
#     '': palette[1],
#     'An. hyrcanus group': palette[2],
#     'An. maculatus group': palette[3],
#     'An. kochi group': palette[4],
#     'An. funestus group': palette[5],
#     'An. barbirostris group': palette[6],
#     'An. dirus complex': palette[7]
# }

# labels = {
#     'An. annularis group': 'An. annularis group',
#     '': 'other',
#     'An. hyrcanus group': 'An. hyrcanus group',
#     'An. maculatus group':'An. maculatus group',
#     'An. kochi group': 'An. kochi group',
#     'An. funestus group': 'An. funestus group',
#     'An. barbirostris group': 'An. barbirostris group',
#     'An. dirus complex':'An. dirus complex'
# }
#use names and palette to define plotting colours

colors = {
    'An. annularis group': palette[0],
    'An. hyrcanus group': palette[3],
    'An. maculatus group': palette[4],
    'An. kochi group': palette[5],
    'An. funestus group': palette[6],
    'An. barbirostris group': palette[7],
    'An. dirus complex': palette[9],
    '': palette[10],
    'An. gambiae': 'r',
    'An. arabiensis':'b',
    'An. funestus': 'k'
}


labels = {
    'An. annularis group': 'An. annularis group',
    'An. hyrcanus group': 'An. hyrcanus group',
    'An. maculatus group':'An. maculatus group',
    'An. kochi group': 'An. kochi group',
    'An. funestus group': 'An. funestus group',
    'An. barbirostris group': 'An. barbirostris group',
    'An. dirus complex':'An. dirus complex',
    '':'other'
}
#fix the df NaNs
df_meta_gud = df_meta_gud.fillna(value='')
#names so we can add labels etc
vbs = np.asarray(df_meta_gud.ROMA_ID)
vsc = np.asarray(df_meta_gud.Species_complex)
vsp = np.asarray(df_meta_gud.Species)
#plot
fig = plt.figure(figsize=(10, 13), )
gs = mpl.gridspec.GridSpec(nrows=1, ncols=2, width_ratios=(2.5, 0.4), wspace=0.22)

#dendro
ax = fig.add_subplot(gs[0, 0])
sns.despine(ax=ax, left=True, right=True, bottom=False)
r = scipy.cluster.hierarchy.dendrogram(
    z, no_labels=False, count_sort=True, 
    color_threshold=0, labels=vbs,
    above_threshold_color='k',
    ax=ax, leaf_font_size=9, orientation='left')
ax.set_xlabel('Distance', size=11)
#legend
handles = [mpl.patches.Patch(color=colors[p], label=labels[p]) for p in list(colors.keys())]
ax.legend(handles=handles, loc='upper right', bbox_to_anchor=(0.5, 1.011), ncol=1, fontsize=12)

#colour bar
ax = fig.add_subplot(gs[0, 1])
sns.despine(ax=ax, left=True, right=True, bottom=True)
#get the clustered order of samples
sp_com = vsc[r['leaves']]
#make list of species complex colours to match the clustered sample order
sp_com_colors = [colors[p] for p in sp_com]
#plot
ax.bar(left=np.full((len(vbs)), 0.5), height=np.full((96), 1), bottom=np.asarray(range(0,96)), color=sp_com_colors, edgecolor='k')
ax.set_xlim(0, 1)
ax.set_xticks([])
ax.set_ylim(0,96)
ax.set_xlabel('Sp. Group', size=11)
ax.set_yticks(np.arange(96)+0.5)
ax.set_yticklabels(vsp[r['leaves']], ha='left', position=(1.1,2), fontsize=10)
ax.tick_params(length=0.0)
ax.text(-0.4, -1.3, 'ID', fontsize=11)
ax.text(1.4, -1.3, 'Sp.', fontsize=11)

fig.savefig('../uniref_nomissing_nocrapsamp.jpeg', jpeg_quality=100, bbox_inches='tight');<jupyter_output><empty_output><jupyter_text>remove the WGA samples too<jupyter_code>#remove the WGA samples too?
isitnotWGA = df_meta_gud.Sample_type == "Genomic DNA"
np.count_nonzero(isitnotWGA)
gt_nowga = gt_nomiss.compress(isitnotWGA, axis=1)
gt_nowga
df_meta_nowga = df_meta_gud[isitnotWGA]
df_meta_nowga.shape
#names so we can add labels etc
vbs = np.asarray(df_meta_nowga.ROMA_ID)
vsc = np.asarray(df_meta_nowga.Species_complex)
vsp = np.asarray(df_meta_nowga.Species)
#measure distance and cluster samples
d = allel.stats.pairwise_distance(gt_nowga.to_n_alt(), metric='cityblock')
z = scipy.cluster.hierarchy.linkage(d, method='average')
#plot
fig = plt.figure(figsize=(10, 13), )
gs = mpl.gridspec.GridSpec(nrows=1, ncols=2, width_ratios=(2.5, 0.4), wspace=0.22)

#dendro
ax = fig.add_subplot(gs[0, 0])
sns.despine(ax=ax, left=True, right=True, bottom=False)
r = scipy.cluster.hierarchy.dendrogram(
    z, no_labels=False, count_sort=True, 
    color_threshold=0, labels=vbs,
    above_threshold_color='k',
    ax=ax, leaf_font_size=9, orientation='left')
ax.set_xlabel('Distance', size=11)
#legend
handles = [mpl.patches.Patch(color=colors[p], label=labels[p]) for p in list(colors.keys())]
ax.legend(handles=handles, loc='upper right', bbox_to_anchor=(0.5, 1.011), ncol=1, fontsize=12)

#colour bar
ax = fig.add_subplot(gs[0, 1])
sns.despine(ax=ax, left=True, right=True, bottom=True)
#get the clustered order of samples
sp_com = vsc[r['leaves']]
#make list of species complex colours to match the clustered sample order
sp_com_colors = [colors[p] for p in sp_com]
#plot
ax.bar(left=np.full((len(vbs)), 0.5), height=np.full((len(vbs)), 1), bottom=np.asarray(range(0,len(vbs))), color=sp_com_colors, edgecolor='k')
ax.set_xlim(0, 1)
ax.set_xticks([])
ax.set_ylim(0,len(vbs))
ax.set_xlabel('Sp. Group', size=11)
ax.set_yticks(np.arange(len(vbs))+0.5)
ax.set_yticklabels(vsp[r['leaves']], ha='left', position=(1.1,2), fontsize=10)
ax.tick_params(length=0.0)
ax.text(-0.4, -1.3, 'ID', fontsize=11)
ax.text(1.4, -1.3, 'Sp.', fontsize=11)

fig.savefig('../uniref_nomissing_nocrapsamp_noWGA.jpeg', jpeg_quality=100, bbox_inches='tight');
#use different metrics

#these seem to generate the cleanest clusters, using method=average is less good
d = allel.stats.pairwise_distance(gt_nowga.to_n_alt(), metric='cityblock')
z = scipy.cluster.hierarchy.linkage(d, method='complete')
#plot
fig = plt.figure(figsize=(10, 13), )
gs = mpl.gridspec.GridSpec(nrows=1, ncols=2, width_ratios=(2.5, 0.4), wspace=0.22)

#dendro
ax = fig.add_subplot(gs[0, 0])
sns.despine(ax=ax, left=True, right=True, bottom=False)
r = scipy.cluster.hierarchy.dendrogram(
    z, no_labels=False, count_sort=True, 
    color_threshold=0, labels=vbs,
    above_threshold_color='k',
    ax=ax, leaf_font_size=9, orientation='left')
ax.set_xlabel('Distance', size=11)
#legend
handles = [mpl.patches.Patch(color=colors[p], label=labels[p]) for p in list(colors.keys())]
ax.legend(handles=handles, loc='upper right', bbox_to_anchor=(0.5, 1.011), ncol=1, fontsize=12)

#colour bar
ax = fig.add_subplot(gs[0, 1])
sns.despine(ax=ax, left=True, right=True, bottom=True)
#get the clustered order of samples
sp_com = vsc[r['leaves']]
#make list of species complex colours to match the clustered sample order
sp_com_colors = [colors[p] for p in sp_com]
#plot
ax.bar(left=np.full((len(vbs)), 0.5), height=np.full((len(vbs)), 1), bottom=np.asarray(range(0,len(vbs))), color=sp_com_colors, edgecolor='k')
ax.set_xlim(0, 1)
ax.set_xticks([])
ax.set_ylim(0,len(vbs))
ax.set_xlabel('Sp. Group', size=11)
ax.set_yticks(np.arange(len(vbs))+0.5)
ax.set_yticklabels(vsp[r['leaves']], ha='left', position=(1.1,2), fontsize=10)
ax.tick_params(length=0.0)
ax.text(-0.4, -1.3, 'ID', fontsize=11)
ax.text(1.4, -1.3, 'Sp.', fontsize=11)

fig.savefig('../uniref_nomissing_nocrapsamp_noWGA_spam.jpeg', jpeg_quality=100, bbox_inches='tight');
#interestingly is does appear to affect the clustering...<jupyter_output><empty_output><jupyter_text>use all the variants - what happens to missing data?<jupyter_code>#take the all var-gud sample genotype, remove the WGA samples then cluster
gt_gudsam
gt_al = gt_gudsam.compress(isitnotWGA, axis=1)
gt_al
#measure distance and cluster samples
d = allel.stats.pairwise_distance(gt_al.to_n_alt(), metric='cityblock')
z = scipy.cluster.hierarchy.linkage(d, method='complete')
#plot
fig = plt.figure(figsize=(10, 13), )
gs = mpl.gridspec.GridSpec(nrows=1, ncols=2, width_ratios=(2.5, 0.4), wspace=0.22)

#dendro
ax = fig.add_subplot(gs[0, 0])
sns.despine(ax=ax, left=True, right=True, bottom=False)
r = scipy.cluster.hierarchy.dendrogram(
    z, no_labels=False, count_sort=True, 
    color_threshold=0, labels=vbs,
    above_threshold_color='k',
    ax=ax, leaf_font_size=9, orientation='left')
ax.set_xlabel('Distance', size=11)
#legend
handles = [mpl.patches.Patch(color=colors[p], label=labels[p]) for p in list(colors.keys())]
ax.legend(handles=handles, loc='upper right', bbox_to_anchor=(0.5, 1.011), ncol=1, fontsize=12)

#colour bar
ax = fig.add_subplot(gs[0, 1])
sns.despine(ax=ax, left=True, right=True, bottom=True)
#get the clustered order of samples
sp_com = vsc[r['leaves']]
#make list of species complex colours to match the clustered sample order
sp_com_colors = [colors[p] for p in sp_com]
#plot
ax.bar(left=np.full((len(vbs)), 0.5), height=np.full((len(vbs)), 1), bottom=np.asarray(range(0,len(vbs))), color=sp_com_colors, edgecolor='k')
ax.set_xlim(0, 1)
ax.set_xticks([])
ax.set_ylim(0,len(vbs))
ax.set_xlabel('Sp. Group', size=11)
ax.set_yticks(np.arange(len(vbs))+0.5)
ax.set_yticklabels(vsp[r['leaves']], ha='left', position=(1.1,2), fontsize=10)
ax.tick_params(length=0.0)
ax.text(-0.4, -1.3, 'ID', fontsize=11)
ax.text(1.4, -1.3, 'Sp.', fontsize=11)

fig.savefig('../uniref_nomissing_nocrapsamp_noWGA_allvar.jpeg', jpeg_quality=100, bbox_inches='tight');
#looks almost identical to no-missingness!<jupyter_output><empty_output><jupyter_text>can we annotate anophele/celia?<jupyter_code>df_meta_nowga
#looks like we need that info from Brandy<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16058276001618776, 'fraction_numerical': 0.027519222986645085, 'mean_word_length': 5.845429362880886, 'pattern_counts': {'":': 0, '<': 35, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '665046', 'license': 'MIT', 'n_tokens_mistral': 5202, 'n_tokens_neox': 4964, 'n_words': 1004, 'path': 'notebooks/20180706-uniref-distance.ipynb', 'repo_name': 'malariagen/vobs-sea-pilot-paper', 'size': '12355'}