{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"**Background**\n- geneactiviy (from [Episcanpy](https://episcanpy.readthedocs.io/en/latest/)) is defined as summing the number of open features (windows, peaks, etc) overlapping genes (gene bodies + 5kb upstream of the TSS).\n- geneactivity is probably enough to predict mRNA expression from chromatin accessibility information, and we can substantially reduce features from ~250,000 as open chromatin to ~20,000 as gene.\n- codes for making geneactivity matrix are shared in [Train](https://www.kaggle.com/code/masato114/chromatin-gene-conversion-train) and [Test](https://www.kaggle.com/code/masato114/chromatin-gene-conversion-test). I am sorry for having splitted Train/Test notebooks because of long running time.\n<br></br>\n- in this notebook, I check the usefulness of geneactivity matrix by replacing open-chromatin inputs data while referencing [ambrosm's great notebook](https://www.kaggle.com/code/ambrosm/msci-multiome-quickstart)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"markdown","source":"# Imports","metadata":{}},{"cell_type":"code","source":"!pip install -q tables","metadata":{"execution":{"iopub.status.busy":"2022-08-31T05:47:46.602259Z","iopub.execute_input":"2022-08-31T05:47:46.602741Z","iopub.status.idle":"2022-08-31T05:48:01.181446Z","shell.execute_reply.started":"2022-08-31T05:47:46.602637Z","shell.execute_reply":"2022-08-31T05:48:01.179907Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nimport gc\nfrom tqdm.auto import tqdm\n\nimport numpy as np\nimport pandas as pd\nfrom sklearn.base import BaseEstimator, TransformerMixin\nfrom sklearn.model_selection import KFold\nfrom sklearn.preprocessing import StandardScaler, scale\nfrom sklearn.decomposition import PCA\nfrom sklearn.dummy import DummyRegressor\nfrom sklearn.pipeline import make_pipeline, Pipeline\nfrom sklearn.linear_model import Ridge, LinearRegression\nfrom sklearn.metrics import mean_squared_error\n%matplotlib inline\nimport matplotlib.pyplot as plt\nfrom matplotlib.ticker import MaxNLocator\nfrom matplotlib_venn import venn2\nfrom colorama import Fore, Back, Style","metadata":{"execution":{"iopub.status.busy":"2022-08-31T05:49:57.343683Z","iopub.execute_input":"2022-08-31T05:49:57.344772Z","iopub.status.idle":"2022-08-31T05:49:57.369450Z","shell.execute_reply.started":"2022-08-31T05:49:57.344717Z","shell.execute_reply":"2022-08-31T05:49:57.368431Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Let's check geneactivity matrix","metadata":{}},{"cell_type":"code","source":"train_geneactivity = pd.read_hdf('../input/open-problems-train-geneactivity/train_multi_geneactivity.h5', start=0, stop=5)\ntrain_geneactivity.head()","metadata":{"execution":{"iopub.status.busy":"2022-08-31T05:49:06.108095Z","iopub.execute_input":"2022-08-31T05:49:06.108728Z","iopub.status.idle":"2022-08-31T05:49:06.278375Z","shell.execute_reply.started":"2022-08-31T05:49:06.108652Z","shell.execute_reply":"2022-08-31T05:49:06.277204Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# in order to avoid OOM, I use only 6000 train data\ntrain_targets = pd.read_hdf('../input/open-problems-multimodal/train_multi_targets.h5', start=0, stop=6000)\ntrain_targets.head()","metadata":{"execution":{"iopub.status.busy":"2022-08-31T05:51:33.041716Z","iopub.execute_input":"2022-08-31T05:51:33.042158Z","iopub.status.idle":"2022-08-31T05:51:38.210372Z","shell.execute_reply.started":"2022-08-31T05:51:33.042103Z","shell.execute_reply":"2022-08-31T05:51:38.209202Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# check gene overlaps\ndiff_a = set(train_geneactivity.columns) - set(train_targets.columns)\ndiff_b = set(train_targets.columns) - set(train_geneactivity.columns)\ninter = set(train_geneactivity.columns) & set(train_targets.columns)\nvenn2(subsets = (len(diff_a), len(diff_b), len(inter)), set_labels = ('Genes in Geneactivity', 'Genes in Targets'))\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-08-31T05:53:57.803592Z","iopub.execute_input":"2022-08-31T05:53:57.804533Z","iopub.status.idle":"2022-08-31T05:53:57.996754Z","shell.execute_reply.started":"2022-08-31T05:53:57.804485Z","shell.execute_reply":"2022-08-31T05:53:57.995471Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"diff_b_mean = train_targets[list(diff_b)].mean(axis=0).values\ninter_mean = train_targets[list(inter)].mean(axis=0).values\nprint(f'mean expression of overlap genes: {np.mean(inter_mean):.4f} and genes only in targets: {np.mean(diff_b_mean):.4f}')\nfig = plt.figure(figsize=(6,6))\nax = fig.add_subplot()\nax.boxplot([inter_mean, diff_b_mean], labels=['overlap genes', 'genes only in targets'])\nax.set_ylabel('mean gene expression')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-08-31T05:58:45.537592Z","iopub.execute_input":"2022-08-31T05:58:45.538068Z","iopub.status.idle":"2022-08-31T05:58:46.281341Z","shell.execute_reply.started":"2022-08-31T05:58:45.538021Z","shell.execute_reply":"2022-08-31T05:58:46.280061Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"1. There are 23418 target genes, which is roughly equal to all known genes, I think. So there should exist some genes not expressed in hematopoietic cells, and the corresponding chromatin should be closed in biology.\n2. Contrary to that consideration(1), the facts were different; genes whose corresponding chromatins were not detected as open state were substantially expressed at least in given atac-seq data.\n3. So we have to predict gene(mRNA) expression not only by the corresponding open chromatin but also by other open chromatins or other gene expressions.","metadata":{}},{"cell_type":"code","source":"del train_geneactivity, train_targets, diff_a, diff_b, inter, diff_b_mean, inter_mean\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-08-31T06:08:03.128744Z","iopub.execute_input":"2022-08-31T06:08:03.129210Z","iopub.status.idle":"2022-08-31T06:08:03.726294Z","shell.execute_reply.started":"2022-08-31T06:08:03.129161Z","shell.execute_reply":"2022-08-31T06:08:03.725176Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Preprocess","metadata":{}},{"cell_type":"code","source":"#%%time\n# Preprocessing\n\nclass PreprocessMultiome(BaseEstimator, TransformerMixin):\n#     columns_to_use = slice(10000, 14000)\n    \n#     @staticmethod\n#     def take_column_subset(X):\n#         return X[:,PreprocessMultiome.columns_to_use]\n    def __init__(self, n_components=4):\n        super().__init__()\n        self.n_components = n_components\n    \n    def transform(self, X):\n        print(X.shape)\n        X = X[:,~self.all_zero_columns]\n        print(X.shape)\n#         X = PreprocessMultiome.take_column_subset(X) # use only a part of the columns\n#         print(X.shape)\n        gc.collect()\n\n        X = self.pca.transform(X)\n        print(X.shape)\n        return X\n\n    def fit_transform(self, X):\n        print(X.shape)\n        self.all_zero_columns = (X == 0).all(axis=0)\n        X = X[:,~self.all_zero_columns]\n        print(X.shape)\n#         X = PreprocessMultiome.take_column_subset(X) # use only a part of the columns\n#         print(X.shape)\n        gc.collect()\n\n        self.pca = PCA(n_components=self.n_components, copy=False, random_state=1)\n        X = self.pca.fit_transform(X)\n        plt.plot(self.pca.explained_variance_ratio_.cumsum())\n        plt.title(\"Cumulative explained variance ratio\")\n        plt.gca().xaxis.set_major_locator(MaxNLocator(integer=True))\n        plt.xlabel('PCA component')\n        plt.ylabel('Cumulative explained variance ratio')\n        plt.show()\n        print(X.shape)\n        return X\n\npreprocessor = PreprocessMultiome(n_components=4)\n\nmulti_train_x = None\nstart, stop = 0, 6000\nmulti_train_x = preprocessor.fit_transform(pd.read_hdf('../input/open-problems-train-geneactivity/train_multi_geneactivity.h5', start=start, stop=stop).values)\n\nmulti_train_y = pd.read_hdf('../input/open-problems-multimodal/train_multi_targets.h5', start=start, stop=stop)\ny_columns = multi_train_y.columns\nmulti_train_y = multi_train_y.values\nprint(multi_train_y.shape)","metadata":{"execution":{"iopub.status.busy":"2022-08-31T06:10:11.286429Z","iopub.execute_input":"2022-08-31T06:10:11.286862Z","iopub.status.idle":"2022-08-31T06:10:24.298849Z","shell.execute_reply.started":"2022-08-31T06:10:11.286825Z","shell.execute_reply":"2022-08-31T06:10:24.297531Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"I found the cumulative explained variance ratio at PC4 was higher in the geneactivity case than in the open-chromatin case (0.006 at PC4).","metadata":{}},{"cell_type":"markdown","source":"# Linear Regression","metadata":{}},{"cell_type":"code","source":"# Metrics\n# https://www.kaggle.com/code/ambrosm/msci-multiome-quickstart\ndef correlation_score(y_true, y_pred):\n    \"\"\"Scores the predictions according to the competition rules. \n    \n    It is assumed that the predictions are not constant.\n    \n    Returns the average of each sample's Pearson correlation coefficient\"\"\"\n    if type(y_true) == pd.DataFrame: y_true = y_true.values\n    if type(y_pred) == pd.DataFrame: y_pred = y_pred.values\n    if y_true.shape != y_pred.shape: raise ValueError(\"Shapes are different.\")\n    corrsum = 0\n    for i in range(len(y_true)):\n        corrsum += np.corrcoef(y_true[i], y_pred[i])[1, 0]\n    return corrsum / len(y_true)","metadata":{"execution":{"iopub.status.busy":"2022-08-31T06:14:19.653647Z","iopub.execute_input":"2022-08-31T06:14:19.654138Z","iopub.status.idle":"2022-08-31T06:14:19.662459Z","shell.execute_reply.started":"2022-08-31T06:14:19.654043Z","shell.execute_reply":"2022-08-31T06:14:19.661153Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# Cross-validation\n\nkf = KFold(n_splits=5, shuffle=True, random_state=1)\nscore_list = []\nfor fold, (idx_tr, idx_va) in enumerate(kf.split(multi_train_x)):\n    model = None\n    gc.collect()\n    X_tr = multi_train_x[idx_tr] # creates a copy, https://numpy.org/doc/stable/user/basics.copies.html\n    y_tr = multi_train_y[idx_tr]\n    del idx_tr\n\n    model = Ridge(copy_X=False)\n    model.fit(X_tr, y_tr)\n    del X_tr, y_tr\n    gc.collect()\n\n    # We validate the model\n    X_va = multi_train_x[idx_va]\n    y_va = multi_train_y[idx_va]\n    del idx_va\n    y_va_pred = model.predict(X_va)\n    mse = mean_squared_error(y_va, y_va_pred)\n    corrscore = correlation_score(y_va, y_va_pred)\n    del X_va, y_va\n\n    print(f\"Fold {fold}: mse = {mse:.5f}, corr =  {corrscore:.3f}\")\n    score_list.append((mse, corrscore))\n\n# Show overall score\nresult_df = pd.DataFrame(score_list, columns=['mse', 'corrscore'])\nprint(f\"{Fore.GREEN}{Style.BRIGHT}{multi_train_x.shape} Average  mse = {result_df.mse.mean():.5f}; corr = {result_df.corrscore.mean():.3f}{Style.RESET_ALL}\")","metadata":{"execution":{"iopub.status.busy":"2022-08-31T06:14:20.366914Z","iopub.execute_input":"2022-08-31T06:14:20.367887Z","iopub.status.idle":"2022-08-31T06:14:34.661324Z","shell.execute_reply.started":"2022-08-31T06:14:20.367841Z","shell.execute_reply":"2022-08-31T06:14:34.660352Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"got slightly better result than the original","metadata":{}},{"cell_type":"markdown","source":"# Retraining","metadata":{}},{"cell_type":"code","source":"# We retrain the model and then delete the training data, which is no longer needed\nmodel, score_list, result_df = None, None, None # free the RAM occupied by the old model\ngc.collect()\nmodel = Ridge(copy_X=False) # we overwrite the training data\nmodel.fit(multi_train_x, multi_train_y)\ndel multi_train_x, multi_train_y # free the RAM\n_ = gc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-08-31T06:15:40.975086Z","iopub.execute_input":"2022-08-31T06:15:40.975507Z","iopub.status.idle":"2022-08-31T06:15:42.604469Z","shell.execute_reply.started":"2022-08-31T06:15:40.975476Z","shell.execute_reply":"2022-08-31T06:15:42.603293Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Inference","metadata":{}},{"cell_type":"code","source":"%%time\n# Read the table of rows and columns required for submission\neval_ids = pd.read_csv('../input/open-problems-multimodal/evaluation_ids.csv', index_col='row_id')\n\n# Convert the string columns to more efficient categorical types\n#eval_ids.cell_id = eval_ids.cell_id.apply(lambda s: int(s, base=16))\neval_ids.cell_id = eval_ids.cell_id.astype(pd.CategoricalDtype())\neval_ids.gene_id = eval_ids.gene_id.astype(pd.CategoricalDtype())\ndisplay(eval_ids)\n\n# Create the set of needed cell_ids\ncell_id_set = set(eval_ids.cell_id)\n\n# Convert the string gene_ids to a more efficient categorical dtype\ny_columns = pd.CategoricalIndex(y_columns, dtype=eval_ids.gene_id.dtype, name='gene_id')","metadata":{"execution":{"iopub.status.busy":"2022-08-31T06:16:10.514144Z","iopub.execute_input":"2022-08-31T06:16:10.515309Z","iopub.status.idle":"2022-08-31T06:18:20.944955Z","shell.execute_reply.started":"2022-08-31T06:16:10.515254Z","shell.execute_reply":"2022-08-31T06:18:20.943852Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Prepare an empty series which will be filled with predictions\nsubmission = pd.Series(name='target',\n                       index=pd.MultiIndex.from_frame(eval_ids), \n                       dtype=np.float32)\nsubmission","metadata":{"execution":{"iopub.status.busy":"2022-08-31T06:18:20.946857Z","iopub.execute_input":"2022-08-31T06:18:20.947204Z","iopub.status.idle":"2022-08-31T06:18:21.122183Z","shell.execute_reply.started":"2022-08-31T06:18:20.947173Z","shell.execute_reply":"2022-08-31T06:18:21.121238Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# Process the test data in chunks of 5000 rows\n\nstart = 0\nchunksize = 5000\ntotal_rows = 0\nwhile True:\n    multi_test_x = None # Free the memory if necessary\n    gc.collect()\n    # Read the 5000 rows and select the 30 % subset which is needed for the submission\n    multi_test_x = pd.read_hdf('../input/open-problems-test-geneactivity/test_multi_geneactivity.h5', start=start, stop=start+chunksize)\n    rows_read = len(multi_test_x)\n    needed_row_mask = multi_test_x.index.isin(cell_id_set)\n    multi_test_x = multi_test_x.loc[needed_row_mask]\n    \n    # Keep the index (the cell_ids) for later\n    multi_test_index = multi_test_x.index\n    \n    # Predict\n    multi_test_x = multi_test_x.values\n    multi_test_x = preprocessor.transform(multi_test_x)\n    test_pred = model.predict(multi_test_x)\n    \n    # Convert the predictions to a dataframe so that they can be matched with eval_ids\n    test_pred = pd.DataFrame(test_pred,\n                             index=pd.CategoricalIndex(multi_test_index,\n                                                       dtype=eval_ids.cell_id.dtype,\n                                                       name='cell_id'),\n                             columns=y_columns)\n    gc.collect()\n    \n    # Fill the predictions into the submission series row by row\n    for i, (index, row) in enumerate(test_pred.iterrows()):\n        row = row.reindex(eval_ids.gene_id[eval_ids.cell_id == index])\n        submission.loc[index] = row.values\n    print('na:', submission.isna().sum())\n\n    #test_pred_list.append(test_pred)\n    total_rows += len(multi_test_x)\n    print(total_rows)\n    if rows_read < chunksize: break # this was the last chunk\n    start += chunksize\n    \ndel multi_test_x, multi_test_index, needed_row_mask","metadata":{"execution":{"iopub.status.busy":"2022-08-31T06:18:21.125010Z","iopub.execute_input":"2022-08-31T06:18:21.125346Z","iopub.status.idle":"2022-08-31T06:19:45.958650Z","shell.execute_reply.started":"2022-08-31T06:18:21.125315Z","shell.execute_reply":"2022-08-31T06:19:45.957544Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Submission\nIn this notebook, I did not train CITE-seq data. I temporarily use the public predictions from [vuonglam's great notebook](https://www.kaggle.com/code/vuonglam/lgbm-baseline-optuna-drop-constant-cite-task).","metadata":{}},{"cell_type":"code","source":"submission.reset_index(drop=True, inplace=True)\nsubmission.index.name = 'row_id'\n\ncite_submission = pd.read_csv(\"../input/lgbm-baseline-optuna-drop-constant-cite-task/submission.csv\")\ncite_submission = cite_submission.set_index(\"row_id\")\ncite_submission = cite_submission[\"target\"]\nsubmission[submission.isnull()] = cite_submission[submission.isnull()]\nsubmission.to_csv(\"submission.csv\")\nsubmission","metadata":{},"execution_count":null,"outputs":[]}]}