{"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":"# What is about ?\n\nEDA for Kaggle competition: Multimodal Single-Cell Integration Competition\n\n(c) Alexander Chervov\n\nWork in progress. We will move towards standard ( +- variations)  bionformatics analysis for single cell data. \n\n(For single cell RNA sequencing part - see e.g. several Scanpy tutorials: \nhttps://scanpy.readthedocs.io/en/stable/tutorials.html , some of them can be found on Kaggle: \nhttps://www.kaggle.com/datasets/alexandervc/scanpy-python-package-for-scrnaseq-analysis .\nIf you are \"R\"-user - google for \"Seurat\") \n\n#### Further works:\n\nContinuation on same part of data: https://www.kaggle.com/alexandervc/mmscel-eda-targets-citeseq-02\n\n\n\n#### Versions:\n\n#### 2,3,4,5, 6 cosmetic changes\n\n#### 1 CiteSeq targets analysis\n\nLook on only CiteSeq targets: Proteomics  (cell surface protein markers labeled mainly CD**). CD means \"Cluster of Differentiation\" - see https://en.wikipedia.org/wiki/Cluster_of_differentiation . \nAs expected we see many of them quite correlate with the cell types. \n\nSome correlated targets are: ['CD71' ,'CD115', 'CD88'] and ['CD155', 'CD112', 'CD47', 'HLA-A-B-C', 'CD45RA', 'CD31', 'CD11a', 'CD13', 'CD29',   'CD81', 'CD18', 'CD45', 'CD49d', 'CD162']\n","metadata":{"execution":{"iopub.status.busy":"2022-08-11T21:21:49.972077Z","iopub.execute_input":"2022-08-11T21:21:49.972474Z","iopub.status.idle":"2022-08-11T21:21:49.977181Z","shell.execute_reply.started":"2022-08-11T21:21:49.972443Z","shell.execute_reply":"2022-08-11T21:21:49.975948Z"}}},{"cell_type":"markdown","source":"#  Install/Import packages\n\n","metadata":{}},{"cell_type":"code","source":"#If you see a urllib warning running this cell, go to \"Settings\" on the right hand side, \n#and turn on internet. Note, you need to be phone verified.\n!pip install --quiet tables","metadata":{"execution":{"iopub.status.busy":"2022-09-23T13:36:11.112012Z","iopub.execute_input":"2022-09-23T13:36:11.112549Z","iopub.status.idle":"2022-09-23T13:36:25.579577Z","shell.execute_reply.started":"2022-09-23T13:36:11.112442Z","shell.execute_reply":"2022-09-23T13:36:25.578236Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import time \nt0start = time.time()\n\nimport os\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport seaborn as sns","metadata":{"execution":{"iopub.status.busy":"2022-09-23T16:14:17.945663Z","iopub.execute_input":"2022-09-23T16:14:17.946407Z","iopub.status.idle":"2022-09-23T16:14:17.952436Z","shell.execute_reply.started":"2022-09-23T16:14:17.946364Z","shell.execute_reply":"2022-09-23T16:14:17.950970Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 1.2. Set filepaths","metadata":{"execution":{"iopub.status.busy":"2022-08-11T21:22:04.888845Z","iopub.execute_input":"2022-08-11T21:22:04.889239Z","iopub.status.idle":"2022-08-11T21:22:04.893408Z","shell.execute_reply.started":"2022-08-11T21:22:04.889209Z","shell.execute_reply":"2022-08-11T21:22:04.892357Z"}}},{"cell_type":"code","source":"os.listdir(\"/kaggle/input/open-problems-multimodal/\")","metadata":{"execution":{"iopub.status.busy":"2022-09-23T13:36:26.286327Z","iopub.execute_input":"2022-09-23T13:36:26.287124Z","iopub.status.idle":"2022-09-23T13:36:26.298429Z","shell.execute_reply.started":"2022-09-23T13:36:26.287039Z","shell.execute_reply":"2022-09-23T13:36:26.297178Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DATA_DIR = \"/kaggle/input/open-problems-multimodal/\"\nFP_CELL_METADATA = os.path.join(DATA_DIR,\"metadata.csv\")\n\nFP_CITE_TRAIN_INPUTS = os.path.join(DATA_DIR,\"train_cite_inputs.h5\")\nFP_CITE_TRAIN_TARGETS = os.path.join(DATA_DIR,\"train_cite_targets.h5\")\nFP_CITE_TEST_INPUTS = os.path.join(DATA_DIR,\"test_cite_inputs.h5\")\n\nFP_MULTIOME_TRAIN_INPUTS = os.path.join(DATA_DIR,\"train_multi_inputs.h5\")\nFP_MULTIOME_TRAIN_TARGETS = os.path.join(DATA_DIR,\"train_multi_targets.h5\")\nFP_MULTIOME_TEST_INPUTS = os.path.join(DATA_DIR,\"test_multi_inputs.h5\")\n\nFP_SUBMISSION = os.path.join(DATA_DIR,\"sample_submission.csv\")\nFP_EVALUATION_IDS = os.path.join(DATA_DIR,\"evaluation_ids.csv\")\n\ndf_cell = pd.read_csv(FP_CELL_METADATA)\ndf_cell","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-09-23T13:40:19.289203Z","iopub.execute_input":"2022-09-23T13:40:19.290757Z","iopub.status.idle":"2022-09-23T13:40:19.778524Z","shell.execute_reply.started":"2022-09-23T13:40:19.290674Z","shell.execute_reply":"2022-09-23T13:40:19.777184Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Load Data","metadata":{"execution":{"iopub.status.busy":"2022-08-11T21:26:56.381326Z","iopub.execute_input":"2022-08-11T21:26:56.381727Z","iopub.status.idle":"2022-08-11T21:26:56.386286Z","shell.execute_reply.started":"2022-08-11T21:26:56.381695Z","shell.execute_reply":"2022-08-11T21:26:56.385262Z"}}},{"cell_type":"code","source":"str_data_inf = ' CITEseq Targets '","metadata":{"execution":{"iopub.status.busy":"2022-09-23T14:07:57.231050Z","iopub.execute_input":"2022-09-23T14:07:57.231919Z","iopub.status.idle":"2022-09-23T14:07:57.237655Z","shell.execute_reply.started":"2022-09-23T14:07:57.231875Z","shell.execute_reply":"2022-09-23T14:07:57.235973Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndf_cite_train_y = pd.read_hdf('../input/open-problems-multimodal/train_cite_targets.h5')\ndisplay(df_cite_train_y.head(2))\n\ndf_cell = pd.read_csv(FP_CELL_METADATA, index_col = 0)\nd = df_cell.join(df_cite_train_y, how = 'right')#.join()\nprint(d.shape, df_cite_train_y.shape )\nprint(df_cell['cell_type'].value_counts())\ndisplay(df_cell.head(2))\n\n","metadata":{"execution":{"iopub.status.busy":"2022-09-23T13:41:39.236202Z","iopub.execute_input":"2022-09-23T13:41:39.236667Z","iopub.status.idle":"2022-09-23T13:41:40.530453Z","shell.execute_reply.started":"2022-09-23T13:41:39.236632Z","shell.execute_reply":"2022-09-23T13:41:40.528928Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Cell types: \n\n    MasP = Mast Cell Progenitor\n    MkP = Megakaryocyte Progenitor\n    NeuP = Neutrophil Progenitor\n    MoP = Monocyte Progenitor\n    EryP = Erythrocyte Progenitor\n    HSC = Hematoploetic Stem Cell\n    BP = B-Cell Progenitor","metadata":{}},{"cell_type":"markdown","source":"# Dimensional reduction with PCA","metadata":{}},{"cell_type":"code","source":"X = df_cite_train_y.values\nlist_X_column_names = list(df_cite_train_y.columns)\n","metadata":{"execution":{"iopub.status.busy":"2022-09-23T16:15:21.572628Z","iopub.execute_input":"2022-09-23T16:15:21.573727Z","iopub.status.idle":"2022-09-23T16:15:21.579680Z","shell.execute_reply.started":"2022-09-23T16:15:21.573682Z","shell.execute_reply":"2022-09-23T16:15:21.578452Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nfrom sklearn.decomposition import PCA\nimport time \nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\npca = PCA(n_components=50)\nt0 = time.time()\nr = pca.fit_transform(X)\nprint(time.time()-t0, 'secs passed for PCA')\n\nfig = plt.figure(figsize = (15,7) )\nc = 0\n# for f in ['cp_dose', 'cp_type','cp_time', 'y_sum']:\n#     c+=1; fig.add_subplot(1, 4 , c) \nsns.scatterplot(x=r[:,0], y=r[:,1])# , hue = df[f]  )\n#    plt.title('Colored by '+f)\nplt.show()\n\nfig = plt.figure(figsize = (15,7) )\nfig.add_subplot(1, 2, 1) \nplt.plot(pca.singular_values_,'o-')\nplt.title('Singular values')\nfig.add_subplot(1, 2, 2) \nplt.plot(pca.explained_variance_ratio_,'o-')\nplt.title('explained variance')\n","metadata":{"execution":{"iopub.status.busy":"2022-08-21T21:03:08.321842Z","iopub.execute_input":"2022-08-21T21:03:08.322698Z","iopub.status.idle":"2022-08-21T21:03:12.141282Z","shell.execute_reply.started":"2022-08-21T21:03:08.322645Z","shell.execute_reply":"2022-08-21T21:03:12.139514Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize = (20,10) )\nc = 0\nfor f in ['day', 'donor', 'cell_type']:# , 'technology']:\n    c+=1; fig.add_subplot(1, 3 , c) \n    sns.scatterplot(x=r[:,0], y=r[:,1] , hue = d[f]  )\n    plt.title('Colored by '+f)\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2022-08-21T21:03:12.154049Z","iopub.execute_input":"2022-08-21T21:03:12.154412Z","iopub.status.idle":"2022-08-21T21:03:21.163295Z","shell.execute_reply.started":"2022-08-21T21:03:12.154381Z","shell.execute_reply":"2022-08-21T21:03:21.162017Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Correlation analysis","metadata":{}},{"cell_type":"code","source":"t0 = time.time()\ncorr_matr = np.corrcoef(X.T) # Hint - use numpy , pandas is MUCH SLOWER   (df.corr() )\nprint(time.time() - t0, 'seconds passed')\nprint(np.min(corr_matr ), 'minimal correlation' )\ncorr_matr_abs = np.abs( corr_matr )\nprint(np.mean(corr_matr_abs ), 'average absolute correlation' )\nprint(np.median(corr_matr_abs), 'median absolute correlation' )\nprint(np.min(corr_matr_abs ), 'min absolute correlation' )\nprint(np.std(corr_matr_abs ), 'std absolute correlation' )\n\nv = corr_matr.flatten()\nplt.figure(figsize=(14,8))\nt0 = time.time()\nplt.hist(v, bins = 50)\nplt.title('correlation coefficients distribution')\nplt.show()\nprint(time.time() - t0, 'seconds passed')\n\n\nv.shape","metadata":{"execution":{"iopub.status.busy":"2022-08-21T21:03:21.167458Z","iopub.execute_input":"2022-08-21T21:03:21.167859Z","iopub.status.idle":"2022-08-21T21:03:21.767805Z","shell.execute_reply.started":"2022-08-21T21:03:21.167810Z","shell.execute_reply":"2022-08-21T21:03:21.766423Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_corr_matr = pd.DataFrame(corr_matr, index = list_X_column_names, columns = list_X_column_names )\n\n#clustermap\nimport seaborn as sns\nt0 = time.time()\nsns.clustermap(np.abs(df_corr_matr),cmap='vlag');\nprint( np.round(time.time()- t0,1), ' seconds passed.')","metadata":{"execution":{"iopub.status.busy":"2022-08-21T21:03:21.769284Z","iopub.execute_input":"2022-08-21T21:03:21.769646Z","iopub.status.idle":"2022-08-21T21:03:22.913528Z","shell.execute_reply.started":"2022-08-21T21:03:21.769614Z","shell.execute_reply":"2022-08-21T21:03:22.912458Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import igraph\n","metadata":{"execution":{"iopub.status.busy":"2022-08-21T21:03:22.914791Z","iopub.execute_input":"2022-08-21T21:03:22.915839Z","iopub.status.idle":"2022-08-21T21:03:22.920665Z","shell.execute_reply.started":"2022-08-21T21:03:22.915780Z","shell.execute_reply":"2022-08-21T21:03:22.919448Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"verbose = 0\ndf_stat = pd.DataFrame() # dict_save_largest_component_size = {} \ni = 0\nfor correlation_threshold in [0.8, 0.7, 0.6, 0.5, 0.4] :\n    t0 = time.time()\n    print()\n    print(correlation_threshold , 'correlation_threshold ')\n    corr_matr_abs_bool = corr_matr_abs > correlation_threshold\n    corr_matr_abs_bool = corr_matr_abs_bool # Restrict to  genes part \n    corr_matr_abs_bool = np.triu(corr_matr_abs_bool,1) # Take upper triangular part \n    g = igraph.Graph().Adjacency(corr_matr_abs_bool.tolist())\n    g.to_undirected(mode = 'collapse')\n    if verbose >= 10:\n        print( corr_matr_abs_bool.astype(int) )\n        print('Number of nodes ', g.vcount())\n        print('Number of edges ', g.ecount() )\n        print('Number of weakly connected compoenents', len( g.clusters(mode='WEAK')))\n\n\n    list_clusters_nodes_lists = list( g.clusters(mode='WEAK') )\n    list_clusers_size = [len(t) for t in list_clusters_nodes_lists ]\n    list_clusers_size = np.sort(list_clusers_size)[::-1]\n    print('Top 5 cluster sizes:', list_clusers_size[:5] , 'seconds passed:', np.round(time.time()-t0 , 2))\n    #dict_save_largest_component_size[correlation_threshold ] = list_clusers_size[0]\n    for t  in list_clusters_nodes_lists:\n        if len(t) == list_clusers_size[0]:\n            print('50 Genes in largest correlated group:')\n            print(np.array(list_X_column_names)[t[:50]])\n    i += 1\n    df_stat.loc[i,'correlation threshold'] = correlation_threshold\n    df_stat.loc[i,'Largest Component Size'] = list_clusers_size[0]\n    df_stat.loc[i,'Second Component Size'] = list_clusers_size[1]\n    \ndf_stat","metadata":{"execution":{"iopub.status.busy":"2022-08-21T21:03:22.921751Z","iopub.execute_input":"2022-08-21T21:03:22.922120Z","iopub.status.idle":"2022-08-21T21:03:22.966397Z","shell.execute_reply.started":"2022-08-21T21:03:22.922085Z","shell.execute_reply":"2022-08-21T21:03:22.964916Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Dimensional reduction with UMAP","metadata":{}},{"cell_type":"code","source":"import umap\n\nt0 = time.time()\nr = umap.UMAP().fit_transform(X)\nprint(time.time()-t0, 'secs passed for umap')\n\n\nfig = plt.figure(figsize = (15,7) )\n# c = 0\n# for f in ['cp_dose', 'cp_type','cp_time','y_sum']:\n#     c+=1; fig.add_subplot(1, 4 , c) \nsns.scatterplot(x=r[:,0], y=r[:,1])#  , hue = df[f]  )\n# plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-08-21T21:03:22.968452Z","iopub.execute_input":"2022-08-21T21:03:22.969708Z","iopub.status.idle":"2022-08-21T21:04:21.908302Z","shell.execute_reply.started":"2022-08-21T21:03:22.969645Z","shell.execute_reply":"2022-08-21T21:04:21.906948Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize = (20,10) )\nc = 0\nfor f in ['day', 'donor', 'cell_type']:# , 'technology']:\n    c+=1; fig.add_subplot(1, 3 , c) \n    sns.scatterplot(x=r[:,0], y=r[:,1] , hue = d[f]  )\n    plt.title('Colored by '+f)\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2022-08-21T21:04:21.910992Z","iopub.execute_input":"2022-08-21T21:04:21.911482Z","iopub.status.idle":"2022-08-21T21:04:28.375371Z","shell.execute_reply.started":"2022-08-21T21:04:21.911445Z","shell.execute_reply":"2022-08-21T21:04:28.374324Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"    MasP = Mast Cell Progenitor\n    MkP = Megakaryocyte Progenitor\n    NeuP = Neutrophil Progenitor\n    MoP = Monocyte Progenitor\n    EryP = Erythrocyte Progenitor\n    HSC = Hematoploetic Stem Cell\n    BP = B-Cell Progenitor","metadata":{}},{"cell_type":"code","source":"# 0.7 correlation_threshold \n# Top 5 cluster sizes: [3 2 1 1 1] seconds passed: 0.0\n# 50 Genes in largest correlated group:\n# ['CD71' 'CD115' 'CD88']\n\n# 0.6 correlation_threshold \n# Top 5 cluster sizes: [14  5  2  1  1] seconds passed: 0.0\n# 50 Genes in largest correlated group:\n# ['CD155' 'CD112' 'CD47' 'HLA-A-B-C' 'CD45RA' 'CD31' 'CD11a' 'CD13' 'CD29'\n#  'CD81' 'CD18' 'CD45' 'CD49d' 'CD162']\n\nimport numbers\n\nn_x_subplots = 8\npalette = 'rainbow'#'viridis'\n\nt0 = time.time()\n\nc = 0\nfor f in ['CD71', 'CD115' ,'CD88'] + ['CD155', 'CD112', 'CD47', 'HLA-A-B-C' ,'CD45RA' ,'CD31' ,'CD11a', 'CD13', 'CD29',\n  'CD81', 'CD18', 'CD45', 'CD49d', 'CD162']: # ['day', 'donor', 'cell_type']:# , 'technology']:\n    if c % n_x_subplots == 0:\n        if c > 0:\n            plt.show()\n        fig = plt.figure(figsize = (20,5) ); c = 0\n        #plt.suptitle(str(k),fontsize = 20 )# str_data_inf + ' n_cells: ' +str(mask.sum()) + ' RED > median expression, BLUE <= median ' )# +' ' + cell_type +' ' + drug )\n        \n    c += 1; fig.add_subplot(1,n_x_subplots ,c)    \n    \n    v = d[f]\n    if isinstance(v[0], numbers.Number):\n        v = np.clip(v,np.percentile(v,5), np.percentile(v,95) )\n        # v = (v > np.median(v)).astype(int)\n    ax = sns.scatterplot(x=r[:,0], y=r[:,1] , hue =  v , palette = palette )\n    plt.setp(ax.get_legend().get_texts(), fontsize=10) # for legend text\n    plt.setp(ax.get_legend().get_title(), fontsize=10) # for legend title    \n    plt.title('Colored by '+f, fontsize = 10)\n    \nplt.show()\nprint(time.time()-t0, 'secs passed')\n","metadata":{"execution":{"iopub.status.busy":"2022-08-21T21:04:28.376950Z","iopub.execute_input":"2022-08-21T21:04:28.378063Z","iopub.status.idle":"2022-08-21T21:05:44.029469Z","shell.execute_reply.started":"2022-08-21T21:04:28.378021Z","shell.execute_reply":"2022-08-21T21:05:44.028031Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numbers\n\nn_x_subplots = 8\npalette = 'rainbow'#'viridis'\n\nt0 = time.time()\n\nc = 0\nfor f in d.columns: # ['day', 'donor', 'cell_type']:# , 'technology']:\n    if c % n_x_subplots == 0:\n        if c > 0:\n            plt.show()\n        fig = plt.figure(figsize = (20,5) ); c = 0\n        #plt.suptitle(str(k),fontsize = 20 )# str_data_inf + ' n_cells: ' +str(mask.sum()) + ' RED > median expression, BLUE <= median ' )# +' ' + cell_type +' ' + drug )\n        \n    c += 1; fig.add_subplot(1,n_x_subplots ,c)    \n    \n    v = d[f]\n    if isinstance(v[0], numbers.Number):\n        v = np.clip(v,np.percentile(v,5), np.percentile(v,95) )\n        # v = (v > np.median(v)).astype(int)\n    ax = sns.scatterplot(x=r[:,0], y=r[:,1] , hue =  v , palette = palette )\n    plt.setp(ax.get_legend().get_texts(), fontsize=10) # for legend text\n    plt.setp(ax.get_legend().get_title(), fontsize=10) # for legend title    \n    plt.title('Colored by '+f, fontsize = 10)\n    \nplt.show()\nprint(time.time()-t0, 'secs passed')\n","metadata":{"execution":{"iopub.status.busy":"2022-08-21T21:05:44.031205Z","iopub.execute_input":"2022-08-21T21:05:44.031875Z","iopub.status.idle":"2022-08-21T21:16:31.217284Z","shell.execute_reply.started":"2022-08-21T21:05:44.031805Z","shell.execute_reply":"2022-08-21T21:16:31.215696Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Other dimensional reductions (from Sklearn)","metadata":{}},{"cell_type":"code","source":"# Based on: \n# https://scikit-learn.org/stable/auto_examples/manifold/plot_compare_methods.html#sphx-glr-auto-examples-manifold-plot-compare-methods-py\n# See also:\n# https://scikit-learn.org/stable/auto_examples/manifold/plot_lle_digits.html#\n\n\n\n# To speed-up reduce dimensions by PCA first\nX_save = X.copy( )\n#r = pca.fit_transform(X)\n#X = r[:1000,:20]\n\n\n\nimport umap \nfrom sklearn import manifold\nfrom sklearn.decomposition import PCA\nfrom sklearn.decomposition import FactorAnalysis\nfrom sklearn.decomposition import NMF\nfrom sklearn.decomposition import FastICA\nfrom sklearn.decomposition import FactorAnalysis\nfrom sklearn.decomposition import LatentDirichletAllocation\nfrom sklearn.ensemble import RandomTreesEmbedding\nfrom sklearn.random_projection import SparseRandomProjection\nfrom sklearn.discriminant_analysis import LinearDiscriminantAnalysis\n\nfrom sklearn.pipeline import make_pipeline\nfrom sklearn.decomposition import TruncatedSVD\n\n\nfrom collections import OrderedDict\nfrom functools import partial\nfrom matplotlib.ticker import NullFormatter\n\n\nn_neighbors = 10\nn_components = 2\n# Set-up manifold methods\nLLE = partial(manifold.LocallyLinearEmbedding,\n              n_neighbors, n_components, eigen_solver='auto')\n\nmethods = OrderedDict()\nmethods['PCA'] = PCA()\nmethods['umap'] = umap.UMAP(n_components = n_components)\nmethods['t-SNE'] = manifold.TSNE(n_components=n_components, init='pca', random_state=0)\nmethods['ICA'] = FastICA(n_components=n_components,         random_state=0)\nmethods['FA'] = FactorAnalysis(n_components=n_components, random_state=0)\n#methods['LLE'] = LLE(method='standard')\n#methods['Modified LLE'] = LLE(method='modified')\n#methods['Isomap'] = manifold.Isomap(n_neighbors, n_components)\nmethods['MDS'] = manifold.MDS(n_components, max_iter=100, n_init=1)\nmethods['SE'] = manifold.SpectralEmbedding(n_components=n_components,\n                                           n_neighbors=n_neighbors)\nmethods['NMF'] = NMF(n_components=n_components,  init='random', random_state=0) \nmethods['RandProj'] = SparseRandomProjection(n_components=n_components, random_state=42)\n\nrand_trees_embed = make_pipeline(RandomTreesEmbedding(n_estimators=200, random_state=0, max_depth=5), TruncatedSVD(n_components=n_components) )\nmethods['RandTrees'] = rand_trees_embed\nmethods['LatDirAll'] = LatentDirichletAllocation(n_components=n_components,  random_state=0)\n#methods['LTSA'] = LLE(method='ltsa') \n#methods['Hessian LLE'] = LLE(method='hessian') \n\nlist_fast_methods = ['PCA','umap','FA', 'NMF','RandProj','RandTrees'] # 'ICA',\nlist_slow_methods = ['t-SNE','LLE','Modified LLE','Isomap','MDS','SE','LatDirAll','LTSA','Hessian LLE']\n\n# transformer = NeighborhoodComponentsAnalysis(init='random',  n_components=2, random_state=0) # Cannot be applied since supervised - requires y \n# methods['LinDisA'] = LinearDiscriminantAnalysis(n_components=n_components)# Cannot be applied since supervised - requires y \n\n\n# Create figure\nfig = plt.figure(figsize=(25, 16))\n\n# Plot results\nc = 0\nfor i, (label, method) in enumerate(methods.items()):\n    if label not in  list_fast_methods :\n        continue\n        \n    t0 = time.time()\n    try:\n        r = method.fit_transform(X)\n    except:\n        print('Got Exception', label )\n        continue \n    t1 = time.time()\n    print(\"%s: %.2g sec\" % (label, t1 - t0))\n    c+=1\n    fig.add_subplot(2, 3 , c) \n    sns.scatterplot(x=r[:,0], y=r[:,1] , hue =  d['cell_type'])\n    plt.title(label )\n    plt.legend('')\n\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2022-08-21T21:16:31.219417Z","iopub.execute_input":"2022-08-21T21:16:31.219806Z","iopub.status.idle":"2022-08-21T21:18:00.048455Z","shell.execute_reply.started":"2022-08-21T21:16:31.219770Z","shell.execute_reply":"2022-08-21T21:18:00.047126Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Trimap - yet another anlogue of tsne/umap ","metadata":{}},{"cell_type":"code","source":"!pip install trimap \nimport trimap\n\n","metadata":{"execution":{"iopub.status.busy":"2022-08-21T21:18:00.051067Z","iopub.execute_input":"2022-08-21T21:18:00.051487Z","iopub.status.idle":"2022-08-21T21:18:11.467388Z","shell.execute_reply.started":"2022-08-21T21:18:00.051448Z","shell.execute_reply":"2022-08-21T21:18:11.465969Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"column4color = 'cell_type'","metadata":{"execution":{"iopub.status.busy":"2022-09-23T16:16:41.767865Z","iopub.execute_input":"2022-09-23T16:16:41.768341Z","iopub.status.idle":"2022-09-23T16:16:41.774529Z","shell.execute_reply.started":"2022-09-23T16:16:41.768305Z","shell.execute_reply":"2022-09-23T16:16:41.773343Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nreducer =  trimap.TRIMAP() #  umap.UMAP()\n\nt0 = time.time()\nif isinstance(X, pd.DataFrame):\n    r = reducer.fit_transform(X.values)\nelse:\n    r = reducer.fit_transform(X)\nprint(time.time()-t0, 'secs passed for dimension reduction')\n\n\nfig = plt.figure(figsize = (15,7) )\n# c = 0\n# for f in ['cp_dose', 'cp_type','cp_time','y_sum']:\n#     c+=1; fig.add_subplot(1, 4 , c) \nsns.scatterplot(x=r[:,0], y=r[:,1]  , hue = d[column4color]  )\nplt.title(str(reducer)[:6], fontsize = 20 )\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-08-21T21:18:11.479871Z","iopub.execute_input":"2022-08-21T21:18:11.480360Z","iopub.status.idle":"2022-08-21T21:19:49.350167Z","shell.execute_reply.started":"2022-08-21T21:18:11.480321Z","shell.execute_reply":"2022-08-21T21:19:49.348895Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# NCVis - similar to UMAP, bust faster (from Skolkovo team)","metadata":{}},{"cell_type":"code","source":"!pip install ncvis\nimport ncvis","metadata":{"execution":{"iopub.status.busy":"2022-09-23T16:05:17.727743Z","iopub.execute_input":"2022-09-23T16:05:17.728259Z","iopub.status.idle":"2022-09-23T16:05:53.987777Z","shell.execute_reply.started":"2022-09-23T16:05:17.728222Z","shell.execute_reply":"2022-09-23T16:05:53.986155Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"reducer =  ncvis.NCVis() # trimap.TRIMAP() #  umap.UMAP()\n\nt0 = time.time()\nr = reducer.fit_transform(X)\nprint(time.time()-t0, 'secs passed for dimension reduction')\n\n\nfig = plt.figure(figsize = (15,7) )\n# c = 0\n# for f in ['cp_dose', 'cp_type','cp_time','y_sum']:\n#     c+=1; fig.add_subplot(1, 4 , c) \nsns.scatterplot(x=r[:,0], y=r[:,1]  , hue = d[column4color]  )\nplt.title(str(reducer), fontsize = 20 )\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-09-23T16:16:47.537386Z","iopub.execute_input":"2022-09-23T16:16:47.538818Z","iopub.status.idle":"2022-09-23T16:17:30.408729Z","shell.execute_reply.started":"2022-09-23T16:16:47.538769Z","shell.execute_reply":"2022-09-23T16:17:30.407762Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"r.shape\nd2 = df_cite_train_y.copy()\nd2[1] = r[:,0]\nd2[2] = r[:,1]\ncm2 = d2.corr()\n#pd.DataFrame(np.corrcoef(r.T, df_cite_train_y.T) )","metadata":{"execution":{"iopub.status.busy":"2022-09-23T16:21:23.707329Z","iopub.execute_input":"2022-09-23T16:21:23.707779Z","iopub.status.idle":"2022-09-23T16:21:27.748629Z","shell.execute_reply.started":"2022-09-23T16:21:23.707742Z","shell.execute_reply":"2022-09-23T16:21:27.747352Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cm2[2].plot()\nv2 = cm2[2].sort_values()\nprint(v2.head(10) .index)\nprint(v2.tail(10) .index)\nv2","metadata":{"execution":{"iopub.status.busy":"2022-09-23T16:23:38.912027Z","iopub.execute_input":"2022-09-23T16:23:38.913170Z","iopub.status.idle":"2022-09-23T16:23:39.136635Z","shell.execute_reply.started":"2022-09-23T16:23:38.913108Z","shell.execute_reply":"2022-09-23T16:23:39.135690Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize = (20,10) )\nc = 0\nfor f in ['day', 'donor', 'cell_type']:# , 'technology']:\n    c+=1; fig.add_subplot(1, 3 , c) \n    sns.scatterplot(x=r[:,0], y=r[:,1] , hue = d[f]  )\n    plt.title('Colored by '+f)\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2022-08-21T21:20:42.560726Z","iopub.execute_input":"2022-08-21T21:20:42.561114Z","iopub.status.idle":"2022-08-21T21:20:51.458336Z","shell.execute_reply.started":"2022-08-21T21:20:42.561081Z","shell.execute_reply":"2022-08-21T21:20:51.457060Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"    MasP = Mast Cell Progenitor\n    MkP = Megakaryocyte Progenitor\n    NeuP = Neutrophil Progenitor\n    MoP = Monocyte Progenitor\n    EryP = Erythrocyte Progenitor\n    HSC = Hematoploetic Stem Cell\n    BP = B-Cell Progenitor","metadata":{}},{"cell_type":"code","source":"l = ['CD88', 'CD32', 'CD36', 'CD71', 'CD41', 'CD115', 'FceRIa', 'CD26',\n       'CD272', 'CD82'] + [     'CD45', 'HLA-A-B-C',     'CD244',     'CD123',      'CD13',\n            'CD44',      'CD31',     'CD49b',     'CD62L',           ]\nl1 = ['CD88', 'CD32', 'CD36', 'CD71', 'CD41', 'CD115', 'FceRIa', 'CD26', 'CD272', 'CD82']\nl2 = [     'CD45', 'HLA-A-B-C',     'CD244',     'CD123',      'CD13', 'CD44',      'CD31',     'CD49b',     'CD62L' ]\n\nlc1 = ['CD71', 'CD115', 'CD88']\nlc2 =  ['CD155', 'CD112', 'CD47', 'HLA-A-B-C' ,'CD45RA' ,'CD31' ,'CD11a', 'CD13', 'CD29',\n  'CD81', 'CD18', 'CD45', 'CD49d']\nprint(set(l1) & set(lc1) )\nprint(set(l1) & set(lc2) )\n\nprint(set(l2) & set(lc1) )\nprint(set(l2) & set(lc2) )\n\n","metadata":{"execution":{"iopub.status.busy":"2022-09-23T16:27:46.402584Z","iopub.execute_input":"2022-09-23T16:27:46.403053Z","iopub.status.idle":"2022-09-23T16:27:46.414461Z","shell.execute_reply.started":"2022-09-23T16:27:46.403018Z","shell.execute_reply":"2022-09-23T16:27:46.413209Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Top correlated with ncivis2 - which seems to reflect the cell type\nl1 = ['CD88', 'CD32', 'CD36', 'CD71', 'CD41', 'CD115', 'FceRIa', 'CD26', 'CD272', 'CD82']\nl2 = [     'CD45', 'HLA-A-B-C',     'CD244',     'CD123',      'CD13', 'CD44',      'CD31',     'CD49b',     'CD62L' ]\n\nimport numbers\n\nn_x_subplots = 5\npalette = 'rainbow'#'viridis'\n\nt0 = time.time()\n\nc = 0\nfor f in ['cell_type', 'day', 'donor', ] + l1 +l2: # ['day', 'donor', 'cell_type']:# , 'technology']:\n    if c % n_x_subplots == 0:\n        if c > 0:\n            plt.show()\n        fig = plt.figure(figsize = (20,7) ); c = 0\n        #plt.suptitle(str(k),fontsize = 20 )# str_data_inf + ' n_cells: ' +str(mask.sum()) + ' RED > median expression, BLUE <= median ' )# +' ' + cell_type +' ' + drug )\n        \n    c += 1; fig.add_subplot(1,n_x_subplots ,c)    \n    \n    v = d[f]\n    if isinstance(v[0], numbers.Number):\n        v = np.clip(v,np.percentile(v,5), np.percentile(v,95) )\n        # v = (v > np.median(v)).astype(int)\n    ax = sns.scatterplot(x=r[:,0], y=r[:,1] , hue =  v , palette = palette )\n    plt.setp(ax.get_legend().get_texts(), fontsize=12) # for legend text\n    plt.setp(ax.get_legend().get_title(), fontsize=12) # for legend title    \n    plt.title('Colored by '+f, fontsize = 10)\n    \nplt.show()\nprint(time.time()-t0, 'secs passed')\n","metadata":{"execution":{"iopub.status.busy":"2022-09-23T16:33:59.667158Z","iopub.execute_input":"2022-09-23T16:33:59.667882Z","iopub.status.idle":"2022-09-23T16:35:35.158406Z","shell.execute_reply.started":"2022-09-23T16:33:59.667843Z","shell.execute_reply":"2022-09-23T16:35:35.156942Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 0.7 correlation_threshold \n# Top 5 cluster sizes: [3 2 1 1 1] seconds passed: 0.0\n# 50 Genes in largest correlated group:\n# ['CD71' 'CD115' 'CD88']\n\n# 0.6 correlation_threshold \n# Top 5 cluster sizes: [14  5  2  1  1] seconds passed: 0.0\n# 50 Genes in largest correlated group:\n# ['CD155' 'CD112' 'CD47' 'HLA-A-B-C' 'CD45RA' 'CD31' 'CD11a' 'CD13' 'CD29'\n#  'CD81' 'CD18' 'CD45' 'CD49d' 'CD162']\n\nimport numbers\n\nn_x_subplots = 8\npalette = 'rainbow'#'viridis'\n\nt0 = time.time()\n\nc = 0\nfor f in ['CD71', 'CD115' ,'CD88'] + ['CD155', 'CD112', 'CD47', 'HLA-A-B-C' ,'CD45RA' ,'CD31' ,'CD11a', 'CD13', 'CD29',\n  'CD81', 'CD18', 'CD45', 'CD49d']:# , 'CD162']: # ['day', 'donor', 'cell_type']:# , 'technology']:\n    if c % n_x_subplots == 0:\n        if c > 0:\n            plt.show()\n        fig = plt.figure(figsize = (20,5) ); c = 0\n        #plt.suptitle(str(k),fontsize = 20 )# str_data_inf + ' n_cells: ' +str(mask.sum()) + ' RED > median expression, BLUE <= median ' )# +' ' + cell_type +' ' + drug )\n        \n    c += 1; fig.add_subplot(1,n_x_subplots ,c)    \n    \n    v = d[f]\n    if isinstance(v[0], numbers.Number):\n        v = np.clip(v,np.percentile(v,5), np.percentile(v,95) )\n        # v = (v > np.median(v)).astype(int)\n    ax = sns.scatterplot(x=r[:,0], y=r[:,1] , hue =  v , palette = palette )\n    plt.setp(ax.get_legend().get_texts(), fontsize=10) # for legend text\n    plt.setp(ax.get_legend().get_title(), fontsize=10) # for legend title    \n    plt.title('Colored by '+f, fontsize = 10)\n    \nplt.show()\nprint(time.time()-t0, 'secs passed')\n","metadata":{"execution":{"iopub.status.busy":"2022-08-21T21:20:51.459976Z","iopub.execute_input":"2022-08-21T21:20:51.461073Z","iopub.status.idle":"2022-08-21T21:22:03.703949Z","shell.execute_reply.started":"2022-08-21T21:20:51.461032Z","shell.execute_reply":"2022-08-21T21:22:03.702928Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n# 0.7 correlation_threshold \n# Top 5 cluster sizes: [3 2 1 1 1] seconds passed: 0.0\n# 50 Genes in largest correlated group:\n# ['CD71' 'CD115' 'CD88']\n\n# 0.6 correlation_threshold \n# Top 5 cluster sizes: [14  5  2  1  1] seconds passed: 0.0\n# 50 Genes in largest correlated group:\n# ['CD155' 'CD112' 'CD47' 'HLA-A-B-C' 'CD45RA' 'CD31' 'CD11a' 'CD13' 'CD29'\n#  'CD81' 'CD18' 'CD45' 'CD49d' 'CD162']\n\n# CD69 экспрессируется в Т клетках после активации. После этого они пролиферируют. Так же могут быть интересны CD25 и CD71\n\n# https://link.springer.com/article/10.1007/BF01305907\n    \nimport numbers\n\nn_x_subplots = 4\npalette = 'rainbow'#'viridis'\n\nt0 = time.time()\n\nc = 0\nfor f in ['CD69', 'CD25' ,'CD71','cell_type']:# , 'CD162']: # ['day', 'donor', 'cell_type']:# , 'technology']:\n    if c % n_x_subplots == 0:\n        if c > 0:\n            plt.show()\n        fig = plt.figure(figsize = (20,5) ); c = 0\n        #plt.suptitle(str(k),fontsize = 20 )# str_data_inf + ' n_cells: ' +str(mask.sum()) + ' RED > median expression, BLUE <= median ' )# +' ' + cell_type +' ' + drug )\n        \n    c += 1; fig.add_subplot(1,n_x_subplots ,c)    \n    \n    v = d[f]\n    if isinstance(v[0], numbers.Number):\n        v = np.clip(v,np.percentile(v,5), np.percentile(v,95) )\n        # v = (v > np.median(v)).astype(int)\n    ax = sns.scatterplot(x=r[:,0], y=r[:,1] , hue =  v , palette = palette )\n    plt.setp(ax.get_legend().get_texts(), fontsize=10) # for legend text\n    plt.setp(ax.get_legend().get_title(), fontsize=10) # for legend title    \n    plt.title('Colored by '+f, fontsize = 10)\n    \nplt.show()\nprint(time.time()-t0, 'secs passed')\n","metadata":{"execution":{"iopub.status.busy":"2022-08-21T21:30:54.754724Z","iopub.execute_input":"2022-08-21T21:30:54.755328Z","iopub.status.idle":"2022-08-21T21:31:11.192082Z","shell.execute_reply.started":"2022-08-21T21:30:54.755283Z","shell.execute_reply":"2022-08-21T21:31:11.190786Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"    MasP = Mast Cell Progenitor\n    MkP = Megakaryocyte Progenitor\n    NeuP = Neutrophil Progenitor\n    MoP = Monocyte Progenitor\n    EryP = Erythrocyte Progenitor\n    HSC = Hematoploetic Stem Cell\n    BP = B-Cell Progenitor","metadata":{}},{"cell_type":"code","source":"\n# 0.7 correlation_threshold \n# Top 5 cluster sizes: [3 2 1 1 1] seconds passed: 0.0\n# 50 Genes in largest correlated group:\n# ['CD71' 'CD115' 'CD88']\n\n# 0.6 correlation_threshold \n# Top 5 cluster sizes: [14  5  2  1  1] seconds passed: 0.0\n# 50 Genes in largest correlated group:\n# ['CD155' 'CD112' 'CD47' 'HLA-A-B-C' 'CD45RA' 'CD31' 'CD11a' 'CD13' 'CD29'\n#  'CD81' 'CD18' 'CD45' 'CD49d' 'CD162']\n\n# CD69 экспрессируется в Т клетках после активации. После этого они пролиферируют. Так же могут быть интересны CD25 и CD71\n\n# https://link.springer.com/article/10.1007/BF01305907\n    \nimport numbers\n\nn_x_subplots = 5\npalette = 'rainbow'#'viridis'\n\nt0 = time.time()\n\nc = 0\nfor f in ['CD3', 'CD4', 'CD8', 'CD25','cell_type']: #  ['CD69', 'CD25' ,'CD71','cell_type']:# , 'CD162']: # ['day', 'donor', 'cell_type']:# , 'technology']:\n    if f not in d.columns: continue \n    if c % n_x_subplots == 0:\n        if c > 0:\n            plt.show()\n        fig = plt.figure(figsize = (20,5) ); c = 0\n        #plt.suptitle(str(k),fontsize = 20 )# str_data_inf + ' n_cells: ' +str(mask.sum()) + ' RED > median expression, BLUE <= median ' )# +' ' + cell_type +' ' + drug )\n        \n    c += 1; fig.add_subplot(1,n_x_subplots ,c)    \n    \n    v = d[f]\n    if isinstance(v[0], numbers.Number):\n        v = np.clip(v,np.percentile(v,5), np.percentile(v,95) )\n        # v = (v > np.median(v)).astype(int)\n    ax = sns.scatterplot(x=r[:,0], y=r[:,1] , hue =  v , palette = palette )\n    plt.setp(ax.get_legend().get_texts(), fontsize=10) # for legend text\n    plt.setp(ax.get_legend().get_title(), fontsize=10) # for legend title    \n    plt.title('Colored by '+f, fontsize = 10)\n    \nplt.show()\nprint(time.time()-t0, 'secs passed')\n","metadata":{"execution":{"iopub.status.busy":"2022-08-21T21:38:27.080792Z","iopub.execute_input":"2022-08-21T21:38:27.081291Z","iopub.status.idle":"2022-08-21T21:38:46.996404Z","shell.execute_reply.started":"2022-08-21T21:38:27.081252Z","shell.execute_reply":"2022-08-21T21:38:46.995434Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Top correlated with ncivis2 - which seems to reflect the cell type\nl1 = ['CD88', 'CD32', 'CD36', 'CD71', 'CD41', 'CD115', 'FceRIa', 'CD26', 'CD272', 'CD82']\nl2 = [     'CD45', 'HLA-A-B-C',     'CD244',     'CD123',      'CD13', 'CD44',      'CD31',     'CD49b',     'CD62L' ]\n\n\nreducer =  ncvis.NCVis() # trimap.TRIMAP() #  umap.UMAP()\n\nt0 = time.time()\nr2 = reducer.fit_transform(df_cite_train_y[l1+l2])\nprint(time.time()-t0, 'secs passed for dimension reduction')\n\n\nfig = plt.figure(figsize = (15,7) )\n# c = 0\n# for f in ['cp_dose', 'cp_type','cp_time','y_sum']:\n#     c+=1; fig.add_subplot(1, 4 , c) \nsns.scatterplot(x=r2[:,0], y=r2[:,1]  , hue = d[column4color]  )\nplt.title('NCCIS on selected genes potentially cell type related ', fontsize = 20 )\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-09-23T16:57:20.828241Z","iopub.execute_input":"2022-09-23T16:57:20.828760Z","iopub.status.idle":"2022-09-23T16:57:51.906919Z","shell.execute_reply.started":"2022-09-23T16:57:20.828721Z","shell.execute_reply":"2022-09-23T16:57:51.905747Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\n# Top correlated with ncivis2 - which seems to reflect the cell type\nl1 = ['CD88', 'CD32', 'CD36', 'CD71', 'CD41', 'CD115', 'FceRIa', 'CD26', 'CD272', 'CD82']\nl2 = [     'CD45', 'HLA-A-B-C',     'CD244',     'CD123',      'CD13', 'CD44',      'CD31',     'CD49b',     'CD62L' ]\n\nX2 = df_cite_train_y[l1+l2].copy()\n\nv1 = X2.quantile(0.95, axis=0)\nv2 = X2.quantile(0.05, axis=0)\n\n# cut off outliers:\nfor col in X2.columns:\n    X2[col] = np.clip(X2[col], a_min = v2[col], a_max = v1[col])\n\n# Take top  most variable \nv = X2.std(axis = 0) / X2.abs().mean(axis=0) \nIX = v.sort_values(ascending = False).iloc[:100]\nX2 = X2[IX.index]\n\n\nfrom sklearn.preprocessing import StandardScaler\nscaler = StandardScaler()\nX2 = scaler.fit_transform(X2)\n\n\nreducer =  ncvis.NCVis() # trimap.TRIMAP() #  umap.UMAP()\nt0 = time.time()\nr2 = reducer.fit_transform(X2)\nprint(time.time()-t0, 'secs passed for dimension reduction')\nprint(X2.shape)\n\nfig = plt.figure(figsize = (15,7) )\n# c = 0\n# for f in ['cp_dose', 'cp_type','cp_time','y_sum']:\n#     c+=1; fig.add_subplot(1, 4 , c) \nsns.scatterplot(x=r2[:,0], y=r2[:,1]  , hue = d[column4color]  )\nplt.title('NCCIS on preprocessed selected genes potentially cell type related ', fontsize = 20 )\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-09-23T17:12:27.916120Z","iopub.execute_input":"2022-09-23T17:12:27.916648Z","iopub.status.idle":"2022-09-23T17:12:59.028226Z","shell.execute_reply.started":"2022-09-23T17:12:27.916610Z","shell.execute_reply":"2022-09-23T17:12:59.026573Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Top correlated with ncivis2 - which seems to reflect the cell type\nl1 = ['CD88', 'CD32', 'CD36', 'CD71', 'CD41', 'CD115', 'FceRIa', 'CD26', 'CD272', 'CD82']\nl2 = [     'CD45', 'HLA-A-B-C',     'CD244',     'CD123',      'CD13', 'CD44',      'CD31',     'CD49b',     'CD62L' ]\n\nimport numbers\n\nn_x_subplots = 5\npalette = 'rainbow'#'viridis'\n\nt0 = time.time()\n\nc = 0\nfor f in ['cell_type', 'day', 'donor', ] + l1 +l2: # ['day', 'donor', 'cell_type']:# , 'technology']:\n    if c % n_x_subplots == 0:\n        if c > 0:\n            plt.show()\n        fig = plt.figure(figsize = (20,7) ); c = 0\n        #plt.suptitle(str(k),fontsize = 20 )# str_data_inf + ' n_cells: ' +str(mask.sum()) + ' RED > median expression, BLUE <= median ' )# +' ' + cell_type +' ' + drug )\n        \n    c += 1; fig.add_subplot(1,n_x_subplots ,c)    \n    \n    v = d[f]\n    if isinstance(v[0], numbers.Number):\n        v = np.clip(v,np.percentile(v,5), np.percentile(v,95) )\n        # v = (v > np.median(v)).astype(int)\n    ax = sns.scatterplot(x=r2[:,0], y=r2[:,1] , hue =  v , palette = palette )\n    plt.setp(ax.get_legend().get_texts(), fontsize=12) # for legend text\n    plt.setp(ax.get_legend().get_title(), fontsize=12) # for legend title    \n    plt.title('Colored by '+f, fontsize = 10)\n    \nplt.show()\nprint(time.time()-t0, 'secs passed')\n","metadata":{"execution":{"iopub.status.busy":"2022-09-23T17:17:42.902680Z","iopub.execute_input":"2022-09-23T17:17:42.903877Z","iopub.status.idle":"2022-09-23T17:19:13.196936Z","shell.execute_reply.started":"2022-09-23T17:17:42.903830Z","shell.execute_reply":"2022-09-23T17:19:13.195816Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Day/Donor dependence - that is very puzzling ","metadata":{"execution":{"iopub.status.busy":"2022-09-23T13:55:06.356764Z","iopub.execute_input":"2022-09-23T13:55:06.357668Z","iopub.status.idle":"2022-09-23T13:55:06.385992Z","shell.execute_reply.started":"2022-09-23T13:55:06.357615Z","shell.execute_reply":"2022-09-23T13:55:06.384286Z"}}},{"cell_type":"code","source":"list_data_cols = df_cite_train_y.columns\n\nst = d.groupby('day').mean()\n#display(d.head(1))\nst = st[list_data_cols]#.iloc[:,1:]\nst.head(1)\nst.T.plot(figsize = (20,3))\nplt.show()\n\nv = st.mean(axis = 0)\nw = st/v\nw\n#plt.figure(figsize = (20,3))\nw.T.plot(figsize = (20,3))\nplt.show()\nt = w.mean(axis = 1)\nt.iat[1]/t.iat[0], t.iat[1]/t.iat[2], t.iat[2]/t.iat[0], ","metadata":{"execution":{"iopub.status.busy":"2022-09-23T14:41:45.883352Z","iopub.execute_input":"2022-09-23T14:41:45.883829Z","iopub.status.idle":"2022-09-23T14:41:46.522654Z","shell.execute_reply.started":"2022-09-23T14:41:45.883793Z","shell.execute_reply":"2022-09-23T14:41:46.521477Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"list_data_cols = df_cite_train_y.columns\nst = d.groupby('donor').mean()\n#display(d.head(1))\nst = st[list_data_cols]# .iloc[:,1:]\n#display(st.head(1))\n\nv = st.mean(axis = 0)\nw = st/v\n# display(w)\n#plt.figure(figsize = (20,3))\nw.T.plot(figsize = (20,3))\nplt.show()\nt = w.mean(axis = 1)\nt.iat[1]/t.iat[0], t.iat[1]/t.iat[2], t.iat[2]/t.iat[0], ","metadata":{"execution":{"iopub.status.busy":"2022-09-23T14:44:15.447097Z","iopub.execute_input":"2022-09-23T14:44:15.448722Z","iopub.status.idle":"2022-09-23T14:44:15.833155Z","shell.execute_reply.started":"2022-09-23T14:44:15.448657Z","shell.execute_reply":"2022-09-23T14:44:15.831972Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"list_data_cols = df_cite_train_y.columns\nst = d.groupby(['donor', 'day',]).mean()\n#display(d.head(1))\nst = st[list_data_cols]# .iloc[:,1:]\n#display(st.head(1))\n\nv = st.mean(axis = 0)\nw = st/v\n# display(w)\n#plt.figure(figsize = (20,3))\nw.T.iloc[:,:3].plot(figsize = (20,5))\nplt.show()\nw.T.iloc[:,3:6].plot(figsize = (20,5))\nplt.show()\nw.T.iloc[:,6:9].plot(figsize = (20,5))\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-09-23T14:29:33.722054Z","iopub.execute_input":"2022-09-23T14:29:33.723119Z","iopub.status.idle":"2022-09-23T14:29:34.626662Z","shell.execute_reply.started":"2022-09-23T14:29:33.723053Z","shell.execute_reply":"2022-09-23T14:29:34.625320Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df2 = df_cite_train_y.copy()\nfor c in df2.columns:\n    v = df2[c]\n    a,b = v.quantile(0.05),v.quantile(0.95),\n    df2[c] = np.clip(v, a, b)\n    \nv = df2.std()\nv.sort_values().plot(  figsize = (20,4) ) \nplt.show()\ndisplay( v.sort_values().tail(10) )\nv = df2.abs().mean()\nv.sort_values().plot(  figsize = (20,4) ) \n\ndisplay( v.sort_values().tail(10) )\n","metadata":{"execution":{"iopub.status.busy":"2022-09-23T15:15:35.851301Z","iopub.execute_input":"2022-09-23T15:15:35.851767Z","iopub.status.idle":"2022-09-23T15:15:37.236226Z","shell.execute_reply.started":"2022-09-23T15:15:35.851734Z","shell.execute_reply":"2022-09-23T15:15:37.235042Z"},"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":"","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) )","metadata":{"execution":{"iopub.status.busy":"2022-09-23T16:14:09.391714Z","iopub.execute_input":"2022-09-23T16:14:09.392511Z","iopub.status.idle":"2022-09-23T16:14:09.417234Z","shell.execute_reply.started":"2022-09-23T16:14:09.392469Z","shell.execute_reply":"2022-09-23T16:14:09.415440Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}