{"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":"code","source":"! pip install scanpy","metadata":{"execution":{"iopub.status.busy":"2022-10-30T18:35:15.943702Z","iopub.execute_input":"2022-10-30T18:35:15.944988Z","iopub.status.idle":"2022-10-30T18:35:28.236862Z","shell.execute_reply.started":"2022-10-30T18:35:15.944926Z","shell.execute_reply":"2022-10-30T18:35:28.235294Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nimport time\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport gseapy as gp # enrichment analysis\nimport scanpy as sc # single cell\n\nfrom sklearn.linear_model import Lasso\nfrom sklearn.model_selection import train_test_split, GridSearchCV, RepeatedKFold\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.linear_model import LassoCV\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-10-30T18:35:56.578392Z","iopub.execute_input":"2022-10-30T18:35:56.579515Z","iopub.status.idle":"2022-10-30T18:35:58.459203Z","shell.execute_reply.started":"2022-10-30T18:35:56.579347Z","shell.execute_reply":"2022-10-30T18:35:58.457628Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DATA_DIR = \"/kaggle/input/open-problems-multimodal/\"\nCITE_SEQ_GSE = \"/kaggle/input/citeseqgse148127\"\n\nCITE_SEQ_GSE148127 = os.path.join(CITE_SEQ_GSE,\"GSE148127_TotalA_Ref.csv\")\nMETA = os.path.join(CITE_SEQ_GSE,\"GSE148127_metadata.csv\")\nCOUNT = os.path.join(CITE_SEQ_GSE,\"GSE148127_ADT.counts.csv\")\nRNA_DATA = os.path.join(CITE_SEQ_GSE,\"GSE148127_SCT.normalized.RNA.counts.csv\")","metadata":{"execution":{"iopub.status.busy":"2022-10-30T18:35:58.460956Z","iopub.execute_input":"2022-10-30T18:35:58.461850Z","iopub.status.idle":"2022-10-30T18:35:58.469048Z","shell.execute_reply.started":"2022-10-30T18:35:58.461815Z","shell.execute_reply":"2022-10-30T18:35:58.467315Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Genes extraction for further enrichment analysis","metadata":{}},{"cell_type":"code","source":"total = pd.read_csv(CITE_SEQ_GSE148127, index_col=0)\nmeta = pd.read_csv(META, index_col=0)\n\ndf = pd.read_csv(COUNT, index_col=0)\ndf_cite = df.T\n\nRNA = pd.read_csv(RNA_DATA, index_col=0)\ndf_RNA = RNA.T","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Building one joined DataFrame\ndf_new = df_cite.join(df_RNA)\n\n# Making list of features\nfeatures_list = df_cite.columns.to_list()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Lasso mastering","metadata":{}},{"cell_type":"code","source":"# Just an example feature df\nexample_feature = features_list[0]\ntrain_df = df_new.drop(columns=features_list[1:])\n\n# Preparing data for Lasso\ny = train_df[example_feature]\nX = train_df.drop(columns=[example_feature])\n\nX_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, random_state=42)\n\n# Defining model\nmodel = Lasso(alpha=0.01)\n\n# Running train\nmodel.fit(X=X_train, y=y_train)\n\n# Squared tests\nprint('R squared training set', round(model.score(X_train, y_train)*100, 2))\nprint('R squared test set', round(model.score(X_test, y_test)*100, 2))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Running prediction\npred = model.predict(X)\nprint(\n    'Min:', np.min(pred), '\\nMax:', np.max(pred), '\\nMedian:', np.median(pred), \n      '\\n90%:', np.quantile(pred, 0.9)\n)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Finding target genes\ngenes = list(zip(pred, X))\nlist_genes = []\nfor gene in genes:\n    if gene[0] > 3:\n        list_genes.append(gene[1])\nprint(len(list_genes), '\\n', list_genes)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## LassoCV","metadata":{}},{"cell_type":"code","source":"# Lasso with 5 fold cross-validation\nmodel_cross = LassoCV(cv=5, random_state=0, max_iter=10000)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Fit model\nmodel_cross.fit(X_train, y_train)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model_cross.alpha_","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"alpha = model_cross.alpha_","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Model\nmodel = Lasso()\n# Evaluation method\ncv = RepeatedKFold(n_splits=10, n_repeats=3, random_state=1)\n\n# Grid\ngrid = dict()\ngrid['alpha'] = np.arange(0, 1, 0.01)\n\nsearch = GridSearchCV(model, grid, scoring='neg_mean_absolute_error', cv=cv, n_jobs=-1)\nresults = search.fit(X, y)\n\n# Summarize\nprint('MAE: %.3f' % results.best_score_)\nprint('Config: %s' % results.best_params_)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Run lasso","metadata":{}},{"cell_type":"code","source":"def lasso_run(feature, idx):\n    train_df = df_new.drop(columns=features_list[idx+1:])\n    \n    y = train_df[feature]\n    X = train_df.drop(columns=[feature])\n    \n    #X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, random_state=42)\n    model = Lasso(alpha=0.04838482756585978)  \n    model.fit(X=X, y=y)\n    # Squared tests\n    print(f'Lasso for {feature}:')\n    print('R squared training set', round(model.score(X, y)*100, 2))\n    print('R squared test set', round(model.score(X, y)*100, 2))\n    \n    pred = model.predict(X)\n    print(\n        'Min:', np.min(pred), '\\nMax:', np.max(pred), '\\nMedian:', np.median(pred), \n        '\\n90%:', np.quantile(pred, 0.9)\n    )\n    \n    genes = list(zip(pred, X))\n    list_genes = []\n    k = np.quantile(pred, 0.99)\n    for gene in genes:\n        if gene[0] > k:\n            list_genes.append(tuple([gene[1], gene[0]]))\n    \n    return np.min(pred), np.max(pred), np.median(pred), np.quantile(pred, 0.9), list_genes","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfinal = pd.DataFrame(columns=['Feature', 'Min.score', 'Max.score', 'Median.score', '0.9.score', 'Genes'])\n\nfor idx, feature in enumerate(features_list):\n    minimum, maximum, median, quantile, list_genes = lasso_run(feature, idx)\n    feature_df = pd.DataFrame([[feature, minimum, maximum, median, quantile, list_genes]],\n                              columns=['Feature', 'Min.score', 'Max.score', 'Median.score', '0.9.score', 'Genes'])\n    final = pd.concat([final, feature_df])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"final.to_csv('GSE148127_genes.csv')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Enrichment analysis","metadata":{}},{"cell_type":"code","source":"GSE_FEATURES = \"/kaggle/input/gse-features\"\nLASSO_OUTPUT = os.path.join(GSE_FEATURES, \"gse_features.csv\")\n\ngse_features = pd.read_csv(LASSO_OUTPUT, index_col=0)","metadata":{"execution":{"iopub.status.busy":"2022-10-30T18:36:07.133529Z","iopub.execute_input":"2022-10-30T18:36:07.134911Z","iopub.status.idle":"2022-10-30T18:36:07.168014Z","shell.execute_reply.started":"2022-10-30T18:36:07.134851Z","shell.execute_reply":"2022-10-30T18:36:07.167253Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gp.get_library_name(organism='Human')","metadata":{"execution":{"iopub.status.busy":"2022-10-30T19:53:56.809118Z","iopub.execute_input":"2022-10-30T19:53:56.809542Z","iopub.status.idle":"2022-10-30T19:53:57.067746Z","shell.execute_reply.started":"2022-10-30T19:53:56.809506Z","shell.execute_reply":"2022-10-30T19:53:57.066485Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_db_info(df: pd.DataFrame, database: str) -> pd.DataFrame:\n    \"\"\"Enrichment Analysis.\"\"\"\n    \n    for col in list(df.columns):\n        cd_genes = df[col].to_list()\n        \n        enrichment_results = gp.enrichr(\n            gene_list=cd_genes,\n            gene_sets=[database],\n            organism=\"Human\", \n            cutoff=0.5\n        )\n        \n        if not os.path.exists(database):\n            os.mkdir(database)\n            print(f\"The new directory {database} is created!\")\n        \n        enrichment_results_df = enrichment_results.results\n        enrichment_results_df = enrichment_results_df.loc[\n            enrichment_results_df[\"Adjusted P-value\"] < 0.05\n        ].reset_index(drop=True)\n        \n        if not enrichment_results_df.empty:\n            file_name = col + \".csv\"\n            enrichment_results_df.to_csv(\n                os.path.join(f\"/kaggle/working/{database}\", file_name), index=False\n            )","metadata":{"execution":{"iopub.status.busy":"2022-10-30T19:16:57.334990Z","iopub.execute_input":"2022-10-30T19:16:57.335381Z","iopub.status.idle":"2022-10-30T19:16:57.343570Z","shell.execute_reply.started":"2022-10-30T19:16:57.335349Z","shell.execute_reply":"2022-10-30T19:16:57.342581Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# get_db_info(gse_features, \"GO_Molecular_Function_2021\")\nget_db_info(gse_features, \"GO_Biological_Process_2021\")\nget_db_info(gse_features, \"KEGG_2021_Human\")","metadata":{"execution":{"iopub.status.busy":"2022-10-30T19:40:08.173105Z","iopub.execute_input":"2022-10-30T19:40:08.173475Z","iopub.status.idle":"2022-10-30T19:42:06.879296Z","shell.execute_reply.started":"2022-10-30T19:40:08.173445Z","shell.execute_reply":"2022-10-30T19:42:06.878387Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}