{"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 ATAC-seq data. Mainly we do dimensional reduction and visualizations.\n\nMain surprise - cell types can be clearly seen via ATAC-seq data, seems even better than with scRNA-seq data.\nDespite typically cell types are typically studied via scRNA-seq ! \n\n### Conclusions:\n    \n#### 1) Cell types - can be seen; days somewhat less clear.  But donors - seems NO - i.e. not \"batch\" effect from donors.\n    The cell types can be seen via ATAC-seq and seems to provide the main variation. Some variation due to  day can also be seen. Seems to match effect from donors.  \n    \n    \n**Remark1. Cell types** The plots in PCA3,4, 5 is SURPRISINGLY clearly \"separate\" different cell types. That seems works even better than for scRNA-seq data. Biostars https://www.biostars.org/p/9538656/ commented that it might be due to gene have many regulatary regions and thus (due to that \"many\") we see clearly pictures. \nErythroid subtype seems to have alongated structure going away from the main part (PCA2, PCA4) - that probably correspond to RNA-seq visualizatin seen in S.Nestorova et.al. paper/dataset  https://www.kaggle.com/datasets/alexandervc/single-cell-rna-seq-nestorova2016-mouse-hspc\n\n\n**Remark 2. Days** days can be seen by straightforward \"imshow\"-like visualization for sparse matrices:\nhttps://www.kaggle.com/code/alexandervc/heatmap-of-atac-seq-multiome-data (based on NORIFUMI IRIE https://www.kaggle.com/code/norifumiirie/heatmap-of-multiome-data )\n\n**Remark 3. Donors** - \"absence\" of \"batch\" effect from donors has been reported before:\nhttps://www.kaggle.com/competitions/open-problems-multimodal/discussion/350953 . But no code was provided.\nI.e. data are quite similar for different donors. \n\n#### 2) PCA1 - correspond to sum of all values in each row. \n    We can also consider sums over each chromosome - but they seems to be quite correlated to simple sum.\n    \n    Is there any biological sense in sum ? \n    \n\n#### 3) TruncatedSVD is fast enough to make experiments\n    We use TruncatedSVD which can work for sparse matrices. \n    It is reasonably fast, especially when considired NOT for ALL samples (e.g. fix donor, or day etc).  \n    E.g. 7264 cells x 228942 features -> 10 components - 16 seconds. \n    33701 cells x 228942 features -> 10 components - 80 seconds. ; 23911  x 228942 features -> 10 components - 56 seconds. \n    PS\n    TruncatedSVD with FULL data is quite long: # 105942x x 228942 features -> 32 components - 9min 23s \n    \n#### 4) Beautiful \"batterfly\" plots \n    Coloring by sum - creates nice vice visualizations. See e.g. PCA2-PCA4. \n    \n#### Questions.\n    PCA1 for ATAC seq data should probably have some known interpretation - what is it ?  probably quality of reads or something else,     probably - not biological variation \n    \n#### Technicalities: \nATAC-seq data is huge - (105942, 228942) , but sparse - so we use sparse fmt created by SBUNZINI https://www.kaggle.com/datasets/sbunzini/open-problems-msci-multiome-sparse-matrices\n\n\nMetadata. Dataframe with metadata for mutliome train part is prepared. (Using h5py loading packages from .h5 files). \n","metadata":{}},{"cell_type":"markdown","source":"# Install/Import","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 numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\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":"2022-09-16T15:30:41.012571Z","iopub.execute_input":"2022-09-16T15:30:41.013066Z","iopub.status.idle":"2022-09-16T15:30:41.052822Z","shell.execute_reply.started":"2022-09-16T15:30:41.012968Z","shell.execute_reply":"2022-09-16T15:30:41.051927Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install scanpy\nimport scanpy as sc\nimport anndata\n\nimport time\nt0start = time.time()\n\nimport pandas as pd\nimport numpy as np\nimport os\nimport sys\n\nimport matplotlib.pyplot as plt\n#plt.style.use('dark_background')\nimport seaborn as sns\n\n#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-16T15:30:41.054406Z","iopub.execute_input":"2022-09-16T15:30:41.054934Z","iopub.status.idle":"2022-09-16T15:31:11.756406Z","shell.execute_reply.started":"2022-09-16T15:30:41.054902Z","shell.execute_reply":"2022-09-16T15:31:11.755174Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import h5py\n!pip install hdf5plugin~=2.0 # https://forum.hdfgroup.org/t/cant-open-directory-usr-local-hdf5-lib-plugin/9738/4\nimport hdf5plugin\n\nprint()\nprint('Look on TARGETs: ')\nfilename = '/kaggle/input/open-problems-multimodal/train_multi_targets.h5'\nf = h5py.File(filename,'r')#, mode)\nprint('fragment of a DNA Ids:')\nprint( f['train_multi_targets']['axis0'][:5] )\nprint();  print('Cell Ids (seems to me):')\nprint( f['train_multi_targets']['axis1'].shape)\nprint(); print('First:')\nprint( f['train_multi_targets']['axis1'][:15] )\nprint(); print('Last:')\nprint( f['train_multi_targets']['axis1'][-15:] )\n","metadata":{"execution":{"iopub.status.busy":"2022-09-16T15:31:11.758389Z","iopub.execute_input":"2022-09-16T15:31:11.758838Z","iopub.status.idle":"2022-09-16T15:31:23.600650Z","shell.execute_reply.started":"2022-09-16T15:31:11.758797Z","shell.execute_reply":"2022-09-16T15:31:23.599018Z"},"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":{"execution":{"iopub.status.busy":"2022-09-16T15:31:23.604526Z","iopub.execute_input":"2022-09-16T15:31:23.605069Z","iopub.status.idle":"2022-09-16T15:31:24.066216Z","shell.execute_reply.started":"2022-09-16T15:31:23.605020Z","shell.execute_reply":"2022-09-16T15:31:24.064896Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Create metadata dataframe - cell types, donors, etc","metadata":{}},{"cell_type":"code","source":"%%time\nbarcodes = [t.decode() for t in f['train_multi_targets']['axis1'] ]\nprint(barcodes[:15])\nprint(barcodes[-15:])\n\ndf_meta_multiome_train = pd.DataFrame(index = barcodes)\ndf_meta_multiome_train['Index in Multiome'] = range(len(df_meta_multiome_train))\nprint(df_meta_multiome_train.shape)\n#obs = pd.DataFrame(index = barcodes)\nd = df_cell.reset_index().set_index('cell_id') \nd.columns = ['Index in df_cell'] + list(d.columns[1:])\ndf_meta_multiome_train = df_meta_multiome_train.join(d, how = 'left', )\ndf_meta_multiome_train['str_donor'] = df_meta_multiome_train['donor'].apply(lambda x: str(x))  # will be useful for visualizations\nprint(df_meta_multiome_train.shape)\ndisplay(df_meta_multiome_train)","metadata":{"execution":{"iopub.status.busy":"2022-09-16T15:31:24.068000Z","iopub.execute_input":"2022-09-16T15:31:24.068358Z","iopub.status.idle":"2022-09-16T15:31:33.834629Z","shell.execute_reply.started":"2022-09-16T15:31:24.068326Z","shell.execute_reply":"2022-09-16T15:31:33.833426Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nprint('Sanity check - we see same indexes when load in a way as orgs proposed:')\nSTART = 0\nSTOP = 5\ndf_multi_train_y = pd.read_hdf(FP_MULTIOME_TRAIN_TARGETS, start=START, stop=STOP)\n#display(df_multi_train_y.head())\nprint(df_multi_train_y.index)\nSTART = 105942-5 # Five last \nSTOP = 105942\ndf_multi_train_y = pd.read_hdf(FP_MULTIOME_TRAIN_TARGETS, start=START, stop=STOP)\n#display(df_multi_train_y.head())\nprint(df_multi_train_y.index)\n","metadata":{"execution":{"iopub.status.busy":"2022-09-16T15:31:33.836479Z","iopub.execute_input":"2022-09-16T15:31:33.837245Z","iopub.status.idle":"2022-09-16T15:31:34.076465Z","shell.execute_reply.started":"2022-09-16T15:31:33.837200Z","shell.execute_reply":"2022-09-16T15:31:34.075195Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load ATAC-seq data in sparse format (created by  SBUNZINI )","metadata":{}},{"cell_type":"code","source":"%%time\nimport scipy.sparse \nfn = '/kaggle/input/open-problems-msci-multiome-sparse-matrices/train_multiome_input_sparse.npz'\ntrain_inputs = scipy.sparse.load_npz(fn)# \"../input/multimodal-single-cell-as-sparse-matrix/train_multi_inputs_values.sparse.npz\")","metadata":{"execution":{"iopub.status.busy":"2022-09-16T15:31:34.077925Z","iopub.execute_input":"2022-09-16T15:31:34.078291Z","iopub.status.idle":"2022-09-16T15:32:37.659691Z","shell.execute_reply.started":"2022-09-16T15:31:34.078259Z","shell.execute_reply":"2022-09-16T15:32:37.658440Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_inputs","metadata":{"execution":{"iopub.status.busy":"2022-09-16T15:32:37.661275Z","iopub.execute_input":"2022-09-16T15:32:37.661726Z","iopub.status.idle":"2022-09-16T15:32:37.671173Z","shell.execute_reply.started":"2022-09-16T15:32:37.661680Z","shell.execute_reply":"2022-09-16T15:32:37.669847Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('8.7G consumed')","metadata":{"execution":{"iopub.status.busy":"2022-09-16T15:32:37.675977Z","iopub.execute_input":"2022-09-16T15:32:37.676763Z","iopub.status.idle":"2022-09-16T15:32:37.683170Z","shell.execute_reply.started":"2022-09-16T15:32:37.676724Z","shell.execute_reply":"2022-09-16T15:32:37.682024Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Calculate additonal features - just sum in row, and sums over each chromosome. (They are quite correlated and correlated with total PCA1.  ","metadata":{}},{"cell_type":"code","source":"%%time\nv = train_inputs.sum(axis = 1)\ndf_meta_multiome_train['Sum'] = np.asarray(v).ravel()\ndf_meta_multiome_train['Sum'].plot()","metadata":{"execution":{"iopub.status.busy":"2022-09-16T15:32:37.684740Z","iopub.execute_input":"2022-09-16T15:32:37.685984Z","iopub.status.idle":"2022-09-16T15:32:38.820784Z","shell.execute_reply.started":"2022-09-16T15:32:37.685911Z","shell.execute_reply":"2022-09-16T15:32:38.819342Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nSTART = int(0)\nSTOP = START+1# 4_000 # Bigger chunk will crash 16G memory at plt.imshow()\n\ndf_multi_train_x = pd.read_hdf(FP_MULTIOME_TRAIN_INPUTS,start=START,stop=STOP)\ndisplay( df_multi_train_x.head() )\n\nlist_col_names = list( df_multi_train_x.columns)\nprint(list_col_names[:5], list_col_names[-5:])","metadata":{"execution":{"iopub.status.busy":"2022-09-16T15:32:38.822458Z","iopub.execute_input":"2022-09-16T15:32:38.822839Z","iopub.status.idle":"2022-09-16T15:32:39.592047Z","shell.execute_reply.started":"2022-09-16T15:32:38.822803Z","shell.execute_reply":"2022-09-16T15:32:39.590704Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nlist_chromosomes_ids = ['chr'+str(i) for i in range(1,23)] + ['chrY' , 'chrX', 'GL0',  'KI2']\nfor c in list_chromosomes_ids:\n    l = [ c in x for x in list_col_names ]\n    I = np.where(l)[0] \n    print(c,np.sum(l))\n    v = train_inputs[:,I].sum(axis = 1)\n    df_meta_multiome_train['Sum '+c] = np.asarray(v).ravel()\n    \nfor c in list_chromosomes_ids:\n    df_meta_multiome_train['Sum '+c].plot()\n\n","metadata":{"execution":{"iopub.status.busy":"2022-09-16T15:32:39.593858Z","iopub.execute_input":"2022-09-16T15:32:39.594217Z","iopub.status.idle":"2022-09-16T15:33:36.574773Z","shell.execute_reply.started":"2022-09-16T15:32:39.594185Z","shell.execute_reply":"2022-09-16T15:33:36.572911Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ncorr_m = df_meta_multiome_train.corr()\ndisplay(corr_m)","metadata":{"execution":{"iopub.status.busy":"2022-09-16T15:33:36.577298Z","iopub.execute_input":"2022-09-16T15:33:36.577737Z","iopub.status.idle":"2022-09-16T15:33:36.959932Z","shell.execute_reply.started":"2022-09-16T15:33:36.577691Z","shell.execute_reply":"2022-09-16T15:33:36.958647Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize = (20,12))\nsns.heatmap(corr_m);","metadata":{"execution":{"iopub.status.busy":"2022-09-16T15:33:36.961448Z","iopub.execute_input":"2022-09-16T15:33:36.961910Z","iopub.status.idle":"2022-09-16T15:33:37.877233Z","shell.execute_reply.started":"2022-09-16T15:33:36.961876Z","shell.execute_reply":"2022-09-16T15:33:37.876067Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Dimensional reduction with TruncatedSVD\n\nTruncatedSVD - similar to PCA, but can work with sparse matrices. \n","metadata":{}},{"cell_type":"code","source":"from sklearn.decomposition import TruncatedSVD","metadata":{"execution":{"iopub.status.busy":"2022-09-16T15:33:37.879509Z","iopub.execute_input":"2022-09-16T15:33:37.879996Z","iopub.status.idle":"2022-09-16T15:33:38.086293Z","shell.execute_reply.started":"2022-09-16T15:33:37.879933Z","shell.execute_reply":"2022-09-16T15:33:38.084946Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Vusuzalizations to show: PCA 1 seems to be quite related to total sum . And other beautilful pictures. \n\nThat sometimes happens. ","metadata":{}},{"cell_type":"code","source":"%%time\nfor donor in [32606]:# df_meta_multiome_train['donor'].unique():\n    #donor = 32606\n    #day = 2\n    str_inf = ''\n    str_inf = ' donor= ' + str(donor) \n    #str_inf += ' day= ' + str(day)\n    #mask = (df_meta_multiome_train['day'] == day) # & (df_meta_multiome_train['donor'] == donor)\n    mask = (df_meta_multiome_train['donor'] == donor)\n\n    print(mask.sum(), str_inf )\n    tsvd = TruncatedSVD(n_components = 10,random_state = 4223)\n    t0 = time.time()\n    r = tsvd.fit_transform(train_inputs[mask,:])\n    print('%.1f'%(time.time()-t0) + ' seconds passed. TruncatedSVD finished.')\n\n    \n    for (i,j) in [ (0,1),(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4)]:\n        plt.figure(figsize = (20,10))\n        ax = sns.scatterplot(r[:,i],r[:,j], hue = df_meta_multiome_train['Sum'][mask], palette = 'rainbow', alpha = 0.9)\n        if 1:\n            plt.setp(ax.get_legend().get_texts(), fontsize=20) # for legend text\n            plt.setp(ax.get_legend().get_title(), fontsize=20) # for legend title        \n        plt.title('ATAC-seq. Train. PCA (TruncSVD). ' + str_inf  +   ' n_cells='+str(len(r)), fontsize = 20)\n        plt.xlabel('PCA'+str(i+1), fontsize = 20)\n        plt.ylabel('PCA'+str(j+1), fontsize = 20)\n        plt.show()    \n","metadata":{"execution":{"iopub.status.busy":"2022-09-16T13:25:28.394275Z","iopub.execute_input":"2022-09-16T13:25:28.395697Z","iopub.status.idle":"2022-09-16T13:27:04.135774Z","shell.execute_reply.started":"2022-09-16T13:25:28.395653Z","shell.execute_reply":"2022-09-16T13:27:04.134386Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Suprisingly clearly seen cell types in PCA3,PCA4, PCA5. For all donors. (All days included). ","metadata":{}},{"cell_type":"code","source":"%%time\nfor donor in df_meta_multiome_train['donor'].unique():\n    #donor = 32606\n    #day = 2\n    str_inf = ''\n    str_inf = ' donor= ' + str(donor) \n    #str_inf += ' day= ' + str(day)\n    #mask = (df_meta_multiome_train['day'] == day) # & (df_meta_multiome_train['donor'] == donor)\n    mask = (df_meta_multiome_train['donor'] == donor)\n\n    print(mask.sum(), str_inf )\n    tsvd = TruncatedSVD(n_components = 10,random_state = 4223)\n    t0 = time.time()\n    r = tsvd.fit_transform(train_inputs[mask,:])\n    print('%.1f'%(time.time()-t0) + ' seconds passed. TruncatedSVD finished.')\n\n    \n    for (i,j) in [(2,3),(2,4), (3,4), ]:# (0,1),(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4)]:\n        plt.figure(figsize = (20,10))\n        ax = sns.scatterplot(r[:,i],r[:,j], hue = df_meta_multiome_train['cell_type'][mask], palette = 'rainbow', alpha = 0.9)\n        if 1:\n            plt.setp(ax.get_legend().get_texts(), fontsize=20) # for legend text\n            plt.setp(ax.get_legend().get_title(), fontsize=20) # for legend title        \n        plt.title('ATAC-seq. Train. PCA (TruncSVD). ' + str_inf  +   ' n_cells='+str(len(r)), fontsize = 20)\n        plt.xlabel('PCA'+str(i+1), fontsize = 20)\n        plt.ylabel('PCA'+str(j+1), fontsize = 20)\n        plt.show()    \n","metadata":{"execution":{"iopub.status.busy":"2022-09-16T08:35:07.695222Z","iopub.execute_input":"2022-09-16T08:35:07.696117Z","iopub.status.idle":"2022-09-16T08:39:03.449745Z","shell.execute_reply.started":"2022-09-16T08:35:07.696067Z","shell.execute_reply":"2022-09-16T08:39:03.448536Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# TruncatedSVD, visualizations fixed day and donor , look on  cell types\n","metadata":{}},{"cell_type":"code","source":"%%time\ndonor = 32606\nday = 2\nstr_inf = ' donor= ' + str(donor) + ' day= ' + str(day)\nmask = (df_meta_multiome_train['day'] == day) & (df_meta_multiome_train['donor'] == donor)\nprint(mask.sum())\ntsvd = TruncatedSVD(n_components = 32,random_state = 4223)\n\nr = tsvd.fit_transform(train_inputs[mask,:])","metadata":{"execution":{"iopub.status.busy":"2022-09-16T08:39:03.451335Z","iopub.execute_input":"2022-09-16T08:39:03.451698Z","iopub.status.idle":"2022-09-16T08:39:35.520254Z","shell.execute_reply.started":"2022-09-16T08:39:03.451665Z","shell.execute_reply":"2022-09-16T08:39:35.517576Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_meta_multiome_train[mask].groupby('cell_type')['Index in Multiome'].mean()","metadata":{"execution":{"iopub.status.busy":"2022-09-16T08:39:35.524034Z","iopub.execute_input":"2022-09-16T08:39:35.524547Z","iopub.status.idle":"2022-09-16T08:39:35.573707Z","shell.execute_reply.started":"2022-09-16T08:39:35.524493Z","shell.execute_reply":"2022-09-16T08:39:35.572538Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for (i,j) in [(0,1),(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4)]:\n    plt.figure(figsize = (20,10))\n    ax = sns.scatterplot(r[:,i],r[:,j], hue = df_meta_multiome_train['cell_type'][mask], palette = 'rainbow')\n    if 1:\n        plt.setp(ax.get_legend().get_texts(), fontsize=20) # for legend text\n        plt.setp(ax.get_legend().get_title(), fontsize=20) # for legend title        \n    plt.title('ATAC-seq. Train. PCA (TruncSVD). ' + str_inf  +   ' n_cells='+str(len(r)), fontsize = 20)\n    plt.xlabel('PCA'+str(i+1), fontsize = 20)\n    plt.ylabel('PCA'+str(j+1), fontsize = 20)\n    plt.show()    ","metadata":{"execution":{"iopub.status.busy":"2022-09-16T08:39:35.575642Z","iopub.execute_input":"2022-09-16T08:39:35.576011Z","iopub.status.idle":"2022-09-16T08:39:42.617467Z","shell.execute_reply.started":"2022-09-16T08:39:35.575976Z","shell.execute_reply":"2022-09-16T08:39:42.616451Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# UMAP after TrancatedSVD . ","metadata":{}},{"cell_type":"code","source":"%%time\nprint(r.shape)\nimport umap\nreducer = umap.UMAP()\nr2 = reducer.fit_transform(r)\n\n","metadata":{"execution":{"iopub.status.busy":"2022-09-16T08:39:42.618785Z","iopub.execute_input":"2022-09-16T08:39:42.619909Z","iopub.status.idle":"2022-09-16T08:40:35.910316Z","shell.execute_reply.started":"2022-09-16T08:39:42.619854Z","shell.execute_reply":"2022-09-16T08:40:35.908925Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for (i,j) in [(0,1)]:# ,(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4)]:\n    plt.figure(figsize = (20,10))\n    ax = sns.scatterplot(r2[:,i],r2[:,j], hue = df_meta_multiome_train['cell_type'][mask], palette = 'rainbow')\n    if 1:\n        plt.setp(ax.get_legend().get_texts(), fontsize=20) # for legend text\n        plt.setp(ax.get_legend().get_title(), fontsize=20) # for legend title        \n    plt.title('ATAC-seq. Train. UMAP after PCA (TruncSVD) 32. ' + str_inf  +   ' n_cells='+str(len(r)), fontsize = 20)\n    plt.xlabel('UMAP'+str(i+1), fontsize = 20)\n    plt.ylabel('UMAP'+str(j+1), fontsize = 20)\n    plt.show()    ","metadata":{"execution":{"iopub.status.busy":"2022-09-16T08:40:35.912507Z","iopub.execute_input":"2022-09-16T08:40:35.913923Z","iopub.status.idle":"2022-09-16T08:40:37.058497Z","shell.execute_reply.started":"2022-09-16T08:40:35.913860Z","shell.execute_reply":"2022-09-16T08:40:37.057412Z"},"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\n","metadata":{"execution":{"iopub.status.busy":"2022-09-16T08:40:37.059966Z","iopub.execute_input":"2022-09-16T08:40:37.060326Z","iopub.status.idle":"2022-09-16T08:41:12.064825Z","shell.execute_reply.started":"2022-09-16T08:40:37.060294Z","shell.execute_reply":"2022-09-16T08:41:12.063047Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n%%time\nprint(r.shape)\nimport ncvis\nreducer = ncvis.NCVis()  # umap.UMAP()\nr2 = reducer.fit_transform(r)","metadata":{"execution":{"iopub.status.busy":"2022-09-16T08:41:12.067063Z","iopub.execute_input":"2022-09-16T08:41:12.067530Z","iopub.status.idle":"2022-09-16T08:41:15.463395Z","shell.execute_reply.started":"2022-09-16T08:41:12.067488Z","shell.execute_reply":"2022-09-16T08:41:15.462147Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for (i,j) in [(0,1)]:# ,(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4)]:\n    plt.figure(figsize = (20,10))\n    ax = sns.scatterplot(r2[:,i],r2[:,j], hue = df_meta_multiome_train['cell_type'][mask], palette = 'rainbow')\n    if 1:\n        plt.setp(ax.get_legend().get_texts(), fontsize=20) # for legend text\n        plt.setp(ax.get_legend().get_title(), fontsize=20) # for legend title        \n    plt.title('ATAC-seq. Train. NCIVIS after PCA (TruncSVD) 32. ' + str_inf  +   ' n_cells='+str(len(r)), fontsize = 20)\n    plt.xlabel('NCIVIS'+str(i+1), fontsize = 20)\n    plt.ylabel('NCIVIS'+str(j+1), fontsize = 20)\n    plt.show()    ","metadata":{"execution":{"iopub.status.busy":"2022-09-16T08:41:15.465758Z","iopub.execute_input":"2022-09-16T08:41:15.466278Z","iopub.status.idle":"2022-09-16T08:41:16.246103Z","shell.execute_reply.started":"2022-09-16T08:41:15.466242Z","shell.execute_reply":"2022-09-16T08:41:16.244903Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# TrunctedSVD and Visualizations for day 2 ","metadata":{}},{"cell_type":"code","source":"%%time\n#donor = 32606\nday = 2\nstr_inf = ''\n#str_inf = ' donor= ' + str(donor) + \nstr_inf += ' day= ' + str(day)\nmask = (df_meta_multiome_train['day'] == day) # & (df_meta_multiome_train['donor'] == donor)\nprint(mask.sum())\ntsvd = TruncatedSVD(n_components = 10,random_state = 4223)\n\nr = tsvd.fit_transform(train_inputs[mask,:])","metadata":{"execution":{"iopub.status.busy":"2022-09-16T08:41:16.247651Z","iopub.execute_input":"2022-09-16T08:41:16.248009Z","iopub.status.idle":"2022-09-16T08:42:05.526870Z","shell.execute_reply.started":"2022-09-16T08:41:16.247975Z","shell.execute_reply":"2022-09-16T08:42:05.525419Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for (i,j) in [(0,1),(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4)]:\n    plt.figure(figsize = (20,10))\n    ax = sns.scatterplot(r[:,i],r[:,j], hue = df_meta_multiome_train['str_donor'][mask], palette = 'rainbow')\n    if 1:\n        plt.setp(ax.get_legend().get_texts(), fontsize=20) # for legend text\n        plt.setp(ax.get_legend().get_title(), fontsize=20) # for legend title        \n    plt.title('ATAC-seq. Train. PCA (TruncSVD). ' + str_inf  +   ' n_cells='+str(len(r)), fontsize = 20)\n    plt.xlabel('PCA'+str(i+1), fontsize = 20)\n    plt.ylabel('PCA'+str(j+1), fontsize = 20)\n    plt.show()    ","metadata":{"execution":{"iopub.status.busy":"2022-09-16T08:42:05.528494Z","iopub.execute_input":"2022-09-16T08:42:05.528891Z","iopub.status.idle":"2022-09-16T08:42:16.701717Z","shell.execute_reply.started":"2022-09-16T08:42:05.528855Z","shell.execute_reply":"2022-09-16T08:42:16.700329Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for (i,j) in [(0,1),(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4)]:\n    plt.figure(figsize = (20,10))\n    ax = sns.scatterplot(r[:,i],r[:,j], hue = df_meta_multiome_train['cell_type'][mask], palette = 'rainbow')\n    if 1:\n        plt.setp(ax.get_legend().get_texts(), fontsize=20) # for legend text\n        plt.setp(ax.get_legend().get_title(), fontsize=20) # for legend title        \n    plt.title('ATAC-seq. Train. PCA (TruncSVD). ' + str_inf  +   ' n_cells='+str(len(r)), fontsize = 20)\n    plt.xlabel('PCA'+str(i+1), fontsize = 20)\n    plt.ylabel('PCA'+str(j+1), fontsize = 20)\n    plt.show()    ","metadata":{"execution":{"iopub.status.busy":"2022-09-16T08:42:16.703287Z","iopub.execute_input":"2022-09-16T08:42:16.704105Z","iopub.status.idle":"2022-09-16T08:42:28.763325Z","shell.execute_reply.started":"2022-09-16T08:42:16.704066Z","shell.execute_reply":"2022-09-16T08:42:28.762282Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Fixed donor, change the days ","metadata":{}},{"cell_type":"code","source":"%%time\ndonor = 32606\n#day = 2\nstr_inf = ''\nstr_inf = ' donor= ' + str(donor) \n#str_inf += ' day= ' + str(day)\n#mask = (df_meta_multiome_train['day'] == day) # & (df_meta_multiome_train['donor'] == donor)\nmask = (df_meta_multiome_train['donor'] == donor)\n\nprint(mask.sum())\ntsvd = TruncatedSVD(n_components = 10,random_state = 4223)\nr = tsvd.fit_transform(train_inputs[mask,:])\n\n","metadata":{"execution":{"iopub.status.busy":"2022-09-16T08:42:28.764873Z","iopub.execute_input":"2022-09-16T08:42:28.766068Z","iopub.status.idle":"2022-09-16T08:43:44.435442Z","shell.execute_reply.started":"2022-09-16T08:42:28.766022Z","shell.execute_reply":"2022-09-16T08:43:44.433507Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for (i,j) in [(0,1),(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4)]:\n    plt.figure(figsize = (20,10))\n    ax = sns.scatterplot(r[:,i],r[:,j], hue = df_meta_multiome_train['day'][mask], palette = 'rainbow', alpha = 0.2)\n    if 1:\n        plt.setp(ax.get_legend().get_texts(), fontsize=20) # for legend text\n        plt.setp(ax.get_legend().get_title(), fontsize=20) # for legend title        \n    plt.title('ATAC-seq. Train. PCA (TruncSVD). ' + str_inf  +   ' n_cells='+str(len(r)), fontsize = 20)\n    plt.xlabel('PCA'+str(i+1), fontsize = 20)\n    plt.ylabel('PCA'+str(j+1), fontsize = 20)\n    plt.show()    \n","metadata":{"execution":{"iopub.status.busy":"2022-09-16T08:43:44.438531Z","iopub.execute_input":"2022-09-16T08:43:44.440293Z","iopub.status.idle":"2022-09-16T08:43:57.933376Z","shell.execute_reply.started":"2022-09-16T08:43:44.440229Z","shell.execute_reply":"2022-09-16T08:43:57.932154Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Visualizations for fixed cell_type","metadata":{}},{"cell_type":"code","source":"df_meta_multiome_train['cell_type'].value_counts()","metadata":{"execution":{"iopub.status.busy":"2022-09-16T08:43:57.935246Z","iopub.execute_input":"2022-09-16T08:43:57.935761Z","iopub.status.idle":"2022-09-16T08:43:57.957890Z","shell.execute_reply.started":"2022-09-16T08:43:57.935710Z","shell.execute_reply":"2022-09-16T08:43:57.956480Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndonor = 32606\n#day = 2\ncell_type = 'EryP'\nstr_inf = ''\nstr_inf = ' cell_type= ' + str(cell_type) \n#str_inf = ' donor= ' + str(donor) \n#str_inf += ' day= ' + str(day)\n#mask = (df_meta_multiome_train['day'] == day) # & (df_meta_multiome_train['donor'] == donor)\n#mask = (df_meta_multiome_train['donor'] == donor)\nmask = (df_meta_multiome_train['cell_type'] == cell_type)\n\nprint(mask.sum())\ntsvd = TruncatedSVD(n_components = 32,random_state = 4223)\nr = tsvd.fit_transform(train_inputs[mask,:])\n\n","metadata":{"execution":{"iopub.status.busy":"2022-09-16T08:43:57.960038Z","iopub.execute_input":"2022-09-16T08:43:57.960449Z","iopub.status.idle":"2022-09-16T08:45:14.970174Z","shell.execute_reply.started":"2022-09-16T08:43:57.960412Z","shell.execute_reply":"2022-09-16T08:45:14.968528Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for (i,j) in [(0,1),(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4)]:\n    plt.figure(figsize = (20,10))\n    ax = sns.scatterplot(r[:,i],r[:,j], hue = df_meta_multiome_train['day'][mask], palette = 'rainbow', alpha = 0.8)\n    if 1:\n        plt.setp(ax.get_legend().get_texts(), fontsize=20) # for legend text\n        plt.setp(ax.get_legend().get_title(), fontsize=20) # for legend title        \n    plt.title('ATAC-seq. Train. PCA (TruncSVD). ' + str_inf  +   ' n_cells='+str(len(r)), fontsize = 20)\n    plt.xlabel('PCA'+str(i+1), fontsize = 20)\n    plt.ylabel('PCA'+str(j+1), fontsize = 20)\n    plt.show() \nfor (i,j) in [(0,1),(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4)]:\n    plt.figure(figsize = (20,10))\n    ax = sns.scatterplot(r[:,i],r[:,j], hue = df_meta_multiome_train['str_donor'][mask], palette = 'rainbow', alpha = 0.8)\n    if 1:\n        plt.setp(ax.get_legend().get_texts(), fontsize=20) # for legend text\n        plt.setp(ax.get_legend().get_title(), fontsize=20) # for legend title        \n    plt.title('ATAC-seq. Train. PCA (TruncSVD). ' + str_inf  +   ' n_cells='+str(len(r)), fontsize = 20)\n    plt.xlabel('PCA'+str(i+1), fontsize = 20)\n    plt.ylabel('PCA'+str(j+1), fontsize = 20)\n    plt.show()   ","metadata":{"execution":{"iopub.status.busy":"2022-09-16T08:45:14.972692Z","iopub.execute_input":"2022-09-16T08:45:14.973087Z","iopub.status.idle":"2022-09-16T08:45:34.522129Z","shell.execute_reply.started":"2022-09-16T08:45:14.973048Z","shell.execute_reply":"2022-09-16T08:45:34.520810Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# UMAP / NSVIS for fixed cell type","metadata":{}},{"cell_type":"code","source":"%%time\nprint(r.shape)\nimport umap\nreducer = umap.UMAP()\nr2 = reducer.fit_transform(r)\n\nfor (i,j) in [(0,1)]:# ,(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4)]:\n    plt.figure(figsize = (20,10))\n    ax = sns.scatterplot(r2[:,i],r2[:,j], hue = df_meta_multiome_train['day'][mask], palette = 'rainbow')\n    if 1:\n        plt.setp(ax.get_legend().get_texts(), fontsize=20) # for legend text\n        plt.setp(ax.get_legend().get_title(), fontsize=20) # for legend title        \n    plt.title('ATAC-seq. Train. UMAP after PCA (TruncSVD) 32. ' + str_inf  +   ' n_cells='+str(len(r)), fontsize = 20)\n    plt.xlabel('UMAP'+str(i+1), fontsize = 20)\n    plt.ylabel('UMAP'+str(j+1), fontsize = 20)\n    plt.show()    \n","metadata":{"execution":{"iopub.status.busy":"2022-09-16T08:45:34.528276Z","iopub.execute_input":"2022-09-16T08:45:34.528649Z","iopub.status.idle":"2022-09-16T08:45:47.392453Z","shell.execute_reply.started":"2022-09-16T08:45:34.528617Z","shell.execute_reply":"2022-09-16T08:45:47.391134Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nprint(r.shape)\nimport ncvis\nreducer = ncvis.NCVis()  # umap.UMAP()\nr2 = reducer.fit_transform(r)\n\nfor (i,j) in [(0,1)]:# ,(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4)]:\n    plt.figure(figsize = (20,10))\n    ax = sns.scatterplot(r2[:,i],r2[:,j], hue = df_meta_multiome_train['day'][mask], palette = 'rainbow')\n    if 1:\n        plt.setp(ax.get_legend().get_texts(), fontsize=20) # for legend text\n        plt.setp(ax.get_legend().get_title(), fontsize=20) # for legend title        \n    plt.title('ATAC-seq. Train. NCIVIS after PCA (TruncSVD) 32. ' + str_inf  +   ' n_cells='+str(len(r)), fontsize = 20)\n    plt.xlabel('NCIVIS'+str(i+1), fontsize = 20)\n    plt.ylabel('NCIVIS'+str(j+1), fontsize = 20)\n    plt.show()  ","metadata":{"execution":{"iopub.status.busy":"2022-09-16T08:45:47.393865Z","iopub.execute_input":"2022-09-16T08:45:47.394212Z","iopub.status.idle":"2022-09-16T08:45:55.418135Z","shell.execute_reply.started":"2022-09-16T08:45:47.394165Z","shell.execute_reply":"2022-09-16T08:45:55.416864Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for (i,j) in [(0,1)]:# ,(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4)]:\n    plt.figure(figsize = (20,10))\n    ax = sns.scatterplot(r2[:,i],r2[:,j], hue = df_meta_multiome_train['str_donor'][mask], palette = 'rainbow')\n    if 1:\n        plt.setp(ax.get_legend().get_texts(), fontsize=20) # for legend text\n        plt.setp(ax.get_legend().get_title(), fontsize=20) # for legend title        \n    plt.title('ATAC-seq. Train. NCIVIS after PCA (TruncSVD) 32. ' + str_inf  +   ' n_cells='+str(len(r)), fontsize = 20)\n    plt.xlabel('NCIVIS'+str(i+1), fontsize = 20)\n    plt.ylabel('NCIVIS'+str(j+1), fontsize = 20)\n    plt.show()  ","metadata":{"execution":{"iopub.status.busy":"2022-09-16T08:45:55.419883Z","iopub.execute_input":"2022-09-16T08:45:55.420613Z","iopub.status.idle":"2022-09-16T08:45:56.411662Z","shell.execute_reply.started":"2022-09-16T08:45:55.420567Z","shell.execute_reply":"2022-09-16T08:45:56.410507Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Loop over  donors - compare PCA, UMAP  - constructed for each of the donors for day 2  ","metadata":{}},{"cell_type":"code","source":"%%time\nday = 2\n\nfor donor in (df_meta_multiome_train['donor'].unique()):\n    #donor = 32606\n    str_inf = ''\n    str_inf = ' donor= ' + str(donor)  \n    str_inf += ' day= ' + str(day)\n    mask = (df_meta_multiome_train['day'] == day)  & (df_meta_multiome_train['donor'] == donor)\n    print(mask.sum(), str_inf )\n    t0 = time.time()\n    tsvd = TruncatedSVD(n_components = 32,random_state = 4223)\n\n    r = tsvd.fit_transform(train_inputs[mask,:])\n    print('%.1f'%(time.time()-t0) + ' seconds passed. TruncatedSVD finished.')\n\n    for (i,j) in [(0,1)]: # ,(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4)]:\n        plt.figure(figsize = (20,10))\n        ax = sns.scatterplot(r[:,i],r[:,j], hue = df_meta_multiome_train['cell_type'][mask], palette = 'rainbow')\n        if 1:\n            plt.setp(ax.get_legend().get_texts(), fontsize=20) # for legend text\n            plt.setp(ax.get_legend().get_title(), fontsize=20) # for legend title        \n        plt.title('ATAC-seq. Train. PCA (TruncSVD). ' + str_inf  +   ' n_cells='+str(len(r)), fontsize = 20)\n        plt.xlabel('PCA'+str(i+1), fontsize = 20)\n        plt.ylabel('PCA'+str(j+1), fontsize = 20)\n        plt.show()    \n\n    import umap\n    reducer = umap.UMAP()\n    r2 = reducer.fit_transform(r)\n    print('%.1f'%(time.time()-t0) + ' seconds passed. UMAP finished.')\n    \n    for (i,j) in [(0,1)]:# ,(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4)]:\n        plt.figure(figsize = (20,10))\n        ax = sns.scatterplot(r2[:,i],r2[:,j], hue = df_meta_multiome_train['cell_type'][mask], palette = 'rainbow')\n        if 1:\n            plt.setp(ax.get_legend().get_texts(), fontsize=20) # for legend text\n            plt.setp(ax.get_legend().get_title(), fontsize=20) # for legend title        \n        plt.title('ATAC-seq. Train. UMAP after PCA (TruncSVD) 32. ' + str_inf  +   ' n_cells='+str(len(r)), fontsize = 20)\n        plt.xlabel('UMAP'+str(i+1), fontsize = 20)\n        plt.ylabel('UMAP'+str(j+1), fontsize = 20)\n        plt.show()               \n              ","metadata":{"execution":{"iopub.status.busy":"2022-09-16T08:45:56.413364Z","iopub.execute_input":"2022-09-16T08:45:56.413721Z","iopub.status.idle":"2022-09-16T08:48:13.657986Z","shell.execute_reply.started":"2022-09-16T08:45:56.413687Z","shell.execute_reply":"2022-09-16T08:48:13.656756Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# At day 2 - TruncatedSVD and UMAP, colored by donors - seems no batch effect from donors. ","metadata":{}},{"cell_type":"code","source":"%%time\nday = 2\n\n#for donor in (df_meta_multiome_train['donor'].unique()):\nif 1:\n    #donor = 32606\n    str_inf = ''\n    #str_inf = ' donor= ' + str(donor)  \n    str_inf += ' day= ' + str(day)\n    mask = (df_meta_multiome_train['day'] == day)#  & (df_meta_multiome_train['donor'] == donor)\n    print(mask.sum(), str_inf )\n    t0 = time.time()\n    tsvd = TruncatedSVD(n_components = 32,random_state = 4223)\n\n    r = tsvd.fit_transform(train_inputs[mask,:])\n    print('%.1f'%(time.time()-t0) + ' seconds passed. TruncatedSVD finished.')\n\n    for (i,j) in [(0,1)]: # ,(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4)]:\n        plt.figure(figsize = (20,10))\n        ax = sns.scatterplot(r[:,i],r[:,j], hue = df_meta_multiome_train['donor'][mask], palette = 'rainbow')\n        if 1:\n            plt.setp(ax.get_legend().get_texts(), fontsize=20) # for legend text\n            plt.setp(ax.get_legend().get_title(), fontsize=20) # for legend title        \n        plt.title('ATAC-seq. Train. PCA (TruncSVD). ' + str_inf  +   ' n_cells='+str(len(r)), fontsize = 20)\n        plt.xlabel('PCA'+str(i+1), fontsize = 20)\n        plt.ylabel('PCA'+str(j+1), fontsize = 20)\n        plt.show()    \n\n    import umap\n    reducer = umap.UMAP()\n    r2 = reducer.fit_transform(r)\n    print('%.1f'%(time.time()-t0) + ' seconds passed. UMAP finished.')\n    \n    for (i,j) in [(0,1)]:# ,(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4)]:\n        plt.figure(figsize = (20,10))\n        ax = sns.scatterplot(r2[:,i],r2[:,j], hue = df_meta_multiome_train['donor'][mask], palette = 'rainbow')\n        if 1:\n            plt.setp(ax.get_legend().get_texts(), fontsize=20) # for legend text\n            plt.setp(ax.get_legend().get_title(), fontsize=20) # for legend title        \n        plt.title('ATAC-seq. Train. UMAP after PCA (TruncSVD) 32. ' + str_inf  +   ' n_cells='+str(len(r)), fontsize = 20)\n        plt.xlabel('UMAP'+str(i+1), fontsize = 20)\n        plt.ylabel('UMAP'+str(j+1), fontsize = 20)\n        plt.show()               \n              ","metadata":{"execution":{"iopub.status.busy":"2022-09-16T08:48:13.659868Z","iopub.execute_input":"2022-09-16T08:48:13.660412Z","iopub.status.idle":"2022-09-16T08:50:02.640748Z","shell.execute_reply.started":"2022-09-16T08:48:13.660372Z","shell.execute_reply":"2022-09-16T08:50:02.639443Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Color by chromosome sums - seems  nothing additional to PCA1 is strongly related to sum. ","metadata":{}},{"cell_type":"code","source":"%%time\nfor donor in [32606]:# df_meta_multiome_train['donor'].unique():\n    #donor = 32606\n    #day = 2\n    str_inf = ''\n    str_inf = ' donor= ' + str(donor) \n    #str_inf += ' day= ' + str(day)\n    #mask = (df_meta_multiome_train['day'] == day) # & (df_meta_multiome_train['donor'] == donor)\n    mask = (df_meta_multiome_train['donor'] == donor)\n\n    print(mask.sum(), str_inf )\n    tsvd = TruncatedSVD(n_components = 10,random_state = 4223)\n    t0 = time.time()\n    r = tsvd.fit_transform(train_inputs[mask,:])\n    print('%.1f'%(time.time()-t0) + ' seconds passed. TruncatedSVD finished.')\n","metadata":{"execution":{"iopub.status.busy":"2022-09-16T15:34:16.019714Z","iopub.execute_input":"2022-09-16T15:34:16.021096Z","iopub.status.idle":"2022-09-16T15:35:28.195499Z","shell.execute_reply.started":"2022-09-16T15:34:16.021027Z","shell.execute_reply":"2022-09-16T15:35:28.193922Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"c = 0\nn_x_subplots = 4\nfor col in df_meta_multiome_train.columns:\n    if 'Sum' not in col: continue \n    v = df_meta_multiome_train[mask][col]\n    for (i,j) in [ (0,1)] : # ,(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4)]:\n        #plt.figure(figsize = (20,10))\n\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        c += 1; fig.add_subplot(1,n_x_subplots ,c)\n        \n        ax = sns.scatterplot(x = r[:,i],y = r[:,j], hue = v, palette = 'rainbow', alpha = 0.9)\n        if 1:\n            plt.setp(ax.get_legend().get_texts(), fontsize=20) # for legend text\n            plt.setp(ax.get_legend().get_title(), fontsize=20) # for legend title        \n        plt.suptitle('ATAC-seq. Train. PCA (TruncSVD). ' + str_inf  +   ' n_cells='+str(len(r)), fontsize = 20)\n        plt.xlabel('PCA'+str(i+1), fontsize = 20)\n        plt.ylabel('PCA'+str(j+1), fontsize = 20)\nplt.show()    \n","metadata":{"execution":{"iopub.status.busy":"2022-09-16T15:42:25.612666Z","iopub.execute_input":"2022-09-16T15:42:25.613079Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_meta_multiome_train.columns","metadata":{"execution":{"iopub.status.busy":"2022-09-16T15:48:57.942896Z","iopub.execute_input":"2022-09-16T15:48:57.943506Z","iopub.status.idle":"2022-09-16T15:48:57.952805Z","shell.execute_reply.started":"2022-09-16T15:48:57.943461Z","shell.execute_reply":"2022-09-16T15:48:57.951817Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"c = 0\nn_x_subplots = 4\nfor col in ['Sum chr1', 'Sum chr13',  'Sum chrY', 'Sum chrX' , 'Sum GL0', 'Sum KI2']: # df_meta_multiome_train.columns:\n    if 'Sum' not in col: continue \n    v = df_meta_multiome_train[mask][col]\n    for (i,j) in [ (0,1),(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4)]:\n        #plt.figure(figsize = (20,10))\n\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        c += 1; fig.add_subplot(1,n_x_subplots ,c)\n        \n        ax = sns.scatterplot(x = r[:,i],y = r[:,j], hue = v, palette = 'rainbow', alpha = 0.9)\n        if 1:\n            plt.setp(ax.get_legend().get_texts(), fontsize=20) # for legend text\n            plt.setp(ax.get_legend().get_title(), fontsize=20) # for legend title        \n        plt.suptitle('ATAC-seq. Train. PCA (TruncSVD). ' + str_inf  +   ' n_cells='+str(len(r)), fontsize = 20)\n        plt.xlabel('PCA'+str(i+1), fontsize = 20)\n        plt.ylabel('PCA'+str(j+1), fontsize = 20)\nplt.show()    \n","metadata":{"execution":{"iopub.status.busy":"2022-09-16T15:50:23.041125Z","iopub.execute_input":"2022-09-16T15:50:23.041689Z","iopub.status.idle":"2022-09-16T15:52:25.739943Z","shell.execute_reply.started":"2022-09-16T15:50:23.041651Z","shell.execute_reply":"2022-09-16T15:52:25.738404Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('%.1f seconds passed total '%(time.time()-t0start) )\n","metadata":{},"execution_count":null,"outputs":[]}]}