<jupyter_start><jupyter_text>define plotting + statistical functions<jupyter_code>def hierarchical_clustering(data, label, linkage, ax):
    import scipy.cluster.hierarchy as shc

    dend = shc.dendrogram(shc.linkage(data, method=linkage), 
                      ax= ax, labels=label,
                      leaf_font_size=8, leaf_rotation=90)
    
    from sklearn.cluster import AgglomerativeClustering
    cluster = AgglomerativeClustering(n_clusters=7, 
                                      affinity = 'euclidean', linkage = linkage)
    cluster_labels = cluster.fit_predict(data)
    
    return dend, cluster_labels
def do_pearson(x, y):
    from scipy.stats import pearsonr
    corr, p_val = pearsonr(x, y)
    return corr, p_val
def capping(input_df):
    # cap values above 1 and below 0 to 20% above/below
    df_capped = input_df.copy()

    for i in range(len(input_df.columns)-1):
        i+=1
        for j in range(len(input_df)):
            if input_df.iloc[j, i] > 1.2:
                df_capped.iloc[j, i] = 1.2
            elif input_df.iloc[j, i] < -0.2:
                df_capped.iloc[j, i ] = -0.2
    return df_capped
def plot_heatmap(data, columns, labels, cmap):
    fig, ax = plt.subplots(figsize=(10,15), dpi=100)

    im = ax.imshow(data, cmap=cmap, aspect='equal')

    # show all ticks...
    ax.set_xticks(np.arange(len(columns)))
    ax.set_yticks(np.arange(len(data)))
    ax.tick_params(top=True, bottom=False,
                   labeltop=True, labelbottom=False)

    # label ticks
    ax.set_xticklabels(columns, fontsize=8, rotation=70, 
                       rotation_mode='anchor', ha='left', va='center')
    ax.set_yticklabels(labels, fontsize=8)

    # set colorbar
    cbar = fig.colorbar(im, cax = fig.add_axes([0.6, 0.5, 0.01, 0.1]), 
                    ax=ax, orientation='vertical', fraction=.1)
    cbar.ax.tick_params(labelsize=6)
    
    # hatch the background axes patch
    #ax.patch.set(hatch='xxx', edgecolor='black')
    ax.patch.set(color = 'grey', facecolor = 'grey')
    
    return fig
def scatter_plotting(labels, label_xval, label_yval, data, x_name, y_name):
    colors = {0:'cornflowerblue', 1:'darkblue', 2:'darkorange',
          3:'gold', 4:'lightsteelblue', 5:'blue',
          6:'blueviolet', 7:'cyan', 8:'black'}
    
    fig = plt.figure(figsize=(10,10), dpi=100)
    ax = fig.add_subplot(111)
    
    for i in range(len(labels)): # plot each point + it's index as text above
        ax.text(label_xval[i], label_yval[i], 
                labels[i], size=8, color='k')

    for j in data['cluster_labels'].unique():
        df_subset = data[data['cluster_labels']==j]
        ax.scatter(df_subset[x_name], df_subset[y_name], 
                   alpha=0.6, s=80, c=colors[j])

    ax.set_xlabel(x_name)
    ax.set_ylabel(y_name)

    return fig, ax
def correlate(df_to_plot, x_name, y_name):

    # plot scatter
    fig, ax = scatter_plotting(df_to_plot['variant'], df_to_plot[x_name], df_to_plot[y_name], df_to_plot, x_name, y_name)
    
    return fig, ax
def save_fig(fig, title=''):
    # output figure
    mpl.rcParams['pdf.fonttype'] = 42
    plt.savefig(f'{title}.pdf', transparent=True)
    return fig<jupyter_output><empty_output><jupyter_text>read in datasets<jupyter_code>df_wing = pd.read_csv('fly_wing_size.csv', index_col=False)
df_eclosion = pd.read_csv('fly_eclosion.csv', index_col=False)
df_yeast = pd.read_csv('PTEN_Mighell_2018.csv', index_col=False)
df_vac14_hek = pd.read_csv('PTEN_VAC14_HEK.csv', index_col=False)
df_mcf10a = pd.read_csv('PTEN_MCF10A.csv', index_col=False)<jupyter_output><empty_output><jupyter_text>inspect datasets and make adjustments<jupyter_code>''' update 4/23/20: Payel requested to drop variants: 
eclosion - E256K, F56C, G132D, I135V, M198I, N340D, N340H, N356H, Q396R, R130X
wing size - N340H, N340D, C136M fxS44'''

