{"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"},{"sourceId":144604808,"sourceType":"kernelVersion"}],"dockerImageVersionId":30626,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# What is about ?\n\nApply batch correction algorithm to single cell RNA seq part of the data \n\nPS\n\nBasic analysis can be found at: \nhttps://www.kaggle.com/code/alexandervc/op2-rna-seq-data-scanpy-adata-cell-cycle","metadata":{}},{"cell_type":"markdown","source":"# Preliminaries ","metadata":{}},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport time\nt0start = time.time() \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-12-27T21:13:00.441800Z","iopub.execute_input":"2023-12-27T21:13:00.442307Z","iopub.status.idle":"2023-12-27T21:13:00.456306Z","shell.execute_reply.started":"2023-12-27T21:13:00.442268Z","shell.execute_reply":"2023-12-27T21:13:00.455133Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install scanpy\n!pip install symphonypy\n!pip install bbknn\n!pip install leidenalg\n\nimport numpy as np\nimport pandas as pd\nimport scanpy as sc\nimport matplotlib.pyplot as plt\nimport gc\nimport symphonypy as sp\n\nsc.settings.verbosity = 3             # verbosity: errors (0), warnings (1), info (2), hints (3)\nsc.logging.print_header()\nsc.settings.set_figure_params(dpi=80, facecolor='white')\n","metadata":{"execution":{"iopub.status.busy":"2023-12-27T21:09:04.246884Z","iopub.execute_input":"2023-12-27T21:09:04.248587Z","iopub.status.idle":"2023-12-27T21:10:54.822119Z","shell.execute_reply.started":"2023-12-27T21:09:04.248530Z","shell.execute_reply":"2023-12-27T21:10:54.820772Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Get adata","metadata":{}},{"cell_type":"code","source":"%%time\nfrom scipy.sparse import load_npz\n\nfn = '/kaggle/input/op2-rna-seq-data-to-sparse-matrix/rna_seq_counts_sparse_matrix.npz'\n# Load the sparse matrix from the saved npz file\nX = load_npz(fn)\nprint(type(X))\n# X = X.tocsr()\nprint(X.shape)\nX[:3,:2].toarray()\n\n# %%time\n# Create an AnnData object with the sparse matrix\nadata = sc.AnnData(X=X)\n\n# %%time\nfn = '/kaggle/input/op2-rna-seq-data-to-sparse-matrix/rna_seq_obs_id.csv'\ndf_obs = pd.read_csv(fn,index_col = 0)\nprint(df_obs.shape)\ndf_obs.columns = ['index cell']\n# df_obs = df_obs.set_index('cell_id')\n# df_obs['index'] = range(len(df_obs))\n#display(df_obs)\n\nfn = '/kaggle/input/open-problems-single-cell-perturbations/adata_obs_meta.csv'\ndf_obs2 = pd.read_csv(fn, index_col = 0)\nprint(df_obs2['cell_type'].unique() )\n#df_obs2\ndf_obs = df_obs.join(df_obs2, how = 'left')\ndisplay(df_obs)\n\n\n# %%time\nfn = '/kaggle/input/op2-rna-seq-data-to-sparse-matrix/rna_seq_gene_symbols.csv'\n\ndf_var = pd.read_csv(fn,index_col = 0)\nprint(df_var.shape)\ndf_var.columns = ['index gene']\n# df_var = df_var.set_index('cell_id')\n# df_var['index gene'] = range(len(df_var))\ndf_var\n\n# %%time\nadata.obs = df_obs\nadata.var = df_var\nadata\n\n\n","metadata":{"execution":{"iopub.status.busy":"2023-12-27T21:13:09.379935Z","iopub.execute_input":"2023-12-27T21:13:09.380452Z","iopub.status.idle":"2023-12-27T21:13:34.977863Z","shell.execute_reply.started":"2023-12-27T21:13:09.380414Z","shell.execute_reply":"2023-12-27T21:13:34.976495Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Preprocess - filter, etc","metadata":{}},{"cell_type":"code","source":"%%time\nprint(adata.shape)\nsc.pp.filter_cells(adata, min_genes=200)\nsc.pp.filter_genes(adata, min_cells=3)\nprint(adata.shape)\n\n# %%time\nadata.var['mt'] = adata.var_names.str.startswith('MT-')  # annotate the group of mitochondrial genes as 'mt'\nsc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)\n\n# %%time\nsc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],\n             jitter=0.4, multi_panel=True)\n\n# %%time\nsc.pl.scatter(adata, x='total_counts', y='pct_counts_mt')\nsc.pl.scatter(adata, x='total_counts', y='n_genes_by_counts')\n\n\n# %%time\n# Changed to more relaxed , defuault values too restrictive\nprint(adata.shape)\nadata = adata[adata.obs.n_genes_by_counts < 10_000, :]\nadata = adata[adata.obs.pct_counts_mt < 20, :]\nprint(adata.shape)\n","metadata":{"execution":{"iopub.status.busy":"2023-12-27T21:16:47.027794Z","iopub.execute_input":"2023-12-27T21:16:47.028313Z","iopub.status.idle":"2023-12-27T21:17:28.790947Z","shell.execute_reply.started":"2023-12-27T21:16:47.028275Z","shell.execute_reply":"2023-12-27T21:17:28.789583Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Norm & Log - preprocessing","metadata":{}},{"cell_type":"code","source":"sc.pp.normalize_total(adata, target_sum=1e4)\nsc.pp.log1p(adata)","metadata":{"execution":{"iopub.status.busy":"2023-12-27T21:17:28.793741Z","iopub.execute_input":"2023-12-27T21:17:28.794734Z","iopub.status.idle":"2023-12-27T21:17:43.407642Z","shell.execute_reply.started":"2023-12-27T21:17:28.794690Z","shell.execute_reply":"2023-12-27T21:17:43.406507Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# PCA ","metadata":{}},{"cell_type":"code","source":"%%time\nsc.tl.pca(adata, svd_solver='arpack')\nsc.pl.pca_variance_ratio(adata, log=True)\n# %%time\nsc.pl.pca(adata, color='CST3')\n","metadata":{"execution":{"iopub.status.busy":"2023-12-27T21:29:26.926467Z","iopub.execute_input":"2023-12-27T21:29:26.927053Z","iopub.status.idle":"2023-12-27T21:35:05.137562Z","shell.execute_reply.started":"2023-12-27T21:29:26.927019Z","shell.execute_reply":"2023-12-27T21:35:05.136315Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n!pip install leidenalg\nimport leidenalg","metadata":{"execution":{"iopub.status.busy":"2023-12-27T21:36:05.372435Z","iopub.execute_input":"2023-12-27T21:36:05.373542Z","iopub.status.idle":"2023-12-27T21:36:20.323201Z","shell.execute_reply.started":"2023-12-27T21:36:05.373497Z","shell.execute_reply":"2023-12-27T21:36:20.321348Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Harmony ","metadata":{}},{"cell_type":"code","source":"def cluster_harmony_batch_data(adata, resolution=1, batch_key = 'donor_id', k = 2 ):\n  adata = adata.copy()\n  sc.pp.highly_variable_genes(adata,min_mean=0.0125, max_mean=3, min_disp=0.5, batch_key= batch_key)\n  #Select genes that are variable in at least k samples \n  var_select = adata.var.highly_variable_nbatches >= k\n  var_genes = var_select.index[var_select]\n  adata.var['var_genes'] = var_select\n  adata = adata[:, adata.var.var_genes]\n  sc.pp.scale(adata, max_value=10)\n  sc.pp.pca(adata, n_comps=30)\n  sp.pp.harmony_integrate(adata, key=batch_key)\n  sc.pp.neighbors(adata, n_neighbors=50, use_rep=\"X_pca_harmony\")\n  sc.tl.leiden(adata)\n  sc.tl.umap(adata)\n  return adata\n","metadata":{"execution":{"iopub.status.busy":"2023-12-27T21:38:56.668389Z","iopub.execute_input":"2023-12-27T21:38:56.668879Z","iopub.status.idle":"2023-12-27T21:38:56.678212Z","shell.execute_reply.started":"2023-12-27T21:38:56.668841Z","shell.execute_reply":"2023-12-27T21:38:56.677125Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nadata2 = cluster_harmony_batch_data(adata, resolution=1, batch_key = 'donor_id', k = 2 )\nadata2","metadata":{"execution":{"iopub.status.busy":"2023-12-27T21:38:57.971102Z","iopub.execute_input":"2023-12-27T21:38:57.971559Z","iopub.status.idle":"2023-12-27T22:40:11.777636Z","shell.execute_reply.started":"2023-12-27T21:38:57.971521Z","shell.execute_reply":"2023-12-27T22:40:11.774515Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nsc.pl.umap(adata2, color=\"donor_id\", title=\"UMAP ...\", frameon=False)\n","metadata":{"execution":{"iopub.status.busy":"2023-12-27T22:40:11.779628Z","iopub.execute_input":"2023-12-27T22:40:11.780008Z","iopub.status.idle":"2023-12-27T22:40:13.819586Z","shell.execute_reply.started":"2023-12-27T22:40:11.779976Z","shell.execute_reply":"2023-12-27T22:40:13.818567Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nsc.pl.umap(adata2, color=\"cell_type\")# , title=\"UMAP ...\", frameon=False)","metadata":{"execution":{"iopub.status.busy":"2023-12-27T22:44:52.829700Z","iopub.execute_input":"2023-12-27T22:44:52.830189Z","iopub.status.idle":"2023-12-27T22:44:54.997959Z","shell.execute_reply.started":"2023-12-27T22:44:52.830151Z","shell.execute_reply":"2023-12-27T22:44:54.996680Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"adata_save = adata.copy()","metadata":{"execution":{"iopub.status.busy":"2023-12-27T22:46:26.485844Z","iopub.execute_input":"2023-12-27T22:46:26.486420Z","iopub.status.idle":"2023-12-27T22:46:29.847141Z","shell.execute_reply.started":"2023-12-27T22:46:26.486348Z","shell.execute_reply":"2023-12-27T22:46:29.845810Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"adata = adata2","metadata":{"execution":{"iopub.status.busy":"2023-12-27T22:46:33.533237Z","iopub.execute_input":"2023-12-27T22:46:33.533686Z","iopub.status.idle":"2023-12-27T22:46:33.539232Z","shell.execute_reply.started":"2023-12-27T22:46:33.533649Z","shell.execute_reply":"2023-12-27T22:46:33.537986Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\nprint( adata.obs['cell_type'].value_counts() )\n\nadata.obs['cell_type cat']= adata.obs['cell_type'].astype('category')\nadata.obs['cell_type T regulatory cells']= (adata.obs['cell_type'] == 'T regulatory cells').astype(float) # .astype('category')\nprint(adata.obs['cell_type T regulatory cells'].sum() )\nadata.obs['cell_type T cells CD4+']= (adata.obs['cell_type'] == 'T cells CD4+').astype(float) # .astype('category')\nprint(adata.obs['cell_type T cells CD4+'].sum() )\nadata.obs['cell_type T cells CD8+']= (adata.obs['cell_type'] == 'T cells CD8+').astype(float) # .astype('category')\nprint(adata.obs['cell_type T cells CD8+'].sum() )\nadata.obs['cell_type NK cells']= (adata.obs['cell_type'] == 'NK cells').astype(float) # .astype('category')\nprint(adata.obs['cell_type NK cells'].sum() )\nadata.obs['cell_type B cells']= (adata.obs['cell_type'] == 'B cells').astype(float) # .astype('category')\nprint(adata.obs['cell_type B cells'].sum() )\nadata.obs['cell_type Myeloid cells']= (adata.obs['cell_type'] == 'Myeloid cells').astype(float) # .astype('category')\nprint(adata.obs['cell_type Myeloid cells'].sum() )\n\nimport warnings\n\nwarnings.simplefilter(\"ignore\")\n\nsc.pl.umap(adata, color=[ 'leiden', 'donor_id', 'cell_type'  ]) # n_genes', 'n_genes_by_counts', 'total_counts'])\nsc.pl.umap(adata, color=['cell_type T regulatory cells','cell_type T cells CD4+','cell_type T cells CD8+']) # n_genes', 'n_genes_by_counts', 'total_counts'])\nsc.pl.umap(adata, color=['cell_type NK cells','cell_type B cells','cell_type Myeloid cells']) # n_genes', 'n_genes_by_counts', 'total_counts'])\n\n# sc.pl.umap(adata, color=['cell_type cat','total_counts_mt','sm_name']) # n_genes', 'n_genes_by_counts', 'total_counts'])\nsc.pl.umap(adata, color=['n_genes', 'n_genes_by_counts', 'total_counts'])\n\nsc.pl.umap(adata, color=['library_id', 'plate_name', 'well', 'row', 'col',]) # n_genes', 'n_genes_by_counts', 'total_counts'])\n\n","metadata":{"execution":{"iopub.status.busy":"2023-12-27T22:46:35.619098Z","iopub.execute_input":"2023-12-27T22:46:35.620272Z","iopub.status.idle":"2023-12-27T22:47:07.034672Z","shell.execute_reply.started":"2023-12-27T22:46:35.620230Z","shell.execute_reply":"2023-12-27T22:47:07.033279Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# AmbrosM found drugs (misclassified cell types) and other drugs - UMAP\n\n    'CGM-097', 'LY2090314','Ganetespib (STA-9090)','IN1451', 'Oprozomib (ONX 0912)','MLN 2238','CEP-18770 (Delanzomib)'\n    \nindicated by AmbrosM as potentially leading to misclassified cell types - see\nhttps://www.kaggle.com/competitions/open-problems-single-cell-perturbations/discussion/458661\n","metadata":{}},{"cell_type":"code","source":"%%time\nprint(adata.obs['sm_name'].value_counts().head(10))\nprint(adata.obs['sm_name'].value_counts().tail(10))\n\nprint(adata.obs['sm_name'].value_counts().head(10).index.to_list() )\n\nprint(adata.obs['sm_name'].unique().to_list() )\n\nlist_selected_drugs = ['Dimethyl Sulfoxide', 'Dabrafenib', 'Belinostat', 'CGM-097', 'LY2090314','Ganetespib (STA-9090)','IN1451', 'Oprozomib (ONX 0912)','MLN 2238','CEP-18770 (Delanzomib)',]\nfor drug in  list_selected_drugs:\n    adata.obs[drug] = (adata.obs['sm_name'] == drug).astype('category') \n    print(drug, (adata.obs['sm_name'] == drug).sum() )","metadata":{"execution":{"iopub.status.busy":"2023-12-27T22:47:46.402227Z","iopub.execute_input":"2023-12-27T22:47:46.402651Z","iopub.status.idle":"2023-12-27T22:47:46.465532Z","shell.execute_reply.started":"2023-12-27T22:47:46.402620Z","shell.execute_reply":"2023-12-27T22:47:46.464204Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nsc.pl.umap(adata, color=['donor_id', 'cell_type'  ]) # n_genes', 'n_genes_by_counts', 'total_counts'])\n\nsc.pl.umap(adata, color=list_selected_drugs[:3]) # n_genes', 'n_genes_by_counts', 'total_counts'])\nsc.pl.umap(adata, color=list_selected_drugs[3:7]) # n_genes', 'n_genes_by_counts', 'total_counts'])\nsc.pl.umap(adata, color=list_selected_drugs[7:]) # n_genes', 'n_genes_by_counts', 'total_counts'])\n\nsc.pl.umap(adata, color=['donor_id', 'cell_type','cell_type T cells CD8+' ,'cell_type T cells CD4+'  ]) # n_genes', 'n_genes_by_counts', 'total_counts'])\n","metadata":{"execution":{"iopub.status.busy":"2023-12-27T22:48:04.016500Z","iopub.execute_input":"2023-12-27T22:48:04.016927Z","iopub.status.idle":"2023-12-27T22:48:35.567388Z","shell.execute_reply.started":"2023-12-27T22:48:04.016895Z","shell.execute_reply":"2023-12-27T22:48:35.565972Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nsc.pl.umap(adata, color=['donor_id', 'cell_type'  ]) # n_genes', 'n_genes_by_counts', 'total_counts'])\nsc.pl.umap(adata, color=['cell_type T regulatory cells','cell_type T cells CD4+','cell_type T cells CD8+']) # n_genes', 'n_genes_by_counts', 'total_counts'])\nsc.pl.umap(adata, color=['cell_type NK cells','cell_type B cells','cell_type Myeloid cells']) # n_genes', 'n_genes_by_counts', 'total_counts'])\n\n# sc.pl.umap(adata, color=['cell_type cat','total_counts_mt','sm_name']) # n_genes', 'n_genes_by_counts', 'total_counts'])\nsc.pl.umap(adata, color=['n_genes', 'n_genes_by_counts', 'total_counts'])\n\nsc.pl.umap(adata, color=['library_id', 'plate_name', 'well', 'row', 'col',]) # n_genes', 'n_genes_by_counts', 'total_counts'])\n\n","metadata":{"execution":{"iopub.status.busy":"2023-12-27T22:48:35.569985Z","iopub.execute_input":"2023-12-27T22:48:35.570404Z","iopub.status.idle":"2023-12-27T22:49:04.610613Z","shell.execute_reply.started":"2023-12-27T22:48:35.570372Z","shell.execute_reply":"2023-12-27T22:49:04.609356Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n# adata_sc = sc.concat([adata_sc1,adata_sc2], axis =0 )\n\n# # Preprocessing \n# def pp_data(adata):\n#   adata = adata.copy()\n#   sc.pp.filter_cells(adata, min_genes=500)  \n#   #sc.pp.filter_genes(adata, min_cells=3) # genes being expressed in at least 'min_cells' cells \n#   adata.var['mt'] = adata.var_names.str.startswith('MT-')  # mitochondrial genes \n#   sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)\n#   return adata\n\n# adata_sc = pp_data(adata_sc)\n\n# # Preprocessing \n# def filter_data(adata, resolution=1):\n#   adata = adata.copy()\n#   adata = adata[adata.obs.pct_counts_mt < 20, :] \n#   #adata = adata[adata.obs.n_genes_by_counts < 10000, :] \n#   sc.pp.normalize_total(adata, target_sum=1e4)\n#   sc.pp.log1p(adata)\n  \n#   return adata\n\n# adata_sc = filter_data(adata_sc)\n# adata_sc.raw = adata_sc\n\n# # No batch correction\n# def cluster_nobatch_data(adata, resolution=1):\n#   adata = adata.copy()\n\n# #   HVG regardless batches\n# #   sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)\n# #   adata = adata[:, adata.var.highly_variable]\n\n# #  Or find HVG in each sample and select genes that are variable in at least k samples\n#   sc.pp.highly_variable_genes(adata,min_mean=0.0125, max_mean=3, min_disp=0.5, batch_key='sample_name')\n#   k = 5\n#   var_select = adata.var.highly_variable_nbatches >= k\n#   var_genes = var_select.index[var_select]\n#   adata.var['var_genes'] = var_select\n#   adata = adata[:, adata.var.var_genes]\n\n#   sc.pp.scale(adata, max_value=10)\n#   sc.pp.pca(adata, n_comps=30)\n#   sc.pp.neighbors(adata)\n#   sc.tl.leiden(adata)\n#   sc.tl.umap(adata)\n#   return adata\n\n\n# adata_sc2 = adata_sc.copy()\n# adata_sc2 = cluster_nobatch_data(adata_sc2)\n# sc.pl.umap(adata_sc2, color=\"sample_name\", title=\"UMAP ...\", frameon=False)\n\n# # Harmony batch correction\n# def cluster_harmony_batch_data(adata, resolution=1):\n#   adata = adata.copy()\n#   sc.pp.highly_variable_genes(adata,min_mean=0.0125, max_mean=3, min_disp=0.5, batch_key='sample_name')\n#   #Select genes that are variable in at least k samples \n#   var_select = adata.var.highly_variable_nbatches >= k\n#   var_genes = var_select.index[var_select]\n#   adata.var['var_genes'] = var_select\n#   adata = adata[:, adata.var.var_genes]\n#   sc.pp.scale(adata, max_value=10)\n#   sc.pp.pca(adata, n_comps=30)\n#   sp.pp.harmony_integrate(adata, key='sample_name')\n#   sc.pp.neighbors(adata, n_neighbors=50, use_rep=\"X_pca_harmony\")\n#   sc.tl.leiden(adata)\n#   sc.tl.umap(adata)\n#   return adata\n\n# adata_sc3 = adata_sc.copy()\n# adata_sc3 = cluster_harmony_batch_data(adata_sc3)\n# sc.pl.umap(adata_sc3, color=\"sample_name\", title=\"UMAP ...\", frameon=False)\n\n# # BBKNN batch correction\n# def cluster_bbknn_batch_data(adata, resolution=1):\n#   adata = adata.copy()\n#   sc.pp.highly_variable_genes(adata,min_mean=0.0125, max_mean=3, min_disp=0.5, batch_key='samples')\n#   var_select = adata.var.highly_variable_nbatches > k\n#   var_genes = var_select.index[var_select]\n#   adata.var['var_genes'] = var_select\n#   adata = adata[:, adata.var.var_genes]\n#   sc.pp.scale(adata, max_value=10)\n#   sc.pp.pca(adata, n_comps=30)\n#   sc.external.pp.bbknn(adata, batch_key='sample_name')  # running bbknn 1.3.6\n#   sc.tl.leiden(adata)\n#   sc.tl.umap(adata)\n#   return adata\n# ​\n# ​\n# # Combat batch correction\n# def cluster_combat_batch_data(adata, resolution=1):\n#   adata = adata.copy()\n#   sc.pp.highly_variable_genes(adata,min_mean=0.0125, max_mean=3, min_disp=0.5, batch_key='samples')\n#   var_select = adata.var.highly_variable_nbatches > k\n#   var_genes = var_select.index[var_select]\n#   adata.var['var_genes'] = var_select\n#   adata = adata[:, adata.var.var_genes]\n#   sc.pp.scale(adata, max_value=10)\n#   sc.pp.combat(adata, key='sample_name')\n#   sc.pp.pca(adata, n_comps=30)\n#   sc.pp.neighbors(adata)\n#   sc.tl.leiden(adata)\n#   sc.tl.umap(adata)\n#   return adata\n\n\n\n\n\n\n\n\n","metadata":{"execution":{"iopub.status.busy":"2023-12-27T22:44:43.910303Z","iopub.execute_input":"2023-12-27T22:44:43.910976Z","iopub.status.idle":"2023-12-27T22:44:43.921501Z","shell.execute_reply.started":"2023-12-27T22:44:43.910941Z","shell.execute_reply":"2023-12-27T22:44:43.919924Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('%.1f seconds passed total '%(time.time()-t0start) )\nprint('%.1f minutes passed total '%( (time.time()-t0start)/60)  )\nprint('%.2f hours passed total '%( (time.time()-t0start)/3600)  )","metadata":{},"execution_count":null,"outputs":[]}]}