{"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":7298890,"sourceType":"datasetVersion","datasetId":4234070},{"sourceId":144604808,"sourceType":"kernelVersion"}],"dockerImageVersionId":30558,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# What is about ?\n\nLoad scRNA-seq data and first look. Use basic scanpy pipeline and brief look on cell (proliferation) cycle. \n\nWe load data prepared in the notebook: https://www.kaggle.com/alexandervc/op2-rna-seq-data-to-sparse-matrix\n    \nMake first look on them following scanpy tutorial: https://scanpy-tutorials.readthedocs.io/en/latest/pbmc3k.html\n        \n        \nAlso brief look on cell (prolifiration) cycle - almost no prolifirating cells is seen \n\n\nVersions:\n\n    10 save created adata: \"adata_OPSCP_rawcounts_sparse_240090_21255.h5ad\", future versions will just load it\n    7 major revision - cell types, drugs, analysis added , analysis of found by AmbrosM drugs leads to interesting observations \n    \n    3,4,5,6 cosmetic changes\n    2 use second version of prepared data - CSR matrix sparse format, and other small changes.\n        Relax filtering by MT and n_genes_by_counts\n    \n    1 data from the first version of the notebook above, we need to update since we changed the format in version 2 of the data prepration nb.  \n        default thresholds on MT and n_genes_by_counts are too restrictive - will change ","metadata":{}},{"cell_type":"markdown","source":"# Preliminaries / Scanpy","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-28T21:53:59.429677Z","iopub.execute_input":"2023-12-28T21:53:59.430201Z","iopub.status.idle":"2023-12-28T21:54:00.993463Z","shell.execute_reply.started":"2023-12-28T21:53:59.430157Z","shell.execute_reply":"2023-12-28T21:54:00.992201Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install scanpy","metadata":{"execution":{"iopub.status.busy":"2023-12-28T21:54:00.995820Z","iopub.execute_input":"2023-12-28T21:54:00.996830Z","iopub.status.idle":"2023-12-28T21:54:22.842234Z","shell.execute_reply.started":"2023-12-28T21:54:00.996790Z","shell.execute_reply":"2023-12-28T21:54:22.840312Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import scanpy as sc","metadata":{"execution":{"iopub.status.busy":"2023-12-28T21:54:22.845630Z","iopub.execute_input":"2023-12-28T21:54:22.846155Z","iopub.status.idle":"2023-12-28T21:54:29.881481Z","shell.execute_reply.started":"2023-12-28T21:54:22.846112Z","shell.execute_reply":"2023-12-28T21:54:29.879868Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load / Create adata\n\n\"adata\" is \"dataframe on steroids\" - the standard format used in single-cell RNA seq data processing in Python.\n\nThe idea is that both \"index\" and \"columns\" of the \"adata(steroid dataframe)\" are pandas dataframes themselves - so allow to store lots of information.\n\nWhile the \"main-meat\" is stored in adata.X - analogue of df.values, which can support \"X\" to be sparse matrix, pay attention it can only store numeric values in contrast to pandas dataframe (all non-numeric should be moved to \"index = adata.obs\").\n\nVocabluary:\n\n    df.values <-> adata.X (can be sparse or numpy but only supports numeric values)\n    df.index <-> adata.obs - it is a kind of \"extended index\" which is dataframe by itself\n    df.columns <-> adata.var - it is a kind of \"extened columns\" which is dataframe by itself \n    \n    \nSee docs: https://scanpy.readthedocs.io/en/stable/usage-principles.html    ","metadata":{}},{"cell_type":"code","source":"%%time\n\nif 1: # Load prepared adata\n    fn = '/kaggle/input/open-problems-single-cell-perturbations-data/adata_OPSCP_rawcounts_sparse_240090_21255.h5ad'\n    print(fn)\n    adata = sc.read_h5ad(fn)\n    print(adata)    \n\nif 0: # Here is code to prepare adata. We run it in versions 1-10, later saved it to Kaggle dataset and now use prepared file \n    from scipy.sparse import load_npz\n\n    fn = '/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\n    X = load_npz(fn)\n    print(type(X))\n    # X = X.tocsr()\n    print(X.shape)\n    X[:3,:2].toarray()\n\n    # %%time\n    # Create an AnnData object with the sparse matrix\n    adata = sc.AnnData(X=X)\n\n    # %%time\n    fn = '/kaggle/input/op2-rna-seq-data-to-sparse-matrix/rna_seq_obs_id.csv'\n    df_obs = pd.read_csv(fn,index_col = 0)\n    print(df_obs.shape)\n    df_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\n    fn = '/kaggle/input/open-problems-single-cell-perturbations/adata_obs_meta.csv'\n    df_obs2 = pd.read_csv(fn, index_col = 0)\n    print(df_obs2['cell_type'].unique() )\n    #df_obs2\n    df_obs = df_obs.join(df_obs2, how = 'left')\n    display(df_obs)\n\n\n    # %%time\n    m = df_obs['sm_name'] == 'Dimethyl Sulfoxide'\n    df_obs['control negative'] = (m).astype(float)\n    print( df_obs['control negative'].sum() )\n\n    m = df_obs['sm_name'] == 'Dabrafenib'\n    m = m | ( df_obs['sm_name'] == 'Belinostat' )\n    df_obs['control positive'] = (m).astype(float)\n    print( df_obs['control positive'].sum() )\n\n\n    # %%time\n    fn = '/kaggle/input/op2-rna-seq-data-to-sparse-matrix/rna_seq_gene_symbols.csv'\n    df_var = pd.read_csv(fn,index_col = 0)\n    print(df_var.shape)\n    df_var.columns = ['index gene']\n    # df_var = df_var.set_index('cell_id')\n    # df_var['index gene'] = range(len(df_var))\n    df_var\n\n    # %%time\n    adata.obs = df_obs\n    adata.var = df_var\n    adata\n\n    # %%time\n    if 0:\n        adata.write('adata_OPSCP_rawcounts_sparse_240090_21255.h5ad', compression='gzip')\n    # %%time\n    if 0:\n        loaded_adata = sc.read_h5ad('adata_OPSCP_rawcounts_sparse_240090_21255.h5ad')\n        print(loaded_adata)    \n        adata = loaded_adata\n","metadata":{"execution":{"iopub.status.busy":"2023-12-28T21:54:29.884124Z","iopub.execute_input":"2023-12-28T21:54:29.885104Z","iopub.status.idle":"2023-12-28T21:54:55.991199Z","shell.execute_reply.started":"2023-12-28T21:54:29.885060Z","shell.execute_reply":"2023-12-28T21:54:55.989922Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"adata.obs","metadata":{"execution":{"iopub.status.busy":"2023-12-28T21:54:55.994508Z","iopub.execute_input":"2023-12-28T21:54:55.994987Z","iopub.status.idle":"2023-12-28T21:54:56.056984Z","shell.execute_reply.started":"2023-12-28T21:54:55.994943Z","shell.execute_reply":"2023-12-28T21:54:56.055342Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"adata.var","metadata":{"execution":{"iopub.status.busy":"2023-12-28T21:54:56.061071Z","iopub.execute_input":"2023-12-28T21:54:56.062088Z","iopub.status.idle":"2023-12-28T21:54:56.077852Z","shell.execute_reply.started":"2023-12-28T21:54:56.062015Z","shell.execute_reply":"2023-12-28T21:54:56.076088Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"m = adata.obs['control'] == True\nprint(m.sum())\nprint( adata.obs[m]['sm_name'].value_counts() )\n\nprint( list(adata.obs['sm_name'].unique() ) )\nadata.obs['sm_name'].value_counts()\n ","metadata":{"execution":{"iopub.status.busy":"2023-12-28T21:54:56.079105Z","iopub.execute_input":"2023-12-28T21:54:56.079502Z","iopub.status.idle":"2023-12-28T21:54:56.127500Z","shell.execute_reply.started":"2023-12-28T21:54:56.079457Z","shell.execute_reply":"2023-12-28T21:54:56.126291Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Some examples to operate with adata","metadata":{}},{"cell_type":"code","source":"adata[:,'E2F1'].X.toarray().sum()","metadata":{"execution":{"iopub.status.busy":"2023-12-28T21:54:56.129170Z","iopub.execute_input":"2023-12-28T21:54:56.130671Z","iopub.status.idle":"2023-12-28T21:54:57.011967Z","shell.execute_reply.started":"2023-12-28T21:54:56.130619Z","shell.execute_reply":"2023-12-28T21:54:57.010285Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"adata[:,['E2F1','E2F2']].X.toarray().sum(axis = 0)","metadata":{"execution":{"iopub.status.busy":"2023-12-28T21:54:57.014203Z","iopub.execute_input":"2023-12-28T21:54:57.014670Z","iopub.status.idle":"2023-12-28T21:54:58.362514Z","shell.execute_reply.started":"2023-12-28T21:54:57.014636Z","shell.execute_reply":"2023-12-28T21:54:58.361111Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"adata[:100,:100].to_df()","metadata":{"execution":{"iopub.status.busy":"2023-12-28T21:54:58.364329Z","iopub.execute_input":"2023-12-28T21:54:58.364759Z","iopub.status.idle":"2023-12-28T21:54:58.401375Z","shell.execute_reply.started":"2023-12-28T21:54:58.364724Z","shell.execute_reply":"2023-12-28T21:54:58.399441Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Genes not exrpessed in certain cell types ","metadata":{}},{"cell_type":"code","source":"%%time\nlist_cell_types = adata.obs['cell_type'].unique()\nfor cell_type in list_cell_types:\n    m0 = adata.obs['cell_type'] == cell_type\n    print(cell_type,'count cells:',  m0.sum() )\n    v2 = (adata[m0,:].X!=0).sum(axis = 0 )# .\n    k = 0\n    m = np.asarray(v2).ravel() == k\n    print(m.sum(), 'genes expressed in  exactly ', k, 'cells for cell type', cell_type  )      \n    print('50 genes:',  list(adata.var[m].index)[:50] )\n    print()\n    ","metadata":{"execution":{"iopub.status.busy":"2023-12-28T17:36:41.195367Z","iopub.execute_input":"2023-12-28T17:36:41.196617Z","iopub.status.idle":"2023-12-28T17:36:46.163438Z","shell.execute_reply.started":"2023-12-28T17:36:41.196569Z","shell.execute_reply":"2023-12-28T17:36:46.162365Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Least and top expressed genes","metadata":{}},{"cell_type":"code","source":"%%time\nv2 = (adata.X!=0).sum(axis = 0 )\nv2 = np.asarray(v2).ravel()\nprint('Statistics: number of cells , where gene is non-zero - distribution over all genes')\nprint(pd.Series(v2).describe() )\nplt.hist(v2 , bins = 100)\nplt.title('number of cells , where gene is non-zero - distribution over all genes',fontsize = 17)\nplt.show()\nprint('Top frequent - numbers - how many cells gene is expressed in:')\nprint( pd.Series(v2).value_counts().head(20) )\nprint(np.sort(v2)[:30])\nprint(np.sort(v2)[-30:])\nprint()\n\n\nprint('Least expressed genes (by number of cells)')\nfor k in [1,2,3,4,5,10,100]:\n    m = np.asarray(v2).ravel() == k\n    print(m.sum(), 'genes expressed in  exactly ', k, 'cells ')      \n    print('50 genes:',  list(adata.var[m].index)[:50] )\n    print()\n\nprint(); print()    \nprint('Top expressed genes (by number of cells)')\nfor k in [ 240_000, 230_000]:    \n    m = np.asarray(v2).ravel() >= k\n    print(m.sum(), 'genes expressed in  >=  ', k, 'cells ')      \n    print('50 genes:',  list(adata.var[m].index)[:50] )\n    print()\n","metadata":{"execution":{"iopub.status.busy":"2023-12-28T21:58:59.917228Z","iopub.execute_input":"2023-12-28T21:58:59.918021Z","iopub.status.idle":"2023-12-28T21:59:04.837490Z","shell.execute_reply.started":"2023-12-28T21:58:59.917957Z","shell.execute_reply":"2023-12-28T21:59:04.834949Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Start analysis (standard pipeline)\n\nStandard single cell rna seq data analysis piplene steps, see scanpy tutorial: https://scanpy-tutorials.readthedocs.io/en/latest/pbmc3k.html\n","metadata":{}},{"cell_type":"code","source":"sc.pl.highest_expr_genes(adata, n_top=40, )","metadata":{"execution":{"iopub.status.busy":"2023-12-28T17:37:03.379127Z","iopub.execute_input":"2023-12-28T17:37:03.379510Z","iopub.status.idle":"2023-12-28T17:37:16.598598Z","shell.execute_reply.started":"2023-12-28T17:37:03.379480Z","shell.execute_reply":"2023-12-28T17:37:16.595551Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Filtering","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","metadata":{"execution":{"iopub.status.busy":"2023-12-28T21:59:38.015108Z","iopub.execute_input":"2023-12-28T21:59:38.015783Z","iopub.status.idle":"2023-12-28T21:59:57.197478Z","shell.execute_reply.started":"2023-12-28T21:59:38.015621Z","shell.execute_reply":"2023-12-28T21:59:57.195794Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%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)","metadata":{"execution":{"iopub.status.busy":"2023-12-28T21:59:57.200602Z","iopub.execute_input":"2023-12-28T21:59:57.201438Z","iopub.status.idle":"2023-12-28T22:00:19.621471Z","shell.execute_reply.started":"2023-12-28T21:59:57.201392Z","shell.execute_reply":"2023-12-28T22:00:19.618149Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nsc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],\n             jitter=0.4, multi_panel=True)","metadata":{"execution":{"iopub.status.busy":"2023-12-28T22:00:25.057721Z","iopub.execute_input":"2023-12-28T22:00:25.058281Z","iopub.status.idle":"2023-12-28T22:00:31.194446Z","shell.execute_reply.started":"2023-12-28T22:00:25.058242Z","shell.execute_reply":"2023-12-28T22:00:31.193302Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%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","metadata":{"execution":{"iopub.status.busy":"2023-12-28T17:38:30.501577Z","iopub.execute_input":"2023-12-28T17:38:30.502012Z","iopub.status.idle":"2023-12-28T17:38:32.593391Z","shell.execute_reply.started":"2023-12-28T17:38:30.501977Z","shell.execute_reply":"2023-12-28T17:38:32.592361Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%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-28T22:02:28.211016Z","iopub.execute_input":"2023-12-28T22:02:28.211624Z","iopub.status.idle":"2023-12-28T22:02:28.599080Z","shell.execute_reply.started":"2023-12-28T22:02:28.211573Z","shell.execute_reply":"2023-12-28T22:02:28.597961Z"},"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)","metadata":{"execution":{"iopub.status.busy":"2023-12-28T22:02:32.337472Z","iopub.execute_input":"2023-12-28T22:02:32.337981Z","iopub.status.idle":"2023-12-28T22:02:43.813339Z","shell.execute_reply.started":"2023-12-28T22:02:32.337943Z","shell.execute_reply":"2023-12-28T22:02:43.811915Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sc.pp.log1p(adata)","metadata":{"execution":{"iopub.status.busy":"2023-12-28T22:02:43.816267Z","iopub.execute_input":"2023-12-28T22:02:43.816816Z","iopub.status.idle":"2023-12-28T22:02:51.024616Z","shell.execute_reply.started":"2023-12-28T22:02:43.816770Z","shell.execute_reply":"2023-12-28T22:02:51.022749Z"},"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","metadata":{"execution":{"iopub.status.busy":"2023-12-28T17:39:08.580124Z","iopub.execute_input":"2023-12-28T17:39:08.581028Z","iopub.status.idle":"2023-12-28T17:44:39.560038Z","shell.execute_reply.started":"2023-12-28T17:39:08.580983Z","shell.execute_reply":"2023-12-28T17:44:39.558939Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nsc.pl.pca(adata, color=['CST3','MALAT1','B2M'] )\nsc.pl.pca(adata, color=['cell_type', 'control negative', 'control positive' ])\nsc.pl.pca(adata, color=['n_genes_by_counts', 'total_counts', 'pct_counts_mt'])\n\n","metadata":{"execution":{"iopub.status.busy":"2023-12-28T17:48:10.283713Z","iopub.execute_input":"2023-12-28T17:48:10.284211Z","iopub.status.idle":"2023-12-28T17:48:22.629751Z","shell.execute_reply.started":"2023-12-28T17:48:10.284176Z","shell.execute_reply":"2023-12-28T17:48:22.628733Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# UMAP","metadata":{}},{"cell_type":"code","source":"%%time\n!pip install leidenalg","metadata":{"execution":{"iopub.status.busy":"2023-12-28T22:02:51.026621Z","iopub.execute_input":"2023-12-28T22:02:51.027030Z","iopub.status.idle":"2023-12-28T22:03:08.582887Z","shell.execute_reply.started":"2023-12-28T22:02:51.026998Z","shell.execute_reply":"2023-12-28T22:03:08.581433Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import leidenalg","metadata":{"execution":{"iopub.status.busy":"2023-12-28T22:03:08.586496Z","iopub.execute_input":"2023-12-28T22:03:08.587468Z","iopub.status.idle":"2023-12-28T22:03:08.817534Z","shell.execute_reply.started":"2023-12-28T22:03:08.587409Z","shell.execute_reply":"2023-12-28T22:03:08.814981Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nsc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)\nsc.tl.leiden(adata)\nsc.tl.paga(adata)\nsc.pl.paga(adata, plot= False)  # remove `plot=False` if you want to see the coarse-grained graph\nsc.tl.umap(adata, init_pos='paga')","metadata":{"execution":{"iopub.status.busy":"2023-12-28T22:03:19.666981Z","iopub.execute_input":"2023-12-28T22:03:19.667463Z","iopub.status.idle":"2023-12-28T22:31:25.095337Z","shell.execute_reply.started":"2023-12-28T22:03:19.667427Z","shell.execute_reply":"2023-12-28T22:31:25.093663Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nsc.tl.umap(adata)\nsc.pl.umap(adata, color=['CST3', 'NKG7', 'PPBP'])","metadata":{"execution":{"iopub.status.busy":"2023-12-28T22:33:16.877870Z","iopub.execute_input":"2023-12-28T22:33:16.878455Z","iopub.status.idle":"2023-12-28T22:39:55.603185Z","shell.execute_reply.started":"2023-12-28T22:33:16.878415Z","shell.execute_reply":"2023-12-28T22:39:55.602104Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"adata","metadata":{"execution":{"iopub.status.busy":"2023-12-28T22:43:49.132608Z","iopub.execute_input":"2023-12-28T22:43:49.133096Z","iopub.status.idle":"2023-12-28T22:43:49.142591Z","shell.execute_reply.started":"2023-12-28T22:43:49.133062Z","shell.execute_reply":"2023-12-28T22:43:49.141129Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nif 0:\n    adata.write('adata_OPSCP_FilteredNormalized_sparse_'+str(adata.shape[0])+'_'+str(adata.shape[1])+'.h5ad', compression='gzip')\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nsc.pl.umap(adata, color=['CST3','MALAT1','B2M'] )\nsc.pl.umap(adata, color=['cell_type', 'control negative', 'control positive' ])\nsc.pl.umap(adata, color=['n_genes_by_counts', 'total_counts', 'pct_counts_mt'])\n","metadata":{"execution":{"iopub.status.busy":"2023-12-28T18:45:07.233202Z","iopub.execute_input":"2023-12-28T18:45:07.233681Z","iopub.status.idle":"2023-12-28T18:45:22.099433Z","shell.execute_reply.started":"2023-12-28T18:45:07.233645Z","shell.execute_reply":"2023-12-28T18:45:22.098297Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# UMAP Cell types, etc","metadata":{}},{"cell_type":"code","source":"print(adata.obs.columns)","metadata":{"execution":{"iopub.status.busy":"2023-12-27T18:59:18.855536Z","iopub.execute_input":"2023-12-27T18:59:18.856449Z","iopub.status.idle":"2023-12-27T18:59:18.863072Z","shell.execute_reply.started":"2023-12-27T18:59:18.856410Z","shell.execute_reply":"2023-12-27T18:59:18.861893Z"},"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-27T19:21:22.937699Z","iopub.execute_input":"2023-12-27T19:21:22.938265Z","iopub.status.idle":"2023-12-27T19:21:55.685141Z","shell.execute_reply.started":"2023-12-27T19:21:22.938225Z","shell.execute_reply":"2023-12-27T19:21:55.684253Z"},"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-27T19:12:03.256315Z","iopub.execute_input":"2023-12-27T19:12:03.256892Z","iopub.status.idle":"2023-12-27T19:12:03.317996Z","shell.execute_reply.started":"2023-12-27T19:12:03.256851Z","shell.execute_reply":"2023-12-27T19:12:03.316805Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import warnings\n\nwarnings.simplefilter(\"ignore\")","metadata":{"execution":{"iopub.status.busy":"2023-12-27T19:14:39.355687Z","iopub.execute_input":"2023-12-27T19:14:39.357111Z","iopub.status.idle":"2023-12-27T19:14:39.363020Z","shell.execute_reply.started":"2023-12-27T19:14:39.357057Z","shell.execute_reply":"2023-12-27T19:14:39.361994Z"},"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-27T19:26:28.385819Z","iopub.execute_input":"2023-12-27T19:26:28.386285Z","iopub.status.idle":"2023-12-27T19:26:58.562954Z","shell.execute_reply.started":"2023-12-27T19:26:28.386244Z","shell.execute_reply":"2023-12-27T19:26:58.561687Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Seaborn scatterplots","metadata":{}},{"cell_type":"code","source":"%%time\nx = adata.obsm['X_umap'][:,0]\ny = adata.obsm['X_umap'][:,1]\nsns.scatterplot(x=x,y=y, hue = adata.obs['cell_type'])\nplt.show()\n# %%time\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() )\n    \nx = adata.obsm['X_umap'][:,0]\ny = adata.obsm['X_umap'][:,1]\nfor col in list_selected_drugs:\n    sns.scatterplot(x=x,y=y, hue = adata.obs[col])\n    plt.show()    ","metadata":{"execution":{"iopub.status.busy":"2023-12-28T22:57:28.015422Z","iopub.execute_input":"2023-12-28T22:57:28.015965Z","iopub.status.idle":"2023-12-28T22:59:32.236747Z","shell.execute_reply.started":"2023-12-28T22:57:28.015931Z","shell.execute_reply":"2023-12-28T22:59:32.235360Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Other drugs - most cannot see on umap - since not so many cells and effect is not so strong to produce separate cluster","metadata":{}},{"cell_type":"code","source":"%%time\nimport warnings\nwarnings.simplefilter(\"ignore\")\n\nprint(adata.obs['sm_name'].value_counts().head(10))\nprint(adata.obs['sm_name'].value_counts().tail(10))\n\n#print(adata.obs['sm_name'].unique().to_list() )\n\nlist_selected_drugs = adata.obs['sm_name'].value_counts().index.to_list() # ['Dimethyl Sulfoxide', 'Dabrafenib', 'Belinostat', 'CGM-097', 'LY2090314','Ganetespib (STA-9090)','IN1451', 'Oprozomib (ONX 0912)','MLN 2238','CEP-18770 (Delanzomib)',]\n\nprint( list_selected_drugs )\n\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() )\n    \nsc.pl.umap(adata, color= list_selected_drugs ) # n_genes', 'n_genes_by_counts', 'total_counts'])\n    ","metadata":{"execution":{"iopub.status.busy":"2023-12-28T19:02:59.310349Z","iopub.execute_input":"2023-12-28T19:02:59.310853Z","iopub.status.idle":"2023-12-28T19:06:33.197955Z","shell.execute_reply.started":"2023-12-28T19:02:59.310815Z","shell.execute_reply":"2023-12-28T19:06:33.196372Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Seaborn scatterplots","metadata":{}},{"cell_type":"code","source":"%%time\nimport warnings\nwarnings.simplefilter(\"ignore\")\n\nx = adata.obsm['X_umap'][:,0]\ny = adata.obsm['X_umap'][:,1]\nsns.scatterplot(x=x,y=y, hue = adata.obs['cell_type'])\nplt.show()\n# %%time\n# list_selected_drugs = ['Dimethyl Sulfoxide', 'Dabrafenib', 'Belinostat', 'CGM-097', 'LY2090314','Ganetespib (STA-9090)','IN1451', 'Oprozomib (ONX 0912)','MLN 2238','CEP-18770 (Delanzomib)',]\nlist_selected_drugs = adata.obs['sm_name'].value_counts().head(12).index.to_list() # ['Dimethyl Sulfoxide', 'Dabrafenib', 'Belinostat', 'CGM-097', 'LY2090314','Ganetespib (STA-9090)','IN1451', 'Oprozomib (ONX 0912)','MLN 2238','CEP-18770 (Delanzomib)',]\n\nx = adata.obsm['X_umap'][:,0]\ny = adata.obsm['X_umap'][:,1]\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() )\n    sns.scatterplot(x=x,y=y, hue = adata.obs[drug])\n    plt.show()        ","metadata":{"execution":{"iopub.status.busy":"2023-12-28T23:10:09.007482Z","iopub.execute_input":"2023-12-28T23:10:09.007944Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Some drug groups","metadata":{}},{"cell_type":"code","source":"l = [t[-3:] for t in adata.obs['sm_name'].unique() ]\npd.Series(l).value_counts().head(10)","metadata":{"execution":{"iopub.status.busy":"2023-12-28T19:12:10.934505Z","iopub.execute_input":"2023-12-28T19:12:10.934917Z","iopub.status.idle":"2023-12-28T19:12:10.945738Z","shell.execute_reply.started":"2023-12-28T19:12:10.934880Z","shell.execute_reply":"2023-12-28T19:12:10.944647Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"[t for t in adata.obs['sm_name'].unique() if t[-3:] == 'one' ], [t for t in adata.obs['sm_name'].unique() if t[-3:] == 'ine' ],  [t for t in adata.obs['sm_name'].unique() if t[-3:] == 'tat'], [t for t in adata.obs['sm_name'].unique() if t[-3:] == 'sib'], ","metadata":{"execution":{"iopub.status.busy":"2023-12-28T19:28:37.115302Z","iopub.execute_input":"2023-12-28T19:28:37.115737Z","iopub.status.idle":"2023-12-28T19:28:37.132159Z","shell.execute_reply.started":"2023-12-28T19:28:37.115706Z","shell.execute_reply":"2023-12-28T19:28:37.131112Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndict_drug_groups = {}\ndict_drug_groups['tyrosine kinase inhibitors'] = [t for t in adata.obs['sm_name'].unique() if t[-3:] == 'nib' ]\nprint(dict_drug_groups['tyrosine kinase inhibitors'])\n# k = 'statins(Lower cholesterol in blood. Inhibit enzyme HMG-CoA reductase)' - mistake of chatGPT\nk  = 'Belinostat Vorinostat Resminostat Ricolinostat Tosedostat'\ndict_drug_groups[k] = [t for t in adata.obs['sm_name'].unique() if t[-4:] == 'stat' ]\nprint(dict_drug_groups[k])\nk  = 'Vorinostat Resminostat Ricolinostat Tosedostat'\ndict_drug_groups[k] = ['Vorinostat','Resminostat','Ricolinostat','Tosedostat' ]\nprint(dict_drug_groups[k])\nk  = 'Idelalisib Dactolisib (PI3K inhibitors)'\ndict_drug_groups[k] = [t for t in adata.obs['sm_name'].unique() if t[-3:] == 'sib' ]\nprint(dict_drug_groups[k])\n\n# k = 'tyrosine kinase inhibitors'\nfor k in dict_drug_groups.keys(): \n    print(k)\n    l = dict_drug_groups[k]\n    print(l)\n    v = adata.obs['sm_name'].isin(l).astype('category')\n    adata.obs[k] = v\n    sc.pl.umap(adata, color= k ) # n_genes', 'n_genes_by_counts', 'total_counts'])\n    plt.show()\n\n","metadata":{"execution":{"iopub.status.busy":"2023-12-28T19:59:43.246773Z","iopub.execute_input":"2023-12-28T19:59:43.247214Z","iopub.status.idle":"2023-12-28T19:59:47.913613Z","shell.execute_reply.started":"2023-12-28T19:59:43.247181Z","shell.execute_reply":"2023-12-28T19:59:47.912644Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Cell (prolifiration) Cycle G1S/G2M plot  - small number of proliferating cells\n\nProliferation cycle or cell cycle - key process for many cells. Many drugs here are anti-cancer - so should affect the ability of cells to proliferate. So might be interesting to look on drugs affect on these particular genes.\n\nTirosh et.al. proposed list of about 100 cell-cycle genes which are the most effectively seen by single cell data. It is good starting point.\n\nIn general there are much more genes related to the cell cycle, with various degree of \"relation\". Some list are cell cycle genes may contain thousands genes, but in fact in such huge lists most of the genes are related to cell cycle very weakly or these genes not good captured by single cell sequencing, while Tirosh list contains genes very strongly related to cell cycle and well captured by single cell technology. Moreover, many genes are associated to various biological processes in the cell, but Tirosh genes are mostly associated with cell cycle, not with the other processes - another argument why they are convenient to work with.\n\nOne may look at Computational challenges of cell cycle analysis using single cell transcriptomics Alexander Chervov, Andrei Zinovyev https://arxiv.org/abs/2208.05229\n\nAnd hundreds Kaggle notebooks/datasets related to that work e.g. that discussion: https://www.kaggle.com/competitions/open-problems-multimodal/discussion/350314 , or that notebook: https://www.kaggle.com/code/alexandervc/mmscel-cell-cycle-03b-daybydaychange-allcelltypes/notebook\n\nPS\n\nOther cell cycle genes sets e.g. by Tom Freeman\n\nSee some comparison e.g. here: https://www.kaggle.com/code/alexandervc/tabmuris-cell-cycle-1-data-one-by-one#Cell-cycle So list is bigger but kind of more \"dirty\", that means containing more non-cell cycle effects, and G1S - G2M split is less prominent.","metadata":{}},{"cell_type":"code","source":"G1S_genes_Tirosh = ['MCM5', 'PCNA', 'TYMS', 'FEN1', 'MCM2', 'MCM4', 'RRM1', 'UNG', 'GINS2', 'MCM6', 'CDCA7', 'DTL', 'PRIM1', 'UHRF1', 'MLF1IP', 'HELLS', 'RFC2', 'RPA2', 'NASP', 'RAD51AP1', 'GMNN', 'WDR76', 'SLBP', 'CCNE2', 'UBR7', 'POLD3', 'MSH2', 'ATAD2', 'RAD51', 'RRM2', 'CDC45', 'CDC6', 'EXO1', 'TIPIN', 'DSCC1', 'BLM', 'CASP8AP2', 'USP1', 'CLSPN', 'POLA1', 'CHAF1B', 'BRIP1', 'E2F8']\nG2M_genes_Tirosh = ['HMGB2', 'CDK1', 'NUSAP1', 'UBE2C', 'BIRC5', 'TPX2', 'TOP2A', 'NDC80', 'CKS2', 'NUF2', 'CKS1B', 'MKI67', 'TMPO', 'CENPF', 'TACC3', 'FAM64A', 'SMC4', 'CCNB2', 'CKAP2L', 'CKAP2', 'AURKB', 'BUB1', 'KIF11', 'ANP32E', 'TUBB4B', 'GTSE1', 'KIF20B', 'HJURP', 'CDCA3', 'HN1', 'CDC20', 'TTK', 'CDC25C', 'KIF2C', 'RANGAP1', 'NCAPD2', 'DLGAP5', 'CDCA2', 'CDCA8', 'ECT2', 'KIF23', 'HMMR', 'AURKA', 'PSRC1', 'ANLN', 'LBR', 'CKAP5', 'CENPE', 'CTCF', 'NEK2', 'G2E3', 'GAS2L3', 'CBX5', 'CENPA']\ngenes_Tirosh = G1S_genes_Tirosh + G2M_genes_Tirosh\n\n# Subset of Tirosh genes to capture \"fast\" cell cycle pattern = see https://arxiv.org/abs/2208.05229\nlist_genes_fastCCsign = ['CDK1', 'UBE2C', 'TOP2A', 'TMPO', 'HJURP', 'RRM1', 'RAD51AP1', 'RRM2', 'CDC45', 'BLM', 'BRIP1', 'E2F8', 'HIST2H2AC']\n\nG1S_genes_Freeman = ['ADAMTS1', 'ASF1B', 'ATAD2', 'BARD1', 'BLM', 'BRCA1', 'BRIP1', 'C17orf75', 'C9orf40', 'CACYBP', 'CASP8AP2', 'CCDC15', 'CCNE1', 'CCNE2', 'CCP110', 'CDC25A', 'CDC45', 'CDC6', 'CDC7', 'CDK2', 'CDT1', 'CENPJ', 'CENPQ', 'CENPU', 'CEP57', 'CHAF1A', 'CHAF1B', 'CHEK1', 'CLSPN', 'CREBZF', 'CRYL1', 'CSE1L', 'DCLRE1B', 'DCTPP1', 'DEK', 'DERA', 'DHFR', 'DNA2', 'DNAJC9', 'DNMT1', 'DONSON', 'DSCC1', 'DSN1', 'DTL', 'E2F8', 'EED', 'EFCAB11', 'ENDOD1', 'ETAA1', 'EXO1', 'EYA2', 'EZH2', 'FAM111A', 'FANCE', 'FANCG', 'FANCI', 'FANCL', 'FBXO5', 'FEN1', 'GGH', 'GINS1', 'GINS2', 'GINS3', 'GLMN', 'GMNN', 'GMPS', 'GPD2', 'HADH', 'HELLS', 'HSF2', 'ITGB3BP', 'KIAA0101', 'KNTC1', 'LIG1', 'MCM10', 'MCM2', 'MCM3', 'MCM4', 'MCM5', 'MCM6', 'MCM7', 'MCMBP', 'METTL9', 'MMD', 'MNS1', 'MPP1', 'MRE11A', 'MSH2', 'MSH6', 'MYO19', 'NASP', 'NPAT', 'NSMCE4A', 'ORC1', 'OSGEPL1', 'PAK1', 'PAQR4', 'PARP2', 'PASK', 'PAXIP1', 'PBX3', 'PCNA', 'PKMYT1', 'PMS1', 'POLA1', 'POLA2', 'POLD3', 'POLE2', 'PRIM1', 'PRPS2', 'PSMC3IP', 'RAB23', 'RAD51', 'RAD51AP1', 'RAD54L', 'RBBP8', 'RBL1', 'RDX', 'RFC2', 'RFC3', 'RFC4', 'RMI1', 'RNASEH2A', 'RPA1', 'RRM1', 'RRM2', 'SLBP', 'SLC25A40', 'SMC2', 'SMC3', 'SSX2IP', 'SUPT16H', 'TEX30', 'TFDP1', 'THAP10', 'THEM6', 'TIMELESS', 'TIPIN', 'TMEM106C', 'TMEM38B', 'TRIM45', 'TRIP13', 'TSPYL4', 'TTI1', 'TUBGCP5', 'TYMS', 'UBR7', 'UNG', 'USP1', 'WDHD1', 'WDR76', 'WRB', 'YEATS4', 'ZBTB14', 'ZWINT']\nG2M_genes_Freeman = ['ADGRE5', 'ARHGAP11A', 'ARHGDIB', 'ARL6IP1', 'ASPM', 'AURKA', 'AURKB', 'BIRC5', 'BORA', 'BRD8', 'BUB1', 'BUB1B', 'BUB3', 'CCNA2', 'CCNB1', 'CCNB2', 'CCNF', 'CDC20', 'CDC25B', 'CDC25C', 'CDC27', 'CDCA3', 'CDCA8', 'CDK1', 'CDKN1B', 'CDKN3', 'CENPE', 'CENPF', 'CENPI', 'CENPN', 'CEP55', 'CEP70', 'CEP85', 'CKAP2', 'CKAP5', 'CKS1B', 'CKS2', 'CTCF', 'DBF4', 'DBF4B', 'DCAF7', 'DEPDC1', 'DLGAP5', 'ECT2', 'ERCC6L', 'ESPL1', 'FAM64A', 'FOXM1', 'FZD2', 'FZD7', 'FZR1', 'GPSM2', 'GTF2E1', 'GTSE1', 'H2AFX', 'HJURP', 'HMGB2', 'HMGB3', 'HMMR', 'HN1', 'INCENP', 'JADE2', 'KIF11', 'KIF14', 'KIF15', 'KIF18A', 'KIF18B', 'KIF20A', 'KIF20B', 'KIF22', 'KIF23', 'KIF2C', 'KIF4A', 'KIF5B', 'KIFC1', 'KPNA2', 'LBR', 'LMNB2', 'MAD2L1', 'MELK', 'MET', 'METTL4', 'MIS18BP1', 'MKI67', 'MPHOSPH9', 'MTMR6', 'NCAPD2', 'NCAPG', 'NCAPG2', 'NCAPH', 'NDC1', 'NDC80', 'NDE1', 'NEIL3', 'NEK2', 'NRF1', 'NUSAP1', 'OIP5', 'PAFAH2', 'PARPBP', 'PBK', 'PLEKHG3', 'PLK1', 'PLK4', 'PRC1', 'PRR11', 'PSRC1', 'PTTG1', 'PTTG3P', 'RACGAP1', 'RAD21', 'RASSF1', 'REEP4', 'SAP30', 'SHCBP1', 'SKA1', 'SLCO1B3', 'SOGA1', 'SPA17', 'SPAG5', 'SPC25', 'SPDL1', 'STIL', 'STK17B', 'TACC3', 'TAF5', 'TBC1D2', 'TBC1D31', 'TMPO', 'TOP2A', 'TPX2', 'TROAP', 'TTF2', 'TTK', 'TUBB4B', 'TUBD1', 'UBE2C', 'UBE2S', 'VANGL1', 'WEE1', 'WHSC1', 'XPO1', 'ZMYM1']\n","metadata":{"execution":{"iopub.status.busy":"2023-12-28T20:04:55.220211Z","iopub.execute_input":"2023-12-28T20:04:55.220585Z","iopub.status.idle":"2023-12-28T20:04:55.248563Z","shell.execute_reply.started":"2023-12-28T20:04:55.220556Z","shell.execute_reply":"2023-12-28T20:04:55.247418Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## G1S-G2M plot - Tirosh signatures - standard ones ","metadata":{}},{"cell_type":"code","source":"%%time\nG1S = np.asarray(adata[:,list(set(G1S_genes_Tirosh)&set(adata.var.index))].X.mean(axis = 1)).ravel()\nG2M = np.asarray(adata[:,list(set(G2M_genes_Tirosh)&set(adata.var.index))].X.mean(axis = 1)).ravel()\n\nplt.figure(figsize = (15,6) )\nsns.scatterplot(x= G1S, y=G2M, hue = adata.obs['cell_type'])#  , palette = 'rainbow' )\nplt.title('G1S/G2M plot. Tirosh  signatures - standard ones.', fontsize= 20)\nplt.xlabel('G1S', fontsize= 20)\nplt.ylabel('G2M', fontsize= 20)\nplt.grid()\nplt.show()\n\nc = np.round( np.corrcoef(G1S,G2M)[0,1] ,2 )\nprint('Correlation G1S-G2M', c )\n\nfor col in ['control', 'control negative', 'control positive', 'sm_name']:\n    plt.figure(figsize = (15,6) )\n    sns.scatterplot(x= G1S, y=G2M, hue = adata.obs[col])#  , palette = 'rainbow' )\n    plt.title('G1S/G2M plot. Tirosh  signatures - standard ones.', fontsize= 20)\n    plt.xlabel('G1S', fontsize= 20)\n    plt.ylabel('G2M', fontsize= 20)\n    plt.grid()\n    plt.show()\n\n# plt.figure(figsize = (15,6) )\n# sns.scatterplot(x= G1S, y=G2M, hue = adata.obs['sm_name'])#  , palette = 'rainbow' )\n# plt.title('G1S/G2M plot. Tirosh  signatures - standard ones.', fontsize= 20)\n# plt.xlabel('G1S', fontsize= 20)\n# plt.ylabel('G2M', fontsize= 20)\n# plt.grid()\n# plt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-12-28T20:04:56.202985Z","iopub.execute_input":"2023-12-28T20:04:56.203367Z","iopub.status.idle":"2023-12-28T20:05:46.358088Z","shell.execute_reply.started":"2023-12-28T20:04:56.203338Z","shell.execute_reply":"2023-12-28T20:05:46.357043Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Freeman G1S-G2M signatures\n\nNote: typically Freeman and Tirosh singatures give similar results,\n\nBut Tirosh is in a sense smaller, more restrictive and cleaner set of genes which capture mainly cell cycle effects.\n\nFreeman is more wider set of genes capturing not only cell cycle effects (which is not good for the task), but for very noisy data we one may hope averaging more genes might help to reduce noise. Thus Tirosh is the first (and the standard) choice, especially for the high quality data, Freeman can be checked for very noisy data.\n\nIn general Freeman G1S and G2M are more correlated than Tirosh - which is not good for the G1S/G2M plot, but it is a price to pay including more genes which might be helpful for noisy data.","metadata":{}},{"cell_type":"code","source":"%%time\nG1S = np.asarray(adata[:,list(set(G1S_genes_Freeman)&set(adata.var.index))].X.mean(axis = 1)).ravel()\nG2M = np.asarray(adata[:,list(set(G2M_genes_Freeman)&set(adata.var.index))].X.mean(axis = 1)).ravel()\nprint( G1S.shape, G2M.shape)\n\nplt.figure(figsize = (15,6) )\nsns.scatterplot(x= G1S, y=G2M)\nplt.title('G1S/G2M plot. Freeman signatures', fontsize= 20)\nplt.xlabel('G1S', fontsize= 20)\nplt.ylabel('G2M', fontsize= 20)\nplt.grid()\nplt.show()\n\nc = np.round( np.corrcoef(G1S,G2M)[0,1] ,2 )\nprint('Correlation G1S-G2M', c )","metadata":{"execution":{"iopub.status.busy":"2023-12-27T17:08:06.485465Z","iopub.execute_input":"2023-12-27T17:08:06.486094Z","iopub.status.idle":"2023-12-27T17:08:10.609881Z","shell.execute_reply.started":"2023-12-27T17:08:06.486048Z","shell.execute_reply":"2023-12-27T17:08:10.608759Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Fast signature plot\n\nUseful to visualize \"fast\" cell cycle (it is not standard pattern of the cell cycle appearing for embryonic stem cells and some cancer cells (in particular many tp53-mutated ) ).  See https://arxiv.org/abs/2208.05229","metadata":{}},{"cell_type":"code","source":"# Subset of Tirosh genes to capture \"fast\" cell cycle pattern = see https://arxiv.org/abs/2208.05229\nlist_genes_fastCCsign = ['CDK1', 'UBE2C', 'TOP2A', 'TMPO', 'HJURP', 'RRM1', 'RAD51AP1', 'RRM2', 'CDC45', 'BLM', 'BRIP1', 'E2F8', 'HIST2H2AC']\n\n# G1S = np.asarray(adata[:,list(set(G1S_genes_Tirosh)&set(adata.var.index))].X.mean(axis = 1)).ravel()\nG2M = np.asarray(adata[:,list(set(G2M_genes_Tirosh)&set(adata.var.index))].X.mean(axis = 1)).ravel()\nfastCC = np.asarray(adata[:,list(set(list_genes_fastCCsign)&set(adata.var.index))].X.mean(axis = 1)).ravel()\nprint( G1S.shape, fastCC.shape)\n\nplt.figure(figsize = (15,6) )\nsns.scatterplot(x= fastCC, y=G2M, hue = adata.obs['cell_type'])\nplt.title('fastCC-G2M signatures', fontsize= 20)\nplt.xlabel('fastCC', fontsize= 20)\nplt.ylabel('G2M', fontsize= 20)\nplt.grid()\nplt.show()\n\nc = np.round( np.corrcoef(G1S,fastCC)[0,1] ,2 )\nprint('Correlation G1S-fastCC', c )","metadata":{"execution":{"iopub.status.busy":"2023-12-27T18:04:05.691294Z","iopub.execute_input":"2023-12-27T18:04:05.692368Z","iopub.status.idle":"2023-12-27T18:04:17.991732Z","shell.execute_reply.started":"2023-12-27T18:04:05.692328Z","shell.execute_reply":"2023-12-27T18:04:17.990340Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## UMAP Colored by cell cycle signatures \n\nVery few cells are proliferating -  the cluster is very small. \n","metadata":{}},{"cell_type":"code","source":"G1S = np.asarray(adata[:,list(set(G1S_genes_Tirosh)&set(adata.var.index))].X.mean(axis = 1)).ravel()\nG2M = np.asarray(adata[:,list(set(G2M_genes_Tirosh)&set(adata.var.index))].X.mean(axis = 1)).ravel()\nadata.obs['G1S'] = G1S\nadata.obs['G2M'] = G2M\nadata.obs['G1S_bin'] = (G1S > 0.4)\nadata.obs['G1S_bin'] = adata.obs['G1S_bin'].apply(lambda x: str(x)).astype('category')\n\n\nsc.pl.umap(adata, color=['G1S_bin', 'G1S', 'G2M',])","metadata":{"execution":{"iopub.status.busy":"2023-12-27T18:26:03.051227Z","iopub.execute_input":"2023-12-27T18:26:03.051911Z","iopub.status.idle":"2023-12-27T18:26:10.536993Z","shell.execute_reply.started":"2023-12-27T18:26:03.051871Z","shell.execute_reply":"2023-12-27T18:26:10.535743Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nx = adata.obsm['X_umap'][:,0]\ny = adata.obsm['X_umap'][:,1]\nsns.scatterplot(x=x,y=y, hue = adata.obs['G1S_bin'])\nplt.show()","metadata":{"trusted":true},"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":{"iopub.status.busy":"2023-12-27T17:08:14.464383Z","iopub.execute_input":"2023-12-27T17:08:14.464801Z","iopub.status.idle":"2023-12-27T17:08:14.472892Z","shell.execute_reply.started":"2023-12-27T17:08:14.464758Z","shell.execute_reply":"2023-12-27T17:08:14.471418Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}