{"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\nCollect some information on genes  for RNA-seq part of data. And pack it into the pandas dataframe. \n\nPython package \"mygene\" is used.\n\n**Brief info** Among information attached: gene \"symbol\" - i.e. typically 4 letter gene name - typically used by human experts,\nchromosome where gene is placed, the position on chromosome, strand. \nAlso:  n_pubmed - number of publications in pubmbed - some soft measure how gene is \"important\",\nand similar n_pathways, n_GO_BP etc... . Also: type of gene - protein coding, etc., summary, name - some information about the gene, etc. \n\n**Remarks**  Out of around 23K genes there are  7363 which are not really well studied - no pubmed, no GO, no pathways (at least according to \"mygene\"). Other sources may give slighly different numbers - but seems to be the general picture would be the same - there are many low expressed and not well studied genes which we may expect would be diffictult to predict in currect challenge. \n\n**Possible use in competition** - for Multiome task one may expect that only \"nearby\" regions of DNA affect the expression,\n(see discussion: https://www.kaggle.com/competitions/open-problems-multimodal/discussion/350863 )\nthat means - regions on the SAME chromosome, and not more than around 1 million bp far from gene position.\nSo one may try to check that hypothesis using position of the gene given here.\nHowever there can be some subtleties - due to different versions of the genome - and it might \nbe better to use the code by organizers: https://www.kaggle.com/competitions/open-problems-multimodal/discussion/349559,\nor code by https://www.kaggle.com/code/masato114/msci-multiome-using-geneactivity \n\n**\"IM\"possible use in competition** - that hardly will work for score improving, but interesting for research and especially for interpretation.  The idea: DNA regions (ATAC-seq part of data) are influencing gene expressions typically because\nspecial proteins \"transcription factors\" bind at those regions - https://en.wikipedia.org/wiki/Transcription_factor.\nFor some transcription factors there are known \"motifs\" e.g. short sequences of nucleotids where they bind.\nSo essentially  one has some approximate guesses on ATAC-seq part of data where transcription factors bind.\nOn the other hand the there is some knowledge on genes which are \"targets\" of  transcription factors.\nThus using domain knowledge one may build some \"solution\" to the Multiome task - connecting ATAC-seq -> transcription factors -> genes. \nThe problem is that domain knowledge is far from being perfect (motif does not guarantee that tf actually binds, targets depend on cell types and conditions, etc)  and most probably any even simplest ML model will be better than such solution.\nBut what  is interesting - to get some interpretation of ML solutions along these lines. \n\nSo we attach some transciption factor info to the pandas dataframe. TF info taken from: http://humantfs.ccbr.utoronto.ca/\n\n\n\n**Technnicalities using \"mygene\":**\n\nUnfortunate case revealed - mygene may return duplicated list by some unclear reason. \nSo we need to drop duplicates. Sometimes \"ensembl\" field is not returned, despite we query by \"ensembl\" and query seems to work correctly, so we finally substitute \"ensembl\" by \"query\". The query on 23K genes took 118 seconds - not bad. \nIn general mygene and other  similar packages query data from databases via internet and that might not be always stable,\nso one should be careful and check results.\n(For information - query returns much of more information, but we cut it to something more managable). \n\n\n\n","metadata":{}},{"cell_type":"markdown","source":"# Top genes in different categories\n\nSo just for some informal info - top 12 genes in different categories - like length, number of publications (here TP53 is the champion - \"guardian of the genome\" - most mutated genes in cancers https://en.wikipedia.org/wiki/P53 ), and other :\n\n    len\n    ['RBFOX1', 'CNTNAP2', 'PTPRD', 'DMD', 'DLG2', 'CSMD1', 'MACROD2', 'EYS', 'LRP1B', 'CTNNA3', 'PCDH15', 'ROBO2']\n\n    n_pubmed\n    ['TP53', 'TNF', 'APOE', 'VEGFA', 'TGFB1', 'MTHFR', 'IL10', 'HIF1A', 'ESR1', 'IL1B', 'ACE', 'STAT3']\n\n    n_pathways\n    ['MAPK1', 'PLCB1', 'HRAS', 'ALDOA', 'PIP5K1A', 'NME1', 'PLA2G4A', 'LDHA', 'SOD1', 'MAP2K2', 'MAP2K1', 'NOS3']\n\n    n_GO_BP\n    ['TNF', 'BMP4', 'NOTCH1', 'CTNNB1', 'WNT5A', 'TP53', 'AKT1', 'TGFB1', 'BCL2', 'SRC', 'SIRT1', 'PRKN']\n\n    n_GO_CC\n    ['LRRK2', 'FMR1', 'ACTB', 'APP', 'PSEN1', 'CTNNB1', 'ITGB1', 'ANXA2', 'HSPA8', 'EZR', 'MAPT', 'SNRPD3']\n\n    n_GO_MF\n    ['TP53', 'EP300', 'DHX9', 'HNRNPU', 'SMAD3', 'RELA', 'LRRK2', 'SRC', 'PARK7', 'ABL1', 'SIRT1', 'HSP90AB1']\n\n    n_transcripts\n    ['PCBP1-AS1', 'MIR99AHG', 'MAPK10', 'PVT1', 'TEX41', 'FANCL', 'MIR663AHG', 'CCDC26', 'SNHG14', 'ANK2', 'LOC100506207', 'SOX2-OT']\n    \n    \nGO - means Gene Ontology - https://en.wikipedia.org/wiki/Gene_Ontology, BP, CC, MF - three branches of GO: \n\n    CC - cellular component, the parts of a cell or its extracellular environment;\n    MF - molecular function, the elemental activities of a gene product at the molecular level, such as binding or catalysis;\n    BP - biological process, operations or sets of molecular events with a defined beginning and end, pertinent to the functioning of integrated living units: cells, tissues, organs, and organisms.\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-18T09:01:51.54641Z","iopub.execute_input":"2022-09-18T09:01:51.547723Z","iopub.status.idle":"2022-09-18T09:01:51.562072Z","shell.execute_reply.started":"2022-09-18T09:01:51.547672Z","shell.execute_reply":"2022-09-18T09:01:51.561044Z"},"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\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\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-09-18T09:01:51.624235Z","iopub.execute_input":"2022-09-18T09:01:51.625275Z","iopub.status.idle":"2022-09-18T09:02:25.937772Z","shell.execute_reply.started":"2022-09-18T09:01:51.625213Z","shell.execute_reply":"2022-09-18T09:02:25.936117Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"START = 0\nSTOP  = 10\n\ndf_cite_test_x = pd.read_hdf(FP_CITE_TEST_INPUTS, START = START, STOP = STOP)\ndf_cite_test_x.head()","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:02:25.940858Z","iopub.execute_input":"2022-09-18T09:02:25.941542Z","iopub.status.idle":"2022-09-18T09:03:07.008405Z","shell.execute_reply.started":"2022-09-18T09:02:25.941498Z","shell.execute_reply":"2022-09-18T09:03:07.007532Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_multi_train_y = pd.read_hdf(FP_MULTIOME_TRAIN_TARGETS, start=START, stop=STOP)\ndf_multi_train_y.head()","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:03:07.009752Z","iopub.execute_input":"2022-09-18T09:03:07.010642Z","iopub.status.idle":"2022-09-18T09:03:07.557482Z","shell.execute_reply.started":"2022-09-18T09:03:07.010605Z","shell.execute_reply":"2022-09-18T09:03:07.556652Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Set list of required ensembl ids","metadata":{}},{"cell_type":"code","source":"list_ensembl_ids = list( df_multi_train_y.columns) #[:100] )","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:03:07.559374Z","iopub.execute_input":"2022-09-18T09:03:07.56Z","iopub.status.idle":"2022-09-18T09:03:07.567444Z","shell.execute_reply.started":"2022-09-18T09:03:07.559967Z","shell.execute_reply":"2022-09-18T09:03:07.566031Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# G1S_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']\n# list_ensembl_ids = G1S_genes_Tirosh","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:03:07.569876Z","iopub.execute_input":"2022-09-18T09:03:07.570324Z","iopub.status.idle":"2022-09-18T09:03:07.579919Z","shell.execute_reply.started":"2022-09-18T09:03:07.570279Z","shell.execute_reply":"2022-09-18T09:03:07.578829Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Install/import mygene","metadata":{}},{"cell_type":"code","source":"!pip install mygene\nimport mygene","metadata":{"execution":{"iopub.status.busy":"2022-11-03T22:27:18.802737Z","iopub.execute_input":"2022-11-03T22:27:18.803151Z","iopub.status.idle":"2022-11-03T22:27:32.613914Z","shell.execute_reply.started":"2022-11-03T22:27:18.803118Z","shell.execute_reply":"2022-11-03T22:27:32.612209Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Make a main quiry by mygene - get info on list_ensembl_ids","metadata":{}},{"cell_type":"code","source":"import time\nt0 = time.time()\nmg = mygene.MyGeneInfo()\n# change ensembles\nensembl_ids=['ENSG00000002330', 'ENSG00000030110', 'ENSG00000087088', 'ENSG00000171791', 'ENSG00000142867', 'ENSG00000105327', 'ENSG00000100290', 'ENSG00000165233', 'ENSG00000163098', 'ENSG00000105483', 'ENSG00000137752', 'ENSG00000106144', 'ENSG00000164305', 'ENSG00000196954', 'ENSG00000137757', 'ENSG00000138794', 'ENSG00000165806', 'ENSG00000064012', 'ENSG00000132906', 'ENSG00000003400', 'ENSG00000204403', 'ENSG00000003402', 'ENSG00000169372', 'ENSG00000184047', 'ENSG00000168040', 'ENSG00000125746', 'ENSG00000117560', 'ENSG00000106511', 'ENSG00000145649', 'ENSG00000100453', 'ENSG00000148346', 'ENSG00000186868', 'ENSG00000116688', 'ENSG00000168404', 'ENSG00000249437', 'ENSG00000165030', 'ENSG00000141682', 'ENSG00000123240', 'ENSG00000132646', 'ENSG00000150593', 'ENSG00000156709', 'ENSG00000177595', 'ENSG00000180644', 'ENSG00000111679', 'ENSG00000105327', 'ENSG00000104332', 'ENSG00000105329', 'ENSG00000092969', 'ENSG00000102871', 'ENSG00000121858', 'ENSG00000101966', 'ENSG00000156709', 'ENSG00000142208', 'ENSG00000105221', 'ENSG00000117020', 'ENSG00000120868', 'ENSG00000149311', 'ENSG00000171552', 'ENSG00000187446', 'ENSG00000166869', 'ENSG00000213341', 'ENSG00000100368', 'ENSG00000172115', 'ENSG00000160049', 'ENSG00000169598', 'ENSG00000149218', 'ENSG00000167136', 'ENSG00000157036', 'ENSG00000026103', 'ENSG00000117560', 'ENSG00000104365', 'ENSG00000269335', 'ENSG00000115008', 'ENSG00000125538', 'ENSG00000115594', 'ENSG00000196083', 'ENSG00000164399', 'ENSG00000185291', 'ENSG00000184216', 'ENSG00000134070', 'ENSG00000090376', 'ENSG00000198001', 'ENSG00000006062', 'ENSG00000172936', 'ENSG00000109320', 'ENSG00000100906', 'ENSG00000134259', 'ENSG00000198400', 'ENSG00000121879', 'ENSG00000051382', 'ENSG00000171608', 'ENSG00000105851','ENSG00000145675', 'ENSG00000105647', 'ENSG00000117461', 'ENSG00000141506', 'ENSG00000138814', 'ENSG00000107758', 'ENSG00000120910', 'ENSG00000221823', 'ENSG00000188386', 'ENSG00000072062', 'ENSG00000142875', 'ENSG00000165059', 'ENSG00000108946', 'ENSG00000188191', 'ENSG00000188191', 'ENSG00000005249', 'ENSG00000183943', 'ENSG00000173039', 'ENSG00000137275', 'ENSG00000232810', 'ENSG00000104689', 'ENSG00000120889', 'ENSG00000173535', 'ENSG00000173530', 'ENSG00000067182', 'ENSG00000121858', 'ENSG00000141510', 'ENSG00000127191']\ng = mg.getgenes(ensembl_ids)# [:1000])\n\nprint('%.f seconds passed'%(time.time()-t0))","metadata":{"execution":{"iopub.status.busy":"2022-11-03T22:27:32.616605Z","iopub.execute_input":"2022-11-03T22:27:32.616962Z","iopub.status.idle":"2022-11-03T22:27:35.784651Z","shell.execute_reply.started":"2022-11-03T22:27:32.616929Z","shell.execute_reply":"2022-11-03T22:27:35.783782Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#create new dict\ng2 = {item.get('query'): item.get('symbol') for item in g}\n# print(g2)\nsymbol_for_compare=['CD95', 'BID',  'BIRC2',  'BIRC3',  'CAPN1' ,  'CAPN2',  'CHP1',  'CHP2',  'CHUK',  'CSF2RB',  'CYCS',  'DFFA',  'DFFB',  'ENDOD1' ,  'ENDOG',  'EXOG',  'FAS',  'FASLG',  'IKBKB',  'IKBKG',  'IL1A',  'IL1B',  'IL1R1',  'IL1RAP',  'IL3',  'IL3RA',  'IRAK1',  'IRAK2',  'IRAK3',  'IRAK4',  'MAP3K14',  'MYD88',  'NFKB1',  'NFKBIA',  'NGF',  'NTRK1',  'PIK3CA',  'PIK3CB' , 'PIK3CD',  'PIK3CG' , 'PIK3R1',  'PIK3R2',  'PIK3R3', 'PIK3R5',  'PPP3CA',  'PPP3CB',  'PPP3CC',  'PPP3R1',  'PPP3R2',  'PRKACA', 'PRKACB',  'PRKACG',  'PRKAR1A',  'PRKAR1B',  'PRKAR2A',  'PRKAR2B',  'PRKX' , 'RELA',  'RIPK1',  'TNF' , 'TNFRSF10A',  'TNFRSF10B',  'TNFRSF10C',  'TNFRSF10D',  'TNFRSF1A',  'TNFSF10',  'TP53',  'TRAF2',  'BAD',  'BAK1',  'BAX',  'BCL2',  'BCL10',  'BclXL',  'BIK',  'BINCARD',  'BIRC8',  'CARD8',  'CASP1',  'CASP2',  'CASP3',  'CASP4',  'CASP5',  'CASP6',  'CASP7',  'CASP8',  'CASP9',  'CASP10', 'CASP12',  'CFLAR',  'CRADD',  'DIABLO',  'FADD' , 'EMAP2',  'FASL',  'GAX/MEOX2',  'GZMA',  'GZMB',  'LCN2',  'MAPT',  'MFN2',  'MLKL',  'NAIP1',  'NFIL3',  'NOXA/PMAIP1',  'OPTN' ,'PCNA', 'PDCD4' ,'PDCD8/AIFM1',  'PIDD' , 'PRF1',  'PTPN6',  'PUMA/BBC3',  'SARP2',  'TGFB1',  'TGFB2',  'TRADD',  'TRAIL/TNFSF10',  'XIAP',  'AIFM1' , 'AKT1',  'AKT2',  'AKT3',  'APAF1',  'ATM',  'BCL2L1']","metadata":{"execution":{"iopub.status.busy":"2022-11-03T23:12:51.907075Z","iopub.execute_input":"2022-11-03T23:12:51.907495Z","iopub.status.idle":"2022-11-03T23:12:51.921098Z","shell.execute_reply.started":"2022-11-03T23:12:51.907455Z","shell.execute_reply":"2022-11-03T23:12:51.919284Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(len(ensembl_ids))\nprint(len(symbol_for_compare))","metadata":{"execution":{"iopub.status.busy":"2022-11-03T23:24:41.763051Z","iopub.execute_input":"2022-11-03T23:24:41.763578Z","iopub.status.idle":"2022-11-03T23:24:41.771456Z","shell.execute_reply.started":"2022-11-03T23:24:41.763530Z","shell.execute_reply":"2022-11-03T23:24:41.770068Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#find the identity between 'symbol' from dataset and list\nm=[k for k,v in g2.items() if symbol_for_compare.count(v)]\nidentity=len(m)/len(g)\nprint( \"identity = \", identity)","metadata":{"execution":{"iopub.status.busy":"2022-11-03T22:27:43.856747Z","iopub.execute_input":"2022-11-03T22:27:43.857129Z","iopub.status.idle":"2022-11-03T22:27:43.864624Z","shell.execute_reply.started":"2022-11-03T22:27:43.857098Z","shell.execute_reply":"2022-11-03T22:27:43.863517Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# mismatching\nprint([k for k,v in g2.items() if v not in symbol_for_compare])","metadata":{"execution":{"iopub.status.busy":"2022-11-03T22:27:46.055177Z","iopub.execute_input":"2022-11-03T22:27:46.055613Z","iopub.status.idle":"2022-11-03T22:27:46.063415Z","shell.execute_reply.started":"2022-11-03T22:27:46.055579Z","shell.execute_reply":"2022-11-03T22:27:46.061949Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\n","metadata":{"execution":{"iopub.status.busy":"2022-11-03T23:22:01.563522Z","iopub.execute_input":"2022-11-03T23:22:01.563970Z","iopub.status.idle":"2022-11-03T23:22:01.570048Z","shell.execute_reply.started":"2022-11-03T23:22:01.563939Z","shell.execute_reply":"2022-11-03T23:22:01.568682Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"table = pd.DataFrame({'symbol': pd.Series(g2), \n                      'symbol_for_compare' : pd.array(symbol_for_compare)\n                      \n                      \n                    })\ntable","metadata":{"execution":{"iopub.status.busy":"2022-11-03T23:23:24.157545Z","iopub.execute_input":"2022-11-03T23:23:24.157999Z","iopub.status.idle":"2022-11-03T23:23:24.191957Z","shell.execute_reply.started":"2022-11-03T23:23:24.157965Z","shell.execute_reply":"2022-11-03T23:23:24.189281Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from tabulate import tabulate\nprint(tabulate(g2.items(), headers=['query', 'symbol'], tablefmt=\"grid\"))","metadata":{"execution":{"iopub.status.busy":"2022-11-03T23:01:55.590117Z","iopub.execute_input":"2022-11-03T23:01:55.590580Z","iopub.status.idle":"2022-11-03T23:01:55.604106Z","shell.execute_reply.started":"2022-11-03T23:01:55.590536Z","shell.execute_reply":"2022-11-03T23:01:55.602837Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Some examples what the query returns:","metadata":{}},{"cell_type":"code","source":"g[0].keys()","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:05:13.171903Z","iopub.execute_input":"2022-09-18T09:05:13.172261Z","iopub.status.idle":"2022-09-18T09:05:13.179579Z","shell.execute_reply.started":"2022-09-18T09:05:13.172228Z","shell.execute_reply":"2022-09-18T09:05:13.178366Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"g[0]['query']","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:05:13.181197Z","iopub.execute_input":"2022-09-18T09:05:13.181653Z","iopub.status.idle":"2022-09-18T09:05:13.19322Z","shell.execute_reply.started":"2022-09-18T09:05:13.181611Z","shell.execute_reply":"2022-09-18T09:05:13.192133Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"d = g[0][ 'pathway' ]['reactome']\nlen(d)\nd","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:05:13.197213Z","iopub.execute_input":"2022-09-18T09:05:13.197709Z","iopub.status.idle":"2022-09-18T09:05:13.20639Z","shell.execute_reply.started":"2022-09-18T09:05:13.197675Z","shell.execute_reply":"2022-09-18T09:05:13.205225Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"g[0]['ensembl']","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:05:13.207824Z","iopub.execute_input":"2022-09-18T09:05:13.208258Z","iopub.status.idle":"2022-09-18T09:05:13.217356Z","shell.execute_reply.started":"2022-09-18T09:05:13.208226Z","shell.execute_reply":"2022-09-18T09:05:13.216289Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"g[0]['genomic_pos'],g[0]['genomic_pos_hg19'],g[0]['uniprot'],g[0]['wikipedia'],g[0]['other_names'],g[0]['entrezgene'],","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:05:13.219021Z","iopub.execute_input":"2022-09-18T09:05:13.219375Z","iopub.status.idle":"2022-09-18T09:05:13.230841Z","shell.execute_reply.started":"2022-09-18T09:05:13.219335Z","shell.execute_reply":"2022-09-18T09:05:13.22959Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"g[0]['alias'], g[0]['map_location'], g[0]['name'],g[0]['name'],g[0]['name'],g[0]['name'],g[0]['name'], ","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:05:13.233Z","iopub.execute_input":"2022-09-18T09:05:13.233513Z","iopub.status.idle":"2022-09-18T09:05:13.243173Z","shell.execute_reply.started":"2022-09-18T09:05:13.23347Z","shell.execute_reply":"2022-09-18T09:05:13.242136Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Main reformation obtained information into the table ","metadata":{}},{"cell_type":"code","source":"%%time\ndf = pd.DataFrame(index = range(len(g)))\nfor IX in range(len(df)):\n    k = 'query'\n    \n    if k not in g[IX].keys():\n        continue\n    else:\n        df.loc[IX,k] = g[IX][k]\n        \n    ###########  Ensembl ID, we take only ensembl id , but it  \n    k1 = 'ensembl'\n    if  ( ( k1 in g[IX].keys()) and ( isinstance(g[IX][k1], dict) ) ):\n        if 'gene' in g[IX][k1].keys():  df.loc[IX,'ensembl_id'] = g[IX][k1]['gene']\n    \n    k = 'symbol';     \n    if  k in g[IX].keys():        df.loc[IX,k] = g[IX][k]\n        \n    ########## Chromose, start, end, strand, len \n    k1 = 'genomic_pos'\n    if  not ( ( k1 in g[IX].keys()) and ( isinstance(g[IX][k1], dict) ) ):\n        k1 = 'genomic_pos_hg19'\n    if  ( ( k1 in g[IX].keys()) and ( isinstance(g[IX][k1], dict) ) ):\n        fl = 0\n        for k in ['chr','start','end', 'strand' ]:\n            if k in g[IX][k1].keys():   \n                df.loc[IX,k] = (g[IX][k1][k])\n                if k == 'start': fl+=1\n                if k == 'end': fl+=1\n        if fl == 2:\n            df.loc[IX,'len'] = (g[IX][k1]['end'] - g[IX][k1]['start'])\n    \n    # Number of publications:\n    if 'generif' in g[IX].keys():\n        df.loc[IX,'n_pubmed'] = len(g[IX]['generif'])\n    \n    # Count pathways involved\n    k = 'pathway'\n    if (k in g[IX].keys()) and ( isinstance(g[IX][k], dict) ) :\n        df.loc[IX,'n_pathways'] = np.sum(  [    len(g[IX][k])    for t in g[IX][k]]    )\n        df.loc[IX,'pathways_dbs'] = str(g[IX][k].keys()).replace('dict_keys','').replace('(','').replace(')','')\n    \n    # Count gene ontology categories involved\n    k = 'go'\n    if (k in g[IX].keys()) and ( isinstance(g[IX][k], dict) ) :\n        d = g[IX]['go']\n        a = dict()\n        for k in ['BP','CC','MF']:\n            if k in d.keys():\n                if isinstance( d[k] , dict ): a[k] = 1\n                if isinstance( d[k] , list ): a[k] = len( d[k] )\n                df.loc[IX,'n_GO_'+k] = a[k]\n        \n        \n    # Detailed location - chromosome and its subpart\n    k = 'map_location';     \n    if  k in g[IX].keys():        df.loc[IX,k] = g[IX][k]\n\n    # Less important\n    k1 = 'ensembl'\n    if  ( ( k1 in g[IX].keys()) and ( isinstance(g[IX][k1], dict) ) ):\n        if 'transcript' in g[IX][k1].keys():  df.loc[IX,'n_transcripts'] = len( g[IX][k1]['transcript'])\n    \n    k = 'other_names';     \n    if  k in g[IX].keys():        df.loc[IX,k] = str(g[IX][k])\n    k = 'alias'\n    if  k in g[IX].keys():        df.loc[IX,k] = str(g[IX][k])\n    k = 'entrezgene';     \n    if  k in g[IX].keys():        df.loc[IX,k] = g[IX][k]\n        \n    # Type of gene, name (i.e. brief info), summary - somewhat detailed summary \n    for k in ['type_of_gene', 'name','summary'  ]:\n        if  k in g[IX].keys():        df.loc[IX,k] = str( g[IX][k] )\n        \n        \n#     if IX >5000:\n#         break\n \nfor k in df.columns:\n    if 'n_' in k:\n        df[k] = df[k].fillna(0)\n        df[k] = df[k].astype(int)\ndf","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:05:13.244912Z","iopub.execute_input":"2022-09-18T09:05:13.245297Z","iopub.status.idle":"2022-09-18T09:11:34.05474Z","shell.execute_reply.started":"2022-09-18T09:05:13.245212Z","shell.execute_reply":"2022-09-18T09:11:34.05348Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Might be problem - duplicated queries were returned - remove duplicates","metadata":{}},{"cell_type":"code","source":"len(g), len(list_ensembl_ids)","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:11:34.056518Z","iopub.execute_input":"2022-09-18T09:11:34.057747Z","iopub.status.idle":"2022-09-18T09:11:34.065024Z","shell.execute_reply.started":"2022-09-18T09:11:34.0577Z","shell.execute_reply":"2022-09-18T09:11:34.063748Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"m = df.duplicated(subset = 'query', keep=False)# .sum()\ndf[m]","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:11:34.066554Z","iopub.execute_input":"2022-09-18T09:11:34.06711Z","iopub.status.idle":"2022-09-18T09:11:34.134692Z","shell.execute_reply.started":"2022-09-18T09:11:34.067065Z","shell.execute_reply":"2022-09-18T09:11:34.13355Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.drop_duplicates(subset = 'query', keep = 'first', inplace = True)\nm = df.duplicated(subset = 'query', keep=False)# .sum()\nprint(m.sum(), df.shape)\ndf[m]","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:11:34.136139Z","iopub.execute_input":"2022-09-18T09:11:34.13658Z","iopub.status.idle":"2022-09-18T09:11:34.203989Z","shell.execute_reply.started":"2022-09-18T09:11:34.136548Z","shell.execute_reply":"2022-09-18T09:11:34.202859Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Some returned ensembl Ids might be NaN - substitute them by original queried Id","metadata":{}},{"cell_type":"code","source":"m = (df['query'] != df['ensembl_id']) # &( df['ensembl_id'].notna())#  df[]\nprint(m.sum())\ndf[m]","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:11:34.205308Z","iopub.execute_input":"2022-09-18T09:11:34.205646Z","iopub.status.idle":"2022-09-18T09:11:34.253982Z","shell.execute_reply.started":"2022-09-18T09:11:34.205615Z","shell.execute_reply":"2022-09-18T09:11:34.25279Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"m = (df['query'] != df['ensembl_id']) &( df['ensembl_id'].notna())#  df[]\nprint(m.sum())\ndf[m]","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:11:34.25583Z","iopub.execute_input":"2022-09-18T09:11:34.256268Z","iopub.status.idle":"2022-09-18T09:11:34.284055Z","shell.execute_reply.started":"2022-09-18T09:11:34.256224Z","shell.execute_reply":"2022-09-18T09:11:34.283234Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df['ensembl_id'] = df['query']\n","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:11:34.285541Z","iopub.execute_input":"2022-09-18T09:11:34.286141Z","iopub.status.idle":"2022-09-18T09:11:34.294504Z","shell.execute_reply.started":"2022-09-18T09:11:34.286107Z","shell.execute_reply":"2022-09-18T09:11:34.293046Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.drop('query', axis = 1, inplace = True)\ndf","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:11:34.295887Z","iopub.execute_input":"2022-09-18T09:11:34.29641Z","iopub.status.idle":"2022-09-18T09:11:34.3483Z","shell.execute_reply.started":"2022-09-18T09:11:34.296365Z","shell.execute_reply":"2022-09-18T09:11:34.346214Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Check ensembl_ids are exactly the same as orginal","metadata":{"execution":{"iopub.status.busy":"2022-09-17T20:20:41.263617Z","iopub.execute_input":"2022-09-17T20:20:41.264147Z","iopub.status.idle":"2022-09-17T20:20:41.269194Z","shell.execute_reply.started":"2022-09-17T20:20:41.264112Z","shell.execute_reply":"2022-09-17T20:20:41.268213Z"}}},{"cell_type":"code","source":"list( df['ensembl_id'] ) == list(list_ensembl_ids )","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:11:34.350227Z","iopub.execute_input":"2022-09-18T09:11:34.351226Z","iopub.status.idle":"2022-09-18T09:11:34.368147Z","shell.execute_reply.started":"2022-09-18T09:11:34.351188Z","shell.execute_reply":"2022-09-18T09:11:34.366768Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df['ensembl_id'].isna().sum()","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:11:34.369918Z","iopub.execute_input":"2022-09-18T09:11:34.370313Z","iopub.status.idle":"2022-09-18T09:11:34.382885Z","shell.execute_reply.started":"2022-09-18T09:11:34.370283Z","shell.execute_reply":"2022-09-18T09:11:34.381677Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Fill absent symbols by Ensembl Ids","metadata":{}},{"cell_type":"code","source":"m = df['symbol'].isna()\nprint(m.sum())\ndf[m]","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:11:34.384485Z","iopub.execute_input":"2022-09-18T09:11:34.385011Z","iopub.status.idle":"2022-09-18T09:11:34.428583Z","shell.execute_reply.started":"2022-09-18T09:11:34.384976Z","shell.execute_reply":"2022-09-18T09:11:34.427715Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.loc[m, 'symbol'] = df[m]['ensembl_id']\ndf[m]","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:11:34.42993Z","iopub.execute_input":"2022-09-18T09:11:34.430261Z","iopub.status.idle":"2022-09-18T09:11:34.471773Z","shell.execute_reply.started":"2022-09-18T09:11:34.430234Z","shell.execute_reply":"2022-09-18T09:11:34.470573Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df['symbol'].nunique(), df.shape, df['ensembl_id'].nunique()","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:11:34.47383Z","iopub.execute_input":"2022-09-18T09:11:34.474251Z","iopub.status.idle":"2022-09-18T09:11:34.512949Z","shell.execute_reply.started":"2022-09-18T09:11:34.47421Z","shell.execute_reply":"2022-09-18T09:11:34.511695Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"m = df.duplicated(subset = 'symbol', keep=False)# .sum()\ndf[m]","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:11:34.515086Z","iopub.execute_input":"2022-09-18T09:11:34.515802Z","iopub.status.idle":"2022-09-18T09:11:34.573985Z","shell.execute_reply.started":"2022-09-18T09:11:34.515768Z","shell.execute_reply":"2022-09-18T09:11:34.57286Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Value counts for strand, type of gene, chromosome, etc ","metadata":{"execution":{"iopub.status.busy":"2022-09-17T19:30:17.195721Z","iopub.execute_input":"2022-09-17T19:30:17.196348Z","iopub.status.idle":"2022-09-17T19:30:17.205585Z","shell.execute_reply.started":"2022-09-17T19:30:17.196307Z","shell.execute_reply":"2022-09-17T19:30:17.204283Z"}}},{"cell_type":"code","source":"for c in ['strand','type_of_gene', 'chr', 'n_pathways']:\n    print(); print(c)\n    display(df[c].value_counts())","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:11:34.579369Z","iopub.execute_input":"2022-09-18T09:11:34.579797Z","iopub.status.idle":"2022-09-18T09:11:34.611824Z","shell.execute_reply.started":"2022-09-18T09:11:34.579746Z","shell.execute_reply":"2022-09-18T09:11:34.610762Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Show top genes with respect to number of publications, pathways, Gene Ontology categories, etc ","metadata":{}},{"cell_type":"code","source":"pd.set_option('display.max_rows', 500)\npd.set_option('display.max_columns', 500)\npd.set_option('display.width', 1000)","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:11:34.613284Z","iopub.execute_input":"2022-09-18T09:11:34.614031Z","iopub.status.idle":"2022-09-18T09:11:34.620737Z","shell.execute_reply.started":"2022-09-18T09:11:34.613981Z","shell.execute_reply":"2022-09-18T09:11:34.619498Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for k in df.columns:\n    if ('n_' in k) or (k=='len'):\n        print(); print(k)\n        print( list(df.sort_values(k, ascending = False)['symbol'][:12]) )#  .head(5)\n","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:11:34.622221Z","iopub.execute_input":"2022-09-18T09:11:34.622643Z","iopub.status.idle":"2022-09-18T09:11:34.77814Z","shell.execute_reply.started":"2022-09-18T09:11:34.622594Z","shell.execute_reply":"2022-09-18T09:11:34.777034Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for k in df.columns:\n    if ('n_' in k) or (k=='len'):\n        print(); print(); print(k)\n        display( df.sort_values(k, ascending = False).head(3) )\n\n","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:11:34.779828Z","iopub.execute_input":"2022-09-18T09:11:34.780295Z","iopub.status.idle":"2022-09-18T09:11:35.059475Z","shell.execute_reply.started":"2022-09-18T09:11:34.780251Z","shell.execute_reply":"2022-09-18T09:11:35.058272Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Histograms for numeric variables ","metadata":{}},{"cell_type":"code","source":"fig = plt.figure(figsize = (20,10))\nn_x_subplots = 3\nc = 0\nfor col in ['len', 'n_pubmed', 'n_pathways',  'n_GO_BP', 'n_GO_CC', 'n_GO_MF']:\n    c += 1; fig.add_subplot(2,n_x_subplots ,c)\n    plt.hist( df[col], bins = 100 )\n    plt.title(col, fontsize = 20 )\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:11:35.06074Z","iopub.execute_input":"2022-09-18T09:11:35.061066Z","iopub.status.idle":"2022-09-18T09:11:36.786432Z","shell.execute_reply.started":"2022-09-18T09:11:35.061037Z","shell.execute_reply":"2022-09-18T09:11:36.78506Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# There are 7363 genes which are not well studied - no pubmed, no pathways, no GO (at least according to \"mygene\")","metadata":{}},{"cell_type":"code","source":"v = df[[ 'n_pubmed', 'n_pathways',  'n_GO_BP', 'n_GO_CC', 'n_GO_MF']].sum(axis = 1)\nprint(df.shape)\nfor i in range(10):\n    print( i,  (v<= i).sum() )\n","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:11:36.788468Z","iopub.execute_input":"2022-09-18T09:11:36.788933Z","iopub.status.idle":"2022-09-18T09:11:36.802449Z","shell.execute_reply.started":"2022-09-18T09:11:36.788888Z","shell.execute_reply":"2022-09-18T09:11:36.801142Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Attach infornation on trascription factors.\n\nData from http://humantfs.ccbr.utoronto.ca/\n\nSee also https://en.wikipedia.org/wiki/List_of_human_transcription_factors, https://en.wikipedia.org/wiki/Transcription_factor","metadata":{}},{"cell_type":"code","source":"fn = '/kaggle/input/genes-information/TranscriptionFactorsHuman_LambertDatabaseExtract_v_101.csv'\ntf = pd.read_csv(fn, index_col = 0)\ntf.head(2)","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:11:36.803882Z","iopub.execute_input":"2022-09-18T09:11:36.804224Z","iopub.status.idle":"2022-09-18T09:11:36.88904Z","shell.execute_reply.started":"2022-09-18T09:11:36.804192Z","shell.execute_reply":"2022-09-18T09:11:36.887738Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tf.columns, df.columns","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:11:36.890675Z","iopub.execute_input":"2022-09-18T09:11:36.891129Z","iopub.status.idle":"2022-09-18T09:11:36.899819Z","shell.execute_reply.started":"2022-09-18T09:11:36.891084Z","shell.execute_reply":"2022-09-18T09:11:36.898519Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tmp = tf[['Ensembl ID',  'DBD',  'Is TF?', 'TF assessment', 'Binding mode',  'Final Notes', 'Final Comments',\n          'EntrezGene Description','TFCat classification', ]]\ndf2 = pd.merge( left = df, right = tmp, left_on = 'ensembl_id', right_on = 'Ensembl ID', how = 'left',  )\ndf2 = df2.drop( 'Ensembl ID', axis = 1 )\ndf2","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:11:36.901312Z","iopub.execute_input":"2022-09-18T09:11:36.901679Z","iopub.status.idle":"2022-09-18T09:11:37.069011Z","shell.execute_reply.started":"2022-09-18T09:11:36.901648Z","shell.execute_reply":"2022-09-18T09:11:37.067853Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Save data to csv on disk","metadata":{}},{"cell_type":"code","source":"%%time\ndf.to_csv('genes_information_mmscel_challenge_multiome_subtask.csv')\nprint(df.shape)\n","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:11:37.070322Z","iopub.execute_input":"2022-09-18T09:11:37.070797Z","iopub.status.idle":"2022-09-18T09:11:37.520256Z","shell.execute_reply.started":"2022-09-18T09:11:37.070754Z","shell.execute_reply":"2022-09-18T09:11:37.51911Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndf2.to_csv('genes_information_mmscel_challenge_multiome_subtask_tf_info.csv')\nprint(df2.shape)\n","metadata":{"execution":{"iopub.status.busy":"2022-09-18T09:11:37.521833Z","iopub.execute_input":"2022-09-18T09:11:37.522162Z","iopub.status.idle":"2022-09-18T09:11:37.996807Z","shell.execute_reply.started":"2022-09-18T09:11:37.522131Z","shell.execute_reply":"2022-09-18T09:11:37.995607Z"},"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-09-18T09:11:37.998152Z","iopub.execute_input":"2022-09-18T09:11:37.99853Z","iopub.status.idle":"2022-09-18T09:11:38.004696Z","shell.execute_reply.started":"2022-09-18T09:11:37.998496Z","shell.execute_reply":"2022-09-18T09:11:38.003476Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}