{"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\nSome examples on \"genes sets enrichment analysis\" with  nice  package \"gseapy\".\nDocumentation: https://gseapy.readthedocs.io/en/latest/index.html - provides more examples, visualizations, analysis etc.\n\n\n\nGenes sets enrichment analysis is comparaison of some given genes set with known sets from biological databases.\nIt is the standard tool nowdays. It also provides some p-values based on hypergeometric test measuring how random the interesection between given set and known sets seems to be. \n\nOne of the largest databases for genes sets is https://www.gsea-msigdb.org/gsea/index.jsp,\nwhich contains more that 30 000 sets, clustered into hundreds groups, in turn clustered into several higher level \"collections\": \"Hallmark gene sets\", \"Regulatory target gene sets \", etc.. \n\n\nThe most commonly used  genes sets are based on the pathway databases: \"Reactome\", \"KEGG\", etc...; also \"Gene Ontology\" project.  The other noteable sets include collection \"Hallmark gene sets (50 gene sets)\" from msgidb itself, various targets of the transcription factros, etc... \n\nThe package \"gseapy\" provides Python API to work with the \"msigdb\" database (https://www.gsea-msigdb.org ) \n\nOne can do it manually by the web-page: https://www.gsea-msigdb.org/gsea/msigdb/human/annotate.jsp - insert your list to: \"Input Gene Identifiers\" field. \n\n\nThe notebook contains several examples  to work with gseapy using some genes sets obtained as top related to some proteins from the single cell rna+protein  dataset from the Kaggle/NIPS competition: https://www.kaggle.com/competitions/open-problems-multimodal\n\nWe also query all available gseapy singnatures for top100 correlated to cd36 - and save resulting dataframe to csv - there are 189274 gene sets processed. \n\n\n","metadata":{}},{"cell_type":"markdown","source":"# Preparations","metadata":{}},{"cell_type":"markdown","source":"## Install GSEAPY","metadata":{}},{"cell_type":"code","source":"!pip install gseapy\nimport gseapy as gp","metadata":{"execution":{"iopub.status.busy":"2023-01-04T07:56:07.556159Z","iopub.execute_input":"2023-01-04T07:56:07.556742Z","iopub.status.idle":"2023-01-04T07:56:22.156725Z","shell.execute_reply.started":"2023-01-04T07:56:07.556629Z","shell.execute_reply":"2023-01-04T07:56:22.155741Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Load some genes sets ","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":"2023-01-04T07:57:36.629030Z","iopub.execute_input":"2023-01-04T07:57:36.629468Z","iopub.status.idle":"2023-01-04T07:57:36.657676Z","shell.execute_reply.started":"2023-01-04T07:57:36.629425Z","shell.execute_reply":"2023-01-04T07:57:36.656785Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Load genes sets top related to CD36","metadata":{}},{"cell_type":"code","source":"fn = '/kaggle/input/research-project-01-around-multimodal-singlecell/CD36importances.csv'\ndf = pd.read_csv(fn)# , index_col = 0)\ndisplay(df.head(2))\nl = list(df.iloc[0,:])\ndf.columns = l\ndf = df.iloc[1:,:]\ndf\n","metadata":{"execution":{"iopub.status.busy":"2023-01-04T07:57:38.430217Z","iopub.execute_input":"2023-01-04T07:57:38.431417Z","iopub.status.idle":"2023-01-04T07:57:38.493045Z","shell.execute_reply.started":"2023-01-04T07:57:38.431368Z","shell.execute_reply":"2023-01-04T07:57:38.492193Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Define list of genes ","metadata":{}},{"cell_type":"code","source":"col = 'CorrInters_4cellType_alldays_alldonors'\ncol = 'CorrInters_EryP>7donors*days'\ncol = 'CorrTop100All'\nlist_genes = [ t.split('_')[-1] for t in df[col][df[col].notnull()]  ]\nprint(len(list_genes), list_genes[:60])","metadata":{"execution":{"iopub.status.busy":"2023-01-04T08:20:43.819359Z","iopub.execute_input":"2023-01-04T08:20:43.819837Z","iopub.status.idle":"2023-01-04T08:20:43.834640Z","shell.execute_reply.started":"2023-01-04T08:20:43.819802Z","shell.execute_reply":"2023-01-04T08:20:43.832949Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pd.set_option('display.max_rows', 500)\n# pd.set_option('display.max_columns', 500)\n# pd.set_option('display.width', 1000)\n","metadata":{"execution":{"iopub.status.busy":"2023-01-04T08:20:45.736561Z","iopub.execute_input":"2023-01-04T08:20:45.737760Z","iopub.status.idle":"2023-01-04T08:20:45.743455Z","shell.execute_reply.started":"2023-01-04T08:20:45.737709Z","shell.execute_reply":"2023-01-04T08:20:45.741939Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Genes enrichment with most famous genes sets: Reactome, KEGG, Gene Ontology,  etc.. ","metadata":{}},{"cell_type":"code","source":"%%time\nprint(len(list_genes), list_genes[:50])\nenr = gp.enrichr(\n    gene_list=list_genes, \n    gene_sets= ['KEGG_2021_Human', 'Reactome_2022','GO_Molecular_Function_2021'  ], # , 'GO_Biological_Process_2021', 'MSigDB_Oncogenic_Signatures',\n    organism='human',\n    outdir=None,\n)\ncol = 'Adjusted P-value' #P-value # Adjusted P-value - with Bonferoni type correction for checking multiple lists, while just \"P-value\" - is without\ndisplay( enr.results.sort_values(col).head(50) )\n    \n","metadata":{"execution":{"iopub.status.busy":"2023-01-04T08:56:51.228869Z","iopub.execute_input":"2023-01-04T08:56:51.229361Z","iopub.status.idle":"2023-01-04T08:56:55.868568Z","shell.execute_reply.started":"2023-01-04T08:56:51.229322Z","shell.execute_reply":"2023-01-04T08:56:55.867348Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Transcription factor targets related lists","metadata":{}},{"cell_type":"markdown","source":"## Lists based on more \"direct\" evidences (like chip-seq, motives... in contrast to less direct - \"coexpression\" approaches - see below)\n\ngene_sets: 'ChEA_2013',  'ChEA_2015', 'ChEA_2016', 'ChEA_2022' - target genes of transcription factors from published ChIP-chip, ChIP-seq, and other transcription factor binding site profiling studies - see  https://maayanlab.cloud/Harmonizome/dataset/CHEA+Transcription+Factor+Targets\n\n\n\nAll(?)  related to TF lists in gseapy : \n\n    list_tf_related_signatures = [ \n    'ChEA_2013',  'ChEA_2015', 'ChEA_2016', 'ChEA_2022',  \n    'ENCODE_TF_ChIP-seq_2015', 'ENCODE_and_ChEA_Consensus_TFs_from_ChIP-X', \n    'TF-LOF_Expression_from_GEO', 'TF_Perturbations_Followed_by_Expression', \n    'TRANSFAC_and_JASPAR_PWMs',\n    'TRRUST_Transcription_Factors_2019', \n    'ARCHS4_TFs_Coexp', 'Enrichr_Submissions_TF-Gene_Coocurrence',\n    'Transcription_Factor_PPIs',\n    ] \n","metadata":{}},{"cell_type":"code","source":"#Enrichr_Submissions_TF-Gene_Coocurrence\nprint(len(list_genes), list_genes[:50])\nfor sig in ['ChEA_2022',     'ENCODE_TF_ChIP-seq_2015', 'ENCODE_and_ChEA_Consensus_TFs_from_ChIP-X', \n    'TF-LOF_Expression_from_GEO', 'TF_Perturbations_Followed_by_Expression',     'TRANSFAC_and_JASPAR_PWMs',\n    'TRRUST_Transcription_Factors_2019',   ]:\n    \n    print(); print(sig)\n    if 'ChEA' in sig:\n        print(' \"ChEA\" - see: https://maayanlab.cloud/Harmonizome/dataset/CHEA+Transcription+Factor+Targets ' )\n        \n    enr = gp.enrichr(\n        gene_list=list_genes, # 'MSigDB_Oncogenic_Signatures',\n        gene_sets= sig,\n        organism='human',\n        outdir=None,\n    )\n    col = 'Adjusted P-value' #P-value # Adjusted P-value - with Bonferoni type correction for checking multiple lists, while just \"P-value\" - is without\n    m = enr.results[col] < 0.05\n    print(m.sum(), 'out of:', len( enr.results ) )\n    display( enr.results[m].sort_values(col).head(20) )\n#     display( enr.results.sort_values('Adjusted P-value').head(20) )","metadata":{"execution":{"iopub.status.busy":"2023-01-04T09:23:16.020545Z","iopub.execute_input":"2023-01-04T09:23:16.021337Z","iopub.status.idle":"2023-01-04T09:23:29.664236Z","shell.execute_reply.started":"2023-01-04T09:23:16.021296Z","shell.execute_reply":"2023-01-04T09:23:29.662878Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## TF targets lists based related to co-expression or PPIs","metadata":{}},{"cell_type":"code","source":"#Enrichr_Submissions_TF-Gene_Coocurrence\nprint(len(list_genes), list_genes[:50])\nfor sig in [    'ARCHS4_TFs_Coexp', 'Enrichr_Submissions_TF-Gene_Coocurrence',\n    'Transcription_Factor_PPIs'  ]:\n    \n    print(); print(sig)\n    if 'ChEA' in sig:\n        print(' \"ChEA\" - see: https://maayanlab.cloud/Harmonizome/dataset/CHEA+Transcription+Factor+Targets ' )\n        \n    enr = gp.enrichr(\n        gene_list=list_genes, # 'MSigDB_Oncogenic_Signatures',\n        gene_sets= sig,\n        organism='human',\n        outdir=None,\n    )\n    col = 'Adjusted P-value' #P-value # Adjusted P-value - with Bonferoni type correction for checking multiple lists, while just \"P-value\" - is without\n    m = enr.results[col] < 0.05\n    print(m.sum(), 'out of:', len( enr.results ) )\n    display( enr.results[m].sort_values(col).head(20) )\n#     display( enr.results.sort_values('Adjusted P-value').head(20) )","metadata":{"execution":{"iopub.status.busy":"2023-01-04T09:27:09.279222Z","iopub.execute_input":"2023-01-04T09:27:09.279741Z","iopub.status.idle":"2023-01-04T09:27:15.190563Z","shell.execute_reply.started":"2023-01-04T09:27:09.279697Z","shell.execute_reply":"2023-01-04T09:27:15.189395Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Get all gseapy available genes groups \n\nNote: Each group contains many genes sets","metadata":{}},{"cell_type":"code","source":"%%time\nnames = gp.get_library_name()\nprint(type(names), len(names))\nprint('First 10:', names[:10])\nprint()\n","metadata":{"execution":{"iopub.status.busy":"2023-01-04T09:34:12.250532Z","iopub.execute_input":"2023-01-04T09:34:12.251040Z","iopub.status.idle":"2023-01-04T09:34:12.856084Z","shell.execute_reply.started":"2023-01-04T09:34:12.250999Z","shell.execute_reply":"2023-01-04T09:34:12.854815Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"s = ''\nfor nm in names:\n    s += nm + ' '*10\n#print(s)\nimport pprint \n# pp = pprint.PrettyPrinter()\npp = pprint.PrettyPrinter(width=90, compact=True)\npp.pprint( s )\n","metadata":{"execution":{"iopub.status.busy":"2023-01-04T09:37:01.136378Z","iopub.execute_input":"2023-01-04T09:37:01.137670Z","iopub.status.idle":"2023-01-04T09:37:01.151542Z","shell.execute_reply.started":"2023-01-04T09:37:01.137622Z","shell.execute_reply":"2023-01-04T09:37:01.149532Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Get and  show particular genes sets for some groups","metadata":{}},{"cell_type":"code","source":"%%time\nfor nm in names[:5]:\n    lb_name = nm#  'ARCHS4_Kinases_Coexp'\n    ## download library or read a .gmt file\n    lb = gp.get_library(name=lb_name)# , organism='Yeast')\n    #print(go_mf['ATP binding (GO:0005524)'])\n    print()\n    print(nm, type(lb), len(lb))\n    if len(lb) <= 0: continue\n    print('First 10 genes sets:', list(lb.keys())[:10] )\n    k = list(lb.keys())[0]\n    l = list(lb[k])[:10]\n    print('First 10 genes in the first set:', l )\n    ","metadata":{"execution":{"iopub.status.busy":"2023-01-04T09:49:40.716172Z","iopub.execute_input":"2023-01-04T09:49:40.716827Z","iopub.status.idle":"2023-01-04T09:50:26.104054Z","shell.execute_reply.started":"2023-01-04T09:49:40.716770Z","shell.execute_reply":"2023-01-04T09:50:26.102876Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Query all groups - do not do that \n\n419 seconds ","metadata":{"execution":{"iopub.status.busy":"2023-01-04T09:43:24.884884Z","iopub.execute_input":"2023-01-04T09:43:24.885392Z","iopub.status.idle":"2023-01-04T09:43:24.893853Z","shell.execute_reply.started":"2023-01-04T09:43:24.885354Z","shell.execute_reply":"2023-01-04T09:43:24.892487Z"}}},{"cell_type":"code","source":"%%time\nimport time\nt0 = time.time()\nprint(len(list_genes), list_genes[:50])\n\nfor i, sig in enumerate( names ):\n    \n    print(); print(sig)\n    if 'ChEA' in sig:\n        print(' \"ChEA\" - see: https://maayanlab.cloud/Harmonizome/dataset/CHEA+Transcription+Factor+Targets ' )\n        \n    enr = gp.enrichr(\n        gene_list=list_genes, # 'MSigDB_Oncogenic_Signatures',\n        gene_sets= sig,\n        organism='human',\n        outdir=None,\n    )\n    if i == 0:\n        df = enr.results.copy()\n    else: \n        if (len(enr.results) > 0) and (isinstance(enr.results, pd.DataFrame) ):\n            df = pd.concat([df, enr.results ], axis = 0)\n            display(df.tail(1))\n            print('Current result - df.shape:', df.shape, '%.1f seconds passed'%(time.time() - t0))","metadata":{"execution":{"iopub.status.busy":"2023-01-04T09:55:33.323851Z","iopub.execute_input":"2023-01-04T09:55:33.324346Z","iopub.status.idle":"2023-01-04T10:02:32.650028Z","shell.execute_reply.started":"2023-01-04T09:55:33.324311Z","shell.execute_reply":"2023-01-04T10:02:32.649132Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\ndf['n_overlap']=0\nfor i in range(len(df)):\n    s = df['Overlap'].iat[i]\n    if isinstance(s,str ) and ('/' in s):\n        n = int(s.split('/')[0])\n        #print(n); break\n        df['n_overlap'].iat[i] = n\ndf.sort_values( 'n_overlap', ascending = False ).head(20)","metadata":{"execution":{"iopub.status.busy":"2023-01-04T10:29:46.063239Z","iopub.execute_input":"2023-01-04T10:29:46.063795Z","iopub.status.idle":"2023-01-04T10:30:03.412224Z","shell.execute_reply.started":"2023-01-04T10:29:46.063754Z","shell.execute_reply":"2023-01-04T10:30:03.411043Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.sort_values( 'n_overlap', ascending = False ).to_csv('CD36_correlated_gseapy_ALL_sets.csv')","metadata":{"execution":{"iopub.status.busy":"2023-01-04T10:39:52.007661Z","iopub.execute_input":"2023-01-04T10:39:52.008157Z","iopub.status.idle":"2023-01-04T10:39:53.917036Z","shell.execute_reply.started":"2023-01-04T10:39:52.008119Z","shell.execute_reply":"2023-01-04T10:39:53.915720Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Analysis of results obtained from all signatures ","metadata":{}},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nv = df.sort_values( 'n_overlap', ascending = False )['n_overlap']\nplt.plot(v.values[:200])\nplt.show()\nplt.hist(v.values)\nplt.show()\nv.describe()\n\nv = df.sort_values( 'Adjusted P-value', ascending = True )['Adjusted P-value']\nplt.plot(np.log(v.values[:150] ) )\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-01-04T10:41:46.759191Z","iopub.execute_input":"2023-01-04T10:41:46.759704Z","iopub.status.idle":"2023-01-04T10:41:47.267094Z","shell.execute_reply.started":"2023-01-04T10:41:46.759662Z","shell.execute_reply":"2023-01-04T10:41:47.266131Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"m = df['Gene_set'] != 'Jensen_TISSUES'\ndf[m].sort_values( 'n_overlap', ascending = False ).head(100)","metadata":{"execution":{"iopub.status.busy":"2023-01-04T10:34:14.533879Z","iopub.execute_input":"2023-01-04T10:34:14.534396Z","iopub.status.idle":"2023-01-04T10:34:14.758718Z","shell.execute_reply.started":"2023-01-04T10:34:14.534345Z","shell.execute_reply":"2023-01-04T10:34:14.757351Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.sort_values( 'Adjusted P-value', ascending = True ).head(50)\n","metadata":{"execution":{"iopub.status.busy":"2023-01-04T10:42:12.562141Z","iopub.execute_input":"2023-01-04T10:42:12.562913Z","iopub.status.idle":"2023-01-04T10:42:12.651629Z","shell.execute_reply.started":"2023-01-04T10:42:12.562871Z","shell.execute_reply":"2023-01-04T10:42:12.650362Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.shape","metadata":{"execution":{"iopub.status.busy":"2023-01-04T11:07:04.593083Z","iopub.execute_input":"2023-01-04T11:07:04.593594Z","iopub.status.idle":"2023-01-04T11:07:04.601560Z","shell.execute_reply.started":"2023-01-04T11:07:04.593548Z","shell.execute_reply":"2023-01-04T11:07:04.600277Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}