'''update 4/9/20: now wing size data is normalized such that
0 is loss of function, 1 is functional'''<jupyter_output><empty_output><jupyter_text>inspect wing size data<jupyter_code>df_wing.head()
# drop 4A so it doesn't skew data
df_wing.drop(2, axis=0, inplace=True)
print(df_wing.head())<jupyter_output>variant  wing_size
0  control   0.000000
1       WT   1.000000
3    C124S  -0.465176
4     D22E  -0.902782
5     P38H  -0.208456<jupyter_text>inspect eclosion time data<jupyter_code>df_eclosion.head()
df_eclosion.rename(columns={'time':'eclosion_time'}, inplace=True)
df_eclosion.head()<jupyter_output><empty_output><jupyter_text>inspect Mighell (yeast) data<jupyter_code>df_yeast.head()
# rename col name for indexing the merge letter
df_yeast.rename(columns={'Variant (one letter)':'variant','Cum_score':'Mighell_score'}, 
                inplace=True)
# keep only two needed columns
df_yeast = df_yeast[df_yeast.columns[[0, 6]]]

# since Mighell's data is already normalized
# I'll put in WT as 0, KO as -4
df_yeast.loc[0, :] = ['WT', 0]
df_yeast.loc[1, :] = ['control', -4]
df_yeast.head()<jupyter_output><empty_output><jupyter_text>inspect HEK, VAC14, MCF10A dsets<jupyter_code>print(df_vac14_hek.head(), '\n')
print(df_mcf10a.head())<jupyter_output>variant     VAC14       HEK
0   no_PTEN  0.000000  0.000000
1        WT  1.000000  1.000000
2      P38H -0.257287 -0.378378
3     R130L -0.255608 -0.189189
4  C124S-4A -0.227811 -0.568047 

   variant  MCF10A
0  no_PTEN  0.0000
1       WT  1.0000
2      K6E  0.6689
3      K6I  0.5456
4     R15S  0.6130<jupyter_text>merge datasets<jupyter_code># start merging different datasets
df_merged = df_wing.merge(df_eclosion[['variant', 'eclosion_time']], on='variant', how='left') # merge wing+eclosion
df_merged = df_merged.merge(df_yeast, on='variant', how='left') # merge wing/eclosion+yeast

# rename control to no_PTEN
df_merged['variant'].replace({'control':'no_PTEN'}, inplace=True)

df_merged.head()<jupyter_output><empty_output><jupyter_text>normalize + cap eclosion and Mighell for plotting heatmap<jupyter_code># normalized eclosion and Mighell
df_normed = df_merged.copy()

for i in range(2):
    i+=2
    for j in range(len(df_merged)):
        df_normed.iloc[j, i] = \
        ((df_merged.iloc[j, i]-df_merged.iloc[0, i])/ (df_merged.iloc[1, i]-df_merged.iloc[0, i]))

df_normed.head()
#merge in VAC14, HEK and MCF10A dsets

df_all = df_normed.copy()

for i in [df_vac14_hek, df_mcf10a]:
    df_all = df_all.merge(i, how='left', on='variant')

df_all.head()
# cap values above 1 and below 0 to 20% above/below for plotting heatmap
df_capped = capping(df_all)
df_capped<jupyter_output><empty_output><jupyter_text>analyze all data  heatmap for all data<jupyter_code># grab AA positions
AA_position = []
for i in df_all['variant']:
    try:
        AA_position.append(int(i[1:-1]))
    except ValueError:
        AA_position.append(0)

df_all['AA_position'] = AA_position

# sort by aa position
df_all = df_all.sort_values(by=['AA_position'], axis=0).reset_index(drop=True)
# capping values above/blow 1/0 to 1.2/-0.2
df_all_capped = capping(df_all)
df_all_capped
#all_data_scaled = MinMaxScaler().fit_transform(df_all.iloc[:, 1:-1])
heatmap_data = df_all_capped.iloc[:, 1:-1]
labels = df_all_capped['variant'].tolist()

#heatmap_data = all_data_scaled
#labels = df_all['variant'].tolist()


fig = plot_heatmap(heatmap_data, df_all_capped.columns[1:-1], labels, 'Reds_r')

# save fig
mpl.rcParams['ps.fonttype'] = 42
fig.savefig('data_comparison_heatmap_by_aa.pdf', transparent=True)

plt.show()<jupyter_output><empty_output><jupyter_text>clustering all data<jupyter_code>linkage = 'complete'

fig = plt.figure(figsize=(12,5))
ax = fig.add_subplot(111)
# hierarchical_clustering(data, label, linkage, ax)
dend, cluster_labels = hierarchical_clustering(
                                data = df_all.iloc[:, 1:-2].dropna().reset_index(drop=True), 
                                label = df_all.iloc[:, :-2].dropna()['variant'].to_list(), 
                                linkage = linkage,
                                ax = ax)

