{"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":30587,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# XGBoost in Compressed Space for Single-Cell Perturbations \n\n## Context\n\n- Competition context: https://www.kaggle.com/competitions/open-problems-single-cell-perturbations/overview\n\n- Data context: https://www.kaggle.com/competitions/open-problems-single-cell-perturbations/data\n\n# Introduction\n\nRecent advances in single-cell technologies offer new insights into the complex functions of the human body's cells, but leveraging this knowledge for drug development requires linking chemical perturbations and cell state changes. The **\"Open Problems – Single-Cell Perturbations\"** competition proposes accurately predicting changes in gene expression in diverse cell types and treatments. The results obtained in this competition might be a step forward in the fast development of new medicines using Machine Learning.\n\nIn this project, we built an XGBoost model to tackle a multi-output regression challenge. *Noteworthy for its computational efficiency and robustness to dataset noise, our approach operates effectively within a compressed space.*  Our model allows estimation of compound impact on gene expression in a target cell type by leveraging averages from other cells in the compressed space. This approximation can serve as a valuable baseline for training advanced models without prior biological knowledge.\n\n#### Goal \n\nIn this project, we aimed to build a simple and accurate ML model to predict differential gene expression for different cell types affected by various chemical substances (referred to further as \"drugs\"). \n\n\n*Note*: We did not integrate any prior biological knowledge for feature augmentation, i.e. feature engineering was done using only cell type, chemical compound name, and 18211 gene differential expression from `de_train.parquet`. \n\n> Here **differential gene expression (DE)**  is  a value $-\\log_{10}(\\text{p-value}) * \\text{sign} \\left(\\log_2  \\frac{\\text{gene expression (treatment) }} { \\text{gene expression (control)} } \\right)$ per compared conditions.\n\n\n#### Performance \n\nWithout ensembling with other models, our method gives `0.594` for public and `0.777` for private scores (time of computing is around 3 min). Bagging of XGB will give an additional insignificant improvement in the scores. \n\n# Project Structure \n\n* Brief Exploratory Data Analysis\n* Metric\n* Modelling \n    * Building Baseline Models\n    * XGB Demo\n    * XBG Training \n* Submission\n* Conclusions \n* References ","metadata":{}},{"cell_type":"code","source":"import numpy as np \nimport pandas as pd \nimport math\nimport random\nimport time\nimport os\n\nfrom sklearn.model_selection import KFold\nfrom sklearn.decomposition import TruncatedSVD \n\nimport xgboost as xgb\nimport category_encoders as ce\n\nimport scipy.stats\nfrom scipy.stats import norm, multinomial \n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n%matplotlib inline\n%config InlineBackend.figure_formats=['svg']\nfrom IPython.display import display, Markdown\n        \nimport warnings\nwarnings.filterwarnings('ignore')","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:07:19.339793Z","iopub.execute_input":"2023-12-11T18:07:19.340259Z","iopub.status.idle":"2023-12-11T18:07:19.358308Z","shell.execute_reply.started":"2023-12-11T18:07:19.340228Z","shell.execute_reply":"2023-12-11T18:07:19.356778Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Brief Exploratory Data Analysis","metadata":{}},{"cell_type":"markdown","source":"#### Columns in `de_train.parquet` \n\n`cell_type` - names of cell types;\n\n`sm_name` - names of drugs;\n\n`sm_lincs_id` - drug ID;\n\n`SMILES` - drug name with encoded molecular structures;\n\n`control` - if it was affected by technical drug/no drug to introduce a bias (for further demultiplexing);\n\n\n`A1BG`, `A1BG-AS1`, ... `ZZEF1` -  18211 gene names \n\n\nAs it is written in the competition's description, for evaluation performance a model will be provided only 2 features `cell_type` and `sm_name` to make a prediction of the target variables - differential expression of 18211 genes. \n\n#### Global Names\n\n`df` - the train dataset; pd.DataFrame of `de_train.parquet` \n\n`target` - list of all genes\n\n`features` - contains \"cell_type\", \"sm_name\"","metadata":{}},{"cell_type":"code","source":"fn = '/kaggle/input/open-problems-single-cell-perturbations/de_train.parquet'\ndf = pd.read_parquet(fn) \ndf.info()\ndf","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:07:20.096788Z","iopub.execute_input":"2023-12-11T18:07:20.097251Z","iopub.status.idle":"2023-12-11T18:07:22.787701Z","shell.execute_reply.started":"2023-12-11T18:07:20.097216Z","shell.execute_reply":"2023-12-11T18:07:22.786292Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"features = [\"cell_type\", \"sm_name\" ]\ntarget = [value for value in df.columns if value not in [\"sm_lincs_id\", \"SMILES\", \"control\"] + features ]","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:07:22.789792Z","iopub.execute_input":"2023-12-11T18:07:22.790219Z","iopub.status.idle":"2023-12-11T18:07:22.806465Z","shell.execute_reply.started":"2023-12-11T18:07:22.790186Z","shell.execute_reply":"2023-12-11T18:07:22.804886Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.countplot(data=df, x='cell_type')\nplt.title('Distribution of Cell Types')\nplt.xticks(rotation=45)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:07:22.808540Z","iopub.execute_input":"2023-12-11T18:07:22.811161Z","iopub.status.idle":"2023-12-11T18:07:23.189893Z","shell.execute_reply.started":"2023-12-11T18:07:22.811112Z","shell.execute_reply":"2023-12-11T18:07:23.188590Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The dataset contains 6 cell types. **\"B cells\"** and **\"Myeloid cells\"** are underrepresented. ","metadata":{}},{"cell_type":"code","source":"df_control = df[df['control'] == True] \n\nprint(\"Names of control drugs:\", list(df_control[\"sm_name\"].unique()))\nprint(\"Cells affected by control drugs:\")\npd.DataFrame(df_control[\"cell_type\"].value_counts())","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:07:23.193032Z","iopub.execute_input":"2023-12-11T18:07:23.193456Z","iopub.status.idle":"2023-12-11T18:07:23.215108Z","shell.execute_reply.started":"2023-12-11T18:07:23.193419Z","shell.execute_reply":"2023-12-11T18:07:23.213778Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"12 observations are from control group: 2 drugs for 6 cell types.  ","metadata":{}},{"cell_type":"code","source":"print(\"Total number of drugs:\", df[\"sm_name\"].nunique())","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:07:23.216837Z","iopub.execute_input":"2023-12-11T18:07:23.217305Z","iopub.status.idle":"2023-12-11T18:07:23.225063Z","shell.execute_reply.started":"2023-12-11T18:07:23.217263Z","shell.execute_reply":"2023-12-11T18:07:23.223814Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(30,5))\nsns.countplot(data=df, x='sm_name')\nplt.title('Occurrence of Chemical Compounds')\nplt.xticks(rotation=90)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:07:23.226966Z","iopub.execute_input":"2023-12-11T18:07:23.227441Z","iopub.status.idle":"2023-12-11T18:07:25.453861Z","shell.execute_reply.started":"2023-12-11T18:07:23.227399Z","shell.execute_reply":"2023-12-11T18:07:25.449028Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"n_cells = df[\"cell_type\"].nunique()\nn_drugs = df[\"sm_name\"].nunique()\n\nprint(\"Number of possible combinations of cell types and drugs: \", n_drugs*n_cells)","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:07:25.456256Z","iopub.execute_input":"2023-12-11T18:07:25.456677Z","iopub.status.idle":"2023-12-11T18:07:25.463607Z","shell.execute_reply.started":"2023-12-11T18:07:25.456644Z","shell.execute_reply":"2023-12-11T18:07:25.462708Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As we know from the competition's description, the results will be tested on the B cells and Myeloid cell data with chemical substances represented for other cell types.\n\n*Have we the observations of the effect of 146 drugs for the rest 4 cell types?* No, because $146*4 + 17*2 = 618 $, but a number of observations is 614 (which is also seen from \"Distribution of Cell Types\" above ) ","metadata":{}},{"cell_type":"code","source":"for c in [\"NK cells\",\"T cells CD4+\", \"T cells CD8+\", \"T regulatory cells\"] :\n    for d in df[\"sm_name\"].unique():\n        if d not in df[df[\"cell_type\"]==c][\"sm_name\"].to_list():\n            print(f\"Missed compound for {c} is {d}\")","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:07:25.465430Z","iopub.execute_input":"2023-12-11T18:07:25.465781Z","iopub.status.idle":"2023-12-11T18:07:32.149176Z","shell.execute_reply.started":"2023-12-11T18:07:25.465751Z","shell.execute_reply":"2023-12-11T18:07:32.147786Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for col in df.columns:\n    if df[col].isnull().sum()>0:\n        print(f\"Found NaNs in {col}\")","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:07:32.152775Z","iopub.execute_input":"2023-12-11T18:07:32.153201Z","iopub.status.idle":"2023-12-11T18:07:34.836412Z","shell.execute_reply.started":"2023-12-11T18:07:32.153168Z","shell.execute_reply":"2023-12-11T18:07:34.835236Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"`de_train.parquet` doesn't contain NaNs. ","metadata":{}},{"cell_type":"markdown","source":"DE for a gene in the training dataset is Gaussian distributed with around zero mean.","metadata":{}},{"cell_type":"code","source":"np.random.seed(seed=42)\nrandom_genes = np.random.choice(list(df[target].columns), 5 )\n\nfor gene in random_genes:\n    plt.figure()\n    sns.histplot(data=df, x=gene, kde=True)\n    plt.title(f'Distribution of {gene}')\n    plt.xlabel(\"DE\")\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:07:34.838755Z","iopub.execute_input":"2023-12-11T18:07:34.840102Z","iopub.status.idle":"2023-12-11T18:07:39.683838Z","shell.execute_reply.started":"2023-12-11T18:07:34.840050Z","shell.execute_reply":"2023-12-11T18:07:39.680980Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"DE distributions differ for different cell types:","metadata":{}},{"cell_type":"code","source":"for gene in random_genes:\n    plt.figure(figsize=(10,5))\n    sns.boxplot(data=df, x='cell_type', y=gene)\n    plt.title(f'Boxplot for {gene} across cell types')\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:07:39.685609Z","iopub.execute_input":"2023-12-11T18:07:39.686020Z","iopub.status.idle":"2023-12-11T18:07:41.814842Z","shell.execute_reply.started":"2023-12-11T18:07:39.685987Z","shell.execute_reply":"2023-12-11T18:07:41.813182Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Below one can see the distributions of all genes for a specific subset of cells and chemical substances. Most of the distributions have a Gaussian form, but they can differ strongly, as was also noticed and described in [[1]](https://www.kaggle.com/code/ambrosm/scp-eda-which-makes-sense/notebook). ","metadata":{}},{"cell_type":"code","source":"row_num = 0\ncell = df.loc[row_num, [\"cell_type\", \"sm_name\"]].values[0]\ndrug = df.loc[row_num, [\"cell_type\", \"sm_name\"]].values[1]\nplt.figure()\nplt.hist(df.loc[row_num,target], bins=100)\nplt.title(f\"Typical distribution of DE in all genes ({cell} and {drug}) \")\nplt.xlabel(\"DE\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:07:41.816527Z","iopub.execute_input":"2023-12-11T18:07:41.816949Z","iopub.status.idle":"2023-12-11T18:07:42.599887Z","shell.execute_reply.started":"2023-12-11T18:07:41.816915Z","shell.execute_reply":"2023-12-11T18:07:42.598735Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"row_num = 183 \ncell = df.loc[row_num, [\"cell_type\", \"sm_name\"]].values[0]\ndrug = df.loc[row_num, [\"cell_type\", \"sm_name\"]].values[1]\nplt.figure()\nplt.hist(df.loc[row_num,target], bins=100)\nplt.title(f\"Non-typical distribution of DE in all genes ({cell} and {drug})\")\nplt.xlabel(\"DE\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:07:42.601391Z","iopub.execute_input":"2023-12-11T18:07:42.601761Z","iopub.status.idle":"2023-12-11T18:07:43.372974Z","shell.execute_reply.started":"2023-12-11T18:07:42.601729Z","shell.execute_reply":"2023-12-11T18:07:43.371757Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"row_num = 292\ncell = df.loc[row_num, [\"cell_type\", \"sm_name\"]].values[0]\ndrug = df.loc[row_num, [\"cell_type\", \"sm_name\"]].values[1]\nplt.figure()\nplt.hist(df.loc[row_num,target], bins=100)\nplt.title(f\"Non-typical distribution of DE in all genes ({cell} and {drug})\")\nplt.xlabel(\"DE\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:07:43.374348Z","iopub.execute_input":"2023-12-11T18:07:43.374717Z","iopub.status.idle":"2023-12-11T18:07:44.128403Z","shell.execute_reply.started":"2023-12-11T18:07:43.374686Z","shell.execute_reply":"2023-12-11T18:07:44.127489Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Maximum value of DE:\", round(df[target].max().max(),2))\nprint(\"Minimum value of DE:\", round(df[target].min().min(),2))","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:07:44.129576Z","iopub.execute_input":"2023-12-11T18:07:44.130187Z","iopub.status.idle":"2023-12-11T18:07:44.298426Z","shell.execute_reply.started":"2023-12-11T18:07:44.130154Z","shell.execute_reply":"2023-12-11T18:07:44.296891Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As one can see above, usually gene distribution is Gaussian-like with a mean around zero. How does the distribution of genes with min/max DE look like?\n\nIt's shown below, that these extreme values are outliers.","metadata":{}},{"cell_type":"code","source":"top_de_genes = df[target].max().sort_values()[-3:]\n\nfor gene in top_de_genes.index:\n    plt.figure()\n    sns.boxplot(data=df, x=gene)\n    plt.title(f'Distribution of {gene}')\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:07:44.300731Z","iopub.execute_input":"2023-12-11T18:07:44.301280Z","iopub.status.idle":"2023-12-11T18:07:45.304949Z","shell.execute_reply.started":"2023-12-11T18:07:44.301228Z","shell.execute_reply":"2023-12-11T18:07:45.303038Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Another look at the data, which will be used further, is a \"signal\" along rows or columns. It is a way to assess what exactly is going on with a model's prediction for the test data. ","metadata":{}},{"cell_type":"code","source":"for g in random_genes:\n    plt.plot(range(df.shape[0]), df.loc[:,g], label=g)\nplt.xlabel(\"Row Number\")\nplt.ylabel(\"DE\")\nplt.legend()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:07:45.307883Z","iopub.execute_input":"2023-12-11T18:07:45.309082Z","iopub.status.idle":"2023-12-11T18:07:45.735377Z","shell.execute_reply.started":"2023-12-11T18:07:45.309007Z","shell.execute_reply":"2023-12-11T18:07:45.733856Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"first_genes = 500\nrow_ind = 1\n\nplt.plot(range(first_genes), df.loc[row_ind,target[:first_genes]], label=f\"row #{row_ind}\")\nplt.plot(range(first_genes), df.loc[row_ind+1,target[:first_genes]], label=f\"row #{row_ind+1}\")\nplt.xlabel(\"Gene Number\")\nplt.ylabel(\"DE\")\nplt.legend()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:07:45.736770Z","iopub.execute_input":"2023-12-11T18:07:45.737175Z","iopub.status.idle":"2023-12-11T18:07:46.085953Z","shell.execute_reply.started":"2023-12-11T18:07:45.737143Z","shell.execute_reply":"2023-12-11T18:07:46.084617Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Metric\nMetric in this competition is Mean Rowwise Root Mean Squared Error, which is \n\n$$\nMRRMSE=\\frac{1}{R}\\sum_{i=1}^R \\left( \\frac{1}{n} \\sum_{j=1}^n (y_{ij} - \\hat{y}_{ij})^2 \\right)^{1/2}\n$$\n\n\n$R$ is the number of scored rows (number of observations), $n$ is the number of genes.","metadata":{}},{"cell_type":"code","source":"def metric(y, y_hat):\n    \"\"\"\n    Calculates the Mean Row-wise Root Mean Squared Error (MR_RMSE) between two matrices.\n\n    Args:\n    - y (array-like): Ground truth matrix.\n    - y_hat (array-like): Predicted matrix.\n\n    Returns:\n    float: Mean Row-wise Root Mean Squared Error (MR_RMSE) between y and y_hat.\n\n    Raises:\n    ValueError: If the input matrices y and y_hat do not have the same shape.\n\n    \"\"\"\n    y = np.array(y)\n    y_hat = np.array(y_hat)\n    \n    # Check if the input matrices have the same shape\n    if y_hat.shape != y.shape:\n        raise ValueError(\"Input matrices must have the same shape\")\n\n    # Calculate the squared differences element-wise\n    squared_diff = (y - y_hat)**2\n\n    # Calculate the mean row-wise RMSE\n    rowwise_rmse = np.sqrt(np.mean(squared_diff, axis=1))\n\n    # Calculate the mean of row-wise RMSE values\n    mr_rmse = np.mean(rowwise_rmse, axis=0)\n            \n    return mr_rmse","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:07:46.088065Z","iopub.execute_input":"2023-12-11T18:07:46.089217Z","iopub.status.idle":"2023-12-11T18:07:46.100278Z","shell.execute_reply.started":"2023-12-11T18:07:46.089165Z","shell.execute_reply":"2023-12-11T18:07:46.098926Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's look at the behavior of the metric with Gaussian noise in the data.","metadata":{}},{"cell_type":"code","source":"np.random.seed(42)\n\nbatch_sizes = range(35,400,20)\nnum_genes = len(target)\n\nrecords = []\n\nfor size in batch_sizes: \n    rand_err_df = df[target].copy()\n    for i in range(size):\n        rand_err_df.loc[:i] -= 0.055* np.random.normal(size=num_genes)\n        \n    m = metric(df[target].loc[:size], rand_err_df.loc[:size] )\n    records.append(m)\n    \nplt.plot(batch_sizes, records)\nplt.plot([255,255], [min(records), max(records)], '--',label=\"public test dataset size\")\n\nplt.legend()\nplt.xlabel(\"Number of Rows\")\nplt.ylabel(\"Metric\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:07:46.102395Z","iopub.execute_input":"2023-12-11T18:07:46.103719Z","iopub.status.idle":"2023-12-11T18:09:26.910085Z","shell.execute_reply.started":"2023-12-11T18:07:46.103663Z","shell.execute_reply":"2023-12-11T18:09:26.908126Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The main contradiction between the metric and the competition's goal is that the dataset has a lot of outliers, but the metric is very sensitive to them by definition (for example, for 3 points from zero-mean  Gaussian distribution and 1 outlier the error $ RMSE = (0.99-0.98)^2 + (0.01-0.03)^2 + (170 - 168)^2 = 4.0005 $, which is mostly caused by the outlier). So it pushes any ML model to reward noise prediction instead of signal. Even if we clean the dataset, our prediction will be evaluated using the noisy test set. ","metadata":{}},{"cell_type":"markdown","source":"# Modelling \n\n## Build Baseline Models\n\nTo compare our model with something we need to create a heuristic (baseline) model. \n\nBecause gene distributions are Gaussian-like with zero mean,  the first baseline model might be an array with zeros; here it is `prediction_0()`. \n\nFor evaluating performance, we used a simple 6-fold cross-validation. ","metadata":{}},{"cell_type":"code","source":"def prediction_0(df_train, df_test):\n    \"\"\"\n    Generates zero predictions for a given test dataset.\n\n    Args:\n    - df_train (pandas.DataFrame): Training dataset.\n    - df_test (pandas.DataFrame): Test dataset.\n\n    Returns:\n    numpy.ndarray: Zero predictions for each sample in df_test.\n    \"\"\"\n    \n    y_pred = np.zeros((len(df_test), len(target)))\n    \n    return y_pred\n    \n\ndef cross_validation(df, prediction=prediction_0, n_folds=6, print_folds=True):\n    \"\"\"\n    Performs cross-validation on a given dataset using a specified prediction function.\n\n    Args:\n    - df (pandas.DataFrame): Dataset for cross-validation.\n    - prediction (function): Prediction function to use during cross-validation.\n    - n_folds (int): Number of folds for cross-validation.\n    - print_folds (bool): If True, print MRRMSE for each fold.\n\n    Returns:\n    float: Mean of the MRRMSE values across all folds.\n    \"\"\"\n    kfold = KFold(n_splits=n_folds,random_state=42, shuffle=True)\n\n    metric_av = [] \n    counter = 0 \n\n    for train, test in kfold.split(df):\n        \n        df_train = df.iloc[train]\n        df_test = df.iloc[test]\n\n        y_test = df_test[target]\n        y_pred = prediction(df_train, df_test) \n        \n        metric_av.append(metric(y_test,y_pred))\n        \n        if print_folds:\n            print(f\"Fold #{counter} MRRMSE:\", round(metric(y_test,y_pred),4))\n            counter +=1\n\n    return round(np.mean(metric_av),4)","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:09:26.923704Z","iopub.execute_input":"2023-12-11T18:09:26.924735Z","iopub.status.idle":"2023-12-11T18:09:26.941414Z","shell.execute_reply.started":"2023-12-11T18:09:26.924680Z","shell.execute_reply":"2023-12-11T18:09:26.939742Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Baseline model #1 MRRMSE:\", \n      cross_validation(df, prediction=prediction_0,print_folds=True)\n     )","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:09:26.943956Z","iopub.execute_input":"2023-12-11T18:09:26.944521Z","iopub.status.idle":"2023-12-11T18:09:27.843093Z","shell.execute_reply.started":"2023-12-11T18:09:26.944472Z","shell.execute_reply":"2023-12-11T18:09:27.841930Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The second possible model might be an average over different drugs after denoising (compress the space with 18211 genes and decompress it back). \n\n","metadata":{}},{"cell_type":"code","source":"def tsvd_inv_drug_av(df_train, df_test, reducer):\n    \"\"\"\n    Applies dimensionality reduction and inverse transform to denoise a dataset,\n    aggregates drug information, and predicts target values for a test dataset.\n    The idea is from  https://www.kaggle.com/code/alexandervc/op2-eda-baseline-s \n\n    Args:\n    - df_train (pandas.DataFrame): Training dataset.\n    - df_test (pandas.DataFrame): Test dataset.\n    - reducer: Dimensionality reduction model with 'transform' and 'inverse_transform' methods.\n\n    Returns:\n    pandas.Series: Predicted target values for the test dataset.\n    \"\"\"\n    n_components = 36\n    quantile = 0.54\n    aggr_feature = \"sm_name\"\n\n    X = df[target] \n\n    # Apply dimensionality reduction and inverse transform to denoise the dataset\n    Xr = reducer.transform(X)\n    Xr_inv_trans = reducer.inverse_transform(Xr)\n\n    # Create a DataFrame with denoised data and original categorical feature\n    df_denoised = pd.DataFrame(Xr_inv_trans, columns=target)\n    df_denoised[aggr_feature] = df[aggr_feature]\n\n    # Aggregate drug information using quantile and group by the specified feature\n    df_drug_aggr = df_denoised.groupby(aggr_feature).quantile(quantile)\n\n    # Merge aggregated drug information with test dataset to get predictions\n    y_pred = df_test[features].merge(df_drug_aggr, how='left', on=aggr_feature)[target]\n    \n    return y_pred","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:09:27.844622Z","iopub.execute_input":"2023-12-11T18:09:27.845230Z","iopub.status.idle":"2023-12-11T18:09:27.857394Z","shell.execute_reply.started":"2023-12-11T18:09:27.845181Z","shell.execute_reply":"2023-12-11T18:09:27.856092Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def cross_validation_reduced(df, prediction=prediction_0, n_folds=6,  print_folds=True):\n    \"\"\"\n    Performs cross-validation with dimensionality reduction on a given dataset.\n\n    Args:\n    - df (pandas.DataFrame): Dataset for cross-validation.\n    - prediction (function): Prediction function to use during cross-validation.\n    - n_folds (int): Number of folds for cross-validation.\n    - print_folds (bool): If True, print MRRMSE for each fold.\n\n    Returns:\n    float: Mean of the MRRMSE values across all folds.\n    \"\"\"\n    n_components = 36 \n    kfold = KFold(n_splits=n_folds,random_state=42, shuffle=True)\n    reducer = TruncatedSVD(n_components=n_components,  n_iter=7, random_state=42) \n    \n    reducer.fit(df[target])\n\n    metric_av = [] \n    counter = 0 \n\n    for train, test in kfold.split(df):\n        \n        df_train = df.iloc[train]\n        df_test = df.iloc[test]\n        \n        # To be fully honest, we need to train reducer only over df_train, like \n        # reducer.fit(df_train[target])\n        # but in this case we overestimate the error \n        # (make the test on compounds that were not in train data, which is not the case \n        # for both public and private evaluation of the model)\n        \n        y_pred = prediction(df_train, df_test, reducer)\n        y_test = df_test[target]\n        \n        metric_av.append(metric(y_test,y_pred))\n        \n        if print_folds:\n            print(f\"Fold #{counter} MRRMSE:\", round(metric(y_test,y_pred),4))\n            counter +=1\n\n    return round(np.mean(metric_av),4)","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:09:27.859343Z","iopub.execute_input":"2023-12-11T18:09:27.861114Z","iopub.status.idle":"2023-12-11T18:09:27.877546Z","shell.execute_reply.started":"2023-12-11T18:09:27.861044Z","shell.execute_reply":"2023-12-11T18:09:27.876056Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Model #2 MRRMSE:\", \n      cross_validation_reduced(df, prediction=tsvd_inv_drug_av, print_folds=True)\n     )","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:09:27.879834Z","iopub.execute_input":"2023-12-11T18:09:27.880260Z","iopub.status.idle":"2023-12-11T18:09:46.096322Z","shell.execute_reply.started":"2023-12-11T18:09:27.880228Z","shell.execute_reply":"2023-12-11T18:09:46.094695Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The third model is the average value of compressed DE  for each drug (prediction is the same for all cell types). ","metadata":{}},{"cell_type":"code","source":"def tsvd_drug_av_inv(df_train, df_test, reducer):\n    \"\"\"\n    Applies dimensionality reduction to a dataset, computes the average of each drug's reduced representation,\n    inversely transforms the averaged representation, and predicts target values for a test dataset.\n\n    Args:\n    - df_train (pandas.DataFrame): Training dataset.\n    - df_test (pandas.DataFrame): Test dataset.\n    - reducer: Dimensionality reduction model with 'transform' and 'inverse_transform' methods.\n\n    Returns:\n    pandas.Series: Predicted target values for the test dataset.\n    \"\"\"\n    X = df[target] \n    aggr_feature = \"sm_name\"\n    \n    # Apply dimensionality reduction\n    Xr = reducer.transform(X)\n    \n    n_components = Xr.shape[1]\n    \n    # Get unique drugs from the dataset\n    drugs = df[aggr_feature].unique()\n    \n    # Initialize an array to store the average representation for each drug\n    X_av_drug = np.zeros((len(drugs), n_components))\n    \n    # Create a DataFrame to store averaged drug representations\n    df_ = pd.DataFrame(columns=[aggr_feature] + target)\n    df_[aggr_feature] = drugs\n    \n    # Compute the average representation for each drug\n    for i, d in enumerate(drugs):\n        observations = df[df[aggr_feature] == d].index.to_list()\n        X_av_drug[i, :] = np.mean(Xr[observations, :], axis=0)\n\n    # Inversely transform the averaged drug representations\n    X_av_drug_inv = reducer.inverse_transform(X_av_drug)\n    \n    # Update the DataFrame with the averaged and inversely transformed drug representations\n    df_[target] = X_av_drug_inv.copy()\n    \n    # Merge averaged drug representations with the test dataset to get predictions\n    y_pred = df_test[features].merge(df_, how='left', on=aggr_feature)[target]\n    \n    return y_pred","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:09:46.098211Z","iopub.execute_input":"2023-12-11T18:09:46.098694Z","iopub.status.idle":"2023-12-11T18:09:46.111455Z","shell.execute_reply.started":"2023-12-11T18:09:46.098657Z","shell.execute_reply":"2023-12-11T18:09:46.109713Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Model #3 MRRMSE:\", \n      cross_validation_reduced(df, prediction=tsvd_drug_av_inv, print_folds=True)\n     )","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:09:46.113186Z","iopub.execute_input":"2023-12-11T18:09:46.113566Z","iopub.status.idle":"2023-12-11T18:11:02.484209Z","shell.execute_reply.started":"2023-12-11T18:09:46.113535Z","shell.execute_reply":"2023-12-11T18:11:02.482780Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## XGBoost Demo\n\nTo illustrate the logic of the approach, we provide detailed annotated code.","metadata":{}},{"cell_type":"code","source":"# Set the number of components for TruncatedSVD\nn_components = 36\n\n# Initialize TruncatedSVD with specified parameters\nreducer = TruncatedSVD(n_components=n_components, n_iter=7, random_state=42) #20 and 7 are best\n\n# Initialize Target Encoder with smoothing parameter\nenc = ce.TargetEncoder(smoothing=8) \n    \nX = df[target] \n\naggr_feature = \"sm_name\"\n\nXr = reducer.fit_transform(X)\n\n# Get unique values for 'sm_name' and 'cell_type'\ndrugs = list(df[aggr_feature].unique())  \ncells = list(df[\"cell_type\"].unique())\n\n# Create names for the zip (after TruncatedSVD) and one-hot encoded columns\nzip_names = [i for i in range(n_components)]\n\n\n# Initialize DataFrames to store zip-encoded features and TruncatedSVD components\nX_zip = pd.DataFrame(np.zeros((len(df),len(zip_names))),\n                            columns = zip_names)\n\ny_zip = pd.DataFrame(Xr)\n\n\n# Initialize an index variable to fill the DataFrame with features\nind = 0\nfor i,d in enumerate(drugs):\n    \n    # Create a copy of the dataframe for the current drug\n    df_one_drug = df[df[aggr_feature]==d].copy() \n    observations = df_one_drug.index.to_list()\n    \n    # Get unique cell types for the current drug\n    cells_ = list(df_one_drug[\"cell_type\"].unique())\n    \n    for j,c in enumerate(cells_):\n        \n        # Set one-hot encoded values for drug and cell type\n        X_zip.loc[ind, d] = 1 \n        X_zip.loc[ind, [k for k in cells_ if k!= c]] = 1\n        \n        # Get the indexes of the rest of the cell types\n        # cell_index = df_one_drug[df_one_drug[\"cell_type\"] == c].index.to_list()\n        rest_indexes = df_one_drug[df_one_drug[\"cell_type\"] != c].index.to_list()\n        \n        # Get the features as the mean of TruncatedSVD components for the rest of the cell types\n        X_zip.loc[ind, zip_names] = np.mean(Xr[rest_indexes,:], axis=0)\n\n        ind+=1\n        \n        ","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:11:02.486631Z","iopub.execute_input":"2023-12-11T18:11:02.487632Z","iopub.status.idle":"2023-12-11T18:11:10.485043Z","shell.execute_reply.started":"2023-12-11T18:11:02.487585Z","shell.execute_reply":"2023-12-11T18:11:10.483539Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_categ = df[features]\n\n# Initialize an array to store encoded categorical features based on TruncatedSVD components\nX_categ_trained = np.zeros((len(df),2*n_components))\n\n# Loop through each TruncatedSVD component\nfor n in range(n_components):\n    # Fit Target Encoder on the categorical features using the current TruncatedSVD component\n    cat_encoding_one_component = enc.fit_transform(X_categ, Xr[:,n])\n    \n    # Assign the encoded values to the corresponding columns \n    X_categ_trained[:,2*n] = cat_encoding_one_component[\"cell_type\"]\n    X_categ_trained[:,2*n+1] = cat_encoding_one_component[\"sm_name\"]\n\n# Create a DataFrame from the encoded categorical features\nX_categ_ = pd.DataFrame(X_categ_trained, columns=[f\"c_{i}\" for i in range(2*n_components)])\n\n# Concatenate the aggregated by drug features and the encoded categorical features\nX_zip_ = pd.concat([X_zip, X_categ_], axis=1)\n\nX_zip_ ","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:11:10.486851Z","iopub.execute_input":"2023-12-11T18:11:10.487302Z","iopub.status.idle":"2023-12-11T18:11:11.468056Z","shell.execute_reply.started":"2023-12-11T18:11:10.487262Z","shell.execute_reply":"2023-12-11T18:11:11.467129Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.model_selection import train_test_split\n\n\nX_train, X_test, y_train, y_test = train_test_split(\n    X_zip_, y_zip, test_size=0.15, random_state=22 \n)","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:11:11.469371Z","iopub.execute_input":"2023-12-11T18:11:11.470342Z","iopub.status.idle":"2023-12-11T18:11:11.488432Z","shell.execute_reply.started":"2023-12-11T18:11:11.470307Z","shell.execute_reply":"2023-12-11T18:11:11.486899Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"params = {\n'tree_method': 'hist',\n'eval_metric': 'rmse',\n'max_depth': 2, \n'learning_rate': 0.2,\n'silent': 1.0,      \n'n_estimators': 500,\n'random_state': 42}\n\nxgb_regressor = xgb.XGBRegressor(**params)\n\n\nxgb_regressor.fit(X_train,y_train) \n\ny_pred = xgb_regressor.predict(X_test)\n\ny_pred = pd.DataFrame(y_pred, index=X_test.index)","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:11:11.489885Z","iopub.execute_input":"2023-12-11T18:11:11.491076Z","iopub.status.idle":"2023-12-11T18:11:34.241964Z","shell.execute_reply.started":"2023-12-11T18:11:11.491028Z","shell.execute_reply":"2023-12-11T18:11:34.241024Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Display a single observation for average by drugs, xgboost predictions, and test compressed data ","metadata":{}},{"cell_type":"code","source":"j = 6\nj_ = X_test.index[j]\n\nplt.plot(zip_names, X_test.loc[j_,zip_names],'--', label=\"average by drugs\")\nplt.plot(zip_names, y_pred.loc[j_], label=\"prediction\")\nplt.plot(zip_names, y_test.loc[j_] , label=\"test\")\nplt.xlabel(\"Index of Singular Values\")\nplt.ylabel(\"Compressed DE\")\n\nplt.legend()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:11:34.243252Z","iopub.execute_input":"2023-12-11T18:11:34.244302Z","iopub.status.idle":"2023-12-11T18:11:34.593715Z","shell.execute_reply.started":"2023-12-11T18:11:34.244264Z","shell.execute_reply":"2023-12-11T18:11:34.592212Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Display values of a single component for average by drugs, xgboost predictions, and test compressed data. ","metadata":{}},{"cell_type":"code","source":"j = 4\n\nplt.plot(range(len(X_test)), X_test[j],'--', label=\"average by drugs\")\nplt.plot(range(len(X_test)), y_pred[j], label=\"prediction\")\nplt.plot(range(len(X_test)), y_test[j] , label=\"test\")\nplt.xlabel(\"Row Number in Test Set\")\nplt.ylabel(\"Compressed DE\")\nplt.legend()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:11:34.595833Z","iopub.execute_input":"2023-12-11T18:11:34.596455Z","iopub.status.idle":"2023-12-11T18:11:34.883248Z","shell.execute_reply.started":"2023-12-11T18:11:34.596423Z","shell.execute_reply":"2023-12-11T18:11:34.881774Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Test for the decompressed result.","metadata":{}},{"cell_type":"code","source":"\nY_pred = pd.DataFrame(reducer.inverse_transform(y_pred), columns=target, index=X_test.index)\n\nY_test = df.loc[y_test.index, target]\n\nprint(\"MRRMSE:\", round(metric(Y_test, Y_pred),4))","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:11:34.885902Z","iopub.execute_input":"2023-12-11T18:11:34.887028Z","iopub.status.idle":"2023-12-11T18:11:35.067171Z","shell.execute_reply.started":"2023-12-11T18:11:34.886968Z","shell.execute_reply":"2023-12-11T18:11:35.065697Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Display a single observation for xgboost prediction and test data (first 200 genes)\n","metadata":{}},{"cell_type":"code","source":"j = 6\nj_ = X_test.index[j]\n\nplt.plot(range(200), Y_test.loc[j_,target[:200]],'--', label=\"test\")\nplt.plot(range(200), Y_pred.loc[j_,target[:200]], label=\"prediction\")\nplt.xlabel(\"Gene Number\")\nplt.ylabel(\"DE\")\n\nplt.legend()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:11:35.069197Z","iopub.execute_input":"2023-12-11T18:11:35.069732Z","iopub.status.idle":"2023-12-11T18:11:35.486152Z","shell.execute_reply.started":"2023-12-11T18:11:35.069682Z","shell.execute_reply":"2023-12-11T18:11:35.484794Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":" Display values of a single gene for xgboost prediction and test data","metadata":{}},{"cell_type":"code","source":"g = 1\n\nplt.plot(range(len(Y_test)), Y_test.loc[:,target[g]],'--', label=\"test\")\nplt.plot(range(len(Y_test)), Y_pred.loc[:, target[g]], label=\"prediction\")\nplt.xlabel(\"Row Number in Test Set\")\nplt.ylabel(\"DE\")\n\nplt.legend()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:11:35.487833Z","iopub.execute_input":"2023-12-11T18:11:35.488272Z","iopub.status.idle":"2023-12-11T18:11:35.846729Z","shell.execute_reply.started":"2023-12-11T18:11:35.488237Z","shell.execute_reply":"2023-12-11T18:11:35.844904Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"xgb.plot_importance(xgb_regressor, max_num_features=10); # top 10 most important features","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:11:35.849836Z","iopub.execute_input":"2023-12-11T18:11:35.850543Z","iopub.status.idle":"2023-12-11T18:11:36.371651Z","shell.execute_reply.started":"2023-12-11T18:11:35.850482Z","shell.execute_reply":"2023-12-11T18:11:36.369906Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Comparison with the average drugs prediction in the compressed space","metadata":{}},{"cell_type":"code","source":"df_train = df.iloc[X_train.index]\ndf_test = df.iloc[X_test.index]\n\ny_pred_av = tsvd_drug_av_inv(df_train,df_test, reducer)\n\nprint(\"MRRMSE:\", metric(Y_test, y_pred_av))","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:11:36.374874Z","iopub.execute_input":"2023-12-11T18:11:36.375512Z","iopub.status.idle":"2023-12-11T18:11:48.819185Z","shell.execute_reply.started":"2023-12-11T18:11:36.375456Z","shell.execute_reply":"2023-12-11T18:11:48.817669Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"j = 6\nj_ = X_test.index[j]\n\nplt.plot(range(200), Y_test.loc[j_,target[:200]],'--', label=\"test\")\nplt.plot(range(200), y_pred_av.loc[j_,target[:200]], label=\"average by drugs\")\n\nplt.legend()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:11:48.821016Z","iopub.execute_input":"2023-12-11T18:11:48.822321Z","iopub.status.idle":"2023-12-11T18:11:49.150667Z","shell.execute_reply.started":"2023-12-11T18:11:48.822270Z","shell.execute_reply":"2023-12-11T18:11:49.149023Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"g = 1\n\nplt.plot(range(len(Y_test)), Y_test.loc[:,target[g]],'--', label=\"test\")\nplt.plot(range(len(Y_test)), y_pred_av.loc[:, target[g]], label=\"average by drugs\")\n\nplt.legend()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:11:49.152998Z","iopub.execute_input":"2023-12-11T18:11:49.153960Z","iopub.status.idle":"2023-12-11T18:11:49.502901Z","shell.execute_reply.started":"2023-12-11T18:11:49.153891Z","shell.execute_reply":"2023-12-11T18:11:49.499515Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## XGBoost Training\n\nTo summarize your approach, we provide the cross-validation and related functions of it below.\nWe also added one-hot encoding because it gives better CV, public, and private scores.","metadata":{}},{"cell_type":"code","source":"# from https://www.kaggle.com/code/ambrosm/scp-eda-which-makes-sense\ndef de_to_t_score(de):\n    \"\"\"Convert log10pvalues to t-scores\n    \n    Parameter:\n    de: array or DataFrame of log10pvalues\n    \n    Return value:\n    t_score: array or DataFrame of t-scores\n    \"\"\"\n    p_value = 10 ** (-np.abs(de))\n    return - scipy.stats.t.ppf(p_value / 2, df=420) * np.sign(de)\n#     return - norm.ppf(p_value / 2) * np.sign(de)\n\n\n# from https://www.kaggle.com/code/ambrosm/scp-eda-which-makes-sense\ndef t_score_to_de(t_score):\n    \"\"\"Convert t-scores to log10pvalues (inverse of de_to_t_score)\n    \n    Parameter:\n    t_score: array or DataFrame of t-scores\n    \n    Return value:\n    de: array or DataFrame of log10pvalues\n    \"\"\"\n    p_value = scipy.stats.t.cdf(- np.abs(t_score), df=420) * 2\n#     p_value = norm.cdf(- np.abs(t_score)) * 2\n    p_value = p_value.clip(1e-180, None)\n    return - np.log10(p_value) * np.sign(t_score)\n\n\ndef obtain_reduced_data(df, reducer, encoder, t_score = False):\n    \n    \"\"\"\n    Obtains reduced data for a given DataFrame using dimensionality reduction and encoding techniques.\n\n    Args:\n        df (pd.DataFrame): Input DataFrame containing the data.\n        reducer: An instance of a dimensionality reduction model.\n        encoder: An instance of an encoding model for categorical features.\n\n    Returns:\n        pd.DataFrame: Reduced and encoded feature matrix.\n        pd.DataFrame: Original reduced feature matrix.\n    \"\"\"\n    \n    # Extracting the target variable for dimensionality reduction\n    if t_score:\n        X = de_to_t_score(df[target])\n    else: \n        X = df[target] \n\n    # Defining the aggregation feature for grouping\n    aggr_feature = \"sm_name\"\n    \n    # Transforming the data using the provided dimensionality reduction technique\n    Xr = reducer.transform(X)\n    \n    # Extracting the number of components after dimensionality reduction\n    n_components = Xr.shape[1]\n    \n    # Extracting unique drug and cell type names\n    drugs = list(df[aggr_feature].unique())\n    cells = list(df[\"cell_type\"].unique())\n    \n    one_hot_names = drugs + cells \n    \n    # Creating names for zip columns\n    zip_names = [i for i in range(n_components)]\n\n    X_zip = pd.DataFrame(np.zeros((len(df),len(one_hot_names)+len(zip_names))),\n                                columns = one_hot_names + zip_names)\n    \n    # Initializing DataFrames for the reduced and encoded features\n    y_zip = pd.DataFrame(Xr)\n    \n    # Looping through drugs and cells to calculate the mean of the cell types for \n    # each drug in compressed space\n    ind = 0\n    for i,d in enumerate(drugs):\n\n        df_one_drug = df[df[aggr_feature]==d].copy()\n        observations = df_one_drug.index.to_list()\n\n        cells_ = list(df_one_drug[\"cell_type\"].unique())\n\n        for j,c in enumerate(cells_):\n            \n            X_zip.loc[ind, d] = 1 \n            X_zip.loc[ind, [k for k in cells_ if k!= c]] = 1\n   \n            rest_indexes = df_one_drug[df_one_drug[\"cell_type\"] != c].index.to_list()\n\n            X_zip.loc[ind, zip_names] = np.mean(Xr[rest_indexes,:], axis=0)\n\n            ind+=1\n            \n    # Extracting categorical features for encoding\n    X_categ = df[features]\n    \n    # Initializing the array for the encoded categorical features\n    X_categ_trained = np.zeros((len(df),2*n_components))\n    \n    # Looping through components for encoding categorical features\n    for n in range(n_components):\n\n        cat_encoding_one_component = encoder.fit_transform(X_categ, Xr[:,n])\n\n        X_categ_trained[:,2*n] = cat_encoding_one_component[\"cell_type\"]\n        X_categ_trained[:,2*n+1] = cat_encoding_one_component[\"sm_name\"]\n    \n    # Concatenating the reduced and encoded features into a single DataFrame\n    X_zip_ = pd.concat([X_zip ,pd.DataFrame(X_categ_trained, columns=[f\"c_{i}\" for i in range(2*n_components)])], axis=1)\n        \n    return X_zip_, y_zip\n    \ndef xgb_prediction(X_train, y_train, X_test, reducer):\n    \"\"\"\n    Perform XGBoost regression for dimensionality-reduced data.\n\n    Args:\n        X_train (pd.DataFrame): Training data features.\n        y_train (pd.Series): Training data labels.\n        X_test (pd.DataFrame): Test data features.\n        reducer: An instance of a dimensionality reduction model.\n    Note:\n        The function uses XGBoost to predict values in the original feature space\n        for the given test data, assuming dimensionality reduction has been applied\n        using the provided reducer.\n    \"\"\"\n    \n    params = {\n            'tree_method': 'hist',\n            'eval_metric': 'rmse',\n            'max_depth': 2, \n            'learning_rate': 0.2,\n            'silent': 1.0,      \n            'n_estimators': 1000,\n            'random_state': 42}\n\n    xgb_regressor = xgb.XGBRegressor(**params)\n\n    xgb_regressor.fit(X_train,y_train) \n    y_pred = xgb_regressor.predict(X_test)\n    \n    Y_pred = pd.DataFrame(reducer.inverse_transform(y_pred), columns=target, index=X_test.index)\n\n    return Y_pred\n\ndef cross_validation_xgb(df, prediction=xgb_prediction, n_folds=6, n_components = 36, print_folds=True, t_score=False):\n    \"\"\"\n    Perform cross-validation using XGBoost with dimensionality reduction.\n\n    Args:\n    - df (pd.DataFrame): The input DataFrame containing features and target variable.\n    - prediction (function): The prediction function (default is xgb_prediction).\n    - n_folds (int): Number of folds for cross-validation (default is 6).\n    - n_components (int): Number of components for dimensionality reduction (default is 36).\n    - print_folds (bool): Flag to print fold-wise results (default is True).\n\n    Returns:\n    - float: Mean metric value across all folds.\n    \"\"\"\n    # Set up KFold for cross-validation    \n    kfold = KFold(n_splits=n_folds,random_state=42, shuffle=True)\n    \n    # Dimensionality reduction using TruncatedSVD\n    reducer = TruncatedSVD(n_components=n_components, n_iter=7, random_state=42) #20 and 7 are the best on CV\n    \n    # Target encoding using TargetEncoder\n    encoder = ce.TargetEncoder(smoothing=8) # 8 and 12 \n    \n    # Extract target variables (genes) from the DataFrame\n    if t_score:\n        X = de_to_t_score(df[target]) \n    else:\n        X = df[target] \n    \n    # Fit reducer on the target variables\n    reducer.fit(X)\n    \n    # Obtain reduced data using dimensionality reduction and target encoding\n    X_zip, y_zip = obtain_reduced_data(df, reducer, encoder, t_score=t_score)\n\n    # Initialize lists for storing metric values and a counter for fold number\n    metric_av = [] \n    counter = 0 \n    \n    # Loop through each fold in KFold\n    for train, test in kfold.split(X_zip):\n        # Split data into training and testing sets\n        X_train_fold = X_zip.iloc[train] \n        y_train_fold = y_zip.iloc[train]\n        \n        X_test_fold = X_zip.iloc[test]\n        y_test_fold = df[target].iloc[test].copy()\n        \n        # Make predictions using the specified prediction function\n        y_pred = prediction(X_train_fold, y_train_fold, X_test_fold, reducer)\n        \n        if t_score:\n            y_pred = t_score_to_de(y_pred)\n            \n        # Calculate and store the metric value for the fold\n        metric_av.append(metric(y_test_fold,y_pred))\n        \n        # Print fold-wise results if specified\n        if print_folds:\n            print(f\"Fold #{counter} MRRMSE:\", round(metric(y_test_fold,y_pred),4))\n            counter +=1\n    # Return the mean of the metric values across all folds\n    return round(np.mean(metric_av),4)","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:11:49.506042Z","iopub.execute_input":"2023-12-11T18:11:49.507393Z","iopub.status.idle":"2023-12-11T18:11:49.551518Z","shell.execute_reply.started":"2023-12-11T18:11:49.507318Z","shell.execute_reply":"2023-12-11T18:11:49.549970Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"XGB Model MRRMSE:\", \n      cross_validation_xgb(df,  print_folds=True)\n  )","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:11:49.553540Z","iopub.execute_input":"2023-12-11T18:11:49.554444Z","iopub.status.idle":"2023-12-11T18:18:21.171490Z","shell.execute_reply.started":"2023-12-11T18:11:49.554394Z","shell.execute_reply":"2023-12-11T18:18:21.169107Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"*Is it better to use a t-score instead of DE for our model? [1]*. It gives a better score on CV, but worse on the public and private test. ","metadata":{}},{"cell_type":"markdown","source":"### Cell type validation \n\nDoes the accuracy of our model depend on cell types in the dataset? To answer this, we compute the metric for cell types excluding one of 4 cell types from the training set (\"NK cells\", \"T cells CD4+\", \"T cells CD8+\", \"T regulatory cells\"). The score for **T cells CD8+** should be a little bit lower than the ground truth, because we have only 142 observations instead of 146 (see above). \n","metadata":{}},{"cell_type":"code","source":"def validate_cells(df, prediction=xgb_prediction, n_components = 36, print_folds=True):\n    \"\"\"\n    Perform cross-validation of a model using for training dataset without one cell type from a list\n    [\"NK cells\", \"T cells CD4+\", \"T cells CD8+\", \"T regulatory cells\"]\n\n    Args:\n    - df (pd.DataFrame): The input DataFrame containing features and target variable.\n    - prediction (function): The prediction function (default is xgb_prediction).\n    - n_components (int): Number of components for dimensionality reduction (default is 36).\n    - print_folds (bool): Flag to print fold-wise results (default is True).\n\n    Returns:\n    - float: Mean metric value across all folds.\n    \"\"\"\n    # Set up KFold for cross-validation    \n    kfold = KFold(n_splits=5,random_state=42, shuffle=True)\n    \n    # Dimensionality reduction using TruncatedSVD\n    reducer = TruncatedSVD(n_components=n_components, n_iter=7, random_state=42) #20 and 7 are the best on CV\n    \n    # Target encoding using TargetEncoder\n    encoder = ce.TargetEncoder(smoothing=8) # 8 and 12 \n    \n    # Extract target variables (genes) from the DataFrame\n    X = df[target] \n    \n    cell_types = [\"NK cells\", \"T cells CD4+\", \"T cells CD8+\", \"T regulatory cells\"]\n    \n    cell_exclude_indexes = [ df[df[\"cell_type\"] == c].index.to_list() for c in cell_types ]\n    \n    # Fit reducer on the target variables\n    reducer.fit(X)\n    \n    # Obtain reduced data using dimensionality reduction and target encoding\n    X_zip, y_zip = obtain_reduced_data(df, reducer, encoder)\n\n    # Initialize lists for storing metric values and a counter for fold number\n    metric_av = [] \n    \n    counter = 0\n    \n    for c in cell_exclude_indexes:\n        \n        metric_cell_av = []\n        X_exclude_cell = X_zip.drop(c)\n        \n        for train, test in kfold.split(X_exclude_cell):\n            # Split data into training and testing sets\n            X_train_fold = X_zip.iloc[train] \n            y_train_fold = y_zip.iloc[train]\n\n            X_test_fold = X_zip.iloc[test]\n            y_test_fold = df[target].iloc[test].copy()\n\n            # Make predictions using the specified prediction function\n            y_pred = prediction(X_train_fold, y_train_fold, X_test_fold, reducer)\n\n            # Calculate and store the metric value for the fold\n            metric_cell_av.append(metric(y_test_fold,y_pred))\n        \n        metric_av.append(np.mean(metric_cell_av))\n\n        # Print fold-wise results if specified\n            \n        if print_folds:\n            print(f\"{cell_types[counter]} excluded MRRMSE:\", round(metric_av[counter],4))\n            counter +=1\n            \n    # Return the mean of the metric values across all cells\n    return round(np.mean(metric_av),4)","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:18:21.181565Z","iopub.execute_input":"2023-12-11T18:18:21.182706Z","iopub.status.idle":"2023-12-11T18:18:21.199226Z","shell.execute_reply.started":"2023-12-11T18:18:21.182660Z","shell.execute_reply":"2023-12-11T18:18:21.197886Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"XGB Model MRRMSE:\", \n      validate_cells(df,  print_folds=True)\n  )","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:18:21.201234Z","iopub.execute_input":"2023-12-11T18:18:21.202526Z","iopub.status.idle":"2023-12-11T18:38:02.255181Z","shell.execute_reply.started":"2023-12-11T18:18:21.202466Z","shell.execute_reply":"2023-12-11T18:38:02.252728Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Notably, the accuracy of the model differs only at T cells CD8+ exclusion, and not significantly. We believe that it happens because our prediction is based on the average of cell types in compressed space, therefore training on data with more cell types gives a better prediction. ","metadata":{}},{"cell_type":"markdown","source":"# Submission","metadata":{}},{"cell_type":"code","source":"def drug_av_zip(df, idtab, reducer, encoder, t_score=False):\n    \"\"\"\n    Create a feature encoding, applying a drug-specific averaging in compressed space to a given dataframe.\n\n    Parameters:\n    - df (pd.DataFrame): Input training dataframe.\n    - idtab (pd.DataFrame): Reference dataframe for encoding features.\n    - reducer: A dimensionality reduction model used for compression.\n\n    Returns:\n    pd.DataFrame: A new dataframe containing feature encoding.\n    \"\"\"\n    \n    # Extracting the target variable for dimensionality reduction\n    if t_score:\n        X = de_to_t_score(df[target])\n    else:\n        X = df[target] \n    \n    # Defining the aggregation feature for grouping\n    arrg_feature = \"sm_name\"\n    \n    # Transforming the data using the provided dimensionality reduction technique\n    Xr = reducer.transform(X)\n    \n    # Extracting the number of components after dimensionality reduction\n    n_components = Xr.shape[1]\n    \n    # Extracting unique drug names\n    drugs = list(df[arrg_feature].unique())\n    \n    X_av_drug = np.zeros((len(drugs), n_components))\n    \n    # Creating names for zip columns\n    zip_names = [i for i in range(n_components)]\n    \n    df_ = pd.DataFrame( columns= [arrg_feature] + zip_names)\n    \n    df_[arrg_feature] = drugs\n    \n    # Looping through drugs to compute the mean of the cell types for \n    # each drug in compressed space\n    \n    for i,d in enumerate(drugs):\n        observations = df[df[arrg_feature]==d].index.to_list()\n    \n        X_av_drug[i,:] =  np.mean(Xr[observations,:], axis=0)\n\n        \n    df_[zip_names] =  X_av_drug.copy()\n    \n    X_zip = idtab[features].merge(df_, how='left', on=\"sm_name\")\n            \n    #---------- categorical features encoding --------------------------#\n    X_categ = idtab[features]\n\n    X_categ_trained = np.zeros((len(idtab),2*n_components))\n    \n    # Looping through components for encoding categorical features\n    for n in range(n_components):\n\n        enc.fit(df[features], Xr[:,n])\n\n        cat_encoding_one_component = enc.transform(X_categ, X_zip[n])\n\n        X_categ_trained[:,2*n] = cat_encoding_one_component[\"cell_type\"]\n        X_categ_trained[:,2*n+1] = cat_encoding_one_component[\"sm_name\"]\n\n     # Concatenating the reduced and encoded features into a single DataFrame\n    X_zip_ =  pd.concat([X_zip, pd.DataFrame(X_categ_trained,  columns=[f\"c_{i}\" for i in range(2*n_components)])], axis=1) \n    \n    one_hot_names = drugs + cells\n    \n    # This loop iterates over each element in one_hot_names and inserts a new \n    # column with that name filled with zeros into the DataFrame X_zip_. \n    # The \"+2\" in the insert function suggests that the new columns are inserted \n    # starting from the third column.\n    for i in range(len(one_hot_names)):\n        X_zip_.insert(i+2, one_hot_names[i], np.zeros(X_zip_.shape[0]))\n        \n    # One-Hot Encoding for Drugs and Cell Types\n    \n    # It then sets the value in the DataFrame at the intersection of the current\n    # row and the column with the name d to 1, indicating the presence of that drug.\n    # Similarly, it sets the values in columns corresponding to cell types not equal \n    # to c for the same drug d to 1, indicating the presence of that drug for those other cell types.\n    for ind in X_zip_.index:\n        d = X_zip_.loc[ind, \"sm_name\"]\n        X_zip_.loc[ind, d] = 1 \n        \n        c = X_zip_.loc[ind, \"cell_type\"]\n        \n        cells_in_train = df[(df[\"cell_type\"]!=c) & (df[\"sm_name\"]==d)][\"cell_type\"]\n        X_zip_.loc[ind, cells_in_train] = 1 \n        \n    X_zip_.drop(columns=features, inplace=True)\n    \n    return X_zip_","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:38:02.259187Z","iopub.execute_input":"2023-12-11T18:38:02.259956Z","iopub.status.idle":"2023-12-11T18:38:02.303379Z","shell.execute_reply.started":"2023-12-11T18:38:02.259881Z","shell.execute_reply":"2023-12-11T18:38:02.300871Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fn = '/kaggle/input/open-problems-single-cell-perturbations/id_map.csv'\ndf_id_map = pd.read_csv(fn)\n\ndf_id_map","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:38:02.305446Z","iopub.execute_input":"2023-12-11T18:38:02.306261Z","iopub.status.idle":"2023-12-11T18:38:02.355738Z","shell.execute_reply.started":"2023-12-11T18:38:02.306214Z","shell.execute_reply":"2023-12-11T18:38:02.354610Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"n_components = 36\nreducer = TruncatedSVD(n_components=n_components, n_iter=7, random_state=42)\nencoder = ce.TargetEncoder(smoothing=8)\nt_score = False\n\nif t_score:\n    X = de_to_t_score(df[target])\nelse: \n    X = df[target] \n\nreducer.fit(X)\n\nX_test = drug_av_zip(df, df_id_map, reducer, encoder, t_score=t_score)\n\nX_zip, y_zip = obtain_reduced_data(df, reducer, encoder, t_score=t_score)\n ","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:38:02.357554Z","iopub.execute_input":"2023-12-11T18:38:02.358306Z","iopub.status.idle":"2023-12-11T18:38:12.765436Z","shell.execute_reply.started":"2023-12-11T18:38:02.358262Z","shell.execute_reply":"2023-12-11T18:38:12.764090Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\nY_pred = xgb_prediction(X_zip, y_zip, X_test, reducer) \n\nif t_score:\n    Y_pred = t_score_to_de(Y_pred)","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:38:12.767456Z","iopub.execute_input":"2023-12-11T18:38:12.768374Z","iopub.status.idle":"2023-12-11T18:39:15.787759Z","shell.execute_reply.started":"2023-12-11T18:38:12.768336Z","shell.execute_reply":"2023-12-11T18:39:15.786240Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"Y_pred.index.name = 'id'\nY_pred","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:39:15.789926Z","iopub.execute_input":"2023-12-11T18:39:15.791080Z","iopub.status.idle":"2023-12-11T18:39:15.863243Z","shell.execute_reply.started":"2023-12-11T18:39:15.791019Z","shell.execute_reply":"2023-12-11T18:39:15.861737Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y_pred_av = tsvd_drug_av_inv(df, df_id_map[features], reducer)\ny_pred_av = pd.DataFrame(y_pred_av, columns=target)\ny_pred_av.index.name = 'id'","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:39:15.865544Z","iopub.execute_input":"2023-12-11T18:39:15.866415Z","iopub.status.idle":"2023-12-11T18:39:28.240310Z","shell.execute_reply.started":"2023-12-11T18:39:15.866345Z","shell.execute_reply":"2023-12-11T18:39:28.239173Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"j= 0\n\n\nplt.plot(range(200), Y_pred.loc[j,target[:200]], label=\"prediction XGB\")\nplt.plot(range(200), y_pred_av.loc[j,target[:200]],'--', label=\"prediction average\")\n\nplt.legend()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:39:28.241888Z","iopub.execute_input":"2023-12-11T18:39:28.242251Z","iopub.status.idle":"2023-12-11T18:39:28.717835Z","shell.execute_reply.started":"2023-12-11T18:39:28.242220Z","shell.execute_reply":"2023-12-11T18:39:28.716852Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"g = 1\n\nplt.plot(range(len(Y_pred)), Y_pred.loc[:, target[g]], label=\"prediction XGB\")\nplt.plot(range(len(y_pred_av)), y_pred_av.loc[:,target[g]],'--', label=\"prediction average\")\n\nplt.legend()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:39:28.719533Z","iopub.execute_input":"2023-12-11T18:39:28.720219Z","iopub.status.idle":"2023-12-11T18:39:29.003432Z","shell.execute_reply.started":"2023-12-11T18:39:28.720182Z","shell.execute_reply":"2023-12-11T18:39:29.001378Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"Y_pred.to_csv('submission.csv')","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:39:29.006616Z","iopub.execute_input":"2023-12-11T18:39:29.007866Z","iopub.status.idle":"2023-12-11T18:39:42.234688Z","shell.execute_reply.started":"2023-12-11T18:39:29.007784Z","shell.execute_reply":"2023-12-11T18:39:42.233293Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!ls","metadata":{"execution":{"iopub.status.busy":"2023-12-11T18:39:42.236330Z","iopub.execute_input":"2023-12-11T18:39:42.236674Z","iopub.status.idle":"2023-12-11T18:39:43.429556Z","shell.execute_reply.started":"2023-12-11T18:39:42.236646Z","shell.execute_reply":"2023-12-11T18:39:43.427720Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Conclusions \n\nIn this project, we employed an XGBoost model to address a high-dimensional multi-output regression problem to predict the differential expression of 18,211 genes across 6 cell types affected by 146 chemical substances. Our approach stands out for its computational efficiency and robustness to noise in the dataset, owing to its operation within a compressed space. Simultaneously, it is important to note that the accuracy of our method is constrained by information loss due to compression. \n\nAs demonstrated, the impact of compounds on gene expression in a target cell type can be estimated by leveraging the average values for the other cells within the compressed space. This approximation can serve as a baseline model and be used for training advanced models in subsequent experiments, and it doesn't require any prior biological knowledge.\n\n","metadata":{}},{"cell_type":"markdown","source":"### Scores\n\n| Model               | Public Score | Private Score |  CV Score |  \n|---------------------|--------------|---------------|-----------|\n| Baseline model #1   | 0.666        | 0.902         |  1.2818   |\n| Model #2            | 0.624        | 0.818         |  1.0760   |\n| Model #3            | 0.615        | 0.820         |  1.0717   |\n| Model #4            | 0.594        | 0.777         |  0.9939   |\n\n\n\n**Baseline model #1:** `y_pred` is all zeros\n\n**Model #2:** TruncatedSVD (35 components) -> Inverse TruncatedSVD -> Aggregate Drugs -> .quantile(0.54)\n\n**Model #3:** TruncatedSVD (36 components) -> Average for each drug ->  inverse TruncatedSVD \n\n**Model #4:** TruncatedSVD (n_components=36, n_iter=7) ->  TargetEncoder (smooting=8)  -> XGBoost (1000 estimators, feature augmentation with average for each drug) -> inverse TruncatedSVD \n","metadata":{}},{"cell_type":"markdown","source":"## References \n\nWe appreciate the authors of the following notebooks for sharing their work\n\n[1] https://www.kaggle.com/code/ambrosm/scp-eda-which-makes-sense/notebook\n\n[2] https://www.kaggle.com/code/alexandervc/op2-eda-baseline-s\n\n[3] https://www.kaggle.com/code/alexandervc/op2-models-cv-tuning\n\n[4] https://www.kaggle.com/code/ambrosm/scp-quickstart\n\nDescriptions of the functions were created with help of ChatGPT.","metadata":{}}]}