{"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\n**Briefly:** Example how to work with huge .h5 files without loading them in RAM. \nSee also example in https://www.kaggle.com/code/alexandervc/use-h5py-for-huge-h5-file-backing-it-on-disk \n\n**Data** from the competition , but denoising preprocessing by MAGIC package (https://github.com/KrishnaswamyLab/MAGIC) has been done by Liza Geraseva. See discussion: https://www.kaggle.com/competitions/open-problems-multimodal/discussion/350856\n\n**Work principle:** The main data can be accessed as f['0']['block0_values'] - without loading to RAM, but  when one takes a slice: f['0']['block0_values'][:,IX] - it is converted to numpy array (and loaded into memory?). \nThere are some details: one can do constant slices like  f['0']['block0_values'][:10,:10], but to make a  variable slice one need to process in two steps: f['0']['block0_values'][:,IX1][IX0]. \n\n**Details:** We load the data. And as an exercise and sanity check we produce plots for cell proliferation cycle visualization - so-called G1/S-G2/M plots. The results are as expected - the \"cyclic\" structes is better seen comparing as to original files (see version 3 of the present notebook). Remark: one can compare with another (more simple denoising (\"pooling\")) here: https://www.kaggle.com/code/alexandervc/mmscel-cell-cycle-03-daybydaychange-megakaryocyte#G1/S-G2/M-plots-to-analyse-the-cell-cycle\n\n**Context:** See paper: https://arxiv.org/abs/2208.05229 for cell cycle analysis. \"Computational challenges of cell cycle analysis using single cell transcriptomics\" Alexander Chervov, Andrei Zinovyev. All the work for that paper has been done on Kaggle - see datasets and code therein: https://www.kaggle.com/alexandervc/datasets?scroll=true https://www.kaggle.com/andreizinovyev Hundreds of single cell RNA seq datasets were analyzed.\n\n\n\n### Remarks: \n    \n    ( in Vesrstion 1- 6 we look mainly on Megakaryocyte Progenitors - just almost smallest and representative (from previous analysis ) \n\n### Versions:\n\n#### 10 - for comparaison - same data as 9, but NO MAGIC \n\n\n#### 9 - added plots all cell types\n    \n    Same MAGIC CITEseq TEST\n\n#### 7,8 cosmetic changes\n\n\n#### 6 - return to MAGIC - same as version 4 - \"test\" part of CITE-seq (features)\n        \n\n#### 5 - for comparaison  NO(!) MAGIC - difference is STRIKING\n\n    For days 2,4 - we cannot see clear \"cyclic\" structures without denoising.\n    Compare with denoised versions - where \"cyclic\" structure can be clearly seen ! (E.g. version 4) \n    \"test\" part of CITE-seq (features)\n\n\n#### 4 - MAGIC on \"test\" part of CITE-seq (features)\n\n    \"Cyclic\" structure is seen VERY well for many days - may be because only one donor, but for \"train\" part we have 3 donors and they mix and create less clear picture\n\n#### 3 - for comparaison look on original input file\n\n    The \"cyclic\" structure is less clearly seen for original files, comparing to denoised - as expected. \n\n#### 1,2 - work with MAGIC denoised data. \n\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-10-17T19:06:47.889287Z","iopub.execute_input":"2022-10-17T19:06:47.890036Z","iopub.status.idle":"2022-10-17T19:06:47.932578Z","shell.execute_reply.started":"2022-10-17T19:06:47.889942Z","shell.execute_reply":"2022-10-17T19:06:47.931465Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\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\n\n\nimport 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\n# !pip install scanpy\n# import scanpy as sc\n# import anndata\n\nDATA_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-10-17T19:06:47.934536Z","iopub.execute_input":"2022-10-17T19:06:47.934842Z","iopub.status.idle":"2022-10-17T19:07:14.259448Z","shell.execute_reply.started":"2022-10-17T19:06:47.934814Z","shell.execute_reply":"2022-10-17T19:07:14.258250Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# \"Load\" data and work with it ","metadata":{}},{"cell_type":"code","source":"print('Look on Features: ')\nfilename = '/kaggle/input/citeseq-denoised/test_cite_inputs_denoised.h5'\nfilename = '/kaggle/input/citeseq-denoised/train_cite_inputs_denoised.h5'\n#filename = '/kaggle/input/open-problems-multimodal/train_multi_inputs.h5'\n#filename = '/kaggle/input/open-problems-multimodal/train_multi_inputs.h5'\nfilename = '/kaggle/input/open-problems-multimodal/train_cite_inputs.h5'\n\nfilename = '/kaggle/input/citeseq-denoised/test_cite_inputs_denoised.h5'\n\nfilename = '/kaggle/input/open-problems-multimodal/test_cite_inputs.h5'\n\nf2 = h5py.File(filename,'r')#, mode)","metadata":{"execution":{"iopub.status.busy":"2022-10-17T19:07:14.260860Z","iopub.execute_input":"2022-10-17T19:07:14.261206Z","iopub.status.idle":"2022-10-17T19:07:14.273875Z","shell.execute_reply.started":"2022-10-17T19:07:14.261172Z","shell.execute_reply":"2022-10-17T19:07:14.272619Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f2.keys() )","metadata":{"execution":{"iopub.status.busy":"2022-10-17T19:07:14.276871Z","iopub.execute_input":"2022-10-17T19:07:14.277460Z","iopub.status.idle":"2022-10-17T19:07:14.288144Z","shell.execute_reply.started":"2022-10-17T19:07:14.277415Z","shell.execute_reply":"2022-10-17T19:07:14.287061Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"str_data_inf = ''\nif filename == '/kaggle/input/open-problems-multimodal/train_cite_inputs.h5':\n    first_key = 'train_cite_inputs'\n    str_data_inf = ' CITEseq Train '\n# elif filename == '/kaggle/input/open-problems-multimodal/train_multi_inputs.h5':\n#     first_key = 'train_multi_inputs'\n#     str_data_inf = ' Multi Train Features '\nelif filename == '/kaggle/input/open-problems-multimodal/test_cite_inputs.h5':\n    first_key = 'test_cite_inputs'#'train_multi_inputs' \n    str_data_inf = ' CITEseq Test '\nelif filename == '/kaggle/input/citeseq-denoised/test_cite_inputs_denoised.h5':\n    first_key = '0'\n    str_data_inf = ' MAGIC CITEseq Test '\nelif filename == '/kaggle/input/citeseq-denoised/train_cite_inputs_denoised.h5':\n    first_key = '0'\n    str_data_inf = ' MAGIC CITEseq Train '\n# else:\n#     first_key = '0'#'train_multi_inputs' \n    \nprint( f2[first_key].keys() )","metadata":{"execution":{"iopub.status.busy":"2022-10-17T19:07:14.289312Z","iopub.execute_input":"2022-10-17T19:07:14.289612Z","iopub.status.idle":"2022-10-17T19:07:14.302146Z","shell.execute_reply.started":"2022-10-17T19:07:14.289586Z","shell.execute_reply":"2022-10-17T19:07:14.301118Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for k in f2[first_key].keys():\n    print(type(f2[first_key][k]), f2[first_key][k].shape)","metadata":{"execution":{"iopub.status.busy":"2022-10-17T19:07:14.303357Z","iopub.execute_input":"2022-10-17T19:07:14.303637Z","iopub.status.idle":"2022-10-17T19:07:14.324804Z","shell.execute_reply.started":"2022-10-17T19:07:14.303611Z","shell.execute_reply":"2022-10-17T19:07:14.323990Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"f2[first_key]['axis0'][0], f2[first_key]['axis1'][0], f2[first_key]['block0_items'][0],  f2[first_key]['block0_values'][0,0], ","metadata":{"execution":{"iopub.status.busy":"2022-10-17T19:07:14.326167Z","iopub.execute_input":"2022-10-17T19:07:14.326461Z","iopub.status.idle":"2022-10-17T19:07:14.352378Z","shell.execute_reply.started":"2022-10-17T19:07:14.326433Z","shell.execute_reply":"2022-10-17T19:07:14.351184Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"type(f2[first_key]['block0_values'][:10,:10])","metadata":{"execution":{"iopub.status.busy":"2022-10-17T19:07:14.353994Z","iopub.execute_input":"2022-10-17T19:07:14.355179Z","iopub.status.idle":"2022-10-17T19:07:14.369814Z","shell.execute_reply.started":"2022-10-17T19:07:14.355132Z","shell.execute_reply":"2022-10-17T19:07:14.368398Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print( dir(f2[first_key]['block0_values'][:10,:10]) )","metadata":{"execution":{"iopub.status.busy":"2022-10-17T19:07:14.374223Z","iopub.execute_input":"2022-10-17T19:07:14.374576Z","iopub.status.idle":"2022-10-17T19:07:14.383195Z","shell.execute_reply.started":"2022-10-17T19:07:14.374545Z","shell.execute_reply":"2022-10-17T19:07:14.382076Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"list_genes_ensembl_ids = [t.decode().split('_')[0] for t in f2[first_key]['axis0'] ]\nprint( len(list_genes_ensembl_ids), list_genes_ensembl_ids[:10] )\nlist_genes_symbols = [t.decode().split('_')[1] for t in f2[first_key]['axis0'] ]\nprint( len(list_genes_symbols), list_genes_symbols[:10] )","metadata":{"execution":{"iopub.status.busy":"2022-10-17T19:07:14.384538Z","iopub.execute_input":"2022-10-17T19:07:14.385286Z","iopub.status.idle":"2022-10-17T19:07:18.314429Z","shell.execute_reply.started":"2022-10-17T19:07:14.385241Z","shell.execute_reply":"2022-10-17T19:07:18.313366Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"list_cells_barcodes = [t.decode() for t in f2[first_key]['axis1'] ]\nprint(len(list_cells_barcodes) , list_cells_barcodes[:10] )","metadata":{"execution":{"iopub.status.busy":"2022-10-17T19:07:18.315843Z","iopub.execute_input":"2022-10-17T19:07:18.316299Z","iopub.status.idle":"2022-10-17T19:07:22.584743Z","shell.execute_reply.started":"2022-10-17T19:07:18.316256Z","shell.execute_reply":"2022-10-17T19:07:22.583655Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print( type(f2[first_key]['block0_values']), dir(f2[first_key]['block0_values'] ) )","metadata":{"execution":{"iopub.status.busy":"2022-10-17T19:07:22.586128Z","iopub.execute_input":"2022-10-17T19:07:22.586545Z","iopub.status.idle":"2022-10-17T19:07:22.594256Z","shell.execute_reply.started":"2022-10-17T19:07:22.586504Z","shell.execute_reply":"2022-10-17T19:07:22.593080Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# G1/S-G2/M cell proliferation cycle plot for all cells \n\nIt will not be very clear - we need to separate days and cell types - see next subsections.","metadata":{}},{"cell_type":"code","source":"\nG1S_genes_Tirosh = ['ENSMUSG00000005410', 'ENSMUSG00000027342', 'ENSMUSG00000025747', 'ENSMUSG00000024742', 'ENSMUSG00000002870', 'ENSMUSG00000022673', 'ENSMUSG00000030978', 'ENSMUSG00000029591', 'ENSMUSG00000031821', 'ENSMUSG00000026355', 'ENSMUSG00000055612', 'ENSMUSG00000037474', 'ENSMUSG00000025395', 'ENSMUSG00000001228', 'ENSMUSG00000031629', 'ENSMUSG00000025001', 'ENSMUSG00000023104', 'ENSMUSG00000028884', 'ENSMUSG00000028693', 'ENSMUSG00000030346', 'ENSMUSG00000006715', 'ENSMUSG00000027242', 'ENSMUSG00000004642', 'ENSMUSG00000028212', 'ENSMUSG00000041712', 'ENSMUSG00000030726', 'ENSMUSG00000024151', 'ENSMUSG00000022360', 'ENSMUSG00000027323', 'ENSMUSG00000020649', 'ENSMUSG00000000028', 'ENSMUSG00000017499', 'ENSMUSG00000039748', 'ENSMUSG00000032397', 'ENSMUSG00000022422', 'ENSMUSG00000030528', 'ENSMUSG00000028282', 'ENSMUSG00000028560', 'ENSMUSG00000042489', 'ENSMUSG00000006678', 'ENSMUSG00000022945', 'ENSMUSG00000034329', 'ENSMUSG00000046179']\nG2M_genes_Tirosh = ['ENSMUSG00000054717', 'ENSMUSG00000019942', 'ENSMUSG00000027306', 'ENSMUSG00000001403', 'ENSMUSG00000017716', 'ENSMUSG00000027469', 'ENSMUSG00000020914', 'ENSMUSG00000024056', 'ENSMUSG00000062248', 'ENSMUSG00000026683', 'ENSMUSG00000028044', 'ENSMUSG00000031004', 'ENSMUSG00000019961', 'ENSMUSG00000026605', 'ENSMUSG00000037313', 'ENSMUSG00000020808', 'ENSMUSG00000034349', 'ENSMUSG00000032218', 'ENSMUSG00000048327', 'ENSMUSG00000037725', 'ENSMUSG00000020897', 'ENSMUSG00000027379', 'ENSMUSG00000012443', 'ENSMUSG00000015749', 'ENSMUSG00000036752', 'ENSMUSG00000022385', 'ENSMUSG00000024795', 'ENSMUSG00000044783', 'ENSMUSG00000023505', 'ENSMUSG00000020737', 'ENSMUSG00000006398', 'ENSMUSG00000038379', 'ENSMUSG00000044201', 'ENSMUSG00000028678', 'ENSMUSG00000022391', 'ENSMUSG00000038252', 'ENSMUSG00000037544', 'ENSMUSG00000048922', 'ENSMUSG00000028873', 'ENSMUSG00000027699', 'ENSMUSG00000032254', 'ENSMUSG00000020330', 'ENSMUSG00000027496', 'ENSMUSG00000068744', 'ENSMUSG00000036777', 'ENSMUSG00000004880', 'ENSMUSG00000040549', 'ENSMUSG00000045328', 'ENSMUSG00000005698', 'ENSMUSG00000026622', 'ENSMUSG00000035293', 'ENSMUSG00000074802', 'ENSMUSG00000009575', 'ENSMUSG00000029177']\n\nG1S_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']\nlist_genes_fastCCsign = ['CDK1', 'UBE2C', 'TOP2A', 'TMPO', 'HJURP', 'RRM1', 'RAD51AP1', 'RRM2', 'CDC45', 'BLM', 'BRIP1', 'E2F8', 'HIST2H2AC']\n\nG1S_genes_Tirosh = ['ENSG00000100297', 'ENSG00000132646', 'ENSG00000176890', 'ENSG00000168496', 'ENSG00000073111', 'ENSG00000104738', 'ENSG00000167325', 'ENSG00000076248', 'ENSG00000131153', 'ENSG00000076003', 'ENSG00000144354', 'ENSG00000143476', 'ENSG00000198056', 'ENSG00000276043', 'ENSG00000151725', 'ENSG00000119969', 'ENSG00000049541', 'ENSG00000117748', 'ENSG00000132780', 'ENSG00000111247', 'ENSG00000112312', 'ENSG00000092470', 'ENSG00000163950', 'ENSG00000175305', 'ENSG00000012963', 'ENSG00000077514', 'ENSG00000095002', 'ENSG00000156802', 'ENSG00000051180', 'ENSG00000171848', 'ENSG00000093009', 'ENSG00000094804', 'ENSG00000174371', 'ENSG00000075131', 'ENSG00000136982', 'ENSG00000197299', 'ENSG00000118412', 'ENSG00000162607', 'ENSG00000092853', 'ENSG00000101868', 'ENSG00000159259', 'ENSG00000136492', 'ENSG00000129173']\nG2M_genes_Tirosh = ['ENSG00000164104', 'ENSG00000170312', 'ENSG00000137804', 'ENSG00000175063', 'ENSG00000089685', 'ENSG00000088325', 'ENSG00000131747', 'ENSG00000080986', 'ENSG00000123975', 'ENSG00000143228', 'ENSG00000173207', 'ENSG00000148773', 'ENSG00000120802', 'ENSG00000117724', 'ENSG00000013810', 'ENSG00000129195', 'ENSG00000113810', 'ENSG00000157456', 'ENSG00000169607', 'ENSG00000136108', 'ENSG00000178999', 'ENSG00000169679', 'ENSG00000138160', 'ENSG00000143401', 'ENSG00000188229', 'ENSG00000075218', 'ENSG00000138182', 'ENSG00000123485', 'ENSG00000111665', 'ENSG00000189159', 'ENSG00000117399', 'ENSG00000112742', 'ENSG00000158402', 'ENSG00000142945', 'ENSG00000100401', 'ENSG00000010292', 'ENSG00000126787', 'ENSG00000184661', 'ENSG00000134690', 'ENSG00000114346', 'ENSG00000137807', 'ENSG00000072571', 'ENSG00000087586', 'ENSG00000134222', 'ENSG00000011426', 'ENSG00000143815', 'ENSG00000175216', 'ENSG00000138778', 'ENSG00000102974', 'ENSG00000117650', 'ENSG00000092140', 'ENSG00000139354', 'ENSG00000094916', 'ENSG00000115163']\nlist_genes_fastCCsign = ['ENSG00000170312', 'ENSG00000175063', 'ENSG00000131747', 'ENSG00000120802', 'ENSG00000123485', 'ENSG00000167325', 'ENSG00000111247', 'ENSG00000171848', 'ENSG00000093009', 'ENSG00000197299', 'ENSG00000136492', 'ENSG00000129173', 'ENSG00000184260']\n\n\n","metadata":{"execution":{"iopub.status.busy":"2022-10-17T19:07:22.596158Z","iopub.execute_input":"2022-10-17T19:07:22.596761Z","iopub.status.idle":"2022-10-17T19:07:22.618782Z","shell.execute_reply.started":"2022-10-17T19:07:22.596717Z","shell.execute_reply":"2022-10-17T19:07:22.617758Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"IX = np.where(pd.Series(list_genes_ensembl_ids ).isin(G1S_genes_Tirosh) >0) [0]\nlen(IX), len( G1S_genes_Tirosh)","metadata":{"execution":{"iopub.status.busy":"2022-10-17T19:07:22.620268Z","iopub.execute_input":"2022-10-17T19:07:22.620701Z","iopub.status.idle":"2022-10-17T19:07:22.646019Z","shell.execute_reply.started":"2022-10-17T19:07:22.620668Z","shell.execute_reply":"2022-10-17T19:07:22.644885Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nv1 = f2[first_key]['block0_values'][:,IX].sum(axis = 1 )","metadata":{"execution":{"iopub.status.busy":"2022-10-17T19:07:22.647075Z","iopub.execute_input":"2022-10-17T19:07:22.647402Z","iopub.status.idle":"2022-10-17T19:07:51.707896Z","shell.execute_reply.started":"2022-10-17T19:07:22.647373Z","shell.execute_reply":"2022-10-17T19:07:51.706778Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"IX = np.where(pd.Series(list_genes_ensembl_ids ).isin(G2M_genes_Tirosh) >0) [0]\nlen(IX), len( G2M_genes_Tirosh)","metadata":{"execution":{"iopub.status.busy":"2022-10-17T19:07:51.709181Z","iopub.execute_input":"2022-10-17T19:07:51.709508Z","iopub.status.idle":"2022-10-17T19:07:51.721472Z","shell.execute_reply.started":"2022-10-17T19:07:51.709478Z","shell.execute_reply":"2022-10-17T19:07:51.720391Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nv2 = f2[first_key]['block0_values'][:,IX].sum(axis = 1 )\nprint( len(v2), type(v2), v2.shape) ","metadata":{"execution":{"iopub.status.busy":"2022-10-17T19:07:51.722882Z","iopub.execute_input":"2022-10-17T19:07:51.723206Z","iopub.status.idle":"2022-10-17T19:08:04.615983Z","shell.execute_reply.started":"2022-10-17T19:07:51.723177Z","shell.execute_reply":"2022-10-17T19:08:04.614699Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.scatterplot(x = v1, y=v2)\nplt.title(str_data_inf)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-17T19:08:04.617499Z","iopub.execute_input":"2022-10-17T19:08:04.617982Z","iopub.status.idle":"2022-10-17T19:08:04.982886Z","shell.execute_reply.started":"2022-10-17T19:08:04.617934Z","shell.execute_reply":"2022-10-17T19:08:04.981692Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# G1/S-G2/M plot for day=4, cell_type = MkP (Megakaryocyte Progenitors)","metadata":{}},{"cell_type":"code","source":"%%time\n\ntrain_target = h5py.File(FP_CITE_TRAIN_TARGETS,'r')#, mode)\nhue_first_key = 'train_cite_targets'\ngenes = train_target[hue_first_key]['axis0'][:]\n\ndf_cite_train_y = pd.read_hdf('../input/open-problems-multimodal/train_cite_targets.h5')\ndf_cell_hue = df_cell.merge(df_cite_train_y, on = 'cell_id')\n#display(df_cell_hue[df_cell_hue['day']==4][df_cell_hue['cell_type']=='MkP'].head(10))\n\ndisplay( df_cell_hue.head(1) )\nprint(df_cell_hue['cell_type'].unique())\nmask = (df_cell_hue['day'] == 4) & (df_cell_hue['cell_type'] == 'MkP' )\n                                \nprint( mask.sum() )\n\nmask = (df_cell_hue['day'] == 4) & (df_cell_hue['cell_type'] == 'MkP' )                      \nmask.sum()\ncell_ids = list(df_cell_hue['cell_id'][mask])\nprint(len(cell_ids) , cell_ids[:10])\n\n\nIX = np.where(pd.Series(list_genes_ensembl_ids ).isin(G1S_genes_Tirosh) >0) [0]\nprint( len(IX), len( G1S_genes_Tirosh) )\nv1 = (f2[first_key]['block0_values'][:,IX]).sum(axis = 1 )\n# Remark: Error - if write [IX0,IX] - \"TypeError: Only one indexing vector or array is currently allowed for fancy indexing\"\n# so we proceed as above - in two steps\nprint( len(v1), type(v1), v1.shape) \n\nIX = np.where(pd.Series(list_genes_ensembl_ids ).isin(G2M_genes_Tirosh) >0) [0]\nprint( len(IX), len( G2M_genes_Tirosh) )\nv2 = (f2[first_key]['block0_values'][:,IX]).sum(axis = 1 )\nprint( len(v2), type(v2), v2.shape) \n\n\nmask = (df_cell_hue['day'] == 4) & (df_cell_hue['cell_type'] == 'MkP' )                      \nmask.sum()\ncell_ids = list(df_cell_hue['cell_id'][mask])\nprint(len(cell_ids) , cell_ids[:10])\n\nIX0 = np.where( pd.Series(list_cells_barcodes).isin(cell_ids) > 0)[0]\nprint(len(IX0), IX0[:10], len(cell_ids ))\n\nn_x_subplots = 6\npalette = 'rainbow'#'viridis'\n\nc = 0\nfor gene_id in genes: # ['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.xlabel('G1/S score', fontsize = 20 )\n        plt.ylabel('G2/M score', fontsize = 20 )\n        \n    c += 1\n    fig.add_subplot(1,n_x_subplots ,c)    \n    \n    v = df_cell_hue[gene_id.decode('utf-8')]\n\n    ax = sns.scatterplot(x=v1[IX0], y=v2[IX0], hue=v[IX0], 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 '+str(gene_id), fontsize = 10)\n\n\n#plt.figure(figsize = (20,10))\n#sns.scatterplot(x = v1[IX0], y=v2[IX0] )\n#plt.xlabel('G1/S score', fontsize = 20 )\n#plt.ylabel('G2/M score', fontsize = 20 )\nplt.title(str_data_inf + 'Day 4, MkP (Megakaryocyte Progenitors) ', fontsize = 20)\n\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2022-10-17T19:08:04.984569Z","iopub.execute_input":"2022-10-17T19:08:04.985290Z","iopub.status.idle":"2022-10-17T19:08:31.867450Z","shell.execute_reply.started":"2022-10-17T19:08:04.985247Z","shell.execute_reply":"2022-10-17T19:08:31.866273Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# G1/S-G2/M cell proliferation cycle plots for all days for MkP (Megakaryocyte Progenitors)\n\nCompare with: https://www.kaggle.com/code/alexandervc/mmscel-cell-cycle-03-daybydaychange-megakaryocyte?scriptVersionId=105307517&cellId=12\n\n","metadata":{}},{"cell_type":"code","source":"for day in [2,3,4,7]:\n    \n    mask = (df_cell['day'] == day) & (df_cell['cell_type'] == 'MkP' )                      \n    print('Day',day, mask.sum() )\n    cell_ids = list(df_cell['cell_id'][mask])\n    print(len(cell_ids) , cell_ids[:10])\n\n    IX0 = np.where( pd.Series(list_cells_barcodes).isin(cell_ids) > 0)[0]\n    print(len(IX0), IX0[:10], len(cell_ids ))\n\n    plt.figure(figsize = (20,10))\n    sns.scatterplot(x = v1[IX0], y=v2[IX0] )\n    plt.xlabel('G1/S score', fontsize = 20 )\n    plt.ylabel('G2/M score', fontsize = 20 )\n    plt.title(str_data_inf + '  Day '+str(day) + ', MkP (Megakaryocyte Progenitors)  ', fontsize = 20)\n\n    plt.show()\n","metadata":{"execution":{"iopub.status.busy":"2022-10-17T19:08:31.869125Z","iopub.execute_input":"2022-10-17T19:08:31.869747Z","iopub.status.idle":"2022-10-17T19:08:33.143409Z","shell.execute_reply.started":"2022-10-17T19:08:31.869696Z","shell.execute_reply":"2022-10-17T19:08:33.141990Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Loop over all cell types ","metadata":{}},{"cell_type":"code","source":"dict_cell_types = { '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":{"execution":{"iopub.status.busy":"2022-10-17T19:08:33.144817Z","iopub.execute_input":"2022-10-17T19:08:33.145190Z","iopub.status.idle":"2022-10-17T19:08:33.151715Z","shell.execute_reply.started":"2022-10-17T19:08:33.145144Z","shell.execute_reply":"2022-10-17T19:08:33.150722Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_cell.head(1)\nd = pd.DataFrame(index =list_cells_barcodes,  )\nd = d.join(df_cell.set_index('cell_id'), how = 'left')\nd['str_donor'] = d['donor'].apply(lambda x: str(x))\nd","metadata":{"execution":{"iopub.status.busy":"2022-10-17T19:08:33.153209Z","iopub.execute_input":"2022-10-17T19:08:33.153639Z","iopub.status.idle":"2022-10-17T19:08:33.283552Z","shell.execute_reply.started":"2022-10-17T19:08:33.153599Z","shell.execute_reply":"2022-10-17T19:08:33.282449Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for cell_type in ['MkP','HSC',  'EryP', 'NeuP', 'MasP' ,  'MoP' , 'BP' ]:\n    print();print();print()\n    print(cell_type)\n    \n    for day in [2,3,4,7]:\n\n        mask = (df_cell['day'] == day) & (df_cell['cell_type'] == cell_type )                      \n        #print('Day',day, mask.sum() )\n        cell_ids = list(df_cell['cell_id'][mask])\n        #print(len(cell_ids) , cell_ids[:10])\n\n        IX0 = np.where( pd.Series(list_cells_barcodes).isin(cell_ids) > 0)[0]\n        print('Day', day,mask.sum(), len(IX0))# , IX0[:10]) # , len(cell_ids ))\n\n        plt.figure(figsize = (20,10))\n        if d['str_donor'].iloc[IX0].nunique() > 1:\n            sns.scatterplot(x = v1[IX0], y=v2[IX0] , hue = d['str_donor'].iloc[IX0], palette = 'rainbow')#  ['blue'])\n        else:\n            sns.scatterplot(x = v1[IX0], y=v2[IX0] , hue = d['str_donor'].iloc[IX0], palette =  ['blue'])\n            \n        plt.xlabel('G1/S score', fontsize = 20 )\n        plt.ylabel('G2/M score', fontsize = 20 )\n        plt.title(str_data_inf + '  Day '+str(day) + ', '+ cell_type +  ' ' + dict_cell_types[cell_type] + ' n_cells=%d'%len(IX0) , fontsize = 20)\n\n        plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-17T19:08:33.285371Z","iopub.execute_input":"2022-10-17T19:08:33.286067Z","iopub.status.idle":"2022-10-17T19:08:44.560688Z","shell.execute_reply.started":"2022-10-17T19:08:33.286023Z","shell.execute_reply":"2022-10-17T19:08:44.559635Z"},"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) )\n","metadata":{"execution":{"iopub.status.busy":"2022-10-17T19:08:44.564513Z","iopub.execute_input":"2022-10-17T19:08:44.564830Z","iopub.status.idle":"2022-10-17T19:08:44.570616Z","shell.execute_reply.started":"2022-10-17T19:08:44.564796Z","shell.execute_reply":"2022-10-17T19:08:44.569468Z"},"trusted":true},"execution_count":null,"outputs":[]}]}