<jupyter_start><jupyter_text>genotypes<jupyter_code>callset = zarr.open_group('new_uniref_merged.zarr/', mode='r')
callset
gt = allel.GenotypeArray(callset['calldata/GT'][:])
gt
chrom = callset['variants/CHROM']
chrom[:]<jupyter_output><empty_output><jupyter_text>meta data<jupyter_code>df_meta = pd.read_csv('1175-VO-KH-STLAURENT-S01-samples.csv', sep=',')
df_meta.columns, len(df_meta)
#new samples to the meta data - we only need some columms
new_meta = pd.DataFrame()
new_meta['ROMA_ID'] = list(np.concatenate([callset['samples'][:10], list(callset['samples'][-5:])]))
new_meta['Status'] = 'sequenced'
new_meta['Species'] = ['An. arabiensis' for _ in range(5)] + ['An. gambiae' for _ in range(5)] + ['An. funestus' for _ in range(5)]
new_meta['Species_complex'] = ['An. arabiensis' for _ in range(5)] + ['An. gambiae' for _ in range(5)] + ['An. funestus' for _ in range(5)]
new_meta['Sample_type'] = 'Genomic DNA'

df_meta = df_meta.append(new_meta, ignore_index=True)
df_meta
#keep only sequenced
df_meta_seq = df_meta[df_meta.Status == 'sequenced']
len(df_meta_seq)
#sort the meta data so in the same order as the genotypes
df_meta_seq = df_meta_seq.sort_values('ROMA_ID')
df_meta_seq.head()
all_128_samples = np.asarray(df_meta_seq.ROMA_ID)
all_128_samples<jupyter_output><empty_output><jupyter_text>remove pf gene<jupyter_code>gt.shape
pfbool = chrom[:]!='Pf_M76611'
gt_nopf = gt.compress(pfbool)
gt_nopf.shape<jupyter_output><empty_output><jupyter_text>missingness<jupyter_code>#first 10 and last 5 are the new samples, we want to keep these.
callset['samples'][:]
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)
#looks fine, we keep the new outgroups
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>with no missingness<jupyter_code>#if we just use the 0 missingness sites what happens?
gt_nomiss = gt_gudsam.compress(full)
gt_nomiss
#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],
    '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':'g',
    'An. funestus': 'b'
}


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',
    'An. gambiae': 'An. gambiae',
    'An. arabiensis': 'An. arabiensis',
    'An. funestus': 'An. funestus',
    '':'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)
#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);<jupyter_output><empty_output><jupyter_text>all variants?<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('../NEW_uniref_nomissing_nocrapsamp_noWGA_allvar.jpeg', jpeg_quality=100, bbox_inches='tight');<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup