app.py

import base64
import json
import os
import pickle
import re
import uuid
from time import time

import numpy as np
import pandas as pd
import pygad
import streamlit as st
from CGRtools import smiles
from CGRtools.containers import QueryContainer, MoleculeContainer
from CGRtools.utils import grid_depict
from VQGAE.models import VQGAE, OrderingNetwork
from VQGAE.utils import frag_counts_to_inds, restore_order, decode_molecules
from streamlit.components.v1 import html
from PIL import Image

MoleculeContainer.depict_settings(aam=False)


def download_button(object_to_download, download_filename, button_text, pickle_it=False):
    """
    Generates a link to download the given object_to_download.
    Params:
    ------
    object_to_download:  The object to be downloaded.
    download_filename (str): filename and extension of file. e.g. mydata.csv,
    some_txt_output.txt download_link_text (str): Text to display for download
    link.
    button_text (str): Text to display on download button (e.g. 'click here to download file')
    pickle_it (bool): If True, pickle file.
    Returns:
    -------
    (str): the anchor tag to download object_to_download
    Examples:
    --------
    download_link(your_df, 'YOUR_DF.csv', 'Click to download data!')
    download_link(your_str, 'YOUR_STRING.txt', 'Click to download text!')
    """
    if pickle_it:
        try:
            object_to_download = pickle.dumps(object_to_download)
        except pickle.PicklingError as e:
            st.write(e)
            return None

    else:
        if isinstance(object_to_download, bytes):
            pass

        elif isinstance(object_to_download, pd.DataFrame):
            object_to_download = object_to_download.to_csv(index=False).encode('utf-8')

        # Try JSON encode for everything else
        else:
            object_to_download = json.dumps(object_to_download)

    try:
        # some strings <-> bytes conversions necessary here
        b64 = base64.b64encode(object_to_download.encode()).decode()

    except AttributeError as e:
        b64 = base64.b64encode(object_to_download).decode()

    button_uuid = str(uuid.uuid4()).replace('-', '')
    button_id = re.sub('\d+', '', button_uuid)

    custom_css = f""" 
        <style>
            #{button_id} {{
                background-color: rgb(255, 255, 255);
                color: rgb(38, 39, 48);
                padding: 0.25em 0.38em;
                position: relative;
                text-decoration: none;
                border-radius: 4px;
                border-width: 1px;
                border-style: solid;
                border-color: rgb(230, 234, 241);
                border-image: initial;
            }} 
            #{button_id}:hover {{
                border-color: rgb(246, 51, 102);
                color: rgb(246, 51, 102);
            }}
            #{button_id}:active {{
                box-shadow: none;
                background-color: rgb(246, 51, 102);
                color: white;
                }}
        </style> """

    dl_link = custom_css + f'<a download="{download_filename}" id="{button_id}" href="data:file/txt;base64,{b64}">{button_text}</a><br></br>'

    return dl_link


def file_selector(folder_path='.'):
    filenames = os.listdir(folder_path)
    selected_filename = st.selectbox('Select a file', filenames)
    return os.path.join(folder_path, selected_filename)


def render_svg(svg_string):
    """Renders the given svg string."""
    c = st.container()
    with c:
        html(svg_string, height=300, scrolling=True)


# define groups to filter
allene = QueryContainer()
allene.add_atom("C")
allene.add_atom("A")
allene.add_atom("A")
allene.add_bond(1, 2, 2)
allene.add_bond(1, 3, 2)

peroxide_charge = QueryContainer()
peroxide_charge.add_atom("O", charge=-1)
peroxide_charge.add_atom("O")
peroxide_charge.add_bond(1, 2, 1)

peroxide = QueryContainer()
peroxide.add_atom("O")
peroxide.add_atom("O")
peroxide.add_bond(1, 2, 1)