ax.set_title(f'wing_size, eclosion_time, Mighell, VAC14, HEK clustering by "{linkage} linkage"')

# output figure
mpl.rcParams['pdf.fonttype'] = 42
plt.savefig(f'all_data_dendrogram ({linkage}).pdf', transparent=True)

plt.show()
# re-index data according to cluster labels
data_clustered = df_all_capped.iloc[:, :-2].dropna().reset_index(drop=True).reindex(dend['leaves'])
data_clustered.head()
# plot heatmap in order of cluster labels
# plot_heatmap(data, columns, labels, cmap)
fig = plot_heatmap(data_clustered.iloc[:,1:], data_clustered.columns[1:], data_clustered['variant'], 'Reds_r')

# save fig
mpl.rcParams['ps.fonttype'] = 42
plt.savefig(f'all_data_heatmap_clustered ({linkage}).pdf', transparent=True)

plt.show()<jupyter_output><empty_output><jupyter_text>3-way comparisons (wing vs HEK vs Mighell)<jupyter_code>df_all.head()
clustering_data = df_all.drop(columns=['eclosion_time','VAC14','MCF10A']).dropna().reset_index(drop=True)
clustering_data
cols = ['variant', 'wing_size', 'HEK', 'Mighell_score']
clustering_data = clustering_data[cols]
clustering_data.head()
linkage_method = 'complete'

fig = plt.figure(figsize=(10,5), dpi=100)
ax = fig.add_subplot(111)
dend, cluster_labels = hierarchical_clustering(
                              clustering_data.iloc[:, 1:], 
                              clustering_data['variant'].tolist(),
                              linkage_method, ax)

ax.set_title(f"3-way clustering wing/HEK/Mighell (by {linkage_method} linkage)")

#save_fig(fig, title=f'WHM_dendrogram ({linkage_method})')

plt.show()
# re-index data according to cluster labels
heatmap_data = df_capped[['variant', 'wing_size', 'Mighell_score', 'HEK']].dropna().reset_index(drop=True).reindex(dend['leaves'])
heatmap_data.head()
heatmap_data = heatmap_data[cols]
# plot heatmap in order of cluster labels
# plot_heatmap(data, columns, labels, cmap)
fig = plot_heatmap(heatmap_data.iloc[:,1:], heatmap_data.columns[1:], heatmap_data['variant'], 'Reds_r')

save_fig(fig, title=f'WHM_heatmap_clustered ({linkage_method})')

plt.show()<jupyter_output><empty_output><jupyter_text>plot correlation (scatter) plots using cluster assignments<jupyter_code># prep the data for scatter plotting
d = {'variant':clustering_data['variant'], 'cluster_labels':cluster_labels}
df_labels = pd.DataFrame(d)
print(df_labels)
np.unique(cluster_labels)
df_to_plot = df_all.merge(df_labels, on='variant', how='left')
display(df_to_plot)
df_to_plot[df_to_plot['variant']=='I203V']
# fill na with 8 to plot them gray later in the scatterplot
filler = {'cluster_labels':8}
df_to_plot.fillna(value=filler, inplace=True)
x = 'Mighell_score'
y = 'HEK'
df_temp = df_to_plot[['variant', x, y, 'cluster_labels']].dropna().reset_index(drop=True)

# calculate overall Pearson's r
corr, p = do_pearson(df_temp[x], df_temp[y])
pear_r = f"Pearson's r={corr:.4f} (p={p:.4f})"

# plot scatterplot
fig, ax = correlate(df_temp, x, y)
text_x = ax.get_xlim()[0]*0.9
text_y = ax.get_ylim()[1]*0.95
ax.text(text_x, text_y, pear_r)
title = f'{y} vs {x}, colored by W-H-M clustering ({linkage_method} linkage)'
ax.set_title(title)

save_fig(fig, title=f"{y}-{x} (WHM clustering)")

plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14024061770515353, 'fraction_numerical': 0.03609265577302927, 'mean_word_length': 4.2741477272727275, 'pattern_counts': {'":': 0, '<': 52, '<?xml version=': 0, '>': 52, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '123872', 'license': 'MIT', 'n_tokens_mistral': 4580, 'n_tokens_neox': 4181, 'n_words': 812, 'path': 'FlyWing_data_comparisons-Copy1.ipynb', 'repo_name': 'jessecanada/PTEN_Fly', 'size': '11138'}