{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":67356,"databundleVersionId":8006601,"sourceType":"competition"}],"dockerImageVersionId":30673,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# ECFP v.s. Topological FP v.s. Atom pair FP\n### Some of the code below is partly based on Andrew D. Blevins's notebook [Leash Tutorial - ECFPs and Random Forest](https://www.kaggle.com/code/andrewdblevins/leash-tutorial-ecfps-and-random-forest)\n\n## Introduction\n\nVarious types of fingerprint descriptors exist, and the effective one depends on the task. \n\nAmong them, ECFP4 is the most commonly used. ECFP4 is the ECFP with a maximum radius R = 2, which counts substructures up to two bonds away from each atom. Obviously, this descriptor reflects only information about the local substructure of the molecule. One way to mitigate this disadvantage is to set radius to a large value, but this is not an efficient method because a larger length (i.e., the number of bits) is required to reflect sufficient information in the descriptor.\n\nThe binding affinity, which is the subject of this competition, is likely to reflect not only the local structure but also the shape, or topology, of the entire molecule. So, let's compare different FPs. Here, in addition to ECFP4, topological FP (as defined in RDKit) and Atom pair FP. Both take into account substructures larger than ECFP4.","metadata":{}},{"cell_type":"code","source":"!pip install rdkit\n!pip install duckdb","metadata":{"execution":{"iopub.status.busy":"2024-04-06T00:47:08.271747Z","iopub.execute_input":"2024-04-06T00:47:08.272191Z","iopub.status.idle":"2024-04-06T00:47:44.697960Z","shell.execute_reply.started":"2024-04-06T00:47:08.272159Z","shell.execute_reply":"2024-04-06T00:47:44.696758Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Preparation","metadata":{}},{"cell_type":"code","source":"import duckdb\nimport pandas as pd\nfrom rdkit import Chem\nfrom rdkit.Chem import AllChem\nfrom sklearn.ensemble import RandomForestClassifier\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import average_precision_score\nfrom sklearn.preprocessing import OneHotEncoder\nimport copy\n\nn_train = 30000 # the number of train data to sample\ntest_size = 0.2\nrandom_state = 42\nclf = RandomForestClassifier(n_estimators=100, random_state=random_state)","metadata":{"execution":{"iopub.status.busy":"2024-04-06T00:47:44.700884Z","iopub.execute_input":"2024-04-06T00:47:44.701340Z","iopub.status.idle":"2024-04-06T00:47:48.504778Z","shell.execute_reply.started":"2024-04-06T00:47:44.701301Z","shell.execute_reply":"2024-04-06T00:47:48.503552Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\ntrain_path = '/kaggle/input/leash-BELKA/train.parquet'\n\ncon = duckdb.connect()\n\ndf = con.query(f\"\"\"(SELECT *\n                        FROM parquet_scan('{train_path}')\n                        WHERE binds = 0\n                        ORDER BY random()\n                        LIMIT {n_train})\n                        UNION ALL\n                        (SELECT *\n                        FROM parquet_scan('{train_path}')\n                        WHERE binds = 1\n                        ORDER BY random()\n                        LIMIT {n_train})\"\"\").df()\n\ncon.close()","metadata":{"execution":{"iopub.status.busy":"2024-04-06T00:47:48.506297Z","iopub.execute_input":"2024-04-06T00:47:48.506808Z","iopub.status.idle":"2024-04-06T00:48:40.481456Z","shell.execute_reply.started":"2024-04-06T00:47:48.506777Z","shell.execute_reply":"2024-04-06T00:48:40.479945Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Feature Preprocessing","metadata":{}},{"cell_type":"code","source":"# Generate ECFPs\ndef generate_ecfp(molecule, radius=2, bits=1024):\n    if molecule is None:\n        return None\n    return list(AllChem.GetMorganFingerprintAsBitVect(molecule, radius, nBits=bits))\n\n# Generate topological FPs\ndef generate_topologicalfp(molecule, bits=1024):\n    if molecule is None:\n        return None\n    return list(AllChem.RDKFingerprint(molecule, fpSize=bits))\n\n# Generate Atom pair FPs\ndef generate_apfp(molecule, bits=1024):\n    if molecule is None:\n        return None\n    return list(AllChem.GetHashedAtomPairFingerprintAsBitVect(molecule, nBits=bits))\n\n# Generate all FPs\ndef generate_allfp(molecule, bits=1024):\n    if molecule is None:\n        return None\n    ecfp = list(AllChem.GetMorganFingerprintAsBitVect(molecule, 2, nBits=bits))\n    topologicalfp = list(AllChem.RDKFingerprint(molecule, fpSize=bits))\n    apfp = list(AllChem.GetHashedAtomPairFingerprintAsBitVect(molecule, nBits=bits))\n    return ecfp + topologicalfp + apfp","metadata":{"execution":{"iopub.status.busy":"2024-04-06T00:48:40.483927Z","iopub.execute_input":"2024-04-06T00:48:40.484414Z","iopub.status.idle":"2024-04-06T00:48:40.494981Z","shell.execute_reply.started":"2024-04-06T00:48:40.484380Z","shell.execute_reply":"2024-04-06T00:48:40.493834Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Train and Evaluate Model","metadata":{}},{"cell_type":"code","source":"# One-hot encode the protein_name\nonehot_encoder = OneHotEncoder(sparse_output=False)\nprotein_onehot = onehot_encoder.fit_transform(df['protein_name'].values.reshape(-1, 1))\n\n# Convert SMILES to RDKit molecules\ndf['molecule'] = df['molecule_smiles'].apply(Chem.MolFromSmiles)\n\ndef prepare_dataset(generate_fp):\n    # Combine ECFPs and one-hot encoded protein_name\n    _df = df.copy()\n    _df['fp'] = _df['molecule'].apply(generate_fp)\n    X = [fp + protein for fp, protein in zip(_df['fp'].tolist(), protein_onehot.tolist())]\n    y = _df['binds'].tolist()\n\n    # Split the data into train and test sets\n    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=test_size, random_state=random_state)\n    return X_train, X_test, y_train, y_test\n\ndef fit_clf(clf, X_train, y_train):\n    # Create and train the model\n    clf_copy = copy.deepcopy(clf)\n    clf_copy.fit(X_train, y_train)\n    return clf_copy\n\ndef evaluate_map(generate_fp, clf):\n    X_train, X_test, y_train, y_test = prepare_dataset(generate_fp)\n    clf_fitted = fit_clf(clf, X_train, y_train)\n    \n    # Make predictions on the test set\n    y_pred_proba = clf_fitted.predict_proba(X_test)[:, 1]  # Probability of the positive class\n\n    # Calculate the mean average precision\n    map_score = average_precision_score(y_test, y_pred_proba)\n    return map_score, clf_fitted\n","metadata":{"execution":{"iopub.status.busy":"2024-04-06T00:48:40.496576Z","iopub.execute_input":"2024-04-06T00:48:40.496877Z","iopub.status.idle":"2024-04-06T00:49:01.965602Z","shell.execute_reply.started":"2024-04-06T00:48:40.496852Z","shell.execute_reply":"2024-04-06T00:49:01.963951Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nprint(\"Calculating Mean Average Precision (mAP) for each FP...\")\nmap_score_ecfp, _ = evaluate_map(generate_fp=generate_ecfp,clf=clf)\nprint(f\"ECFP: {map_score_ecfp:.2f}\")\nmap_score_topologicalfp, _ = evaluate_map(generate_fp=generate_topologicalfp,clf=clf)\nprint(f\"Topological FP: {map_score_topologicalfp:.2f}\")\nmap_score_apfp, _ = evaluate_map(generate_fp=generate_apfp,clf=clf)\nprint(f\"Atom pair FP: {map_score_apfp:.2f}\")\nmap_score_allfp, clf_allfp = evaluate_map(generate_fp=generate_allfp,clf=clf)\nprint(f\"All FP: {map_score_allfp:.2f}\")","metadata":{"execution":{"iopub.status.busy":"2024-04-06T00:49:01.967610Z","iopub.execute_input":"2024-04-06T00:49:01.968044Z","iopub.status.idle":"2024-04-06T01:06:30.504621Z","shell.execute_reply.started":"2024-04-06T00:49:01.967999Z","shell.execute_reply":"2024-04-06T01:06:30.503073Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The best FP was ECFP. Feature selection or further treatments may be necessary to use topological or atom pair FPs more effectively. \n\nContinuing to predict test dataset with all FPs.","metadata":{}},{"cell_type":"code","source":"%%time\nimport os\n# Process the test.parquet file chunk by chunk\ntest_file = '/kaggle/input/leash-BELKA/test.csv'\noutput_file = 'submission.csv'  # Specify the path and filename for the output file\n\ngenerate_fp = generate_allfp\n# Read the test.parquet file into a pandas DataFrame\nfor df_test in pd.read_csv(test_file, chunksize=100000):\n\n    # Generate ECFPs for the molecule_smiles\n    df_test['molecule'] = df_test['molecule_smiles'].apply(Chem.MolFromSmiles)\n    df_test['fp'] = df_test['molecule'].apply(generate_fp)\n\n    # One-hot encode the protein_name\n    protein_onehot = onehot_encoder.transform(df_test['protein_name'].values.reshape(-1, 1))\n\n    # Combine ECFPs and one-hot encoded protein_name\n    X_test = [fp + protein for fp, protein in zip(df_test['fp'].tolist(), protein_onehot.tolist())]\n\n    # Predict the probabilities\n    probabilities = clf_allfp.predict_proba(X_test)[:, 1]\n\n    # Create a DataFrame with 'id' and 'probability' columns\n    output_df = pd.DataFrame({'id': df_test['id'], 'binds': probabilities})\n\n    # Save the output DataFrame to a CSV file\n    output_df.to_csv(output_file, index=False, mode='a', header=not os.path.exists(output_file))","metadata":{"execution":{"iopub.status.busy":"2024-04-06T01:24:29.934711Z","iopub.execute_input":"2024-04-06T01:24:29.936440Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}