def tanimoto_kernel(x, y):
    """
    "The Tanimoto coefficient is a measure of the similarity between two sets.
    It is defined as the size of the intersection divided by the size of the union of the sample sets."

    The Tanimoto coefficient is also known as the Jaccard index

    Adoppted from https://github.com/cimm-kzn/CIMtools/blob/master/CIMtools/metrics/pairwise.py

    :param x: 2D array of features.
    :param y: 2D array of features.
    :return: The Tanimoto coefficient between the two arrays.
    """
    x_dot = np.dot(x, y.T)

    x2 = (x ** 2).sum(axis=1)
    y2 = (y ** 2).sum(axis=1)

    len_x2 = len(x2)
    len_y2 = len(y2)

    result = x_dot / (np.array([x2] * len_y2).T + np.array([y2] * len_x2) - x_dot)
    result[np.isnan(result)] = 0

    return result


def fitness_func_batch(ga_instance, solutions, solutions_indices):
    frag_counts = np.array(solutions)
    if len(frag_counts.shape) == 1:
        frag_counts = frag_counts[np.newaxis, :]

    # prediction of activity by random forest
    rf_score = rf_model.predict_proba(frag_counts)[:, 1]

    # size penalty if molecule too small
    mol_size = frag_counts.sum(-1).astype(np.int64)
    size_penalty = np.where(mol_size < 18, -1.0, 0.)

    # adding dissimilarity so it generates different solutions
    dissimilarity_score = 1 - tanimoto_kernel(frag_counts, X).max(-1)
    dissimilarity_score += np.where(dissimilarity_score == 0, -5, 0)

    # full fitness function
    fitness = 0.5 * rf_score + 0.3 * dissimilarity_score + size_penalty

    # prediction of ordering score
    if use_ordering_score:
        frag_inds = frag_counts_to_inds(frag_counts, max_atoms=51)
        _, ordering_scores = restore_order(frag_inds, ordering_model)
        ordering_scores = np.array(ordering_scores)
        fitness += 0.2 * ordering_scores

    return fitness.tolist()


def on_generation_progress(ga):
    global ga_progress
    global ga_bar
    ga_progress = ga_progress + 1
    ga_bar.progress(ga_progress / num_generations, text=ga_progress_text)


@st.cache_data
def load_data(batch_size):
    X = np.load("saved_model/tubulin_qsar_class_train_data_vqgae.npz")["x"]
    Y = np.load("saved_model/tubulin_qsar_class_train_data_vqgae.npz")["y"]
    with open("saved_model/rf_class_train_tubulin.pickle", "rb") as inp:
        rf_model = pickle.load(inp)

    vqgae_model = VQGAE.load_from_checkpoint(
        "saved_model/vqgae.ckpt",
        task="decode",
        batch_size=batch_size,
        map_location="cpu"
    )
    vqgae_model = vqgae_model.eval()

    ordering_model = OrderingNetwork.load_from_checkpoint(
        "saved_model/ordering_network.ckpt",
        batch_size=batch_size,
        map_location="cpu"
    )
    ordering_model = ordering_model.eval()
    return X, Y, rf_model, vqgae_model, ordering_model


st.title('Inverse QSAR of Tubulin with VQGAE')

with st.sidebar:
    st.subheader("Optimisation parameters")
    with st.form("ga_options"):
        num_generations = st.slider(
            'Number of generations for GA',
            min_value=2,
            max_value=15,
            value=3
        )

        parent_selection_type = st.selectbox(
            label='Parent selection type',
            options=(
                'Steady-state selection',
                'Roulette wheel selection',
                'Stochastic universal selection',
                'Rank selection',
                'Random selection',
                'Tournament selection'
            ),
            index=1
        )

        parent_selection_translator = {
            "Steady-state selection": "sss",
            "Roulette wheel selection": "rws",
            "Stochastic universal selection": "sus",
            "Rank selection": "rank",
            "Random selection": "random",
            "Tournament selection": "tournament",
        }

        parent_selection_type = parent_selection_translator[parent_selection_type]

        crossover_type = st.selectbox(
            label='Crossover type',
            options=(
                'Single point',
                'Two points',
            ),
            index=0
        )

        crossover_translator = {
            "Single point": "single_point",
            "Two points": "two_points",
        }

        crossover_type = crossover_translator[crossover_type]

        num_parents_mating = int(
            st.slider(
                'Pecentage of parents mating taken from initial population',
                min_value=0,
                max_value=100,
                step=1,
                value=20,
            ) * 603 // 100  # TODO: remove hardcoded number
        )

        keep_parents = int(
            st.slider(
                'Percentage of parents kept taken from number of parents mating',
                min_value=0,
                max_value=100,
                step=1,
                value=50
            ) * num_parents_mating // 100
        )
        # 2/3 of num_parents_mating
        use_ordering_score = st.toggle('Use ordering score', value=False)
        batch_size = int(st.number_input("Batch size", value=250, placeholder="Type a number..."))
        random_seed = int(st.number_input("Random seed", value=44, placeholder="Type a number..."))
        submit_side = st.form_submit_button('Start optimisation')

