{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":59094,"databundleVersionId":7010844,"sourceType":"competition"}],"dockerImageVersionId":30558,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport time\n\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport networkx as nx\n\n\n\nfrom sklearn.cluster import DBSCAN\nfrom sklearn.ensemble import RandomForestClassifier\nfrom sklearn.linear_model import Ridge\nfrom sklearn.manifold import TSNE\nfrom sklearn.decomposition import PCA, FastICA, TruncatedSVD\nfrom sklearn.model_selection import train_test_split, KFold\nfrom sklearn.metrics import r2_score, mean_squared_error\nfrom sklearn.preprocessing import LabelEncoder, OrdinalEncoder, OneHotEncoder\n\nimport os\nimport umap\nimport shap\nimport scipy.cluster.hierarchy as sch\nimport scipy.stats as stats\nfrom scipy.spatial.distance import squareform\n\n\nimport catboost\nfrom catboost import CatBoostRegressor, Pool\nfrom sklearn.multioutput import MultiOutputRegressor","metadata":{"execution":{"iopub.status.busy":"2023-11-23T09:54:12.652801Z","iopub.execute_input":"2023-11-23T09:54:12.653375Z","iopub.status.idle":"2023-11-23T09:54:12.664757Z","shell.execute_reply.started":"2023-11-23T09:54:12.653321Z","shell.execute_reply":"2023-11-23T09:54:12.663199Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-11-23T09:54:12.667618Z","iopub.execute_input":"2023-11-23T09:54:12.668135Z","iopub.status.idle":"2023-11-23T09:54:12.703904Z","shell.execute_reply.started":"2023-11-23T09:54:12.668082Z","shell.execute_reply":"2023-11-23T09:54:12.702860Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_data = pd.read_parquet(\"/kaggle/input/open-problems-single-cell-perturbations/de_train.parquet\")\nid_map = pd.read_csv(\"/kaggle/input/open-problems-single-cell-perturbations/id_map.csv\")\nsample_submission = pd.read_csv(\"/kaggle/input/open-problems-single-cell-perturbations/sample_submission.csv\")","metadata":{"execution":{"iopub.status.busy":"2023-11-23T09:54:12.705266Z","iopub.execute_input":"2023-11-23T09:54:12.705887Z","iopub.status.idle":"2023-11-23T09:54:18.657981Z","shell.execute_reply.started":"2023-11-23T09:54:12.705854Z","shell.execute_reply":"2023-11-23T09:54:18.656727Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fn = '/kaggle/input/open-problems-single-cell-perturbations/de_train.parquet'\ndf_de_train = pd.read_parquet(fn)# , index_col = 0)\nprint(df_de_train.shape)\ndf_de_train","metadata":{"execution":{"iopub.status.busy":"2023-11-23T09:54:18.659288Z","iopub.execute_input":"2023-11-23T09:54:18.659845Z","iopub.status.idle":"2023-11-23T09:54:20.145242Z","shell.execute_reply.started":"2023-11-23T09:54:18.659802Z","shell.execute_reply":"2023-11-23T09:54:20.144027Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"id_map.head()","metadata":{"execution":{"iopub.status.busy":"2023-11-23T09:54:20.149240Z","iopub.execute_input":"2023-11-23T09:54:20.149746Z","iopub.status.idle":"2023-11-23T09:54:20.160822Z","shell.execute_reply.started":"2023-11-23T09:54:20.149712Z","shell.execute_reply":"2023-11-23T09:54:20.159812Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_data.head(5)","metadata":{"execution":{"iopub.status.busy":"2023-11-23T09:54:20.162528Z","iopub.execute_input":"2023-11-23T09:54:20.162907Z","iopub.status.idle":"2023-11-23T09:54:20.195705Z","shell.execute_reply.started":"2023-11-23T09:54:20.162878Z","shell.execute_reply":"2023-11-23T09:54:20.194804Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_data.info()","metadata":{"execution":{"iopub.status.busy":"2023-11-23T09:54:20.196904Z","iopub.execute_input":"2023-11-23T09:54:20.197826Z","iopub.status.idle":"2023-11-23T09:54:21.845982Z","shell.execute_reply.started":"2023-11-23T09:54:20.197793Z","shell.execute_reply":"2023-11-23T09:54:21.844790Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Adjust the sample size to the number of rows in your DataFrame if needed\nsample_size = min(1000, len(train_data))  # Ensuring sample size is not greater than the dataset size\n\n# Perform t-SNE\ntsne = TSNE(n_components=2, perplexity=30, random_state=42)\ntsne_results = tsne.fit_transform(train_data.drop(['cell_type', 'sm_name', 'sm_lincs_id', 'SMILES', 'control'], axis=1).sample(sample_size))\n\n# Visualization\ntrain_data_subset = train_data.sample(sample_size)  # Ensure you're sampling the same rows for coloring by cell_type\n\nplt.figure(figsize=(16,10))\nsns.scatterplot(\n    x=tsne_results[:,0], y=tsne_results[:,1],\n    hue=train_data_subset['cell_type'],\n    palette=sns.color_palette(\"hsv\", len(train_data_subset['cell_type'].unique())),\n    legend=\"full\",\n    alpha=0.7\n)\nplt.title('t-SNE Visualization of Gene Expressions')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-11-23T09:54:21.847785Z","iopub.execute_input":"2023-11-23T09:54:21.848153Z","iopub.status.idle":"2023-11-23T09:54:27.167655Z","shell.execute_reply.started":"2023-11-23T09:54:21.848123Z","shell.execute_reply":"2023-11-23T09:54:27.166411Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# UMAP for dimensionality reduction\nreducer = umap.UMAP(random_state=42)\nembedding = reducer.fit_transform(train_data.drop(['cell_type', 'sm_name', 'sm_lincs_id', 'SMILES', 'control'], axis=1))\n\n# Using DBSCAN for clustering UMAP output\nclustering = DBSCAN(eps=0.5, min_samples=10).fit(embedding)\ntrain_data['umap_cluster'] = clustering.labels_\n\n# Visualizing UMAP output with clusters\nplt.figure(figsize=(12, 10))\nsns.scatterplot(\n    x=embedding[:, 0], \n    y=embedding[:, 1], \n    hue=train_data['umap_cluster'], \n    palette=sns.color_palette(\"hsv\", len(set(clustering.labels_))),\n    legend='full'\n)\nplt.title('UMAP projection of the Dataset', fontsize=20)\nplt.xlabel('UMAP Dimension 1')\nplt.ylabel('UMAP Dimension 2')\nplt.show()\n\n\n# Prepare data for SHAP analysis\nX = train_data.drop(['cell_type', 'sm_name', 'sm_lincs_id', 'SMILES', 'control', 'umap_cluster'], axis=1)\ny = train_data['umap_cluster']\n\n# Train a model for SHAP analysis\nmodel = RandomForestClassifier(random_state=0, n_estimators=100)\nmodel.fit(X, y)\n\n# Explain the model's predictions using SHAP\nexplainer = shap.TreeExplainer(model)\nshap_values = explainer.shap_values(X.sample(100))  # Using a sample for speed\n\n# Visualize the first prediction's explanation\nshap.initjs()\nshap.force_plot(explainer.expected_value[0], shap_values[0][0], X.iloc[0])","metadata":{"execution":{"iopub.status.busy":"2023-11-23T09:54:27.169209Z","iopub.execute_input":"2023-11-23T09:54:27.169647Z","iopub.status.idle":"2023-11-23T09:54:40.676495Z","shell.execute_reply.started":"2023-11-23T09:54:27.169607Z","shell.execute_reply":"2023-11-23T09:54:40.675190Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for gene in ['A1BG', 'A1BG-AS1', 'A2M']:  # Replace with genes of interest\n    result = stats.kruskal(*[group[gene].values for name, group in train_data.groupby('cell_type')])\n    print(f\"Kruskal-Wallis test for {gene}: H-statistic={result.statistic}, p-value={result.pvalue}\")","metadata":{"execution":{"iopub.status.busy":"2023-11-23T09:54:40.678472Z","iopub.execute_input":"2023-11-23T09:54:40.679148Z","iopub.status.idle":"2023-11-23T09:54:40.822534Z","shell.execute_reply.started":"2023-11-23T09:54:40.679106Z","shell.execute_reply":"2023-11-23T09:54:40.821300Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"le = LabelEncoder()\ntrain_data['cell_type_encoded'] = le.fit_transform(train_data['cell_type'])\n\n# Define features and target\nX = train_data.drop(['cell_type', 'cell_type_encoded', 'sm_name', 'sm_lincs_id', 'SMILES', 'control'], axis=1)\ny = train_data['cell_type_encoded']\n\n# Split the data into training and test sets\nX_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)\n\n# Fit the RandomForest model on the training data\nrf = RandomForestClassifier(n_estimators=100, random_state=42)\nrf.fit(X_train, y_train)\n\n# Extract feature importances from the model\nimportances = rf.feature_importances_\n\n# Plot the top N feature importances\ntop_n = 20  # Specify how many top features you'd like to visualize\nindices = np.argsort(importances)[-top_n:]\nplt.figure(figsize=(15, 7))\nplt.title('Top Feature Importances')\nplt.barh(range(len(indices)), importances[indices], color='b', align='center')\nplt.yticks(range(len(indices)), [X_train.columns[i] for i in indices])\nplt.xlabel('Relative Importance')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-11-23T09:54:40.824203Z","iopub.execute_input":"2023-11-23T09:54:40.824555Z","iopub.status.idle":"2023-11-23T09:54:44.529776Z","shell.execute_reply.started":"2023-11-23T09:54:40.824526Z","shell.execute_reply":"2023-11-23T09:54:44.528789Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Model Training","metadata":{}},{"cell_type":"code","source":"\nplt.figure(figsize=(10, 6))\nplt.hist(train_data['A1BG'], bins=50, color='blue', alpha=0.7)\nplt.title(\"Distribution of Gene Expression (A1BG)\")\nplt.xlabel(\"Gene Expression Value\")\nplt.ylabel(\"Frequency\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-11-23T09:54:44.531208Z","iopub.execute_input":"2023-11-23T09:54:44.531621Z","iopub.status.idle":"2023-11-23T09:54:44.905685Z","shell.execute_reply.started":"2023-11-23T09:54:44.531587Z","shell.execute_reply":"2023-11-23T09:54:44.904730Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Training Data Shape:\", train_data.shape)\nprint(\"Id Map Data Shape:\", id_map.shape)","metadata":{"execution":{"iopub.status.busy":"2023-11-23T09:54:44.907144Z","iopub.execute_input":"2023-11-23T09:54:44.907826Z","iopub.status.idle":"2023-11-23T09:54:44.913733Z","shell.execute_reply.started":"2023-11-23T09:54:44.907792Z","shell.execute_reply":"2023-11-23T09:54:44.912563Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time \nn_components = 35\n\npredict_method = 'train_aggregation_by_compounds_with_denoising_TSVD'\n\nif '_pca' in predict_method:\n    str_inf_target_dimred = 'PCA' \n    reducer = PCA(n_components=n_components )\nelif '_ICA' in predict_method:\n    str_inf_target_dimred = 'ICA' \n    reducer = FastICA(n_components=n_components, random_state=0, whiten='unit-variance')\nelif '_TSVD' in predict_method:\n    str_inf_target_dimred = 'TSVD' \n    reducer = TruncatedSVD(n_components=n_components, n_iter=7, random_state=42)\nelse:\n    str_inf_target_dimred = ''\n    \nprint(str_inf_target_dimred, reducer)","metadata":{"execution":{"iopub.status.busy":"2023-11-23T09:54:44.919012Z","iopub.execute_input":"2023-11-23T09:54:44.919447Z","iopub.status.idle":"2023-11-23T09:54:44.933375Z","shell.execute_reply.started":"2023-11-23T09:54:44.919413Z","shell.execute_reply":"2023-11-23T09:54:44.932105Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"Y = df_de_train.iloc[:,5:].values\nYr = reducer.fit_transform(Y)\npd.DataFrame(Yr).corr(method = 'spearman').round(2)","metadata":{"execution":{"iopub.status.busy":"2023-11-23T09:54:44.935330Z","iopub.execute_input":"2023-11-23T09:54:44.935799Z","iopub.status.idle":"2023-11-23T09:54:47.001944Z","shell.execute_reply.started":"2023-11-23T09:54:44.935758Z","shell.execute_reply":"2023-11-23T09:54:46.999229Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_de_train.iloc[:,5:].shape","metadata":{"execution":{"iopub.status.busy":"2023-11-23T09:54:47.004429Z","iopub.execute_input":"2023-11-23T09:54:47.004869Z","iopub.status.idle":"2023-11-23T09:54:47.041768Z","shell.execute_reply.started":"2023-11-23T09:54:47.004833Z","shell.execute_reply":"2023-11-23T09:54:47.040362Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nn_components_for_compound_encoding = 3\nquantile = 0.54\ndf_tmp = pd.DataFrame(Yr[:, :n_components_for_compound_encoding  ], index = df_de_train.index  )\ndf_tmp['column for aggregation'] = df_de_train['sm_name']\ndf_compound_encoded = df_tmp.groupby('column for aggregation').quantile( quantile )\nprint('df_compound_encoded.shape', df_compound_encoded.shape )\ndisplay( df_compound_encoded )\n\nX = pd.DataFrame(index = df_de_train['sm_name']) # \nX['IX'] = np.arange(len(df_de_train))\nX = X.join( df_compound_encoded , how = 'left').sort_values('IX')\nX = X.iloc[:,1:]\nX = X.values \nX.shape","metadata":{"execution":{"iopub.status.busy":"2023-11-23T09:54:47.043383Z","iopub.execute_input":"2023-11-23T09:54:47.043714Z","iopub.status.idle":"2023-11-23T09:54:47.076419Z","shell.execute_reply.started":"2023-11-23T09:54:47.043687Z","shell.execute_reply":"2023-11-23T09:54:47.075214Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"enc = OrdinalEncoder()\nX = enc.fit_transform(df_de_train[ ['cell_type', 'sm_name']])\nX.shape","metadata":{"execution":{"iopub.status.busy":"2023-11-23T09:54:47.077825Z","iopub.execute_input":"2023-11-23T09:54:47.078267Z","iopub.status.idle":"2023-11-23T09:54:47.091801Z","shell.execute_reply.started":"2023-11-23T09:54:47.078228Z","shell.execute_reply":"2023-11-23T09:54:47.090432Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"n_components_for_compound_encoding = 20\nquantile = 0.54\ndf_tmp = pd.DataFrame(Yr[:, :n_components_for_compound_encoding  ], index = df_de_train.index  )\ndf_tmp['column for aggregation'] = df_de_train['sm_name']\ndf_compound_encoded = df_tmp.groupby('column for aggregation').mean() # quantile( quantile )\nprint('df_compound_encoded.shape', df_compound_encoded.shape )\n\nX = pd.DataFrame(index = df_de_train['sm_name']) # \nX['IX'] = np.arange(len(df_de_train))\nX = X.join( df_compound_encoded , how = 'left').sort_values('IX')\nX = X.iloc[:,1:]\nX = X.values \nprint( X.shape )\n\nX_submit = pd.DataFrame(index = id_map['sm_name']) # \nX_submit['IX'] = np.arange(len(X_submit))\nX_submit = X_submit.join( df_compound_encoded , how = 'left').sort_values('IX')\nX_submit = X_submit.iloc[:,1:]\nX_submit = X_submit.values \nprint( X_submit.shape )","metadata":{"execution":{"iopub.status.busy":"2023-11-23T09:54:47.093793Z","iopub.execute_input":"2023-11-23T09:54:47.094247Z","iopub.status.idle":"2023-11-23T09:54:47.119710Z","shell.execute_reply.started":"2023-11-23T09:54:47.094206Z","shell.execute_reply":"2023-11-23T09:54:47.117993Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"n_splits = 3\nkf = KFold(n_splits=n_splits, random_state = 42, shuffle = True )\nIX = pd.DataFrame(); IX['IX'] = range(len(df_de_train))\n\nalpha_regularization_for_linear_models = 100000\nmodel = Ridge(alpha=alpha_regularization_for_linear_models)\n\nY = df_de_train.iloc[:,5:].values\n\nY_reduced_submit = np.zeros( (len(id_map) , Yr.shape[1] )   ); cnt_blend_submit = 0\nY_submit = np.zeros( (len(id_map) , 18211 )   ); cnt_blend_submit = 0\nY_oof = np.zeros( (len(df_de_train) , 18211 )   ); cnt_blend_submit = 0","metadata":{"execution":{"iopub.status.busy":"2023-11-23T09:54:47.121246Z","iopub.execute_input":"2023-11-23T09:54:47.121741Z","iopub.status.idle":"2023-11-23T09:54:47.172973Z","shell.execute_reply.started":"2023-11-23T09:54:47.121700Z","shell.execute_reply":"2023-11-23T09:54:47.171658Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i, (train_index, test_index) in enumerate(kf.split(IX)):\n    print(f\"Fold {i}:\", len(test_index) )\n    \n    Yr_train = reducer.fit_transform(Y[train_index,:])\n    Yr_test = reducer.transform(Y[test_index,:])\n    X_train = X[train_index,:]\n    X_test = X[test_index,:]\n    model.fit(X_train, Yr_train)\n    Yr_pred = model.predict(X_test) # , Yr_train)\n    r2 = r2_score(Yr_test,  Yr_pred )\n    print('r2 test:', r2)\n    r2_test = r2\n    Y_oof[test_index,: ] = reducer.inverse_transform(Yr_pred)\n    \n    Yr_pred = model.predict(X_train) # , Yr_train)\n    r2 = r2_score(Yr_train,  Yr_pred )\n    print('r2 train', r2)\n    Y_reduced_submit =  (Y_reduced_submit * cnt_blend_submit + model.predict(X_submit) ) / ( cnt_blend_submit + 1)\n    Y_submit =  (Y_submit * cnt_blend_submit + reducer.inverse_transform(model.predict(X_submit) ) ) / ( cnt_blend_submit + 1)\n    cnt_blend_submit += 1\n    \ns = mean_squared_error(Y,Y_oof, squared = False)\nprint(s)","metadata":{"execution":{"iopub.status.busy":"2023-11-23T09:54:47.174428Z","iopub.execute_input":"2023-11-23T09:54:47.175142Z","iopub.status.idle":"2023-11-23T09:54:52.422833Z","shell.execute_reply.started":"2023-11-23T09:54:47.175109Z","shell.execute_reply":"2023-11-23T09:54:52.422006Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time \nn_splits = 3\nkf = KFold(n_splits=n_splits, random_state = 42, shuffle = True )\nIX = pd.DataFrame(); IX['IX'] = range(len(df_de_train))\n\nalpha_regularization_for_linear_models = 100000\nmodel = Ridge(alpha=alpha_regularization_for_linear_models)\n\ncategorical_features = ['cell_type','sm_name']\nmodel = MultiOutputRegressor( CatBoostRegressor(cat_features=categorical_features, verbose = 0,  # Categorical features ) ) \n                        iterations=300,  # Number of boosting iterations\n                          depth=5,        # Depth of the trees\n                          learning_rate=0.015,  # Learning rate\n                          loss_function='RMSE'))  # Specify your loss function (e.g., RMSE for regression)\n#                           verbose=0) ) # Set verbose to 0 to suppress output\n\n\nY = df_de_train.iloc[:,5:].values","metadata":{"execution":{"iopub.status.busy":"2023-11-23T09:54:52.423877Z","iopub.execute_input":"2023-11-23T09:54:52.424876Z","iopub.status.idle":"2023-11-23T09:54:52.466184Z","shell.execute_reply.started":"2023-11-23T09:54:52.424844Z","shell.execute_reply":"2023-11-23T09:54:52.465045Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_submit = pd.DataFrame(Y_submit, columns = df_de_train.columns[5:])\ndf_submit.index.name = 'id'\nprint( df_submit.shape )\ndisplay(df_submit)\ndf_submit.to_csv('submission.csv')","metadata":{"execution":{"iopub.status.busy":"2023-11-23T09:54:52.467811Z","iopub.execute_input":"2023-11-23T09:54:52.468298Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}