{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# What is about ?\n\nSimple script to load CITE-seq part of the Kaggle competition dataset\non Single Cell multimodal data: \nhttps://www.kaggle.com/competitions/open-problems-multimodal\n\n\nAnd some simple prediction models.\n\n","metadata":{}},{"cell_type":"markdown","source":"# Preparations","metadata":{}},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport time\nt0start = time.time()\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor 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":{"execution":{"iopub.status.busy":"2022-12-24T18:36:37.257256Z","iopub.execute_input":"2022-12-24T18:36:37.257692Z","iopub.status.idle":"2022-12-24T18:36:37.269986Z","shell.execute_reply.started":"2022-12-24T18:36:37.257656Z","shell.execute_reply":"2022-12-24T18:36:37.268745Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#If you see a urllib warning running this cell, go to \"Settings\" on the right hand side, \n#and turn on internet. Note, you need to be phone verified.\n!pip install --quiet tables","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load data for targets -  140  CD (Cluster of differentiation) proteins \n\nhttps://en.wikipedia.org/wiki/Cluster_of_differentiation\n","metadata":{}},{"cell_type":"code","source":"%%time\ndf_y = pd.read_hdf('/kaggle/input/open-problems-multimodal/train_cite_targets.h5')\ndf_y","metadata":{"execution":{"iopub.status.busy":"2022-12-24T18:19:03.971090Z","iopub.execute_input":"2022-12-24T18:19:03.971500Z","iopub.status.idle":"2022-12-24T18:19:04.814124Z","shell.execute_reply.started":"2022-12-24T18:19:03.971467Z","shell.execute_reply":"2022-12-24T18:19:04.812505Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load features - RNA expression data \n","metadata":{}},{"cell_type":"code","source":"%%time\ndf_rna = pd.read_hdf('/kaggle/input/open-problems-multimodal/train_cite_inputs.h5')\ndisplay(df_rna) \n","metadata":{"execution":{"iopub.status.busy":"2022-12-24T18:21:05.361155Z","iopub.execute_input":"2022-12-24T18:21:05.361618Z","iopub.status.idle":"2022-12-24T18:21:56.914610Z","shell.execute_reply.started":"2022-12-24T18:21:05.361584Z","shell.execute_reply":"2022-12-24T18:21:56.913465Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Auxilliary meta data - not necessary for predictions","metadata":{}},{"cell_type":"code","source":"fn = '/kaggle/input/open-problems-multimodal/metadata.csv'\ndf_meta = pd.read_csv(fn, index_col = 0 )\ndf_meta","metadata":{"execution":{"iopub.status.busy":"2022-12-24T18:24:48.115252Z","iopub.execute_input":"2022-12-24T18:24:48.115761Z","iopub.status.idle":"2022-12-24T18:24:48.455332Z","shell.execute_reply.started":"2022-12-24T18:24:48.115715Z","shell.execute_reply":"2022-12-24T18:24:48.454118Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Keep meta data only relevant for train part of the CITE-seq data \ndf_meta = pd.DataFrame(index = df_y.index ).join(df_meta, how = 'left' )\ndf_meta","metadata":{"execution":{"iopub.status.busy":"2022-12-24T18:26:18.406526Z","iopub.execute_input":"2022-12-24T18:26:18.406941Z","iopub.status.idle":"2022-12-24T18:26:18.517812Z","shell.execute_reply.started":"2022-12-24T18:26:18.406909Z","shell.execute_reply":"2022-12-24T18:26:18.516679Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#  Example of simple modeling  ","metadata":{}},{"cell_type":"code","source":"%%time\nfrom sklearn.model_selection import train_test_split\n\nX = df_rna\ny = df_y['CD36']\nX_train, X_test, y_train, y_test = train_test_split(X, y, shuffle=True, train_size=0.1, random_state=0)\n","metadata":{"execution":{"iopub.status.busy":"2022-12-24T18:30:55.056190Z","iopub.execute_input":"2022-12-24T18:30:55.057001Z","iopub.status.idle":"2022-12-24T18:31:07.718168Z","shell.execute_reply.started":"2022-12-24T18:30:55.056965Z","shell.execute_reply":"2022-12-24T18:31:07.716984Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfrom sklearn.linear_model import Lasso\nfrom sklearn.metrics import r2_score\nfrom sklearn.metrics import mean_squared_error\nimport time\n\ndf_stat = pd.DataFrame()\nfor alpha in [0.005, 0.1, 1, 10, 100]:\n    t0 = time.time()\n    model = Lasso(alpha = alpha )\n    model.fit(X_train,y_train)\n    \n    y_pred = model.predict(X_test)\n    y_true = y_test\n    df_stat.loc['Corr Test', 'Lasso'+str(alpha)] = np.corrcoef( y_true, y_pred  )[0,1]\n    df_stat.loc['r2 Test', 'Lasso'+str(alpha)] = r2_score(y_true, y_pred)\n    df_stat.loc['RMSE Test', 'Lasso'+str(alpha)] = mean_squared_error(y_true,y_pred, squared=False) # squared=False -> RMSE, not MSE \n    \n    y_pred = model.predict(X_train)\n    y_true = y_train\n    df_stat.loc['Corr Train', 'Lasso'+str(alpha)] = np.corrcoef( y_true, y_pred  )[0,1]\n    df_stat.loc['r2 Train', 'Lasso'+str(alpha)] = r2_score(y_true, y_pred)\n    df_stat.loc['RMSE Train', 'Lasso'+str(alpha)] = mean_squared_error(y_true,y_pred, squared=False) # squared=False -> RMSE, not MSE \n    \n    df_stat.loc['n_important features', 'Lasso'+str(alpha)] = np.sum( model.coef_ != 0  )\n    df_stat.loc['Time', 'Lasso'+str(alpha)] = np.round( time.time() - t0,2 )\n    \n    IX = np.argsort( np.abs(model.coef_) )[::-1]\n    for i in range(12):\n        if np.abs(model.coef_[IX][i]) > 0.0000001:\n            df_stat.loc[str(i)+'-th Top feature', 'Lasso'+str(alpha)] = X.columns[IX][i]\n    for i in range(12):\n        df_stat.loc[str(i)+'-th Top coefficient', 'Lasso'+str(alpha)] = np.round( model.coef_[IX][i],4)\n    print('alpha = ',alpha, dict(df_stat['Lasso'+str(alpha)].iloc[:10] ))\n    print('Top 12 features:', X.columns[IX][:12])\n    print('Top 12 Coefs:', np.round( model.coef_[IX][:12],3) ) \n    print()\n        \n    \ndf_stat.head(50)","metadata":{"execution":{"iopub.status.busy":"2022-12-24T19:29:15.128900Z","iopub.execute_input":"2022-12-24T19:29:15.129382Z","iopub.status.idle":"2022-12-24T19:32:08.983537Z","shell.execute_reply.started":"2022-12-24T19:29:15.129348Z","shell.execute_reply":"2022-12-24T19:32:08.981935Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_stat.to_csv('Lasso_Models_predictions_statistics.csv')\ndf_stat","metadata":{"execution":{"iopub.status.busy":"2022-12-24T19:36:04.040863Z","iopub.execute_input":"2022-12-24T19:36:04.042198Z","iopub.status.idle":"2022-12-24T19:36:04.065260Z","shell.execute_reply.started":"2022-12-24T19:36:04.042143Z","shell.execute_reply":"2022-12-24T19:36:04.063914Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_stat_lasso = df_stat.copy()","metadata":{"execution":{"iopub.status.busy":"2022-12-24T19:36:00.813745Z","iopub.execute_input":"2022-12-24T19:36:00.814366Z","iopub.status.idle":"2022-12-24T19:36:00.821287Z","shell.execute_reply.started":"2022-12-24T19:36:00.814326Z","shell.execute_reply":"2022-12-24T19:36:00.820225Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Similar for Ridge","metadata":{}},{"cell_type":"code","source":"%%time\nfrom sklearn.linear_model import Lasso\nfrom sklearn.linear_model import Ridge\nfrom sklearn.metrics import r2_score\nfrom sklearn.metrics import mean_squared_error\nimport time\n\ndf_stat = pd.DataFrame()\nfor alpha in [ 1000,10_000,100_000,1_000_000]:\n    t0 = time.time()\n    model = Ridge(alpha = alpha )\n    model.fit(X_train,y_train)\n    \n    col = 'Ridge'+str(alpha)\n    \n    y_pred = model.predict(X_test)\n    y_true = y_test\n    df_stat.loc['Corr Test', col] = np.corrcoef( y_true, y_pred  )[0,1]\n    df_stat.loc['r2 Test', col] = r2_score(y_true, y_pred)\n    df_stat.loc['RMSE Test', col] = mean_squared_error(y_true,y_pred, squared=False) # squared=False -> RMSE, not MSE \n    \n    y_pred = model.predict(X_train)\n    y_true = y_train\n    df_stat.loc['Corr Train', col] = np.corrcoef( y_true, y_pred  )[0,1]\n    df_stat.loc['r2 Train', col] = r2_score(y_true, y_pred)\n    df_stat.loc['RMSE Train', col] = mean_squared_error(y_true,y_pred, squared=False) # squared=False -> RMSE, not MSE \n    \n    df_stat.loc['n_important features', col] = np.sum( model.coef_ != 0  )\n    df_stat.loc['Time', col] = np.round( time.time() - t0,2 )\n    \n    IX = np.argsort( np.abs(model.coef_) )[::-1]\n    for i in range(12):\n        if np.abs(model.coef_[IX][i]) > 0.0000001:\n            df_stat.loc[str(i)+'-th Top feature', col] = X.columns[IX][i]\n    for i in range(12):\n        df_stat.loc[str(i)+'-th Top coefficient', col] = np.round( model.coef_[IX][i],4)\n    print('alpha = ',alpha, dict(df_stat[col].iloc[:10] ))\n    print('Top 12 features:', X.columns[IX][:12])\n    print('Top 12 Coefs:', np.round( model.coef_[IX][:12],3) ) \n    print()\n        \n    \ndf_stat.head(50)","metadata":{"execution":{"iopub.status.busy":"2022-12-24T19:37:37.386026Z","iopub.execute_input":"2022-12-24T19:37:37.386526Z","iopub.status.idle":"2022-12-24T19:38:30.378495Z","shell.execute_reply.started":"2022-12-24T19:37:37.386490Z","shell.execute_reply":"2022-12-24T19:38:30.374587Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_stat.to_csv('Ridge_Models_predictions_statistics.csv')\ndf_stat","metadata":{"execution":{"iopub.status.busy":"2022-12-24T19:40:10.093999Z","iopub.execute_input":"2022-12-24T19:40:10.094511Z","iopub.status.idle":"2022-12-24T19:40:10.113069Z","shell.execute_reply.started":"2022-12-24T19:40:10.094475Z","shell.execute_reply":"2022-12-24T19:40:10.111826Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"t = time.time()-t0start\nprint('%.1f hours  = %.1f minutes = %.1f seconds passed total '%( t/3600,t/60, t ) )","metadata":{"execution":{"iopub.status.busy":"2022-12-23T20:02:41.854788Z","iopub.execute_input":"2022-12-23T20:02:41.856081Z","iopub.status.idle":"2022-12-23T20:02:41.863063Z","shell.execute_reply.started":"2022-12-23T20:02:41.85603Z","shell.execute_reply":"2022-12-23T20:02:41.861687Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}