import streamlit as st import re import numpy as np import pandas as pd import pickle import sklearn import catboost import shap import plotly.tools as tls from dash import dcc from sklearn.cluster import KMeans from sklearn.decomposition import PCA from sklearn.metrics import silhouette_score import plotly.express as px import matplotlib.pyplot as plt import plotly.graph_objs as go import plotly.graph_objects as go try: import matplotlib.pyplot as pl from matplotlib.colors import LinearSegmentedColormap from matplotlib.ticker import MaxNLocator except ImportError: pass st.set_option('deprecation.showPyplotGlobalUse', False) seed = 0 annotations = pd.read_csv("all_genes_imputed_features.csv") annotations.fillna(0, inplace=True) annotations = annotations.set_index("Gene") model_path = "best_model_fitted.pkl" with open(model_path, 'rb') as file: catboost_model = pickle.load(file) probabilities = catboost_model.predict_proba(annotations) prob_df = pd.DataFrame(probabilities, index=annotations.index, columns=['Probability_Most_Likely', 'Probability_Probable', 'Probability_Least_Likely']) df_total = pd.concat([prob_df, annotations], axis=1) # Create tabs for navigation with st.sidebar: st.sidebar.title("Navigation") tab = st.sidebar.radio("Go to", ("Gene Prioritisation", "Interactive SHAP Plot", "Supervised SHAP Clustering")) st.title('Blood Pressure Gene Prioritisation Post-GWAS') st.markdown("""A machine learning pipeline for predicting disease-causing genes post-genome-wide association study in blood pressure.""") # Define a function to collect genes from input collect_genes = lambda x: [str(i) for i in re.split(",|,\s+|\s+", x) if i != ""] input_gene_list = st.text_input("Input a list of multiple HGNC genes (enter comma separated):") gene_list = collect_genes(input_gene_list) explainer = shap.TreeExplainer(catboost_model) @st.cache_data def convert_df(df): return df.to_csv(index=False).encode('utf-8') probability_columns = ['Probability_Most_Likely', 'Probability_Probable', 'Probability_Least_Likely'] features_list = [column for column in df_total.columns if column not in probability_columns] features = df_total[features_list] # Page 1: Gene Prioritisation if tab == "Gene Prioritisation": if len(gene_list) > 1: df = df_total[df_total.index.isin(gene_list)] df['Gene'] = df.index df.reset_index(drop=True, inplace=True) required_columns = ['Gene'] + probability_columns + [column for column in df.columns if column not in probability_columns and column != 'Gene'] df = df[required_columns] st.dataframe(df) output = df[['Gene'] + probability_columns] csv = convert_df(output) st.download_button("Download Gene Prioritisation", csv, "bp_gene_prioritisation.csv", "text/csv", key='download-csv') df_shap = df.drop(columns=probability_columns + ['Gene']) shap_values = explainer.shap_values(df_shap) col1, col2 = st.columns(2) class_names = ["Most likely", "Probable", "Least likely"] with col1: st.subheader("Global SHAP Summary Plot") shap.summary_plot(shap_values, df_shap, plot_type="bar", class_names=class_names) st.pyplot(bbox_inches='tight', clear_figure=True) with col2: st.subheader(f"{class_names[0]} Gene Prediction") shap.summary_plot(shap_values[0], df_shap) st.pyplot(bbox_inches='tight', clear_figure=True) col3, col4 = st.columns(2) with col3: st.subheader(f"{class_names[1]} Gene Prediction") shap.summary_plot(shap_values[1], df_shap) st.pyplot(bbox_inches='tight', clear_figure=True) with col4: st.subheader(f"{class_names[2]} Gene Prediction") shap.summary_plot(shap_values[2], df_shap) st.pyplot(bbox_inches='tight', clear_figure=True) else: pass input_gene = st.text_input("Input an individual HGNC gene:") if input_gene: df2 = df_total[df_total.index == input_gene] class_names = ["Most likely", "Probable", "Least likely"] if not df2.empty: df2['Gene'] = df2.index df2.reset_index(drop=True, inplace=True) required_columns = ['Gene'] + probability_columns + [col for col in df2.columns if col not in probability_columns and col != 'Gene'] df2 = df2[required_columns] st.dataframe(df2) if ' ' in input_gene or ',' in input_gene: st.write('Input Error: Please input only a single HGNC gene name with no white spaces or commas.') else: df2_shap = df_total.loc[[input_gene], [col for col in df_total.columns if col not in probability_columns + ['Gene']]] print(df2_shap.columns) shap_values = explainer.shap_values(df2_shap) shap.getjs() for i in range(3): st.subheader(f"Force Plot for {class_names[i]} Prediction") force_plot = shap.force_plot( explainer.expected_value[i], shap_values[i], df2_shap, matplotlib=True, show=False ) st.pyplot(fig=force_plot) else: st.write("Gene not found in the dataset.") else: pass url = f"https://astrazeneca-cgr-publications.github.io/DrugnomeAI/geneview.html?gene={input_gene}" markdown_link = f"[{input_gene} druggability in DrugnomeAI]({url})" st.markdown(markdown_link, unsafe_allow_html=True) st.markdown(""" ### Total Gene Prioritisation Results for All Genes: """) df_total_output = df_total df_total_output['Gene'] = df_total_output.index #df_total_output.reset_index(drop=True, inplace=True) st.dataframe(df_total_output) csv = convert_df(df_total_output) st.download_button("Download Gene Prioritisation", csv, "all_genes_bp_prioritisation.csv", "text/csv", key='download-all-csv') elif tab == "Interactive SHAP Plot": st.title("Interactive SHAP Plot") if len(gene_list) > 1: df = df_total[df_total.index.isin(gene_list)] df['Gene'] = df.index df.reset_index(drop=True, inplace=True) required_columns = ['Gene'] + probability_columns + [column for column in df.columns if column not in probability_columns and column != 'Gene'] df = df[required_columns] st.dataframe(df) output = df[['Gene'] + probability_columns] csv = convert_df(output) st.download_button("Download Gene Prioritisation", csv, "bp_gene_prioritisation.csv", "text/csv", key='download-csv') df_shap = df.drop(columns=probability_columns + ['Gene']) shap_values = explainer.shap_values(df_shap) shap_values_first_class = shap_values[0] feature_importance = np.abs(shap_values_first_class).mean(axis=0) top_features_indices = np.argsort(feature_importance)[-20:] features_top = df_shap.columns[top_features_indices][::-1] shap_values_top = shap_values_first_class[:, top_features_indices][..., ::-1] # Prepare data for a single trace x_values = [] y_values = [] hover_texts = [] for i, feature_name in enumerate(features_top): for gene, value in zip(df['Gene'], shap_values_top[:, i]): x_values.append(value) y_values.append(feature_name) hover_texts.append(f'{gene}: {value:.3f}') # Create a single trace for the plot fig = go.Figure(data=go.Scatter( x=x_values, y=y_values, mode='markers', marker=dict( color=x_values, # Set color to SHAP values colorbar=dict(title="SHAP Value"), colorscale=[(0, "blue"), (1, "red")], # Blue to Red color scale ), text=hover_texts, # Set hover text hoverinfo="text+x" # Display hover text and x-value (SHAP value) )) fig.update_layout( title="SHAP Summary Plot - Top 20 Features", xaxis_title="SHAP Value", yaxis=dict(autorange="reversed", title="Feature"), showlegend=False, ) st.plotly_chart(fig, use_container_width=True) st.caption("SHAP Summary Plot of All Input Genes - Top 20 Features") elif tab == "Supervised SHAP Clustering": st.title("Supervised SHAP Clustering") training_genes = pd.read_csv("training_cleaned.csv") training_genes = training_genes[training_genes['BPlabel_encoded'] == 0] training_genes.set_index('Gene', inplace=True) # Calculate SHAP values for the full dataset shap_values_full = explainer.shap_values(annotations) shap_values_full_array = np.array(shap_values_full[0]) # Apply PCA to reduce dimensionality for visualization pca = PCA(n_components=2) shap_values_pca = pca.fit_transform(shap_values_full_array) # Apply clustering on the PCA-reduced SHAP values kmeans = KMeans(n_clusters=3, random_state=0).fit(shap_values_pca) # Get cluster labels for each point in the dataset labels = kmeans.labels_ # Prepare a DataFrame for visualization df_for_plot = pd.DataFrame({ 'PCA_1': shap_values_pca[:, 0], 'PCA_2': shap_values_pca[:, 1], 'Cluster': labels.astype(str), 'Gene': annotations.index, 'Type': 'Clustered Gene' }) # Add a new column for marking the special groups df_for_plot['SpecialGroup'] = 'None' df_for_plot.loc[df_for_plot['Gene'].isin(training_genes.index), 'SpecialGroup'] = 'Most Likely Training Gene' if gene_list: df_for_plot.loc[df_for_plot['Gene'].isin(gene_list), 'SpecialGroup'] = 'User Input Gene' # Initialize an empty figure fig = go.Figure() # Define color mapping for clusters cluster_colors = ['#ADD8E6', '#87CEEB', '#1E90FF'] # Plot clustered genes based on PCA components for i, cluster in enumerate(df_for_plot['Cluster'].unique()): filtered_df = df_for_plot[(df_for_plot['Cluster'] == cluster) & (df_for_plot['SpecialGroup'] == 'None')] fig.add_trace(go.Scatter( x=filtered_df['PCA_1'], y=filtered_df['PCA_2'], mode='markers', name=f'Cluster {cluster}', text=filtered_df['Gene'], marker=dict(color=cluster_colors[i]), hoverinfo="text+x+y", )) # Overlay "Most Likely Training Gene" filtered_df = df_for_plot[df_for_plot['SpecialGroup'] == 'Most Likely Training Gene'] fig.add_trace(go.Scatter( x=filtered_df['PCA_1'], y=filtered_df['PCA_2'], mode='markers', name='Most Likely Training Gene', text=filtered_df['Gene'], marker=dict(color='black'), hoverinfo="text+x+y", )) # Overlay "User Input Gene" filtered_df = df_for_plot[df_for_plot['SpecialGroup'] == 'User Input Gene'] fig.add_trace(go.Scatter( x=filtered_df['PCA_1'], y=filtered_df['PCA_2'], mode='markers', name='User Input Gene', text=filtered_df['Gene'], marker=dict(color='purple'), hoverinfo="text+x+y", )) # Customize layout fig.update_layout( title='Supervised SHAP Clustering with PCA', xaxis_title='First Principal Component', yaxis_title='Second Principal Component', showlegend=True, legend_title_text='Gene Category', ) st.plotly_chart(fig, use_container_width=True)