{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.7.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":38128,"databundleVersionId":4230952,"sourceType":"competition"},{"sourceId":4127126,"sourceType":"datasetVersion","datasetId":2438837}],"dockerImageVersionId":30301,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Notebook to get important multiome features by predicting sex by RNA\n\nY - sex, X - RNA gene expression levels\n\nresults are stored in tables\n\n<..._weights.csv> - table with obtained coeffs after training\n\n<lasso_nonzero.csv> - table with non-zero genes by coeffs\n    \n<lasso_zero.csv> - table with zero genes by coeffs  \n\n<uni_variable_nonzero.csv> - table with univariable auc score for features with non-zero coefficients \n    \n#### Results:\n##### **The first ten largest weights affecting sex obtained for RNA**:\nENSG00000129824_RPS4Y1, ENSG00000229807_XIST, ENSG00000198034_RPS4X, ENSG00000198692_EIF1AY, ENSG00000230202_AL450405.1, ENSG00000164587_RPS14, ENSG00000198502_HLA-DRB5, ENSG00000251562_MALAT1, ENSG00000111640_GAPDH, ENSG00000197728_RPS26\n\n##### **The 10 weights NOT affecting sex obtained for RNA**:\nENSG00000121410_A1BG, ENSG00000229474_PATL2, ENSG00000166889_PATL1, ENSG00000132849_PATJ, ENSG00000115687_PASK, ENSG00000138964_PARVG, ENSG00000188677_PARVB, ENSG00000152931_PART1, ENSG00000162396_PARS2, ENSG00000185480_PARPBP<br> ","metadata":{"execution":{"iopub.status.busy":"2023-01-16T15:28:26.648165Z","iopub.execute_input":"2023-01-16T15:28:26.648692Z","iopub.status.idle":"2023-01-16T15:28:26.658298Z","shell.execute_reply.started":"2023-01-16T15:28:26.648648Z","shell.execute_reply":"2023-01-16T15:28:26.656397Z"}}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nfrom sklearn.linear_model import Lasso, LogisticRegression\nfrom sklearn.model_selection import cross_val_predict, cross_val_score\nfrom sklearn.metrics import r2_score, mean_squared_error, roc_auc_score\nimport matplotlib.colors as mcolors\nimport matplotlib.pyplot as plt\nimport lightgbm as lgb","metadata":{"execution":{"iopub.status.busy":"2023-12-04T12:55:14.053257Z","iopub.execute_input":"2023-12-04T12:55:14.055079Z","iopub.status.idle":"2023-12-04T12:55:16.737624Z","shell.execute_reply.started":"2023-12-04T12:55:14.054895Z","shell.execute_reply":"2023-12-04T12:55:16.736434Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"meta_df = pd.read_csv('/kaggle/input/open-problems-multimodal/metadata.csv', index_col = 0)\ndisplay(meta_df) ","metadata":{"execution":{"iopub.status.busy":"2023-12-04T12:55:16.740304Z","iopub.execute_input":"2023-12-04T12:55:16.741139Z","iopub.status.idle":"2023-12-04T12:55:17.369947Z","shell.execute_reply.started":"2023-12-04T12:55:16.741087Z","shell.execute_reply":"2023-12-04T12:55:17.368004Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X = pd.read_hdf('/kaggle/input/open-problems-multimodal/train_cite_inputs.h5') \n# cite_inputs = lgb.Dataset('/kaggle/input/open-problems-multimodal/test_cite_inputs.h5')\ndisplay(X)","metadata":{"execution":{"iopub.status.busy":"2023-12-04T12:55:17.372393Z","iopub.execute_input":"2023-12-04T12:55:17.373900Z","iopub.status.idle":"2023-12-04T12:56:20.135813Z","shell.execute_reply.started":"2023-12-04T12:55:17.373844Z","shell.execute_reply":"2023-12-04T12:56:20.134783Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_join = pd.merge(meta_df, \n                   X, \n                   on ='cell_id', \n                   how ='inner')","metadata":{"execution":{"iopub.status.busy":"2023-12-04T12:56:20.138438Z","iopub.execute_input":"2023-12-04T12:56:20.138899Z","iopub.status.idle":"2023-12-04T12:56:30.454065Z","shell.execute_reply.started":"2023-12-04T12:56:20.138860Z","shell.execute_reply":"2023-12-04T12:56:30.452844Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_join","metadata":{"execution":{"iopub.status.busy":"2023-12-04T12:56:30.455546Z","iopub.execute_input":"2023-12-04T12:56:30.455980Z","iopub.status.idle":"2023-12-04T12:56:30.547032Z","shell.execute_reply.started":"2023-12-04T12:56:30.455907Z","shell.execute_reply":"2023-12-04T12:56:30.545821Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"meta_df['sex'] = np.where(meta_df['donor'] > 30000, 'Male', 'Female')\ny = meta_df[meta_df.index.isin(X.index)].set_index(X.index)[['sex']]\ny = pd.get_dummies(y)\nprint(y, type(y))","metadata":{"execution":{"iopub.status.busy":"2023-12-04T12:56:30.548823Z","iopub.execute_input":"2023-12-04T12:56:30.549187Z","iopub.status.idle":"2023-12-04T12:56:30.708177Z","shell.execute_reply.started":"2023-12-04T12:56:30.549155Z","shell.execute_reply":"2023-12-04T12:56:30.706620Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rna_index = X.columns\nrna_index","metadata":{"execution":{"iopub.status.busy":"2023-12-04T12:56:30.710586Z","iopub.execute_input":"2023-12-04T12:56:30.711074Z","iopub.status.idle":"2023-12-04T12:56:30.720211Z","shell.execute_reply.started":"2023-12-04T12:56:30.711017Z","shell.execute_reply":"2023-12-04T12:56:30.718774Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y = y['sex_Male'].values\nprint(y.shape, type(y))","metadata":{"execution":{"iopub.status.busy":"2023-12-04T12:56:30.722503Z","iopub.execute_input":"2023-12-04T12:56:30.722915Z","iopub.status.idle":"2023-12-04T12:56:30.731198Z","shell.execute_reply.started":"2023-12-04T12:56:30.722880Z","shell.execute_reply":"2023-12-04T12:56:30.729588Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# X = X.astype('float64')\nX = X.values\nprint(X.shape, type(X))","metadata":{"execution":{"iopub.status.busy":"2023-12-04T12:56:30.733027Z","iopub.execute_input":"2023-12-04T12:56:30.733486Z","iopub.status.idle":"2023-12-04T12:56:30.747031Z","shell.execute_reply.started":"2023-12-04T12:56:30.733448Z","shell.execute_reply":"2023-12-04T12:56:30.745643Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfrom sklearn.preprocessing import StandardScaler\nscaler = StandardScaler()\nX = scaler.fit_transform(X)\nprint(X.shape)","metadata":{"execution":{"iopub.status.busy":"2023-12-04T12:56:30.753380Z","iopub.execute_input":"2023-12-04T12:56:30.753889Z","iopub.status.idle":"2023-12-04T12:56:59.100101Z","shell.execute_reply.started":"2023-12-04T12:56:30.753849Z","shell.execute_reply":"2023-12-04T12:56:59.098826Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import gc\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2023-12-04T12:56:59.104125Z","iopub.execute_input":"2023-12-04T12:56:59.104557Z","iopub.status.idle":"2023-12-04T12:56:59.247436Z","shell.execute_reply.started":"2023-12-04T12:56:59.104500Z","shell.execute_reply":"2023-12-04T12:56:59.245858Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Lasso","metadata":{}},{"cell_type":"code","source":"%%time\n\nimport time \n\ntrain_size = 0.7\nprint('X.shape, y.shape:', X.shape, y.shape, 'train_size:', train_size)\n\n\np = np.random.permutation(len(y))\nN = int(len(y) *  train_size)\nIX_train = np.arange(len(y))[p][:N]\nIX_test =  np.arange(len(y))[p][N:]\nprint(\" len(IX_train), len(IX_test):\",  len(IX_train), len(IX_test) )","metadata":{"execution":{"iopub.status.busy":"2023-12-04T12:56:59.249000Z","iopub.execute_input":"2023-12-04T12:56:59.249395Z","iopub.status.idle":"2023-12-04T12:56:59.265746Z","shell.execute_reply.started":"2023-12-04T12:56:59.249362Z","shell.execute_reply":"2023-12-04T12:56:59.264294Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_models_1 = pd.DataFrame()\n\n\nfor i, alpha in enumerate([1e-4, 1e-3, 2e-3, 5e-3, 1e-2, 1e-1, 1]):\n    t0 = time.time()\n\n    model = Lasso(alpha = alpha)\n\n    model.fit(X[IX_train,:], y[IX_train])\n\n    col = i\n    df_models_1.loc[col,'alpha'] = alpha\n    y_pred = model.predict(X[IX_train])\n    c = np.corrcoef(y[IX_train], y_pred)[0,1]\n    print(alpha, 'Corr Coef Train', c)\n    df_models_1.loc[col,'Corr Coef Train'] = c\n    auc = roc_auc_score(y[IX_train], y_pred)\n    print(alpha, 'ROC AUC Train', auc)\n    df_models_1.loc[col,'ROC AUC Train'] = auc\n\n    y_pred = model.predict(X[IX_test])\n    c = np.corrcoef(y[IX_test], y_pred)[0,1]\n    print(alpha, 'Corr Coef Test', c)\n    df_models_1.loc[col,'Corr Coef Test'] = c\n    auc = roc_auc_score(y[IX_test], y_pred)\n    print(alpha, 'ROC AUC Test', auc)\n    df_models_1.loc[col,'ROC AUC Test'] = auc\n\n    df_models_1.loc[col,'n_nonzeros'] = (model.coef_ != 0 ).sum()\n\n    print('%.1f secs passed'%(time.time()-t0))\n\n    alpha_selected = df_models_1.sort_values('ROC AUC Test', ascending = False)['alpha'].iat[0] \n    print('Best alpha: ', alpha_selected)    \n    display(df_models_1)  \n    \n    \nprint('Best alpha: ', alpha_selected)    ","metadata":{"execution":{"iopub.status.busy":"2023-12-04T12:56:59.267647Z","iopub.execute_input":"2023-12-04T12:56:59.268849Z","iopub.status.idle":"2023-12-04T13:25:50.387431Z","shell.execute_reply.started":"2023-12-04T12:56:59.268808Z","shell.execute_reply":"2023-12-04T13:25:50.384811Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model = Lasso(alpha = alpha_selected)\n\nmodel.fit(X[IX_train,:], y[IX_train])\n\ny_pred = model.predict(X[IX_train])\nc = np.corrcoef(y[IX_train], y_pred)[0,1]\nprint(alpha, 'Corr Coef Train', c)\nauc = roc_auc_score(y[IX_train], y_pred)\nprint(alpha, 'ROC AUC Train', auc)\n\ny_pred = model.predict(X[IX_test])\nc = np.corrcoef(y[IX_test], y_pred)[0,1]\nprint(alpha, 'Corr Coef Test', c)\nauc = roc_auc_score(y[IX_test], y_pred)\nprint(alpha, 'ROC AUC Test', auc)\nweights_dict = {}\nweights_dict['coef'] = model.coef_\nprint('non-zero coeffs:', (model.coef_ != 0 ).sum())","metadata":{"execution":{"iopub.status.busy":"2023-12-04T13:25:50.390944Z","iopub.execute_input":"2023-12-04T13:25:50.393302Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"nonzero_idx_ws = {}\nzero_idx_ws = {}\nzero_idx_ws['sex'] = np.where(weights_dict['coef'] == 0)[0]\nnonzero_idx_ws['sex'] = np.where(weights_dict['coef'] != 0)[0]\nprint(\"zero coef:\", len(zero_idx_ws['sex']), \"non-zero coef:\", len(nonzero_idx_ws['sex']))\n\nlasso_weights = pd.DataFrame(weights_dict, index = rna_index)\nlasso_weights.to_csv('/kaggle/working/lasso_weigths.csv')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"nonzero = lasso_weights['coef'].take(nonzero_idx_ws['sex']).sort_values(ascending = False)\nzero = lasso_weights['coef'].take(zero_idx_ws['sex'])\n\npd.DataFrame(zero).to_csv('/kaggle/working/lasso_zero.csv', header=['coef'], index = True)\npd.DataFrame(nonzero).to_csv('/kaggle/working/lasso_nonzero.csv', header=['coef'], index = True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"nonzero.index[:10].tolist()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"zero.index[:10].tolist()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"str_model_inf = 'Lasso {}'.format(alpha_selected)\nimportance = [abs(ele) for ele in nonzero]\nplt.bar([x for x in range(len(importance))], importance)\nplt.title('Importances '+ str_model_inf, fontsize = 20)\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Univariable for nonzero features","metadata":{}},{"cell_type":"code","source":"auc_scores = []\nfor idx, feature in zip(nonzero_idx_ws['sex'], nonzero.index):\n    # model = LogisticRegression()\n    model = Lasso(alpha = alpha_selected)\n    model.fit(X[IX_train, idx].reshape(-1, 1), y[IX_train])\n    y_pred = model.predict(X[IX_test, idx].reshape(-1, 1))\n    auc = roc_auc_score(y[IX_test], y_pred)\n    auc_scores.append(auc)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = pd.DataFrame({'idx': nonzero_idx_ws['sex'], 'feature': np.asarray(nonzero.index), 'lasso':np.asarray(lasso_weights['coef'].take(nonzero_idx_ws['sex'])), 'roc_auc uni': np.asarray(auc_scores)})\ndf = df.sort_values(by=['roc_auc uni'], ascending = False)\npd.DataFrame(df).to_csv('/kaggle/working/uni_variable_nonzero.csv', index = True)\ndf.head(50)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Plots","metadata":{}},{"cell_type":"code","source":"df_join.head(5)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_join['donor'].unique()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import seaborn as sns\nfor gene in nonzero.index[:20]:\n    sns.boxplot(x = df_join['donor'], y = df_join[gene])\n    plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Comparing with Logistic Regression","metadata":{}},{"cell_type":"code","source":"# https://www.kaggle.com/code/lnalinaf/sex-by-rna-lr-c100\n# https://www.kaggle.com/code/lnalinaf/sex-by-rna-lr-c10\n# https://www.kaggle.com/code/lnalinaf/sex-by-rna-lr-c1\n# https://www.kaggle.com/lnalinaf/sex-by-rna-lr-c0-01\n# https://www.kaggle.com/lnalinaf/sex-by-rna-lr-c0-1\n# https://www.kaggle.com/lnalinaf/sex-by-rna-lr-c0-001","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Comparing by coeficient of regularization (LR):  \n* C = 100: AUC Score:  1.0  \n          non-zero coeffs: 21601\n* C = 10: AUC Score:  0.9999657369971904  \n          non-zero coeffs: 21351 \n* C = 0.001: AUC Score:  0.9990456714383095  \n             non-zero coeffs: 4 \n* C = 0.01: AUC Score:  0.9997940126338918\n            non-zero coeffs: 23\n* C = 1: AUC Score:  0.9999658609859348\n         non-zero coeffs: 15329\n* C = 0.1: AUC Score:  0.999965858654831\n              non-zero coeffs: 500 \n* Lasso (alpha=0.001) ROC AUC Test: 0.9999997532343775\n        non-zero coeffs: 300","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}