intro_text = '''
This is a demo app for de novo design of inhinitors of Tubulin target based on 
[Vector Quantized Graph AutoEncoder](https://github.com/Laboratoire-de-Chemoinformatique/VQGAE) (VQGAE).
The search is done by Genetic Algorithm (GA) implemented in the [PyGAD libary](https://pygad.readthedocs.io/en/latest/).

The general scheme of *de novo* molecular optimisation presented in the figure below.
'''

image = Image.open('GA_optimization.png')
st.markdown(intro_text)
st.image(image, caption='The protocol of Inverse QSAR with VQGAE and Genetic Algorithm')

launch_text = '''
If you want to edit parameters for opimisation including GA parameters, batch size and seed, you can do it in the left
context menu. Another option is to launch calculations with a button below:
'''
st.markdown(launch_text)

submit_main = st.button('Run optimisation with default parameters')

if submit_side or submit_main:
    if submit_main:
        num_generations = 3
        parent_selection_type = "rws"
        crossover_type = "single_point"
        num_parents_mating = 120
        keep_parents = 60
        use_ordering_score = False
        batch_size = 250
        random_seed = 44

    st.toast('Optimisation is started. The progress will be printed below')
    X, Y, rf_model, vqgae_model, ordering_model = load_data(batch_size)
    assert X.shape == (603, 4096)
    start_time = time()
    ga_instance = pygad.GA(
        fitness_func=fitness_func_batch,
        on_generation=on_generation_progress,
        initial_population=X,
        num_genes=X.shape[-1],
        fitness_batch_size=batch_size,
        num_generations=num_generations,
        num_parents_mating=num_parents_mating,
        parent_selection_type=parent_selection_type,
        crossover_type=crossover_type,
        mutation_type="adaptive",
        mutation_percent_genes=[10, 5],
        # https://pygad.readthedocs.io/en/latest/pygad.html#use-adaptive-mutation-in-pygad
        save_best_solutions=False,
        save_solutions=True,
        keep_elitism=0,  # turn it off to make keep_parents work
        keep_parents=keep_parents,
        suppress_warnings=True,
        random_seed=random_seed,
        gene_type=int
    )

    ga_progress = 0
    ga_progress_text = "Genetic optimisation in progress. Please wait."
    ga_bar = st.progress(0, text=ga_progress_text)
    ga_instance.run()
    ga_bar.empty()

    with st.spinner('Getting unique solutions'):
        unique_solutions = list(set(tuple(s) for s in ga_instance.solutions))

    st.success(f'{len(unique_solutions)} solutions were obtained')

    scores = {
        "rf_score": [],
        "similarity_score": []
    }

    if use_ordering_score:
        scores["ordering_score"] = []

    rescoring_progress = 0
    rescoring_progress_text = "Rescoring obtained solutions"
    rescoring_bar = st.progress(0, text=rescoring_progress_text)
    total_rescoring_steps = len(unique_solutions) // batch_size + 1
    for rescoring_step in range(total_rescoring_steps):
        vqgae_latents = unique_solutions[rescoring_step * batch_size: (rescoring_step + 1) * batch_size]
        frag_counts = np.array(vqgae_latents)
        rf_scores = rf_model.predict_proba(frag_counts)[:, 1]
        similarity_scores = tanimoto_kernel(frag_counts, X).max(-1)
        scores["rf_score"].extend(rf_scores.tolist())
        scores["similarity_score"].extend(similarity_scores.tolist())
        if use_ordering_score:
            frag_inds = frag_counts_to_inds(frag_counts, max_atoms=51)
            _, ordering_scores = restore_order(frag_inds, ordering_model)
            scores["ordering_score"].extend(ordering_scores)
        rescoring_bar.progress(rescoring_step / total_rescoring_steps, text=rescoring_progress_text)
    sc_df = pd.DataFrame(scores)
    rescoring_bar.empty()

    if use_ordering_score:
        chosen_gen = sc_df[
            (sc_df["similarity_score"] < 0.95) & (sc_df["rf_score"] > 0.5) & (sc_df["ordering_score"] > 0.7)
            ]
    else:
        chosen_gen = sc_df[
            (sc_df["similarity_score"] < 0.95) & (sc_df["rf_score"] > 0.5)
            ]

    chosen_ids = chosen_gen.index.to_list()
    chosen_solutions = np.array([unique_solutions[ind] for ind in chosen_ids])
    gen_frag_inds = frag_counts_to_inds(chosen_solutions, max_atoms=51)

    st.info(f'The number of chosen solutions is {gen_frag_inds.shape[0]}', icon="ℹ️")

    results = {"smiles": [], "ordering_score": [], "validity": []}
    decoding_progress = 0
    decoding_progress_text = "Decoding chosen solutions"
    decoding_bar = st.progress(0, text=decoding_progress_text)
    total_decoding_steps = gen_frag_inds.shape[0] // batch_size + 1
    for decoding_step in range(total_decoding_steps):
        inputs = gen_frag_inds[decoding_step * batch_size: (decoding_step + 1) * batch_size]
        canon_order_inds, scores = restore_order(
            frag_inds=inputs,
            ordering_model=ordering_model,
        )
        molecules, validity = decode_molecules(
            ordered_frag_inds=canon_order_inds,
            vqgae_model=vqgae_model,
            clean_2d=False
        )
        results["smiles"].extend([str(molecule) for molecule in molecules])
        results["ordering_score"].extend(scores)
        results["validity"].extend([1 if i else 0 for i in validity])
        decoding_bar.progress(decoding_step / total_decoding_steps, text=decoding_progress_text)
    gen_stats = pd.DataFrame(results)
    decoding_bar.empty()
    full_stats = pd.concat([gen_stats, chosen_gen.reset_index()[["similarity_score", "rf_score"]]], axis=1)
    full_stats = full_stats[["smiles", "similarity_score", "rf_score", "ordering_score", "validity"]]
    valid_gen_stats = full_stats[full_stats.validity == 1]

    time_used = time() - start_time
    mins_used = int(round(time_used // 60, 0))
    seconds_used = int(round(time_used % 60, 0))
    st.success(f'{valid_gen_stats.shape[0]} valid solutions were obtained in {mins_used} mins {seconds_used} secs')
    st.balloons()

    st.subheader('Generation results', divider='rainbow')
    st.dataframe(valid_gen_stats)
    download_button_valid = download_button(
        object_to_download=valid_gen_stats,
        download_filename='vqgae_tubulin_inhibitors_valid.csv',
        button_text="Download results as CSV"
    )
    st.markdown(download_button_valid, unsafe_allow_html=True)
    # st.download_button(
    #     label="Download results as CSV",
    #     data=convert_df(filtered_gen_stats),
    #     file_name='vqgae_tubulin_inhibitors_valid.csv',
    #     mime='text/csv',
    # )

    st.subheader('Examples of generated molecules')

    examples_smiles = valid_gen_stats.sort_values(
        by=["similarity_score"],
        ascending=False
    ).iloc[:15:3].smiles.to_list()

    examples = []
    for smi in examples_smiles:
        mol = smiles(smi)
        mol.clean2d()
        examples.append(mol)
    svg = grid_depict(examples, 2)
    render_svg(svg)

    with st.expander("Show full stats"):
        st.dataframe(full_stats)

        download_button_full = download_button(
            object_to_download=full_stats,
            download_filename='vqgae_tubulin_inhibitors_full.csv',
            button_text="Download full results as CSV"
        )
        st.markdown(download_button_full, unsafe_allow_html=True)

        # st.download_button(
        #     label="Download full results as CSV",
        #     data=convert_df(full_stats),
        #     file_name='vqgae_tubulin_inhibitors_full.csv',
        #     mime='text/csv',
        # )

    if st.button("Restart"):
        st.rerun()