{"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\nHere is analysis of the cell cycle based on the single cell RNA sequencing data from Kaggle competition: Multimodal Single-Cell Integration Competition. (Bone marrow cells - mostly progentors. More precisely: mobilized peripheral CD34+ hematopoietic stem and progenitor cells (HSPCs) isolated from four healthy human donors). \n\n\n(с) A.Chervov\n    \nSee paper: [https://arxiv.org/abs/2208.05229](https://arxiv.org/abs/2208.05229) for cell cycle analysis.\n\"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:\n    [https://www.kaggle.com/alexandervc/datasets?scroll=true](https://www.kaggle.com/alexandervc/datasets?scroll=true)  [https://www.kaggle.com/andreizinovyev](https://www.kaggle.com/andreizinovyev)    Hundreds of single cell RNA seq datasets were analyzed. \n    \n    \n#### Remarks:\n\nSee previous also (similar analysis, but only for one cell type - Megakaryocyte ) : \nhttps://www.kaggle.com/alexandervc/mmscel-cell-cycle-03-daybydaychange-megakaryocyte\n\nSee also discussion topic: https://www.kaggle.com/competitions/open-problems-multimodal/discussion/350314\n\nSlides (to be updated): https://docs.google.com/presentation/d/1XxjNe9YJagZ9UayVod2WE_hxIktyeNd6sigOb-TgRbM/edit?usp=sharing\n\n\n#### Conclusions: \n\n    In details - see below for each cell type separately - correspoding to different versions of the script.\n    \n    Overall:\n    \n    Day 2 - most active proliferation - depending on cell type it is classical \"fast\" cell cycle pattern, or more often \"triangle\" but shifted to right (not the standard triangle pattern with G0 part). \n    Only HSC seems to have G0 part for Day 2 other cell types more actively proliferating.\n    Days 3,4 - proliferation activity decreases \n    Day 7 - for some cell types even \"standard\" triangle appears  \n    \n    Seems not batch effect from donors and SigFast plots are quite nice. \n\n    Strange things:\n    Day 4 sometimes seems more active proliferating than day 3.\n    For HSC and Megakaryocyte (and may be others): \n    For day 2 - we see \"triangle\" pattern, but  when we select \"fast\" region - then PCA shows another \"cycle\" !\n    So what do we have ? Such a strange \"double\" cell cycle pattern with \"fast\" and \"very fast\" patterns? \n    Need to investigate further. \n    \n    \n#### Conclusions from previous scripts: \n\nhttps://www.kaggle.com/alexandervc/mmscel-cell-cycle-03-daybydaychange-megakaryocyte/edit\n\n    Megakaryicyte: \n    The cell cycle pattern is changing day by day. \n    Day 2 - active proliferation, no G0 part. \n    Day 3 - G0 part appears\n    Day 4 - G0 part almost DISappears (that is strange, may be incorrect cell type identification )\n    Day 7 - \"standard\"(big triangle with G0) pattern appear, G0 appear, but \"fast\"-like part is still present\n    \n    By the previous analysis scripts  it became clear that for better analysis of the cell cycle for the current dataset \n    we need to 1) split by cell types 2) make denoisings \n    Script (denoising analysis) 2: https://www.kaggle.com/code/alexandervc/mmscel-cell-cycle-02-with-denoising\n    Script (various splits by cell types, donors, days,  etc) 1: https://www.kaggle.com/code/alexandervc/mmscel-cell-cycle-01-preliminary-draft\n\n\nCells are \"mobilized peripheral CD34+ \" - here \"MOBILIZED\" seems to mean that donors were given special drugs for weeks before the experiment such that their bone marrow cells move to the blood. It seems it is usually achieved by the proliferation stimulation. That is why we may expect stronger proliferation - so it might be different from other datasets with similar cells (typically there are datasets WITHOUT\" mobilization\").\n\nSo active proliferation at the begining and fading later is consitent with the information above.\nWhat is strange that on day 3 - there is more fading of prolifetation , than day 4. But it might be technical reason  due to incorrect cell type identification. \n\n\n#### Technicality\n\n    We use preprocessed \"denoised\" scRNA-seq data - we use our own \"pooling\" (simple and fast) with n_neigbour = 20\n    The preprocessing script: https://www.kaggle.com/competitions/open-problems-multimodal\n    It requires huge RAM so the Saturn cloud was used to run it. \n    \n    The huge h5ad would take 12G out 16 Kaggle memory - so can do nothing in downstream, \n    So we back it on disk, but extract part corresponding the particular cell type\n\n### Versions:\n\n#### 8: cosmetic changes\n\n#### 7:  'HSC'  42722 'Hematoploetic Stem Cell'\n\n    Day 2: in contrast to many other cell types G0-like part is present even at day 2.\n    Day 3,4,7: somewhat similar and similar to day 2. Some strange  thing that day 4 seems more active proliferation than day 3, but not very clear. \n\n    The same strange thing as for Megakaryocyte:\n    For day 2 - we see \"triangle\" pattern, but  when we select \"fast\" region - then PCA shows another \"cycle\" !\n    So what do we have ? Such a strange \"double\" cell cycle pattern with \"fast\" and \"very fast\" patterns? \n    \n    As always: seems not batch effect from donors and SigFast plots are quite nice. \n\n#### 6: 'MkP' 10740 'Megakaryocyte Progenitor'\n\n    See also https://www.kaggle.com/alexandervc/mmscel-cell-cycle-03-daybydaychange-megakaryocyte\n\n    Previous script: Megakaryicyte \n    The cell cycle pattern is changing day by day. \n    Day 2 - active proliferation, no G0 part. \n    Day 3 - G0 part appears\n    Day 4 - G0 part almost DISappears (that is strange, may be incorrect cell type identification )\n    Day 7 - \"standard\"(big triangle with G0) pattern appear, G0 appear, but \"fast\"-like part is still present\n    \n    Here we also observe strange thing - need to understand further.\n    For day 2 - we see \"triangle\" pattern, but  when we select \"fast\" region - then PCA shows another \"cycle\" !\n    So what do we have ? Such a strange \"double\" cell cycle pattern with \"fast\" and \"very fast\" patterns? \n    Or it is some mistake in analysis ? \n\n\n#### 5: MoP      1796 'Monocyte Progenitor'\n\n    Day 2: not many cells - not clear what is going on\n    Close look at day 4 seems to reveal double cell cycle - fast + standard\n    Day 3,4,7 - more probably looks like standard cell cycle - (different from other cell types)\n        (even though we cannot fully exclude small fraction of \"fast\" cell cycle\n        \n\n#### 4: EryP    24076     'Erythrocyte Progenitor'\n\n    Day 2: \"fast\" cell cycle \n    Days 3,4,7 - probably double cell cycle with \"fast\" and \"standard\" mixed \n    \n    Again seems NO donor dependence\n    \n#### 3: B-cell Progenitors \n\n    The lowest number of cells - 270 - it might be just mistake of the cell type labeling software.\n    (All the other cells are myeloid type, only these lymphoid).\n    \n    \"Fast\" cell cycle part and part around G0. Seems no much dependence on day(!) or other features. \n\n#### 2: NeuP = Neutrophil Progenitor\n    \n    At that script we look on 'MasP': 'Mast Cell Progenitor' - 21177 cells. 4.4G RAM consumed\n    \n    Again - day by day activity of proliferation is going down.\n    Day 2 - \"fast\" cell cycle pattern, quite clear, almost all cells - \"fast\" \n    Day 3,4 - part of cells -  \"fast\" cell cycle pattern, but other part - seems in G0, or \"standard\"\n    Day 7 - seems \"fast\" disappear - most \"standard\" and G0\n    \n    SigFastScore also shows quite nice figures - the cyclic structure can be seen in all cases. \n    \n    Seems no batch effect from donor. \n    \n#### 1 'MasP': 'Mast Cell Progenitor'\n    Version 1 -  look on   'MasP': 'Mast Cell Progenitor' -   17899 cells. \n    We are afrain of  RAM crash - so do not load all cell types at once\n    \n    Mainly \"fast\" cell cycle pattern. However for day 7 - it is quite different. \n    The days in between: 3,4 - mainly \"fast\", but some fraction of not so active cells.\n    So probably days 3,4 - double cell cycle structure. \n    The \"cyclic\" structure can be seen only making PCA from Tirosh genes, (and addtional filtering G1S+G2M > 8 - for days greater than 2 ) \n    \n    The quite striking - using G1/S-G2/M plots it is impossible to believe that there is \"cyclic\" structure at the \"fast\" part, but when we mask it and make PCA - we clearly see it. (As we usually do - see our paper). \n    \n    SigFastScore also shows quite nice figures - the cyclic structure can be seen in all cases. \n","metadata":{}},{"cell_type":"markdown","source":"# Key param","metadata":{}},{"cell_type":"code","source":"cell_type2consider = 'HSC'# 'MkP' # 'MoP' # 'HSC'#  'EryP' #  'BP'# 'NeuP' # Only one cell type will be considered for each version of the notebook \n# 'MasP' # 'MkP'#  'MoP'\n# 'MkP': 'Megakaryocyte Progenitor'\n","metadata":{"execution":{"iopub.status.busy":"2022-09-11T06:49:30.553689Z","iopub.execute_input":"2022-09-11T06:49:30.554156Z","iopub.status.idle":"2022-09-11T06:49:30.559159Z","shell.execute_reply.started":"2022-09-11T06:49:30.554121Z","shell.execute_reply":"2022-09-11T06:49:30.558098Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Install/import and auxilliary functions ","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":{"execution":{"iopub.status.busy":"2022-09-11T06:49:30.565384Z","iopub.execute_input":"2022-09-11T06:49:30.566145Z","iopub.status.idle":"2022-09-11T06:49:30.577020Z","shell.execute_reply.started":"2022-09-11T06:49:30.566102Z","shell.execute_reply":"2022-09-11T06:49:30.575672Z"},"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\ndict_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'}\n\n\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 = ['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\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\n\ndef phase_plot(adata, list_genes1, list_genes2, mask = None, list_color_by = [None], str_suptitle = None, \n               str_title = None,\n               xfigsize = 20, yfigsize = 10,\n               xlabel = 'G1S score', ylabel = 'G2M score',\n               n_x_subplots = 1, n_y_subplots = 1 #  #len(list_selected_cell_lines)\n              ):\n\n    if mask is None:\n        mask = np.ones( adata.X.shape[0]).astype(bool) #  \n\n    #list_color_by = ['total_counts']# , 'Cell cycle phase']\n\n    \n    total_count = 0\n    for color_by_mode in  list_color_by: # , 'pct_counts_mt']:\n        \n        if n_x_subplots * n_y_subplots <= 1:\n            fig = plt.figure(figsize = (xfigsize,yfigsize) ); c = 0\n            if str_suptitle is not None:\n                plt.suptitle(str_suptitle + ' N_cells= ' +str(mask.sum() ) + ' ' , fontsize = 20 )#' \n        else:\n            if ( total_count % (n_x_subplots * n_y_subplots) ) == 0: \n                fig = plt.figure(figsize = (xfigsize,yfigsize) ); c = 0\n                if str_suptitle is not None:\n                    plt.suptitle(str_suptitle + ' N_cells= ' +str(mask.sum() ) + ' ' , fontsize = 20 )#' \n            c+=1; fig.add_subplot(n_y_subplots,n_x_subplots,c);\n            \n        total_count += 1\n        \n        \n        color_by_field_name = color_by_mode\n        color_by = None\n        if color_by_mode is None:\n            pass\n        elif color_by_mode in adata.var.index:\n            \n            color_by = pd.Series( np.asarray(adata[mask,color_by_mode].X.sum(axis=1) ).ravel() , name = color_by_mode )\n        elif color_by_mode in adata.obs:\n            color_by = (adata[mask].obs[color_by_field_name])\n        elif 'median_binarize_' in color_by_mode:\n            color_by_field_name = color_by_mode[16:]\n            if color_by_field_name in adata.obs:\n                color_by = (adata[mask].obs[color_by_field_name]) > np.median( (adata[mask].obs[color_by_field_name]) )\n        elif 'threshold_binarize_' in color_by_mode:\n            color_by_field_name = color_by_mode.split('_')[3]\n            threshold_binarize = float( color_by_mode.split('_')[2] )\n            if color_by_field_name in adata.obs:\n                color_by = (adata[mask].obs[color_by_field_name]) > threshold_binarize\n        else:\n            color_by = None\n\n        ll = list_genes1\n        I = np.where(adata[mask].var.index.isin(ll) > 0 )[0] # .sum()\n        v1 = np.asarray(adata[mask].X[:,I].mean(axis = 1)).ravel()\n        #print(len(I), v1.mean(), len(v1))    \n        \n        ll = list_genes2#  G2M_genes_Tirosh\n        I = np.where(adata[mask].var.index.isin(ll) > 0 )[0] # .sum()\n        #print(len(I))\n        v2 = np.asarray(adata[mask].X[:,I].mean(axis = 1)).ravel()    \n\n\n        if str_title is not None:\n            plt.title( str_title  , fontsize = 20 )\n        if color_by is None:\n            ax = sns.scatterplot(x=v1, y = v2)# ,  hue= color_by,   alpha = 0.8, marker = '.')#, legend=None)\n        else:\n            #color_by = (adata.obs[color_by_field_name]) > np.median( adata.obs[color_by_field_name].values ) \n            ax = sns.scatterplot(x=v1, y = v2,  hue= color_by, palette = \"rainbow\")# 'viridis\")# sns.color_palette(\"viridis\", as_cmap=True),\n                                #)# ,   alpha = 0.8, marker = '.')#, )#, legend=None)\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\n\n        #plt.title(str_inf)# str_data_inf + ' ' + str_reducer + str_genes_inf + ' ' + col)\n        plt.xlabel(xlabel , fontsize = 20)\n        plt.ylabel(ylabel, fontsize = 20 )\n    plt.show()         \n        ","metadata":{"execution":{"iopub.status.busy":"2022-09-11T06:49:30.645858Z","iopub.execute_input":"2022-09-11T06:49:30.646270Z","iopub.status.idle":"2022-09-11T06:49:43.359444Z","shell.execute_reply.started":"2022-09-11T06:49:30.646235Z","shell.execute_reply":"2022-09-11T06:49:43.357989Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load data and extract only one cell type","metadata":{}},{"cell_type":"code","source":"#fn = '/kaggle/input/scrnaseq-tabula-sapiens-human-part-2/TS_epithelial.h5ad'\nfn = '/kaggle/input/multimodal-singlecell-integration-related-data-01/cite_seq_inputs.h5ad'\nfn = '/kaggle/input/multimodal-singlecell-integration-related-data-01/cite_seq_inputs_denoising_pooling_n_neighbours20.h5ad'\nstr_data_inf = ' CITEseq '\nprint(str_data_inf, fn )\nt0 = time.time()\nadata_full = sc.read(fn,  backed='r' )\nprint(np.round(time.time()-t0,0),'second passed on loading')\nadata_full# ","metadata":{"execution":{"iopub.status.busy":"2022-09-11T06:49:43.362439Z","iopub.execute_input":"2022-09-11T06:49:43.363321Z","iopub.status.idle":"2022-09-11T06:49:44.331966Z","shell.execute_reply.started":"2022-09-11T06:49:43.363276Z","shell.execute_reply":"2022-09-11T06:49:44.330771Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('820 M RAM consumed - quite little')","metadata":{"execution":{"iopub.status.busy":"2022-09-11T06:49:44.333339Z","iopub.execute_input":"2022-09-11T06:49:44.333814Z","iopub.status.idle":"2022-09-11T06:49:44.339885Z","shell.execute_reply.started":"2022-09-11T06:49:44.333779Z","shell.execute_reply":"2022-09-11T06:49:44.338761Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print( adata_full.obs['cell_type'].value_counts() )\nprint( adata_full.obs['Train0Test1'].value_counts() )\nprint( adata_full.obs['day'].value_counts() )\nprint( adata_full.obs['donor'].value_counts() )\nprint( adata_full.obs['Train0Test1'].value_counts() )","metadata":{"execution":{"iopub.status.busy":"2022-09-11T06:49:44.341682Z","iopub.execute_input":"2022-09-11T06:49:44.342359Z","iopub.status.idle":"2022-09-11T06:49:44.363862Z","shell.execute_reply.started":"2022-09-11T06:49:44.342314Z","shell.execute_reply":"2022-09-11T06:49:44.362662Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Select one cell type part from the data","metadata":{}},{"cell_type":"code","source":"import time\nt0 = time.time()\ncol = 'cell_type'\nuv = cell_type2consider # 'NeuP' # 'MasP' # 'MkP'#  'MoP'\nmask = adata_full.obs[col] == uv\nprint(mask.sum())\n\nadata = adata_full[mask,:].to_memory()# Otherwise cannot change like: a.uns['Cell type'] = cell_type_loc\nprint(('%.f seconds passsed'%(time.time()-t0)))\nstr_data_inf = ' CITEseq  ' + uv + ' ' +dict_cell_types[uv]   + ' n_neighbours='+str(20)#n_neighbours)\nprint( str_data_inf )\nadata","metadata":{"execution":{"iopub.status.busy":"2022-09-11T06:49:44.366898Z","iopub.execute_input":"2022-09-11T06:49:44.367244Z","iopub.status.idle":"2022-09-11T06:53:04.850212Z","shell.execute_reply.started":"2022-09-11T06:49:44.367207Z","shell.execute_reply":"2022-09-11T06:53:04.848858Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"adata.obs['str_donor'] = adata.obs['donor'].apply(lambda x: str(x))","metadata":{"execution":{"iopub.status.busy":"2022-09-11T06:53:04.852179Z","iopub.execute_input":"2022-09-11T06:53:04.852550Z","iopub.status.idle":"2022-09-11T06:53:04.877152Z","shell.execute_reply.started":"2022-09-11T06:53:04.852516Z","shell.execute_reply":"2022-09-11T06:53:04.875907Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"adata.obs['day'].value_counts()","metadata":{"execution":{"iopub.status.busy":"2022-09-11T06:53:04.878861Z","iopub.execute_input":"2022-09-11T06:53:04.879216Z","iopub.status.idle":"2022-09-11T06:53:04.888252Z","shell.execute_reply.started":"2022-09-11T06:53:04.879182Z","shell.execute_reply":"2022-09-11T06:53:04.887352Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"adata.obs.groupby(['day','donor'])['day'].count().to_frame()","metadata":{"execution":{"iopub.status.busy":"2022-09-11T06:53:04.889315Z","iopub.execute_input":"2022-09-11T06:53:04.889681Z","iopub.status.idle":"2022-09-11T06:53:04.913185Z","shell.execute_reply.started":"2022-09-11T06:53:04.889648Z","shell.execute_reply":"2022-09-11T06:53:04.911986Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Day by day big figures","metadata":{}},{"cell_type":"code","source":"col = 'day'\nfor uv in adata.obs[col].unique():\n    mask = adata.obs[col] == uv\n    print(mask.sum())\n    phase_plot(adata, list_genes1 =G1S_genes_Tirosh , list_genes2=G2M_genes_Tirosh, \n               mask = mask,\n               str_title =  ' G1/S-G2/M plot.   Tirosh genes ' + col + ' = ' + str(uv)  ,\n               str_suptitle = str_data_inf, #  + ' n_cells = ',# + str( adata.shape[0] ),\n                n_x_subplots = 1, n_y_subplots = 1, \n          list_color_by = [  'str_donor' ],) # 'CCNB1','day',  'n_genes_by_counts', 'total_counts', ,'CCNE1'","metadata":{"execution":{"iopub.status.busy":"2022-09-11T06:53:04.914876Z","iopub.execute_input":"2022-09-11T06:53:04.915242Z","iopub.status.idle":"2022-09-11T06:53:08.658659Z","shell.execute_reply.started":"2022-09-11T06:53:04.915204Z","shell.execute_reply":"2022-09-11T06:53:08.657241Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# PCA from Tirosh Genes ","metadata":{}},{"cell_type":"code","source":"col = 'day'\nfor uv in adata.obs[col].unique():\n    mask = adata.obs[col] == uv\n    \n    print(mask.sum())\n    from sklearn.decomposition import PCA\n    pca = PCA(n_components=5)\n\n    ll = G1S_genes_Tirosh + G2M_genes_Tirosh\n    I = np.where(pd.Series(adata.var.index).isin(ll) > 0 )[0] # .sum()\n    X = adata[mask].X[:,I]# df[mask].iloc[:,I]# np.asarray(df.iloc[:,I].mean(axis = 1)).ravel()\n    print(X.shape, type(X))      #      str(col) + ' = ' + str(uv) + \n    r = pca.fit_transform(X)\n\n    fig = plt.figure(figsize = (20,10) ); c = 0\n\n    ax = sns.scatterplot(x=r[:,0], y = r[:,1] , hue = adata[mask].obs['str_donor'])# pd.Series(v1[mask], name = 'G1S score'), palette = 'rainbow' )#,  hue= color_by, palette = \"rainbow\")# 'viridis\")# sns.color_palette(\"viridis\", as_cmap=True),\n    #ax = sns.scatterplot(x=r[:,0], y = r[:,1])# , hue = pd.Series(v1[mask], name = 'G1S score'), palette = 'rainbow' )#,  hue= color_by, palette = \"rainbow\")# 'viridis\")# sns.color_palette(\"viridis\", as_cmap=True),\n                    #)# ,   alpha = 0.8, marker = '.')#, )#, legend=None)\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\n\n    #plt.title(str_inf)# str_data_inf + ' ' + str_reducer + str_genes_inf + ' ' + col)\n#             plt.title(str_data_inf + '  ' + col1+ ' = ' + str(uv1) +'  '+col3 + ' = '+str(uv3) + ' ' + col2 +' = ' + uv2 + ' ' + str(dict_cell_types[uv2]), fontsize = 30)\n\n    #str_title =  str_inf + ' G1/S-G2/M plot.   Tirosh genes' \n    str_suptitle = str_data_inf # + ' n_cells = ' + str( adata.shape[0] ),\n    plt.suptitle(str_data_inf , fontsize = 20)# str_title)\n    plt.title(str(col) + ' = ' + str(uv) + ' PCA from Tirosh cell cycle genes ' , fontsize = 20 )\n    plt.xlabel('PCA1' , fontsize = 20)\n    plt.ylabel('PCA2', fontsize = 20 )\n    plt.show()             \n","metadata":{"execution":{"iopub.status.busy":"2022-09-11T06:53:08.660426Z","iopub.execute_input":"2022-09-11T06:53:08.660792Z","iopub.status.idle":"2022-09-11T06:53:12.715866Z","shell.execute_reply.started":"2022-09-11T06:53:08.660757Z","shell.execute_reply":"2022-09-11T06:53:12.715007Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# PCA from cells with which are in potential \"fast\" part (i.e. v1+v2>threshold)\n\nThe way to see \"cyclic\" structure for the \"fast\" cell cycle - it is indeed can be seen ! \n","metadata":{}},{"cell_type":"code","source":"\ncol = 'day'\nfor uv in adata.obs[col].unique():\n    \n    mask = np.ones(len(adata), dtype = bool )\n    ll = G1S_genes_Tirosh\n    I = np.where(adata[mask].var.index.isin(ll) > 0 )[0] # .sum()\n    v1 = np.asarray(adata[mask].X[:,I].mean(axis = 1)).ravel()\n    #print(len(I), v1.mean(), len(v1))    \n\n    ll = G2M_genes_Tirosh # list_genes2#  G2M_genes_Tirosh\n    I = np.where(adata[mask].var.index.isin(ll) > 0 )[0] # .sum()\n    #print(len(I))\n    v2 = np.asarray(adata[mask].X[:,I].mean(axis = 1)).ravel()    \n    \n    mask0 = adata.obs[col] == uv\n\n    mask = mask0 & ( (v1 + v2) > 7  ) \n    \n    print(mask.sum())\n    from sklearn.decomposition import PCA\n    pca = PCA(n_components=5)\n\n    ll = G1S_genes_Tirosh + G2M_genes_Tirosh\n    I = np.where(pd.Series(adata.var.index).isin(ll) > 0 )[0] # .sum()\n    X = adata[mask].X[:,I]# df[mask].iloc[:,I]# np.asarray(df.iloc[:,I].mean(axis = 1)).ravel()\n    print(X.shape, type(X))            \n    r = pca.fit_transform(X)\n\n    fig = plt.figure(figsize = (20,10) ); c = 0\n    ax = sns.scatterplot(x=v1[mask0], y = v2[mask0], hue = ( (v1 + v2) > 8  )[mask0] )#  , hue = adata[mask].obs['str_donor'])# pd.Series(v1[mask], name = 'G1S score'), palette = 'rainbow' )#,  hue= color_by, palette = \"rainbow\")# 'viridis\")# sns.color_palette(\"viridis\", as_cmap=True),\n    plt.xlabel('G1/S score', fontsize = 20)\n    plt.ylabel('G2/M score', fontsize = 20)\n    plt.title(str(col) + ' = ' + str(uv) + ' Highlight the area which will be preserved on the next plot. Other part will be dropped off.', fontsize = 20 )\n    \n    str_suptitle = str_data_inf # + ' n_cells = ' + str( adata.shape[0] ),\n    plt.suptitle(str_data_inf , fontsize = 20)# str_title)\n    \n    plt.show()\n\n    \n    fig = plt.figure(figsize = (20,10) ); c = 0\n\n    ax = sns.scatterplot(x=r[:,0], y = r[:,1] , hue = adata[mask].obs['str_donor'])# pd.Series(v1[mask], name = 'G1S score'), palette = 'rainbow' )#,  hue= color_by, palette = \"rainbow\")# 'viridis\")# sns.color_palette(\"viridis\", as_cmap=True),\n    #ax = sns.scatterplot(x=r[:,0], y = r[:,1])# , hue = pd.Series(v1[mask], name = 'G1S score'), palette = 'rainbow' )#,  hue= color_by, palette = \"rainbow\")# 'viridis\")# sns.color_palette(\"viridis\", as_cmap=True),\n                    #)# ,   alpha = 0.8, marker = '.')#, )#, legend=None)\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\n\n    #plt.title(str_inf)# str_data_inf + ' ' + str_reducer + str_genes_inf + ' ' + col)\n#             plt.title(str_data_inf + '  ' + col1+ ' = ' + str(uv1) +'  '+col3 + ' = '+str(uv3) + ' ' + col2 +' = ' + uv2 + ' ' + str(dict_cell_types[uv2]), fontsize = 30)\n\n    #str_title =  str_inf + ' G1/S-G2/M plot.   Tirosh genes' \n    str_suptitle = str_data_inf # + ' n_cells = ' + str( adata.shape[0] ),\n    plt.suptitle(str_data_inf , fontsize = 20)# str_title)\n    plt.title(str(col) + ' = ' + str(uv) + ' PCA from Tirosh cell cycle genes ' , fontsize = 20 )\n    plt.xlabel('PCA1' , fontsize = 20)\n    plt.ylabel('PCA2', fontsize = 20 )\n    plt.show()             \n        ","metadata":{"execution":{"iopub.status.busy":"2022-09-11T06:53:12.717417Z","iopub.execute_input":"2022-09-11T06:53:12.717969Z","iopub.status.idle":"2022-09-11T06:53:23.520833Z","shell.execute_reply.started":"2022-09-11T06:53:12.717931Z","shell.execute_reply":"2022-09-11T06:53:23.519644Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"col = 'day'\nfor uv in adata.obs[col].unique():\n    mask = adata.obs[col] == uv\n    print(mask.sum())\n    phase_plot(adata, list_genes1 =G1S_genes_Tirosh , list_genes2=G2M_genes_Tirosh, \n               mask = mask,\n               str_title =  ' G1/S-G2/M plot.   Tirosh genes ' + col + ' = ' + str(uv)  ,\n               str_suptitle = str_data_inf, #  + ' n_cells = ',# + str( adata.shape[0] ),\n                n_x_subplots = 1, n_y_subplots = 1, \n          list_color_by = [  'str_donor' ],) # 'CCNB1','day',  'n_genes_by_counts', 'total_counts', ,'CCNE1'","metadata":{"execution":{"iopub.status.busy":"2022-09-11T06:53:23.522383Z","iopub.execute_input":"2022-09-11T06:53:23.522776Z","iopub.status.idle":"2022-09-11T06:53:27.236280Z","shell.execute_reply.started":"2022-09-11T06:53:23.522741Z","shell.execute_reply":"2022-09-11T06:53:27.235171Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#  G1/S-G2/M plots to analyse the cell cycle","metadata":{}},{"cell_type":"code","source":"mask = np.ones(len(adata))\nphase_plot(adata, list_genes1 =G1S_genes_Tirosh , list_genes2=G2M_genes_Tirosh, \n           #mask = mask,\n           str_title =  ' G1/S-G2/M plot.   Tirosh genes' ,\n           str_suptitle = str_data_inf , # + ' n_cells = ',# + str( adata.shape[0] ),\n      list_color_by = [ 'CCNB1','day',  'n_genes_by_counts', 'total_counts', 'donor', 'Train0Test1' ],) # ,'CCNE1'","metadata":{"execution":{"iopub.status.busy":"2022-09-11T06:53:27.237892Z","iopub.execute_input":"2022-09-11T06:53:27.239007Z","iopub.status.idle":"2022-09-11T06:53:43.837236Z","shell.execute_reply.started":"2022-09-11T06:53:27.238953Z","shell.execute_reply":"2022-09-11T06:53:43.836039Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# SigFast Plot for all cells ","metadata":{}},{"cell_type":"code","source":"mask = np.ones(len(adata))\nphase_plot(adata, list_genes1 =list_genes_fastCCsign , list_genes2=G2M_genes_Tirosh, \n           #mask = mask, \n           xlabel = 'SigFast score',\n           str_title =  ' SigFast-G2/M plot.   Tirosh genes ' + col + ' = ' + str(uv)  ,\n           str_suptitle = str_data_inf, #  + ' n_cells = ',# + str( adata.shape[0] ),\n            n_x_subplots = 1, n_y_subplots = 1, \n      list_color_by = [ 'CCNB1','day',  'n_genes_by_counts', 'total_counts', 'donor', 'Train0Test1' ]) # ,'CCNE1'","metadata":{"execution":{"iopub.status.busy":"2022-09-11T06:53:43.840888Z","iopub.execute_input":"2022-09-11T06:53:43.841284Z","iopub.status.idle":"2022-09-11T06:54:00.836791Z","shell.execute_reply.started":"2022-09-11T06:53:43.841247Z","shell.execute_reply":"2022-09-11T06:54:00.835559Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Day by day ","metadata":{}},{"cell_type":"code","source":"col = 'day'\nfor uv in adata.obs[col].unique():\n    mask = adata.obs[col] == uv\n    print(mask.sum())\n    phase_plot(adata, list_genes1 =G1S_genes_Tirosh , list_genes2=G2M_genes_Tirosh, \n               mask = mask,\n               str_title =  ' G1/S-G2/M plot.   Tirosh genes ' + col + ' = ' + str(uv)  ,\n               str_suptitle = str_data_inf , # + ' n_cells = ',# + str( adata.shape[0] ),\n                n_x_subplots = 2, n_y_subplots = 1, \n          list_color_by = [  'donor', 'Train0Test1' ],) # 'CCNB1','day',  'n_genes_by_counts', 'total_counts', ,'CCNE1'","metadata":{"execution":{"iopub.status.busy":"2022-09-11T06:54:00.838204Z","iopub.execute_input":"2022-09-11T06:54:00.838550Z","iopub.status.idle":"2022-09-11T06:54:07.460337Z","shell.execute_reply.started":"2022-09-11T06:54:00.838516Z","shell.execute_reply":"2022-09-11T06:54:07.459165Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# SigFast Scores Day by day big figure - quite nice\n","metadata":{}},{"cell_type":"code","source":"# list_genes_fastCCsign = ['CDK1', 'UBE2C', 'TOP2A', 'TMPO', 'HJURP', 'RRM1', 'RAD51AP1', 'RRM2', 'CDC45', 'BLM', 'BRIP1', 'E2F8', 'HIST2H2AC']\n\ncol = 'day'\nfor uv in adata.obs[col].unique():\n    mask = adata.obs[col] == uv\n    print(mask.sum())\n    phase_plot(adata, list_genes1 =list_genes_fastCCsign , list_genes2=G2M_genes_Tirosh, \n               mask = mask, xlabel = 'SigFast score',\n               str_title =  ' SigFast-G2/M plot.   Tirosh genes ' + col + ' = ' + str(uv)  ,\n               str_suptitle = str_data_inf, #  + ' n_cells = ',# + str( adata.shape[0] ),\n                n_x_subplots = 1, n_y_subplots = 1, \n          list_color_by = [  'donor' ],) # 'CCNB1','day',  'n_genes_by_counts', 'total_counts', ,'CCNE1'","metadata":{"execution":{"iopub.status.busy":"2022-09-11T06:54:07.462047Z","iopub.execute_input":"2022-09-11T06:54:07.462389Z","iopub.status.idle":"2022-09-11T06:54:11.157895Z","shell.execute_reply.started":"2022-09-11T06:54:07.462355Z","shell.execute_reply":"2022-09-11T06:54:11.156712Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}