{"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":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"\n### goal of this competition is to predict how small molecules change gene expression in different cell types\n# Single Cell Experiment R Example\n - https://colab.research.google.com/drive/1rWZz9UHW2ZVWieejz9N3EjlgLr-w7xJs?usp=sharing\n# What is about ?  \n\nThe first brief look on the \"Open Problems – Single-Cell Perturbations\" data.\n\nAnd some baselines\n\n\n###  Notes \n\n- number of train samples - 614, test 255 - not much ! ML for such low number of samples - hmmm.... very unclearly how successful it can be ...   \n\n- 18 211 targets (genes)  \n\n- in some sense we have ONLY TWO features - cell type and compound. Both of them are categorical. Cell types - 6, compounds - 146 categories. Thus it is natural to use target encodings. Alternative way - one might think in direction of recsys approach, but with the modification that we have 18211 targets. \n\n\n    V40 EDA : Genes top affected by Clotrimazole, Vorinostat :\n    Clotrimazole\n    Up: ['RPL10', 'RPL13', 'MAPRE3', 'TPT1', 'TRAV6', 'RPL3P4', 'AC107068.1', 'TMEM191C', 'SYNGR1', 'MTATP6P1', 'RPL7A', 'CHD3', 'GAL3ST4', 'SRPK1', 'HIST2H2BE', 'CCM2', 'RPLP1', 'PRKY', 'TOB2', 'SLC25A42']\n    Down: ['MX1', 'MX2', 'MRPS18B', 'RAB18', 'USP18', 'ISG20', 'ISG15', 'C19ORF66', 'SHISA5', 'IFI16', 'IFI44', 'AC090114.2', 'C6ORF62', 'TDRD7', 'LAP3', 'IFI44L', 'TRIM22', 'XAF1', 'AC116407.2', 'HMGN2']\n\n    Vorinostat\n    Up: ['HIST1H1D', 'H1FX', 'METRN', 'HIST1H1C', 'RPS18', 'HIST1H1B', 'STMN1', 'RPLP0', 'RPL7A', 'FXYD7', 'H2AFJ', 'RPS16', 'RPL8', 'HIST1H1E', 'CTNNAL1', 'RPS7', 'MARCKSL1', 'FXYD1', 'PCSK1N', 'RPL13A']\n    Down: ['CORO1A', 'DDX5', 'SRSF5', 'LIMD2', 'CD3E', 'EMP3', 'GPSM3', 'S100A6', 'RGL4', 'LTB', 'AES', 'HCST', 'ADI1', 'TRBC2', 'SRSF7', 'BIN2', 'UFC1', 'GZMM', 'CD37', 'RBM3']\n\n\n    V39 EDA: Genes top affected by 'Topotecan': \n        UP: ['RPS27L', 'HNRNPC', 'VAMP2', 'TMEM14B', 'TTC39C', 'AP3M2', 'RSRP1', 'CXCR4', 'BAX', 'HINT1', 'FTH1', 'HMGB1', 'ITM2A', 'NDUFS5', 'RPL12', 'CREM', 'COX20', 'NSD1', 'RPL28', 'RPLP1']\n        Down: ['ARHGAP15', 'ZBTB20', 'STX8', 'IMMP2L', 'SYNE2', 'PTPRC', 'SMAP1', 'SND1', 'MRPS6', 'DOCK2', 'PRKCB', 'SKAP1', 'ACTB', 'ITPR2', 'TBC1D5', 'LPP', 'COMMD1', 'FAF1', 'SSU72', 'NDUFAF2']\n        \n        \n    V37 0.707 approach 5 - aggr by compound + Ridge\n    \n    V34 LB 0.688 iterations=3,  depth = 6  \n    V33 LB 0.709 Catboost iterations 10,  depth = 6    - try to reduce overfit - but even worse  \n    V32 LB 0.912 - Catboost depth = 2, iteration 100 - try to reduce overfit - but even worse  \n    V30 LB 0.905 ( terrible :)  - Catboost first draft -  depth = 6 , iteration 100\n\n    Versions 23-27 - first modeling approach - target encode cell type and compound and simple Ridge model\n    LB 0.616:  https://www.kaggle.com/code/alexandervc/op2-models-cv-tuning?scriptVersionId=143324657 (V4)\n    It was neccesary to drasticaly increase the regalarization - alpha = 100_000 to get better results. \n    Without that - result were worse than naive approach used before:\n        V27 LB 0.659 LB Ridge nCT1 nCD25 Al10 TSVD35 - even more relaxed : alpha and encoded_compound size \n        V26 LB 0.668 try to relax a model a bit: Ridge nCT1 nCD10 Al100 TSVD35 - results are better near mean/zero submission\n        V25 LB 0.677 - stronger constraining the model: Ridge nCT1 nCD5 Al100 TSVD35, but still we are worse than even predict by mean \n        V24 LB 0.702 - try avoid overfit Ridge nCT3 nCD10 Al100 TSVD35 - better results but still bad \n        V23 LB 0.747 Ridge model in target encoded features Ridge nCT10 nCD35 Al1 TSVD35 - modeling gives worse results than simple approach, might be bug or overfit\n    \n    V22 EDA - cell cycle genes brief analysis \n    V19,20,21 EDA - clustering 10000, 15000, 18211(all) genes by sns.clustermap - 15min,47min,RAMcrash - see two clear clusters in genes - what are they ? \n    V18 LB 0.623 TSVD-35 denoising  quantile 0.54, tsvd is better than pca/ica - similar to previous OP\n    V15 - bug -  LB 0.626 - NO it was not: TSVD-35 denoising  quantile 0.54 - so worse than pca,ICA, at least for these params\n    V14 LB 0.624 - ICA-35 denoising   quantile 0.54 - so ICA is not better than just pca at least for these the same params\n    Fork: LB 0.624 - pca-35 denoising, quantile 0.54 \n    V13 LB 0.626 pca25 denoising \n    V12 LB 0.627 same with pca100 denoising\n            That means: first reduce data to pca100, and  take pca inverse transfrom (denoising)\n            and then apply same as in V9: groupby by compound and quantile(0.6)         \n            so idea is - pca100 reductions hopefully kills some noise\n            In the previous challenge it worked okay, but the number of samples was not 614, but near 100 000\n    V11 - bug - pay attention that column names in submit file should correspond to genes\n    V9 LB 0.638 groupby by compound and quantile(0.6)\n    V6 LB 0.666 quantile(0.7)\n    V5 LB 0.657 quantile(0.6) \n            results are again better, that probably indicates some shift between train and public data \n    V4 LB 0.659 - median instead of mean - results are a bit better, \n            it might mean either a bit of presense of outliers, or  public is somewhat different from train - next experiments suggests second is true \n    V1 LB 0.664 - submission of train means - the simplest baseline ","metadata":{}},{"cell_type":"markdown","source":"# Preliminaries","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() \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\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":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-11-28T11:49:18.750544Z","iopub.execute_input":"2023-11-28T11:49:18.751154Z","iopub.status.idle":"2023-11-28T11:49:21.322742Z","shell.execute_reply.started":"2023-11-28T11:49:18.751088Z","shell.execute_reply":"2023-11-28T11:49:21.321251Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Train data","metadata":{}},{"cell_type":"code","source":"%%time\nfn = '/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-28T11:56:27.453073Z","iopub.execute_input":"2023-11-28T11:56:27.453609Z","iopub.status.idle":"2023-11-28T11:56:30.513748Z","shell.execute_reply.started":"2023-11-28T11:56:27.453574Z","shell.execute_reply":"2023-11-28T11:56:30.512364Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_de_train.iloc[:,5:].head(1)","metadata":{"execution":{"iopub.status.busy":"2023-11-28T11:56:37.873296Z","iopub.execute_input":"2023-11-28T11:56:37.873786Z","iopub.status.idle":"2023-11-28T11:56:37.953409Z","shell.execute_reply.started":"2023-11-28T11:56:37.873749Z","shell.execute_reply":"2023-11-28T11:56:37.951776Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Genes top affected by some compounds (drugs)","metadata":{}},{"cell_type":"code","source":"print( list(df_de_train['sm_name'].unique()) )","metadata":{"execution":{"iopub.status.busy":"2023-11-28T11:56:44.722315Z","iopub.execute_input":"2023-11-28T11:56:44.722774Z","iopub.status.idle":"2023-11-28T11:56:44.732881Z","shell.execute_reply.started":"2023-11-28T11:56:44.722739Z","shell.execute_reply":"2023-11-28T11:56:44.731136Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize = (20,4) )\nv = df_de_train.iloc[:,5:].max(axis = 0 ).sort_values(ascending = False, key = abs )\nplt.plot(v.values,'*-')\nplt.title('Max-abs DE for genes',fontsize = 20 )\nplt.grid()\nplt.show()\n# display(v.head(15))\nt = v\nn =0;d=15;\nfor n in [0,d,2*d,3*d,4*d,5*d,6*d,7*d,8*d,9*d]:\n    t2 = t.iloc[n:n+d].to_frame().reset_index()\n    t2.columns = ['Gene','DE']\n    if n == 0:\n        t3 = t2.copy()\n    else:\n        t3 = pd.concat( (t3,t2),axis = 1 )\n        \ndisplay(t3.round(1))\n","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:02:25.201274Z","iopub.execute_input":"2023-11-28T12:02:25.201843Z","iopub.status.idle":"2023-11-28T12:02:25.748619Z","shell.execute_reply.started":"2023-11-28T12:02:25.201798Z","shell.execute_reply":"2023-11-28T12:02:25.747124Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nprint(' \"Most often\" (compound/cell-types averaged) affected genes ')\nt =  (-(df_de_train.iloc[:,5:].T.abs())).rank().median(axis=1) # Use ranks and average (over pairs - compound/cell-types) them . Minus - to get top gene as rank=1 not as rank=len   \nt = t.sort_values()\nprint(' \"Most often\" (compound/cell-types averaged) affected genes ')\nn =0;d=15;\nfor n in [0,d,2*d,3*d,4*d,5*d,6*d,7*d,8*d,9*d]:\n    t2 = t.iloc[n:n+d].to_frame().reset_index()\n    t2.columns = ['Gene','Rank']\n    if n == 0:\n        t3 = t2.copy()\n    else:\n        t3 = pd.concat( (t3,t2),axis = 1 )\ndisplay(t3)\n\nplt.figure(figsize = (20,4) )\nplt.plot(t.values,'*')\nplt.ylabel('Rank (averaged) ',fontsize = 15)\nplt.grid()\nplt.title( '  \"Most often\" (compound/cell-types averaged) affected genes ', fontsize = 20)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:02:32.382738Z","iopub.execute_input":"2023-11-28T12:02:32.383479Z","iopub.status.idle":"2023-11-28T12:02:36.299759Z","shell.execute_reply.started":"2023-11-28T12:02:32.383433Z","shell.execute_reply":"2023-11-28T12:02:36.298525Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nprint(' Top affected genes at least once ')\nt =  (-(df_de_train.iloc[:,5:].T.abs())).rank().min(axis=1) # Use ranks and average (over pairs - compound/cell-types) them . Minus - to get top gene as rank=1 not as rank=len   \nt = t.sort_values()\nn =0;d=15;\nfor n in [0,d,2*d,3*d,4*d,5*d,6*d,7*d,8*d,9*d]:\n    t2 = t.iloc[n:n+d].to_frame().reset_index()\n    t2.columns = ['Gene','Rank']\n    if n == 0:\n        t3 = t2.copy()\n    else:\n        t3 = pd.concat( (t3,t2),axis = 1 )\n        \nprint('Count top1 at least once:',  (t==1).sum() )\ndisplay(t3)\n\nplt.figure(figsize = (20,4) )\nplt.plot(t.values,'*')\nplt.ylabel('Rank (averaged) ',fontsize = 15)\nplt.grid()\nplt.title( '  Top affected genes at least once ', fontsize = 20)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:02:47.251623Z","iopub.execute_input":"2023-11-28T12:02:47.252305Z","iopub.status.idle":"2023-11-28T12:02:49.850864Z","shell.execute_reply.started":"2023-11-28T12:02:47.252259Z","shell.execute_reply":"2023-11-28T12:02:49.849396Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"t =  (-(df_de_train.iloc[:,5:].T.abs())).rank()\nm = (t == 1)\nv = (m.sum(axis = 1))\nv = v.sort_values(ascending = False)\nfor t in list(range(1,14))[::-1]:\n    print(t, (v == t ).sum() )\nprint(' Genes Most often top1 ')\nt = v\nn =0;d=15;\nfor n in [0,d,2*d,3*d,4*d,5*d,6*d,7*d,8*d,9*d]:\n    t2 = t.iloc[n:n+d].to_frame().reset_index()\n    t2.columns = ['Gene','Rank']\n    if n == 0:\n        t3 = t2.copy()\n    else:\n        t3 = pd.concat( (t3,t2),axis = 1 )\ndisplay(t3)\n\n\nt =  (-(df_de_train.iloc[:,5:].T.abs())).rank()\nm = (t == 1)| (t == 2)\nv = (m.sum(axis = 1))\nv = v.sort_values(ascending = False)\nfor t in list(range(1,14))[::-1]:\n    print(t, (v == t ).sum() )\nprint(' Genes Most often top1 or top2 ')\nt = v\nn =0;d=15;\nfor n in [0,d,2*d,3*d,4*d,5*d,6*d,7*d,8*d,9*d]:\n    t2 = t.iloc[n:n+d].to_frame().reset_index()\n    t2.columns = ['Gene','Rank']\n    if n == 0:\n        t3 = t2.copy()\n    else:\n        t3 = pd.concat( (t3,t2),axis = 1 )\ndisplay(t3)\n\nt =  (-(df_de_train.iloc[:,5:].T.abs())).rank()\nm = (t == 1)| (t == 2) | (t == 3)| (t == 4)  | (t == 5)| (t == 6) \nv = (m.sum(axis = 1))\nv = v.sort_values(ascending = False)\nfor t in list(range(1,14))[::-1]:\n    print(t, (v == t ).sum() )\nprint(' Genes Most often top1 - top5 ')\nt = v\nn =0;d=15;\nfor n in [0,d,2*d,3*d,4*d,5*d,6*d,7*d,8*d,9*d]:\n    t2 = t.iloc[n:n+d].to_frame().reset_index()\n    t2.columns = ['Gene','Rank']\n    if n == 0:\n        t3 = t2.copy()\n    else:\n        t3 = pd.concat( (t3,t2),axis = 1 )\ndisplay(t3)\n","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:02:53.572800Z","iopub.execute_input":"2023-11-28T12:02:53.573635Z","iopub.status.idle":"2023-11-28T12:03:00.475120Z","shell.execute_reply.started":"2023-11-28T12:02:53.573597Z","shell.execute_reply":"2023-11-28T12:03:00.474276Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"drug = 'Topotecan'\nv_drug = (df_de_train['sm_name'] == drug).astype(float)\nv_drug.sum()\ndf = df_de_train.copy()\ndf[drug] = v_drug\n#for gene in df_de_train.columns[5:]:\nd_agg = df.iloc[:,5:].groupby(drug).mean().T\nd_agg = d_agg.sort_values(1)\ndisplay(d_agg.sort_values(1,ascending = False).head(10))\nprint(list(d_agg.sort_values(1,ascending = False).head(20).index))\ndisplay(d_agg.sort_values(1,ascending = True).head(10))\nprint(list(d_agg.sort_values(1,ascending = True).head(20).index))\n# display(d_agg.head(5) )\n# display(d_agg.tail(5) )\nm = d_agg[1] > 4\ndisplay(d_agg[m])\n\nplt.figure(figsize=(20,5))\nplt.plot(d_agg[0].values)\nplt.plot(d_agg[1].values)\nplt.grid()\nplt.title(drug,fontsize=20)\nplt.xlabel('Genes sorted',fontsize = 20 )\nplt.show()\n\nplt.figure(figsize=(20,5))\nd_agg[3] = d_agg[1].values - d_agg[0].values \nplt.plot(d_agg[3].values)\nplt.plot()\nplt.grid()\nplt.title(drug,fontsize=20)\nplt.xlabel('Genes sorted',fontsize = 20 )\nplt.show()\ndisplay( d_agg.sort_values(3).head(20) )\ndisplay( d_agg.sort_values(3).tail(20) )\n","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:03:00.476837Z","iopub.execute_input":"2023-11-28T12:03:00.477772Z","iopub.status.idle":"2023-11-28T12:03:01.637631Z","shell.execute_reply.started":"2023-11-28T12:03:00.477737Z","shell.execute_reply":"2023-11-28T12:03:01.636284Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"drug = 'Vorinostat' # 'Topotecan' 'Belinostat', 'Clotrimazole',\nv_drug = (df_de_train['sm_name'] == drug).astype(float)\nv_drug.sum()\ndf = df_de_train.copy()\ndf[drug] = v_drug\n#for gene in df_de_train.columns[5:]:\nd_agg = df.iloc[:,5:].groupby(drug).mean().T\nd_agg = d_agg.sort_values(1)\ndisplay(d_agg.sort_values(1,ascending = False).head(10))\nprint(list(d_agg.sort_values(1,ascending = False).head(20).index))\ndisplay(d_agg.sort_values(1,ascending = True).head(10))\nprint(list(d_agg.sort_values(1,ascending = True).head(20).index))\n# display(d_agg.head(5) )\n# display(d_agg.tail(5) )\nm = d_agg[1] > 4\ndisplay(d_agg[m])\n\nplt.figure(figsize=(20,5))\nplt.plot(d_agg[0].values)\nplt.plot(d_agg[1].values)\nplt.grid()\nplt.title(drug,fontsize=20)\nplt.xlabel('Genes sorted',fontsize = 20 )\nplt.show()\n\nplt.figure(figsize=(20,5))\nd_agg[3] = d_agg[1].values - d_agg[0].values \nplt.plot(d_agg[3].values)\nplt.plot()\nplt.grid()\nplt.title(drug,fontsize=20)\nplt.xlabel('Genes sorted',fontsize = 20 )\nplt.show()\ndisplay( d_agg.sort_values(3).head(20) )\ndisplay( d_agg.sort_values(3).tail(20) )\n","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:03:10.652055Z","iopub.execute_input":"2023-11-28T12:03:10.653417Z","iopub.status.idle":"2023-11-28T12:03:11.642519Z","shell.execute_reply.started":"2023-11-28T12:03:10.653360Z","shell.execute_reply":"2023-11-28T12:03:11.641154Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"drug = 'Clotrimazole'#  'Vorinostat' # 'Topotecan' 'Belinostat', ,\nv_drug = (df_de_train['sm_name'] == drug).astype(float)\nv_drug.sum()\ndf = df_de_train.copy()\ndf[drug] = v_drug\n#for gene in df_de_train.columns[5:]:\nd_agg = df.iloc[:,5:].groupby(drug).mean().T\nd_agg = d_agg.sort_values(1)\ndisplay(d_agg.sort_values(1,ascending = False).head(10))\nprint(list(d_agg.sort_values(1,ascending = False).head(20).index))\ndisplay(d_agg.sort_values(1,ascending = True).head(10))\nprint(list(d_agg.sort_values(1,ascending = True).head(20).index))\n# display(d_agg.head(5) )\n# display(d_agg.tail(5) )\nm = d_agg[1] > 4\ndisplay(d_agg[m])\n\nplt.figure(figsize=(20,5))\nplt.plot(d_agg[0].values)\nplt.plot(d_agg[1].values)\nplt.grid()\nplt.title(drug,fontsize=20)\nplt.xlabel('Genes sorted',fontsize = 20 )\nplt.show()\n\nplt.figure(figsize=(20,5))\nd_agg[3] = d_agg[1].values - d_agg[0].values \nplt.plot(d_agg[3].values)\nplt.plot()\nplt.grid()\nplt.title(drug,fontsize=20)\nplt.xlabel('Genes sorted',fontsize = 20 )\nplt.show()\ndisplay( d_agg.sort_values(3).head(20) )\ndisplay( d_agg.sort_values(3).tail(20) )\n","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:03:24.101109Z","iopub.execute_input":"2023-11-28T12:03:24.101743Z","iopub.status.idle":"2023-11-28T12:03:25.571109Z","shell.execute_reply.started":"2023-11-28T12:03:24.101699Z","shell.execute_reply":"2023-11-28T12:03:25.569529Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Dimensional reductions (pca, umap,...), visualizations, clustering train target data ","metadata":{}},{"cell_type":"code","source":"X = df_de_train.iloc[:,5:]\nprint(X.shape)","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:03:33.821682Z","iopub.execute_input":"2023-11-28T12:03:33.823292Z","iopub.status.idle":"2023-11-28T12:03:33.883603Z","shell.execute_reply.started":"2023-11-28T12:03:33.823226Z","shell.execute_reply":"2023-11-28T12:03:33.882188Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfrom sklearn.decomposition import PCA\n\nv1_color = df_de_train[  'cell_type']\nv2_color = df_de_train[  'sm_name'].copy()\nv3_color = df_de_train[  'sm_name'].copy()\nl = [t for t in df_de_train[  'sm_name'] if t.endswith('nib') ]\nm = v2_color.isin( l)\nv2_color[~m] = 'non -nib'\nv3_color[m] = '*nib'\n\nv4_color = df_de_train[  'control']#.copy()\n\nlist_top_drugs = ['MLN 2238', 'Resminostat', 'CEP-18770 (Delanzomib)', 'Oprozomib (ONX 0912)', 'Belinostat', 'Vorinostat', 'Ganetespib (STA-9090)', 'Scriptaid', 'Proscillaridin A;Proscillaridin-A', 'Alvocidib', 'IN1451']\nm = df_de_train[  'sm_name'].isin( list_top_drugs)\nv5_color = df_de_train[  'sm_name'].copy()\nv5_color[~m] = np.nan\n\n\n\nlist_cfg = [ ['cell type',v1_color], ['control' , v4_color ] , ['top compounds',v5_color] ]\n#     str_inf1 = ''\n    #X = np.clip(df.iloc[N0:N1,33:137].fillna(0),0, 1)\nstr_inf = 'PCA' \nreducer = PCA(n_components=100 )\nXr = reducer.fit_transform(X)\nfor i,j in [[0,1],[0,2],[1,2],[3,4],[5,6],[7,8]]:\n    plt.figure(figsize = (20,10)); ic=0\n    for str_inf1, v_for_color in list_cfg: # , ['*nib compounds ',v2_color], ['non *nib compounds',v3_color ] ]:\n        ic+=1; plt.subplot(1,len(list_cfg),ic)\n        sns.scatterplot(x= Xr[:,i], y = Xr[:,j], hue =  v_for_color ,s = 100) # df['reads'])\n        plt.xlabel(str_inf+str(i+1), fontsize = 20)\n        plt.ylabel(str_inf+str(j+1), fontsize = 20)\n        plt.title(str_inf1 + ' ', fontsize = 20 )\n\n    plt.show()\n    \n\n","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:03:37.082397Z","iopub.execute_input":"2023-11-28T12:03:37.083066Z","iopub.status.idle":"2023-11-28T12:03:49.748436Z","shell.execute_reply.started":"2023-11-28T12:03:37.083022Z","shell.execute_reply":"2023-11-28T12:03:49.747162Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"d = df_de_train.iloc[:,:5]\nd['PCA1'] = Xr[:,0]\nd['PCA2'] = Xr[:,1]\nd['PCA3'] = Xr[:,2]\nlist_top_drugs = []\ndisplay( d.sort_values('PCA1', ascending = False ).head(8) )\nlist_top_drugs += d.sort_values('PCA1', ascending = False ).head(8)['sm_name'].to_list()\nprint(list_top_drugs)\ndisplay( d.sort_values('PCA2', ascending = False ).head(8) )\nlist_top_drugs += d.sort_values('PCA2', ascending = False ).head(8)['sm_name'].to_list()\ndisplay( d.sort_values('PCA3', ascending = False ).head(8) )\nlist_top_drugs += d.sort_values('PCA3', ascending = False ).head(8)['sm_name'].to_list()\nprint(list(set(list_top_drugs)))","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:04:05.522615Z","iopub.execute_input":"2023-11-28T12:04:05.523189Z","iopub.status.idle":"2023-11-28T12:04:05.594241Z","shell.execute_reply.started":"2023-11-28T12:04:05.523147Z","shell.execute_reply":"2023-11-28T12:04:05.592587Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Clustering  Cell Types","metadata":{}},{"cell_type":"code","source":"%%time\nN = df_de_train.shape[1]# 5000\nprint(N)\nX = df_de_train[ ['cell_type'] + list(df_de_train.columns[5:N]) ].groupby('cell_type').median()\nprint(X.shape)\ncm = np.corrcoef(X)\nprint(cm[:3,:2])\ncm = np.abs(cm)\nl = list(X.index)# [df_de_train['sm_name'].iat[i] +' '+ df_de_train['cell_type'].iat[i]  for i in range(len(df_de_train))] # .columns[5:N]\ncm = pd.DataFrame(cm, index =l , columns = l )\nprint(cm.shape)\nsns.clustermap(cm,  annot=True, fmt=\".2f\", cmap=\"coolwarm\" )\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:04:21.571116Z","iopub.execute_input":"2023-11-28T12:04:21.571583Z","iopub.status.idle":"2023-11-28T12:04:22.657737Z","shell.execute_reply.started":"2023-11-28T12:04:21.571547Z","shell.execute_reply":"2023-11-28T12:04:22.656343Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Clustering compounds","metadata":{}},{"cell_type":"code","source":"%%time\nN = df_de_train.shape[1]# 5000\nprint(N)\nX = df_de_train[ ['sm_name'] + list(df_de_train.columns[5:N]) ].groupby('sm_name').median()\nprint(X.shape)\ncm = np.corrcoef(X)\nprint(cm[:3,:2])\ncm = np.abs(cm)\nl = list(X.index)# [df_de_train['sm_name'].iat[i] +' '+ df_de_train['cell_type'].iat[i]  for i in range(len(df_de_train))] # .columns[5:N]\nl = [t[:20] for t in l] # cut long names\ncm = pd.DataFrame(cm, index =l , columns = l )\nprint(cm.shape)\nclustergrid = sns.clustermap(cm,cmap=\"coolwarm\" )# ,  annot=True, fmt=\".2f\", \nplt.show()\nreordered_columns = clustergrid.dendrogram_col.reordered_ind\nreordered_rows = clustergrid.dendrogram_row.reordered_ind\nprint(len(reordered_rows), len(reordered_columns) )\nprint( list(cm.index[reordered_rows]) )\n# print( list(X.columns[reordered_columns]) )\n\nsns.clustermap(cm,  annot=True, fmt=\".2f\", cmap=\"coolwarm\" )\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:04:27.540355Z","iopub.execute_input":"2023-11-28T12:04:27.540886Z","iopub.status.idle":"2023-11-28T12:05:22.907272Z","shell.execute_reply.started":"2023-11-28T12:04:27.540846Z","shell.execute_reply":"2023-11-28T12:05:22.905480Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Clustering samples (i.e. pairs cell + compound)","metadata":{}},{"cell_type":"code","source":"%%time\nN = df_de_train.shape[1]# 5000\nX = df_de_train.iloc[:,5:N]\nprint(X.shape)\ncm = np.corrcoef(X)\nprint(cm[:3,:2])\ncm = np.abs(cm)\nl = [df_de_train['sm_name'].iat[i] +' '+ df_de_train['cell_type'].iat[i]  for i in range(len(df_de_train))] # .columns[5:N]\ncm = pd.DataFrame(cm, index =l , columns = l )\nprint(cm.shape)\nclustergrid = sns.clustermap(cm)\nplt.show()\nreordered_columns = clustergrid.dendrogram_col.reordered_ind\nreordered_rows = clustergrid.dendrogram_row.reordered_ind\nprint(len(reordered_rows), len(reordered_columns) )\nprint( list(cm.index[reordered_rows]) )\n# print( list(cm.columns[reordered_columns]) )\n","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:06:24.501114Z","iopub.execute_input":"2023-11-28T12:06:24.501988Z","iopub.status.idle":"2023-11-28T12:06:27.119006Z","shell.execute_reply.started":"2023-11-28T12:06:24.501932Z","shell.execute_reply":"2023-11-28T12:06:27.117424Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Look in genes space ","metadata":{}},{"cell_type":"code","source":"%%time\nfrom sklearn.decomposition import PCA\n\nX = df_de_train.iloc[:,5:].T\nprint(X.shape)\n\nv1_color = pd.Series(range(len(X)), name = 'index') # df_de_train[  'cell_type']\n\n\n\nlist_cfg = [ ['Genes',v1_color]]# , ['control' , v4_color ] , ['top compounds',v5_color] ]\n#     str_inf1 = ''\n    #X = np.clip(df.iloc[N0:N1,33:137].fillna(0),0, 1)\nstr_inf = 'PCA' \nreducer = PCA(n_components=10 )\nXr = reducer.fit_transform(X)\nfor i,j in [[0,1],[0,2],[1,2],[3,4],[5,6],[7,8]]:\n    plt.figure(figsize = (20,10)); ic=0\n    for str_inf1, v_for_color in list_cfg: # , ['*nib compounds ',v2_color], ['non *nib compounds',v3_color ] ]:\n        ic+=1; plt.subplot(1,len(list_cfg),ic)\n        sns.scatterplot(x= Xr[:,i], y = Xr[:,j], hue =  v_for_color ,s = 100) # df['reads'])\n        plt.xlabel(str_inf+str(i+1), fontsize = 20)\n        plt.ylabel(str_inf+str(j+1), fontsize = 20)\n        plt.title(str_inf1 + ' ', fontsize = 20 )\n\n    plt.show()\n    \n\n","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:06:37.851847Z","iopub.execute_input":"2023-11-28T12:06:37.852402Z","iopub.status.idle":"2023-11-28T12:06:50.857409Z","shell.execute_reply.started":"2023-11-28T12:06:37.852360Z","shell.execute_reply":"2023-11-28T12:06:50.855997Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfrom sklearn.decomposition import PCA\nimport umap \n\nX = df_de_train.iloc[:,5:].T\nprint(X.shape)\n\nv1_color = pd.Series(range(len(X)), name = 'index') # df_de_train[  'cell_type']\n\n\n\nlist_cfg = [ ['Genes',v1_color]]# , ['control' , v4_color ] , ['top compounds',v5_color] ]\n# str_inf = 'PCA' \n# reducer = PCA(n_components=10 )\nstr_inf = 'UMAP' \nreducer = umap.UMAP()# (n_components=10 )\n\nXr = reducer.fit_transform(X)\nfor i,j in [[0,1] ]:# ,[0,2],[1,2],[3,4],[5,6],[7,8]]:\n    plt.figure(figsize = (20,10)); ic=0\n    for str_inf1, v_for_color in list_cfg: # , ['*nib compounds ',v2_color], ['non *nib compounds',v3_color ] ]:\n        ic+=1; plt.subplot(1,len(list_cfg),ic)\n        sns.scatterplot(x= Xr[:,i], y = Xr[:,j], hue =  v_for_color ,s = 100) # df['reads'])\n        plt.xlabel(str_inf+str(i+1), fontsize = 20)\n        plt.ylabel(str_inf+str(j+1), fontsize = 20)\n        plt.title(str_inf1 + ' ', fontsize = 20 )\n\n    plt.show()\n    \n\n","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:07:04.520641Z","iopub.execute_input":"2023-11-28T12:07:04.521157Z","iopub.status.idle":"2023-11-28T12:08:28.703029Z","shell.execute_reply.started":"2023-11-28T12:07:04.521110Z","shell.execute_reply":"2023-11-28T12:08:28.701873Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Genes clustering ","metadata":{}},{"cell_type":"code","source":"%%time\nN = 1000#  18211 #  15_000# 10000 #\nX = df_de_train.iloc[:,5:N].T\nprint(X.shape)\ncm = np.corrcoef(X)\nprint(cm[:3,:2])\ncm = np.abs(cm)\ncm = pd.DataFrame(cm, index = df_de_train.columns[5:N], columns = df_de_train.columns[5:N] )\nprint(cm.shape)\nsns.clustermap(cm)\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:15:47.361829Z","iopub.execute_input":"2023-11-28T12:15:47.362266Z","iopub.status.idle":"2023-11-28T12:15:51.323324Z","shell.execute_reply.started":"2023-11-28T12:15:47.362233Z","shell.execute_reply":"2023-11-28T12:15:51.322025Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Genes related to proliferation (cell cycle genes )\n\nProliferation cycle or cell cycle - key process for many cells. Many drugs here are anti-cancer - so should affect the ability of cells to proliferate. So might be interesting to look on drugs affect on these particular genes.\n\nTirosh et.al. proposed list of about 100 cell-cycle genes which are the most effectively seen by single cell data. It is good starting point.\n\nIn general there are much more genes related to the cell cycle, with various degree of \"relation\". Some list are cell cycle genes may contain thousands genes, but in fact in such huge lists most of the genes are related to cell cycle very weakly or these genes not good captured by single cell sequencing, while Tirosh list contains genes very strongly related to cell cycle and well captured by single cell technology. Moreover, many genes are associated to various biological processes in the cell, but Tirosh genes are mostly associated with cell cycle, not with the other processes - another argument why they are convenient to work with. \n\nOne may look at Computational challenges of cell cycle analysis using single cell transcriptomics\nAlexander Chervov, Andrei Zinovyev https://arxiv.org/abs/2208.05229\n\nAnd hundreds Kaggle notebooks/datasets related to that work e.g. that discussion: https://www.kaggle.com/competitions/open-problems-multimodal/discussion/350314 , or that notebook: https://www.kaggle.com/code/alexandervc/mmscel-cell-cycle-03b-daybydaychange-allcelltypes/notebook\n\nPS\n\nOther cell cycle genes sets e.g. by Tom Freeman\n\nSee some comparison e.g. here: https://www.kaggle.com/code/alexandervc/tabmuris-cell-cycle-1-data-one-by-one#Cell-cycle\nSo list is bigger but kind of more \"dirty\", that means containing more non-cell cycle effects, and G1S - G2M split is less prominent. \n","metadata":{}},{"cell_type":"code","source":"G1S_genes_Tirosh = ['MCM5', 'PCNA', 'TYMS', 'FEN1', 'MCM2', 'MCM4', 'RRM1', 'UNG', 'GINS2', 'MCM6', 'CDCA7', 'DTL', 'PRIM1', 'UHRF1', 'MLF1IP', 'HELLS', 'RFC2', 'RPA2', 'NASP', 'RAD51AP1', 'GMNN', 'WDR76', 'SLBP', 'CCNE2', 'UBR7', 'POLD3', 'MSH2', 'ATAD2', 'RAD51', 'RRM2', 'CDC45', 'CDC6', 'EXO1', 'TIPIN', 'DSCC1', 'BLM', 'CASP8AP2', 'USP1', 'CLSPN', 'POLA1', 'CHAF1B', 'BRIP1', 'E2F8']\nG2M_genes_Tirosh = ['HMGB2', 'CDK1', 'NUSAP1', 'UBE2C', 'BIRC5', 'TPX2', 'TOP2A', 'NDC80', 'CKS2', 'NUF2', 'CKS1B', 'MKI67', 'TMPO', 'CENPF', 'TACC3', 'FAM64A', 'SMC4', 'CCNB2', 'CKAP2L', 'CKAP2', 'AURKB', 'BUB1', 'KIF11', 'ANP32E', 'TUBB4B', 'GTSE1', 'KIF20B', 'HJURP', 'CDCA3', 'HN1', 'CDC20', 'TTK', 'CDC25C', 'KIF2C', 'RANGAP1', 'NCAPD2', 'DLGAP5', 'CDCA2', 'CDCA8', 'ECT2', 'KIF23', 'HMMR', 'AURKA', 'PSRC1', 'ANLN', 'LBR', 'CKAP5', 'CENPE', 'CTCF', 'NEK2', 'G2E3', 'GAS2L3', 'CBX5', 'CENPA']\ngenes_Tirosh = G1S_genes_Tirosh + G2M_genes_Tirosh\n\n# Subset of Tirosh genes to capture \"fast\" cell cycle pattern = see https://arxiv.org/abs/2208.05229\nlist_genes_fastCCsign = ['CDK1', 'UBE2C', 'TOP2A', 'TMPO', 'HJURP', 'RRM1', 'RAD51AP1', 'RRM2', 'CDC45', 'BLM', 'BRIP1', 'E2F8', 'HIST2H2AC']\n\nG1S_genes_Freeman = ['ADAMTS1', 'ASF1B', 'ATAD2', 'BARD1', 'BLM', 'BRCA1', 'BRIP1', 'C17orf75', 'C9orf40', 'CACYBP', 'CASP8AP2', 'CCDC15', 'CCNE1', 'CCNE2', 'CCP110', 'CDC25A', 'CDC45', 'CDC6', 'CDC7', 'CDK2', 'CDT1', 'CENPJ', 'CENPQ', 'CENPU', 'CEP57', 'CHAF1A', 'CHAF1B', 'CHEK1', 'CLSPN', 'CREBZF', 'CRYL1', 'CSE1L', 'DCLRE1B', 'DCTPP1', 'DEK', 'DERA', 'DHFR', 'DNA2', 'DNAJC9', 'DNMT1', 'DONSON', 'DSCC1', 'DSN1', 'DTL', 'E2F8', 'EED', 'EFCAB11', 'ENDOD1', 'ETAA1', 'EXO1', 'EYA2', 'EZH2', 'FAM111A', 'FANCE', 'FANCG', 'FANCI', 'FANCL', 'FBXO5', 'FEN1', 'GGH', 'GINS1', 'GINS2', 'GINS3', 'GLMN', 'GMNN', 'GMPS', 'GPD2', 'HADH', 'HELLS', 'HSF2', 'ITGB3BP', 'KIAA0101', 'KNTC1', 'LIG1', 'MCM10', 'MCM2', 'MCM3', 'MCM4', 'MCM5', 'MCM6', 'MCM7', 'MCMBP', 'METTL9', 'MMD', 'MNS1', 'MPP1', 'MRE11A', 'MSH2', 'MSH6', 'MYO19', 'NASP', 'NPAT', 'NSMCE4A', 'ORC1', 'OSGEPL1', 'PAK1', 'PAQR4', 'PARP2', 'PASK', 'PAXIP1', 'PBX3', 'PCNA', 'PKMYT1', 'PMS1', 'POLA1', 'POLA2', 'POLD3', 'POLE2', 'PRIM1', 'PRPS2', 'PSMC3IP', 'RAB23', 'RAD51', 'RAD51AP1', 'RAD54L', 'RBBP8', 'RBL1', 'RDX', 'RFC2', 'RFC3', 'RFC4', 'RMI1', 'RNASEH2A', 'RPA1', 'RRM1', 'RRM2', 'SLBP', 'SLC25A40', 'SMC2', 'SMC3', 'SSX2IP', 'SUPT16H', 'TEX30', 'TFDP1', 'THAP10', 'THEM6', 'TIMELESS', 'TIPIN', 'TMEM106C', 'TMEM38B', 'TRIM45', 'TRIP13', 'TSPYL4', 'TTI1', 'TUBGCP5', 'TYMS', 'UBR7', 'UNG', 'USP1', 'WDHD1', 'WDR76', 'WRB', 'YEATS4', 'ZBTB14', 'ZWINT']\nG2M_genes_Freeman = ['ADGRE5', 'ARHGAP11A', 'ARHGDIB', 'ARL6IP1', 'ASPM', 'AURKA', 'AURKB', 'BIRC5', 'BORA', 'BRD8', 'BUB1', 'BUB1B', 'BUB3', 'CCNA2', 'CCNB1', 'CCNB2', 'CCNF', 'CDC20', 'CDC25B', 'CDC25C', 'CDC27', 'CDCA3', 'CDCA8', 'CDK1', 'CDKN1B', 'CDKN3', 'CENPE', 'CENPF', 'CENPI', 'CENPN', 'CEP55', 'CEP70', 'CEP85', 'CKAP2', 'CKAP5', 'CKS1B', 'CKS2', 'CTCF', 'DBF4', 'DBF4B', 'DCAF7', 'DEPDC1', 'DLGAP5', 'ECT2', 'ERCC6L', 'ESPL1', 'FAM64A', 'FOXM1', 'FZD2', 'FZD7', 'FZR1', 'GPSM2', 'GTF2E1', 'GTSE1', 'H2AFX', 'HJURP', 'HMGB2', 'HMGB3', 'HMMR', 'HN1', 'INCENP', 'JADE2', 'KIF11', 'KIF14', 'KIF15', 'KIF18A', 'KIF18B', 'KIF20A', 'KIF20B', 'KIF22', 'KIF23', 'KIF2C', 'KIF4A', 'KIF5B', 'KIFC1', 'KPNA2', 'LBR', 'LMNB2', 'MAD2L1', 'MELK', 'MET', 'METTL4', 'MIS18BP1', 'MKI67', 'MPHOSPH9', 'MTMR6', 'NCAPD2', 'NCAPG', 'NCAPG2', 'NCAPH', 'NDC1', 'NDC80', 'NDE1', 'NEIL3', 'NEK2', 'NRF1', 'NUSAP1', 'OIP5', 'PAFAH2', 'PARPBP', 'PBK', 'PLEKHG3', 'PLK1', 'PLK4', 'PRC1', 'PRR11', 'PSRC1', 'PTTG1', 'PTTG3P', 'RACGAP1', 'RAD21', 'RASSF1', 'REEP4', 'SAP30', 'SHCBP1', 'SKA1', 'SLCO1B3', 'SOGA1', 'SPA17', 'SPAG5', 'SPC25', 'SPDL1', 'STIL', 'STK17B', 'TACC3', 'TAF5', 'TBC1D2', 'TBC1D31', 'TMPO', 'TOP2A', 'TPX2', 'TROAP', 'TTF2', 'TTK', 'TUBB4B', 'TUBD1', 'UBE2C', 'UBE2S', 'VANGL1', 'WEE1', 'WHSC1', 'XPO1', 'ZMYM1']\n\n\n\n","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:16:57.984987Z","iopub.execute_input":"2023-11-28T12:16:57.985609Z","iopub.status.idle":"2023-11-28T12:16:58.013959Z","shell.execute_reply.started":"2023-11-28T12:16:57.985569Z","shell.execute_reply":"2023-11-28T12:16:58.012197Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nN = 1000#  18211 #  15_000# 10000 #\n\nfor l,str_inf in [ [G1S_genes_Tirosh, 'G1S Tirosh'], [G2M_genes_Tirosh, 'G2M Tirosh'],  [G1S_genes_Tirosh + G2M_genes_Tirosh, 'All Tirosh'],\n                  [list_genes_fastCCsign, 'FastCC Signature'],\n                  [ G1S_genes_Freeman, 'G1S Freeman' ],  [ G2M_genes_Freeman, 'G2M Freeman' ], [ G1S_genes_Freeman + G2M_genes_Freeman, 'All Freeman' ],   ]: \n    ll = set(l) & set(df_de_train.columns) \n    ll = list(ll)\n    print(len(ll), str_inf )\n    X = df_de_train[ll].T # .iloc[:,5:N].T\n    print(X.shape)\n    cm = np.corrcoef(X)\n    print(cm[:3,:2])\n    cm = np.abs(cm)\n    cm = pd.DataFrame(cm, index = ll, columns = ll )\n    print(cm.shape)\n    clustergrid = sns.clustermap(cm)\n    plt.title(str_inf, fontsize = 20 )\n    plt.show()\n    reordered_columns = clustergrid.dendrogram_col.reordered_ind\n    reordered_rows = clustergrid.dendrogram_row.reordered_ind\n    print(len(reordered_rows), len(reordered_columns) )\n    print( list(cm.index[reordered_rows]) )\n#     print( list(cm.columns[reordered_columns]) )    \n","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:16:59.981632Z","iopub.execute_input":"2023-11-28T12:16:59.982059Z","iopub.status.idle":"2023-11-28T12:17:08.393264Z","shell.execute_reply.started":"2023-11-28T12:16:59.982026Z","shell.execute_reply":"2023-11-28T12:17:08.392027Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_de_train","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:17:22.960913Z","iopub.execute_input":"2023-11-28T12:17:22.961399Z","iopub.status.idle":"2023-11-28T12:17:23.010497Z","shell.execute_reply.started":"2023-11-28T12:17:22.961363Z","shell.execute_reply":"2023-11-28T12:17:23.008842Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# for l,str_inf in [ [G1S_genes_Tirosh, 'G1S Tirosh'], [G2M_genes_Tirosh, 'G2M Tirosh'],  [G1S_genes_Tirosh + G2M_genes_Tirosh, 'All Tirosh'],\n#                   [list_genes_fastCCsign, 'FastCC Signature'],\n#                   [ G1S_genes_Freeman, 'G1S Freeman' ],  [ G2M_genes_Freeman, 'G2M Freeman' ], [ G1S_genes_Freeman + G2M_genes_Freeman, 'All Freeman' ],   ]: \nfor l,str_inf in [  [G1S_genes_Tirosh + G2M_genes_Tirosh, 'All Tirosh']   ]: \n    ll = set(l) & set(df_de_train.columns) \n    ll = list(ll)\n    print(len(ll), str_inf )\n    X = df_de_train[ ['cell_type'] + list(df_de_train.columns[5:]) ].groupby('cell_type').median()\n    X = X[ll]\n    print(X.shape)\n    clustergrid = sns.clustermap(X)# ,  annot=True, fmt=\".2f\", cmap=\"coolwarm\" )\n    plt.title(str_inf, fontsize = 20 )\n    plt.show()\n    reordered_columns = clustergrid.dendrogram_col.reordered_ind\n    reordered_rows = clustergrid.dendrogram_row.reordered_ind\n    print(len(reordered_rows), len(reordered_columns) )\n    print( list(X.index[reordered_rows]) )\n    print( list(X.columns[reordered_columns]) )\n    \ncol = 'sm_name'    \nfor l,str_inf in [  [G1S_genes_Tirosh + G2M_genes_Tirosh, 'All Tirosh']   ]: \n    ll = set(l) & set(df_de_train.columns) \n    ll = list(ll)\n    print(len(ll), str_inf )\n    X = df_de_train[ [col] + list(df_de_train.columns[5:]) ].groupby(col).median()\n    X = X[ll]\n    print(X.shape)\n    X.index = [t[:20] for t in X.index] # cut too long names\n    clustergrid = sns.clustermap(X)# ,  annot=True, fmt=\".2f\", cmap=\"coolwarm\" )\n    plt.title(str_inf, fontsize = 20 )\n    plt.show()    \n    reordered_columns = clustergrid.dendrogram_col.reordered_ind\n    reordered_rows = clustergrid.dendrogram_row.reordered_ind\n    print(len(reordered_rows), len(reordered_columns) )\n    print( list(X.index[reordered_rows]) )\n    print( list(X.columns[reordered_columns]) )","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:17:29.763236Z","iopub.execute_input":"2023-11-28T12:17:29.763715Z","iopub.status.idle":"2023-11-28T12:17:32.629031Z","shell.execute_reply.started":"2023-11-28T12:17:29.763681Z","shell.execute_reply":"2023-11-28T12:17:32.627630Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Look on compounds ( count = 146  )\n\n15 compounds - 6 times data - only in train \n","metadata":{}},{"cell_type":"code","source":"d = df_de_train[['sm_name','sm_lincs_id','SMILES']].drop_duplicates()\nprint(d.shape)\nd.to_csv('compounds.csv')\ndisplay( d.head(10) )\n\nprint(list(df_de_train['sm_name'].unique() ) )","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:18:36.322587Z","iopub.execute_input":"2023-11-28T12:18:36.323028Z","iopub.status.idle":"2023-11-28T12:18:36.348884Z","shell.execute_reply.started":"2023-11-28T12:18:36.322996Z","shell.execute_reply":"2023-11-28T12:18:36.347159Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"l = [len(s) for s in df_de_train['SMILES']]\nnp.sort(list(set(l)) )","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:18:46.961751Z","iopub.execute_input":"2023-11-28T12:18:46.962282Z","iopub.status.idle":"2023-11-28T12:18:46.973355Z","shell.execute_reply.started":"2023-11-28T12:18:46.962246Z","shell.execute_reply":"2023-11-28T12:18:46.971691Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"display( df_de_train['sm_name'].value_counts().head(20) )\ndisplay( df_de_train['sm_name'].value_counts().tail(10) )\ndf_de_train['sm_name'].value_counts().value_counts()","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:18:50.972745Z","iopub.execute_input":"2023-11-28T12:18:50.973231Z","iopub.status.idle":"2023-11-28T12:18:51.004425Z","shell.execute_reply.started":"2023-11-28T12:18:50.973194Z","shell.execute_reply":"2023-11-28T12:18:51.002569Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Aggregations by compounds, cell_types \n\nIt is used for prediction in early versions of the notebook","metadata":{}},{"cell_type":"code","source":"%%time\ntrain_aggregate_mean_or_median_or_whatever = df_de_train.iloc[:,5:].quantile(0.7)# median()\ntrain_aggregate_mean_or_median_or_whatever\n","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:18:57.652073Z","iopub.execute_input":"2023-11-28T12:18:57.652580Z","iopub.status.idle":"2023-11-28T12:18:58.189756Z","shell.execute_reply.started":"2023-11-28T12:18:57.652546Z","shell.execute_reply":"2023-11-28T12:18:58.188462Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nd = train_aggregate_mean_or_median_or_whatever\nplt.figure(figsize = (20,4) )\nplt.plot(d.values)\nplt.show()\nplt.figure(figsize = (10,4) )\nplt.hist(d.values, bins = 100)\nplt.show()\n\ndisplay( d.describe() )","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:19:01.853480Z","iopub.execute_input":"2023-11-28T12:19:01.853894Z","iopub.status.idle":"2023-11-28T12:19:02.775483Z","shell.execute_reply.started":"2023-11-28T12:19:01.853863Z","shell.execute_reply":"2023-11-28T12:19:02.774130Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load,look on submission data","metadata":{}},{"cell_type":"code","source":"%%time\nfn = '/kaggle/input/open-problems-single-cell-perturbations/id_map.csv'\ndf_id_map = pd.read_csv(fn)\nprint(df_id_map.shape)\ndisplay(df_id_map)\nfn = '/kaggle/input/open-problems-single-cell-perturbations/sample_submission.csv'\ndf = pd.read_csv(fn, index_col = 0)\nprint(df.shape)\ndf","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:19:09.181741Z","iopub.execute_input":"2023-11-28T12:19:09.182536Z","iopub.status.idle":"2023-11-28T12:19:14.325346Z","shell.execute_reply.started":"2023-11-28T12:19:09.182499Z","shell.execute_reply":"2023-11-28T12:19:14.322882Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Baselines. Predictions and submissions. Several approaches.","metadata":{}},{"cell_type":"markdown","source":"# Target mean/median/quantile . Approch 0. \n\nThat is the first baseline for any ML task one needs to do. \n\nConstant prediction. (Different constant for different targets, but it is constant - that means does not depend on sample)\n\n\n    V6 LB 0.666 quantile(0.7)\n    V5 LB 0.657 quantile(0.6) \n            results are again better, that probably indicates some shift between train and public data \n    V4 LB 0.659 - median instead of mean - results are a bit better, \n            it might mean either a bit of presense of outliers, or  public is somewhat different from train - next experiments suggests second is true \n    V1 LB 0.664 - submission of train means - the simplest baseline \n","metadata":{}},{"cell_type":"code","source":"# %%time\n# train_aggregate_mean_or_median_or_whatever = df_de_train.iloc[:,5:].quantile(0.6)# median()\n# train_aggregate_mean_or_median_or_whatever\n\n# # consant for each target submission:\n# for i,col in enumerate( df.columns ):\n#     df[col] = train_aggregate_mean_or_median_or_whatever[col]\n#     if (i%1000) == 0: print(i,col)\n        \n# print(df.shape )        \n# display(df )\n\n# df.to_csv('submission.csv')","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:26:55.813833Z","iopub.execute_input":"2023-11-28T12:26:55.814346Z","iopub.status.idle":"2023-11-28T12:26:55.820939Z","shell.execute_reply.started":"2023-11-28T12:26:55.814312Z","shell.execute_reply":"2023-11-28T12:26:55.819427Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Aggregation by compounds. Approach 1\n\nThat is very similar to the previous approach, but now prediction is not quite constant, but depends on compound, but not on anything else - i.e. not from cell type).  \n\nThe values for each compound are obtained by aggregation of targets over the train. \n\nOne can consider aggregated  mean, median, quantile ... Seems quantile 0.54 is better for LB. \n\n    V9 LB 0.638 groupby by compound and quantile(0.6)\n","metadata":{}},{"cell_type":"code","source":"# %%time\n# train_aggr = df_de_train[ ['sm_name'] + list(df_de_train.columns[5:])  ].groupby('sm_name' ).quantile(0.6)# median()\n# train_aggr\n\n# df = df_id_map.merge(train_aggr, how = 'left', on = 'sm_name')\n# df = df.set_index('id').iloc[:,2:]\n# display( df )\n\n# df.to_csv('submission.csv')","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:22:41.063515Z","iopub.execute_input":"2023-11-28T12:22:41.063893Z","iopub.status.idle":"2023-11-28T12:22:56.762355Z","shell.execute_reply.started":"2023-11-28T12:22:41.063863Z","shell.execute_reply":"2023-11-28T12:22:56.761404Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#  Aggregation by compounds with denoising. Approach 2. \n\nIt is similar to the approach above but we use \"denoising for poors\" - i.e. reduce dimension and then inverse transform it back, the common wisdom (which was e,g, helpful on the previous challenge) that keeping only N main components helps to reduce noise in the data, thus prediction migh become better. We see that it is indeed true for the that challenge.\n\nThat means: first reduce dimensions of the data by say pca100, and  take pca inverse transfrom - we get kind of \"denoised data\".\n\nAnd then apply same as in the previous approach i.e. groupby by compound and aggregation - mean/median/quantile()         \nSo idea is - pca100 reductions hopefully kills some noise. \nIn the previous challenge it worked okay, but the number of samples was not 614, but near 100 000, still any biological data is very noisy and highly correlated - thus it is natural to expect it will help and it does indeed. \n\n    V18 LB 0.623 TSVD-35 denoising  quantile 0.54, tsvd is better than pca/ica - similar to previous OP\n    V15 - bug -  LB 0.626 - NO it was not: TSVD-35 denoising  quantile 0.54 - so worse than pca,ICA, at least for these params\n    V14 LB 0.624 - ICA-35 denoising   quantile 0.54 - so ICA is not better than just pca at least for these the same params\n    Fork: LB 0.624 - pca-35 denoising, quantile 0.54 \n    V13 LB 0.626 pca25 denoising \n    V12 LB 0.627 same with pca100 denoising\n","metadata":{}},{"cell_type":"markdown","source":"## Key params","metadata":{}},{"cell_type":"code","source":"# # predict_method = 'train_aggregation_by_compounds'\n# # predict_method = 'train_aggregation_by_compounds_with_denoising_pca'\n# # predict_method = 'train_aggregation_by_compounds_with_denoising_ICA'\n# predict_method = 'train_aggregation_by_compounds_with_denoising_TSVD'\n# quantile = 0.54\n# n_components = 35\n","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:24:28.295600Z","iopub.execute_input":"2023-11-28T12:24:28.296197Z","iopub.status.idle":"2023-11-28T12:24:28.304072Z","shell.execute_reply.started":"2023-11-28T12:24:28.296153Z","shell.execute_reply":"2023-11-28T12:24:28.302213Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Denoising (by pca/ica/tsvd) and aggregation ","metadata":{}},{"cell_type":"code","source":"# %%time \n# from sklearn.decomposition import PCA\n# from sklearn.decomposition import FastICA\n# from sklearn.decomposition import TruncatedSVD\n\n# Y = df_de_train.iloc[:,5:]\n# print(X.shape)\n# if '_pca' in predict_method:\n#     str_inf_target_dimred = 'PCA' \n#     reducer = PCA(n_components=n_components )\n# elif '_ICA' in predict_method:\n#     str_inf_target_dimred = 'ICA' \n# #     reducer = PCA(n_components=n_components )\n#     reducer = FastICA(n_components=n_components, random_state=0, whiten='unit-variance')\n# elif '_TSVD' in predict_method:\n#     str_inf_target_dimred = 'TSVD' \n# #     reducer = PCA(n_components=n_components )\n# #     reducer = FastICA(n_components=n_components, random_state=0, whiten='unit-variance')\n#     reducer = TruncatedSVD(n_components=n_components, n_iter=7, random_state=42)\n# else:\n#     str_inf_target_dimred = ''\n    \n# print(str_inf_target_dimred , reducer)\n# Yr = reducer.fit_transform(Y)\n# Yr_inv_trans = reducer.inverse_transform(Yr)\n# df_red_inv_trans = pd.DataFrame(Yr_inv_trans, columns = df_de_train.columns[5:])\n# df_red_inv_trans['sm_name'] = df_de_train['sm_name']\n\n# train_aggr_denoised = df_red_inv_trans.groupby('sm_name' ).quantile(quantile)# median()\n# train_aggr_denoised\n","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:24:49.215259Z","iopub.execute_input":"2023-11-28T12:24:49.216841Z","iopub.status.idle":"2023-11-28T12:24:53.840776Z","shell.execute_reply.started":"2023-11-28T12:24:49.216793Z","shell.execute_reply":"2023-11-28T12:24:53.839864Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# # Basic Edit of code above\nimport pandas as pd\nfrom sklearn.decomposition import PCA, FastICA, TruncatedSVD\n\ndef apply_dimension_reduction(data, method='_pca', n_components=10, quantile=0.5):\n    if method == '_pca':\n        reducer = PCA(n_components=n_components)\n        str_inf_target_dimred = 'PCA'\n    elif method == '_ICA':\n        reducer = FastICA(n_components=n_components, random_state=0, whiten='unit-variance')\n        str_inf_target_dimred = 'ICA'\n    elif method == '_TSVD':\n        reducer = TruncatedSVD(n_components=n_components, n_iter=7, random_state=42)\n        str_inf_target_dimred = 'TSVD'\n    else:\n        str_inf_target_dimred = ''\n        reducer = None\n\n    if reducer is not None:\n        Yr = reducer.fit_transform(data)\n        Yr_inv_trans = reducer.inverse_transform(Yr)\n        df_red_inv_trans = pd.DataFrame(Yr_inv_trans, columns=data.columns)\n        df_red_inv_trans['sm_name'] = df_de_train['sm_name']\n\n        train_aggr_denoised = df_red_inv_trans.groupby('sm_name').quantile(quantile)\n        return train_aggr_denoised, str_inf_target_dimred, reducer\n    else:\n        print(f\"Unsupported dimension reduction method: {method}\")\n        return None\n\n# Usage\npredict_method = '_pca'  # Replace with your desired method\nn_components = 5  # Replace with your desired number of components\nquantile = 0.75  # Replace with your desired quantile\n\ntrain_aggr_denoised, str_inf_target_dimred, reducer = apply_dimension_reduction(df_de_train.iloc[:, 5:], method=predict_method, n_components=n_components, quantile=quantile)\n\nif train_aggr_denoised is not None:\n    print(f\"Dimension Reduction Method: {str_inf_target_dimred}\")\n    print(\"Reduced Data:\")\n    print(train_aggr_denoised)\n    print(\"Reducer:\")\n    print(reducer)\n","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:32:46.883922Z","iopub.execute_input":"2023-11-28T12:32:46.884980Z","iopub.status.idle":"2023-11-28T12:32:50.021950Z","shell.execute_reply.started":"2023-11-28T12:32:46.884943Z","shell.execute_reply":"2023-11-28T12:32:50.020503Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Prepapre predictions for submit","metadata":{}},{"cell_type":"code","source":"%%time\n\nif predict_method.startswith('train_aggregation_by_compounds_with_denoising_'):\n    df = df_id_map.merge(train_aggr_denoised, how = 'left', on = 'sm_name')\n    df = df.set_index('id').iloc[:,2:]\nelif predict_method == 'train_aggregation_by_compounds': # Older - approach 1 \n    ## Prepare predictions by direct aggregation of targets by compounds\n    train_aggr_direct = df_de_train[ ['sm_name'] + list(df_de_train.columns[5:])  ].groupby('sm_name' ).quantile(quantile)# median()\n    train_aggr_direct\n    df = df_id_map.merge(train_aggr, how = 'left', on = 'sm_name')\n    df = df.set_index('id').iloc[:,2:]\nelse: # Older - approach 0 \n    # consant for each target submission:\n    for i,col in enumerate( df.columns ):\n        df[col] = train_aggregate_mean_or_median_or_whatever[col]\n        if (i%1000) == 0: print(i,col)\n    \ndf    ","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:33:09.165775Z","iopub.execute_input":"2023-11-28T12:33:09.166264Z","iopub.status.idle":"2023-11-28T12:33:18.563466Z","shell.execute_reply.started":"2023-11-28T12:33:09.166227Z","shell.execute_reply":"2023-11-28T12:33:18.561746Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndf.to_csv('submission.csv')","metadata":{"execution":{"iopub.status.busy":"2023-11-28T12:33:22.813635Z","iopub.execute_input":"2023-11-28T12:33:22.814138Z","iopub.status.idle":"2023-11-28T12:34:17.473844Z","shell.execute_reply.started":"2023-11-28T12:33:22.814087Z","shell.execute_reply":"2023-11-28T12:34:17.472454Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Target endcoding and Ridge. Approach 3. \n\nHere we will try to train the first ML model. \n\nWe will do target encodnging and train Ridge model upon the target encoded cell type and compound - two features which are actually categorical, \n\n    Versions 23-27 - first modeling approach - target encode cell type and compound and simple Ridge model - results currently worse than naive approach used before - without models. \n        V27 LB 0.659 LB Ridge nCT1 nCD25 Al10 TSVD35 - even more relaxed : alpha and encoded_compound size \n        V26 LB 0.668 try to relax a model a bit: Ridge nCT1 nCD10 Al100 TSVD35 - results are better near mean/zero submission\n        V25 LB 0.677 - stronger constraining the model: Ridge nCT1 nCD5 Al100 TSVD35, but still we are worse than even predict by mean \n        V24 LB 0.702 - try avoid overfit Ridge nCT3 nCD10 Al100 TSVD35 - better results but still bad \n        V23 LB 0.747 Ridge model in target encoded features Ridge nCT10 nCD35 Al1 TSVD35 - modeling gives worse results than simple approach, might be bug or overfit\n","metadata":{}},{"cell_type":"markdown","source":"## Key params","metadata":{}},{"cell_type":"code","source":"# n_components_for_cell_type_encoding = 1\n# n_components_for_compound_encoding = 25\n# alpha_regularization_for_linear_models = 10\n\n# # predict_method\n\n# model_type = 'Ridge'\n","metadata":{"execution":{"iopub.status.busy":"2023-11-28T06:01:16.903675Z","iopub.execute_input":"2023-11-28T06:01:16.904070Z","iopub.status.idle":"2023-11-28T06:01:16.909349Z","shell.execute_reply.started":"2023-11-28T06:01:16.904037Z","shell.execute_reply":"2023-11-28T06:01:16.908011Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# str_model_id = model_type\n# str_model_id += ' nCT'+ str(n_components_for_cell_type_encoding)\n# str_model_id += ' nCD'+ str(n_components_for_compound_encoding)\n# str_model_id += ' Al'+ str(alpha_regularization_for_linear_models)\n# str_model_id += ' ' +str_inf_target_dimred+str( n_components )\n\n# print( str_model_id )\n","metadata":{"execution":{"iopub.status.busy":"2023-11-28T06:01:18.864583Z","iopub.execute_input":"2023-11-28T06:01:18.864977Z","iopub.status.idle":"2023-11-28T06:01:18.870815Z","shell.execute_reply.started":"2023-11-28T06:01:18.864946Z","shell.execute_reply":"2023-11-28T06:01:18.869580Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Target encoded features ","metadata":{"execution":{"iopub.status.busy":"2023-09-15T09:29:29.58894Z","iopub.execute_input":"2023-09-15T09:29:29.589366Z","iopub.status.idle":"2023-09-15T09:29:29.623736Z","shell.execute_reply.started":"2023-09-15T09:29:29.589335Z","shell.execute_reply":"2023-09-15T09:29:29.623004Z"}}},{"cell_type":"code","source":"# %%time\n# # Yr = reducer.fit_transform(X)\n# # n_components_for_cell_type_encoding = 10\n# df_tmp = pd.DataFrame(Yr[:, :n_components_for_cell_type_encoding  ], index = df_de_train.index  )\n# df_tmp['column for aggregation'] = df_de_train['cell_type']\n# df_cell_type_encoded = df_tmp.groupby('column for aggregation').quantile( quantile )\n# print('df_cell_type_encoded.shape', df_cell_type_encoded.shape )\n# display( df_cell_type_encoded )\n\n\n# # n_components_for_compound_encoding = 10\n# df_tmp = pd.DataFrame(Yr[:, :n_components_for_compound_encoding  ], index = df_de_train.index  )\n# df_tmp['column for aggregation'] = df_de_train['sm_name']\n# df_compound_encoded = df_tmp.groupby('column for aggregation').quantile( quantile )\n# print('df_compound_encoded.shape', df_compound_encoded.shape )\n# display( df_compound_encoded )\n\n","metadata":{"execution":{"iopub.status.busy":"2023-11-28T06:01:21.404463Z","iopub.execute_input":"2023-11-28T06:01:21.404836Z","iopub.status.idle":"2023-11-28T06:01:21.434564Z","shell.execute_reply.started":"2023-11-28T06:01:21.404808Z","shell.execute_reply":"2023-11-28T06:01:21.433476Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Prepare X_train, X_submit - target encoded cell type and compound features","metadata":{}},{"cell_type":"code","source":"# %%time\n# X_train = np.zeros( (len( df_de_train ) , n_components_for_cell_type_encoding + n_components_for_compound_encoding ))\n\n# for i in range(len( X_train )):\n#     cell_type = df_de_train['cell_type'].iat[i] \n#     X_train[i,:n_components_for_cell_type_encoding] = df_cell_type_encoded.loc[cell_type,:].values  \n#     compound = df_de_train['sm_name'].iat[i] \n#     X_train[i,n_components_for_cell_type_encoding:] = df_compound_encoded.loc[ compound, : ].values\n# print( X_train.shape)     \n# print( X_train)     \n    \n\n# X_submit = np.zeros( (len( df_id_map ) , n_components_for_cell_type_encoding + n_components_for_compound_encoding ))\n# for i in range(len( X_submit )):\n#     cell_type = df_id_map['cell_type'].iat[i] \n#     X_submit[i,:n_components_for_cell_type_encoding] = df_cell_type_encoded.loc[cell_type,:].values  \n#     compound = df_id_map['sm_name'].iat[i] \n#     X_submit[i,n_components_for_cell_type_encoding:] = df_compound_encoded.loc[ compound, : ].values\n    \n    \n# print( X_submit.shape)     \n# print( X_submit)     \n    ","metadata":{"execution":{"iopub.status.busy":"2023-11-28T06:01:28.866911Z","iopub.execute_input":"2023-11-28T06:01:28.867309Z","iopub.status.idle":"2023-11-28T06:01:28.902244Z","shell.execute_reply.started":"2023-11-28T06:01:28.867278Z","shell.execute_reply":"2023-11-28T06:01:28.901081Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Modeling","metadata":{}},{"cell_type":"code","source":"# %%time\n# from sklearn.linear_model import Ridge\n\n# model = Ridge(alpha=alpha_regularization_for_linear_models)\n# print(model)\n# model.fit(X_train, Yr)\n\n# Y_submit = reducer.inverse_transform(   model.predict(X_submit) )\n# print(Y_submit.shape)\n# Y_submit","metadata":{"execution":{"iopub.status.busy":"2023-11-28T06:00:38.424785Z","iopub.execute_input":"2023-11-28T06:00:38.425223Z","iopub.status.idle":"2023-11-28T06:00:38.513543Z","shell.execute_reply.started":"2023-11-28T06:00:38.425181Z","shell.execute_reply":"2023-11-28T06:00:38.512399Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Save submission CSV","metadata":{}},{"cell_type":"code","source":"# %%time\n# df_submit = pd.DataFrame(Y_submit, columns = df_de_train.columns[5:])\n# df_submit.index.name = 'id'\n# print( df_submit.shape )\n# display(df_submit)\n# df_submit.to_csv('submission.csv')","metadata":{"execution":{"iopub.status.busy":"2023-11-28T05:58:38.189154Z","iopub.execute_input":"2023-11-28T05:58:38.189582Z","iopub.status.idle":"2023-11-28T05:58:47.473611Z","shell.execute_reply.started":"2023-11-28T05:58:38.189549Z","shell.execute_reply":"2023-11-28T05:58:47.472390Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Catboost -  approach 4\n\n    V34 LB 0.688 iterations=3,  depth = 6 \n    V33 LB 0.709 Catboost iterations 10,  depth = 6    - try to reduce overfit - but even worse  \n    V32 LB 0.912 - Catboost depth = 2, iteration 100 - try to reduce overfit - but even worse  \n    V30 LB 0.905 ( terrible :)  - Catboost first draft -  depth = 6 , iteration 100\n","metadata":{}},{"cell_type":"code","source":"# import catboost\n# from catboost import CatBoostRegressor, Pool\n# categorical_features = ['cell_type','sm_name']\n# model = CatBoostRegressor(iterations=3,  # Number of boosting iterations\n#                           depth=6,        # Depth of the trees\n#                           learning_rate=0.1,  # Learning rate\n#                           loss_function='RMSE',  # Specify your loss function (e.g., RMSE for regression)\n#                           cat_features=categorical_features,  # Categorical features\n#                           verbose=0)  # Set verbose to 0 to suppress output\n","metadata":{"execution":{"iopub.status.busy":"2023-09-19T20:05:32.875914Z","iopub.execute_input":"2023-09-19T20:05:32.876166Z","iopub.status.idle":"2023-09-19T20:05:33.250242Z","shell.execute_reply.started":"2023-09-19T20:05:32.876142Z","shell.execute_reply":"2023-09-19T20:05:33.249267Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Yr.shape","metadata":{"execution":{"iopub.status.busy":"2023-09-19T20:05:33.25142Z","iopub.execute_input":"2023-09-19T20:05:33.251866Z","iopub.status.idle":"2023-09-19T20:05:33.258857Z","shell.execute_reply.started":"2023-09-19T20:05:33.251842Z","shell.execute_reply":"2023-09-19T20:05:33.257646Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# %%time\n# Y_reduced_submit = np.zeros( (len(df_id_map) , Yr.shape[1] )   )\n# df_train = df_de_train[['cell_type','sm_name']]\n\n# for k in range(Yr.shape[1]):\n#     train_data = Pool(data=df_train, \n#                       label=Yr[:,k],\n#                       cat_features=categorical_features)\n#     test_data = Pool(data=df_id_map[['cell_type','sm_name']], \n#                       cat_features=categorical_features)\n#     model.fit(train_data)    \n#     predicted_value = model.predict(test_data)\n#     Y_reduced_submit[:,k] = predicted_value\n    \n    ","metadata":{"execution":{"iopub.status.busy":"2023-09-19T20:05:33.260251Z","iopub.execute_input":"2023-09-19T20:05:33.260578Z","iopub.status.idle":"2023-09-19T20:05:34.891362Z","shell.execute_reply.started":"2023-09-19T20:05:33.260546Z","shell.execute_reply":"2023-09-19T20:05:34.890333Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# %%time\n# Y_submit = reducer.inverse_transform(  Y_reduced_submit )\n# print(Y_submit.shape)\n# Y_submit","metadata":{"execution":{"iopub.status.busy":"2023-09-19T20:05:34.892304Z","iopub.execute_input":"2023-09-19T20:05:34.892599Z","iopub.status.idle":"2023-09-19T20:05:34.920332Z","shell.execute_reply.started":"2023-09-19T20:05:34.892577Z","shell.execute_reply":"2023-09-19T20:05:34.91958Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# %%time\n# df_submit = pd.DataFrame(Y_submit, columns = df_de_train.columns[5:])\n# df_submit.index.name = 'id'\n# print( df_submit.shape )\n# display(df_submit)\n# df_submit.to_csv('submission.csv')","metadata":{"execution":{"iopub.status.busy":"2023-09-19T20:05:34.921758Z","iopub.execute_input":"2023-09-19T20:05:34.925775Z","iopub.status.idle":"2023-09-19T20:05:40.88626Z","shell.execute_reply.started":"2023-09-19T20:05:34.925736Z","shell.execute_reply":"2023-09-19T20:05:40.885132Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Target Encoded simplified. Approach 5\n\nTarget encoding and Ridge, but here is Ridge is done only on 1 feature and thus it should be quite similar to approach 3 with just denoise+groupby, but a kind of Ridge instead of \"mean/median/quantile\". \n\nIf it were sucessful we will add more features. \n\nBut prelimary result is not good: \n\n    V37 0.707 approach 5 - aggr by compound + Ridge.\n","metadata":{}},{"cell_type":"code","source":"# %%time \n\n# # predict_method = 'train_aggregation_by_compounds'\n# # predict_method = 'train_aggregation_by_compounds_with_denoising_pca'\n# # predict_method = 'train_aggregation_by_compounds_with_denoising_ICA'\n# predict_method = 'train_aggregation_by_compounds_with_denoising_TSVD'\n# quantile = 0.54\n# n_components = 35\n\n# from sklearn.decomposition import PCA\n# from sklearn.decomposition import FastICA\n# from sklearn.decomposition import TruncatedSVD\n\n# Y = df_de_train.iloc[:,5:]\n# print(X.shape)\n# if '_pca' in predict_method:\n#     str_inf_target_dimred = 'PCA' \n#     reducer = PCA(n_components=n_components )\n# elif '_ICA' in predict_method:\n#     str_inf_target_dimred = 'ICA' \n# #     reducer = PCA(n_components=n_components )\n#     reducer = FastICA(n_components=n_components, random_state=0, whiten='unit-variance')\n# elif '_TSVD' in predict_method:\n#     str_inf_target_dimred = 'TSVD' \n# #     reducer = PCA(n_components=n_components )\n# #     reducer = FastICA(n_components=n_components, random_state=0, whiten='unit-variance')\n#     reducer = TruncatedSVD(n_components=n_components, n_iter=7, random_state=42)\n# else:\n#     str_inf_target_dimred = ''\n    \n# print(str_inf_target_dimred , reducer)\n# Yr = reducer.fit_transform(Y)\n# # Yr_inv_trans = reducer.inverse_transform(Yr)\n# df_red = pd.DataFrame(Yr)\n# df_red['sm_name'] = df_de_train['sm_name']\n\n# train_aggr_red = df_red.groupby('sm_name' ).quantile(quantile)# median()\n# train_aggr_red = train_aggr_red.reset_index()\n# train_aggr_red","metadata":{"execution":{"iopub.status.busy":"2023-09-19T20:05:40.887721Z","iopub.execute_input":"2023-09-19T20:05:40.888282Z","iopub.status.idle":"2023-09-19T20:05:42.317576Z","shell.execute_reply.started":"2023-09-19T20:05:40.88825Z","shell.execute_reply":"2023-09-19T20:05:42.31685Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# # Basic Edit based on code above\n# %%time\n# import pandas as pd\n# from sklearn.decomposition import PCA, FastICA, TruncatedSVD\n\n# # Set parameters\n# predict_method = 'train_aggregation_by_compounds_with_denoising_TSVD'\n# quantile = 0.54\n# n_components = 35\n\n# # Extract features from DataFrame\n# Y = df_de_train.iloc[:, 5:]\n\n# # Dimensionality reduction based on predict_method\n# if '_pca' in predict_method:\n#     str_inf_target_dimred = 'PCA'\n#     reducer = PCA(n_components=n_components)\n# elif '_ICA' in predict_method:\n#     str_inf_target_dimred = 'ICA'\n#     reducer = FastICA(n_components=n_components, random_state=0, whiten='unit-variance')\n# elif '_TSVD' in predict_method:\n#     str_inf_target_dimred = 'TSVD'\n#     reducer = TruncatedSVD(n_components=n_components, n_iter=7, random_state=42)\n# else:\n#     str_inf_target_dimred = ''\n#     reducer = None  # or any default reducer\n\n# print(str_inf_target_dimred, reducer)\n\n# # Perform dimensionality reduction\n# if reducer is not None:\n#     Yr = reducer.fit_transform(Y)\n#     df_red = pd.DataFrame(Yr)\n#     df_red['sm_name'] = df_de_train['sm_name']\n\n#     # Group by 'sm_name' and calculate quantile\n#     train_aggr_red = df_red.groupby('sm_name').quantile(quantile).reset_index()\n# else:\n#     # No dimensionality reduction, use original DataFrame\n#     train_aggr_red = df_de_train.copy()\n\n# print(train_aggr_red)\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# %%time\n# df_X =  df_de_train[ ['sm_name']].copy()\n# df_X =  df_X.merge( train_aggr_red, on = 'sm_name' )\n# df_X = df_X.iloc[:,1:]\n# display( df_X )\n \n\n# df_X_submit = df_id_map[['sm_name']].copy()\n# df_X_submit = df_X_submit.merge( train_aggr_red, on = 'sm_name' ).iloc[:,1:]\n# display(df_X_submit)","metadata":{"execution":{"iopub.status.busy":"2023-11-28T06:02:37.495758Z","iopub.execute_input":"2023-11-28T06:02:37.496426Z","iopub.status.idle":"2023-11-28T06:02:37.560571Z","shell.execute_reply.started":"2023-11-28T06:02:37.496364Z","shell.execute_reply":"2023-11-28T06:02:37.559833Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# %%time \n# from sklearn.linear_model import Ridge\n# alpha_regularization_for_linear_models = 1\n# model = Ridge(alpha=alpha_regularization_for_linear_models)\n\n# Yr_submit = np.zeros(   (len(df_id_map) , Yr.shape[1] )   )\n# for i in range(df_X.shape[1]):\n#     model.fit(df_X[[ df_X.columns[i] ]], Yr[:,i])\n#     Yr_submit[:,i] = model.predict(df_X_submit[[ df_X_submit.columns[i] ]])\n    ","metadata":{"execution":{"iopub.status.busy":"2023-11-28T06:02:41.163467Z","iopub.execute_input":"2023-11-28T06:02:41.164204Z","iopub.status.idle":"2023-11-28T06:02:41.326330Z","shell.execute_reply.started":"2023-11-28T06:02:41.164171Z","shell.execute_reply":"2023-11-28T06:02:41.325041Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# %%time\n# Y_submit = reducer.inverse_transform(  Yr_submit )\n# print(Y_submit.shape)\n# Y_submit\n\n# df_submit = pd.DataFrame(Y_submit, columns = df_de_train.columns[5:])\n# df_submit.index.name = 'id'\n# print( df_submit.shape )\n# display(df_submit)\n# df_submit.to_csv('submission.csv')","metadata":{"execution":{"iopub.status.busy":"2023-11-28T06:02:43.554615Z","iopub.execute_input":"2023-11-28T06:02:43.555029Z","iopub.status.idle":"2023-11-28T06:02:52.775786Z","shell.execute_reply.started":"2023-11-28T06:02:43.554996Z","shell.execute_reply":"2023-11-28T06:02:52.774708Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('%.1f seconds passed total '%(time.time()-t0start) )\nprint('%.1f minutes passed total '%( (time.time()-t0start)/60)  )\nprint('%.2f hours passed total '%( (time.time()-t0start)/3600)  )","metadata":{"execution":{"iopub.status.busy":"2023-11-28T06:02:57.503936Z","iopub.execute_input":"2023-11-28T06:02:57.505332Z","iopub.status.idle":"2023-11-28T06:02:57.513306Z","shell.execute_reply.started":"2023-11-28T06:02:57.505284Z","shell.execute_reply":"2023-11-28T06:02:57.512187Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}