{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceType":"competition","sourceId":59094,"databundleVersionId":7010844},{"sourceType":"datasetVersion","sourceId":7068220,"datasetId":3955392,"databundleVersionId":7155770},{"sourceType":"datasetVersion","sourceId":7121955,"datasetId":3948965,"databundleVersionId":7209943},{"sourceType":"datasetVersion","sourceId":6891925,"datasetId":3841755,"databundleVersionId":6977933}],"dockerImageVersionId":30558,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# What is about ?\n\nAnalysis of the average correlation public submissions, and distributions for selected genes/genes groups.\n\n    The least correlated with current top - are candidates to include to the blend.\n    \n    Clusters of higly correlated - should be treated separately - first choose top from the cluster (or blend within cluster), and only then begin to blend from different clusters. \n    \nThe precise suggestions how we can try to improve the blend are here:\n\nhttps://www.kaggle.com/competitions/open-problems-single-cell-perturbations/discussion/453110#2513189    \n    \n\nDatasets with collections of practically ALL available public submits are attached to that notebook. \n\nThe presentation with comments is here: \nhttps://docs.google.com/presentation/d/1wiz0Wmt4D54pqMMsIOyJHuQYMZ3hTBZQQnjbLzwoGYY/edit?usp=sharing\n\n\nSome comments on the collection are here: \nhttps://docs.google.com/spreadsheets/d/1APN63PMaWZygVjYimK9Ivt0RvifdAU5JRYkxiDn4szw/edit?usp=sharing\n    ","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 time\nt0start = time.time() \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\nimport pandas as pd\nimport tensorflow as tf\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport seaborn as sns\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\nc = 0\nprint('Fist 30 files:')\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        c+=1\n        if c<=30:\n            #print(os.path.join(dirname, filename))\n            print(filename) # 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-12-04T15:20:16.857613Z","iopub.execute_input":"2023-12-04T15:20:16.858040Z","iopub.status.idle":"2023-12-04T15:20:31.143957Z","shell.execute_reply.started":"2023-12-04T15:20:16.857974Z","shell.execute_reply":"2023-12-04T15:20:31.142483Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Select files ","metadata":{}},{"cell_type":"code","source":"list_fn = []\nlist_ids = []\n# Current top public: \nlist_fn.append('/kaggle/input/open-problems-2-blends/LB0577_ZacharyScheben_blend_SCPblend058_nbV19.csv')\nlist_ids.append('LB577_TopBlend_Zachary')\nlist_fn.append('/kaggle/input/open-problems-2-blends/LB583_JJLee_blend_EnsembleSubmitionOP_nbV3.csv')\nlist_ids.append('LB583_ExTopBlend_JJLee')\n\n\nwork_mode = 'main_public_submits'\n# work_mode = 'submits_factory'\n# work_mode = 'all'\n\nif (work_mode == 'main_public_submits') or ( work_mode == 'all'):\n    for dirname, _, filenames in os.walk('/kaggle/input'):\n        for filename in filenames:\n            if 'open-problems-2-submits-collection' in dirname:\n            #if 1:# ('LB0608' in dirname):\n                #if  ('simple' in filename) and ('blend' not in filename):\n                #   ('Random' in filename) and ('blend' not in filename) and ('.csv' in filename):\n                #if  ('Random' in filename) and ('priors' in filename) and ('.csv' in filename):\n                if ('.csv' in filename):\n                    f = os.path.join(dirname, filename)\n                    list_fn.append( f )\n                    list_ids.append(f.split('/')[-1][:20]) \n\nif (work_mode == 'submits_factory') or (work_mode == 'all'):\n    for dirname, _, filenames in os.walk('/kaggle/input'):\n        for filename in filenames:\n            if 'open-problems-single-cell-perturbations-submitsetc' in dirname:\n                #if 1:# ('LB0608' in dirname):\n                if  ('Random' in filename) and ('blend' not in filename) and ('.csv' in filename):\n                # if  ('simple' in filename) and ('blend' not in filename) and ('.csv' in filename):\n                #   ('Random' in filename) and ('blend' not in filename) and ('.csv' in filename):\n                #if  ('Random' in filename) and ('priors' in filename) and ('.csv' in filename):\n                    f = os.path.join(dirname, filename)\n                    list_fn.append( f )\n                    list_ids.append(f.split('/')[-2][:30]) \n\nif 1: # Filter:                     \n    list_fn2 = []\n    list_ids2 = []\n    for i in range(len(list_fn)):\n        if ('pyboost' in list_ids[i].lower()) or ( 'CATB' in list_ids[i].upper() ):\n            list_fn2.append( list_fn[i] )\n            list_ids2.append( list_ids[i] )\n    list_fn = list_fn2\n    list_ids = list_ids2\n\n        \nprint('Total files:', len(list_fn))    \nprint(len(list_ids))    \nprint('Short ids:',list_ids)\nprint()\nprint('First 5 filepaths:')\nfor i in range(5):\n    if i < len(list_fn):\n        print(i, list_fn[i])","metadata":{"execution":{"iopub.status.busy":"2023-11-22T21:16:58.454674Z","iopub.execute_input":"2023-11-22T21:16:58.455102Z","iopub.status.idle":"2023-11-22T21:16:58.522125Z","shell.execute_reply.started":"2023-11-22T21:16:58.455071Z","shell.execute_reply":"2023-11-22T21:16:58.520695Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load and StandardScale","metadata":{}},{"cell_type":"code","source":"%%time\n\nflag_add_blend_submits = False\n\n\nprint('Start loading ', len(list_ids), 'files. Might take 5-10 minutes.')\nprint()\nlist_df = []\n\ni_blend = 0\nfor i,fn in enumerate(list_fn): \n#     print(i, fn)\n    df = pd.read_csv(fn, index_col = 'id')\n    print(i,df.shape, list_ids[i] , np.round( [df.iloc[0,0],  df.iloc[0,1],  df.iloc[1,0] ,  df.iloc[1,1] ], 4) , df.columns[0],   df.columns[1]   )\n    #display(df.head(2) )\n    list_df.append(df)\n    if flag_add_blend_submits:\n        if i > 0: # not blend the first one\n            if i_blend == 0:\n                df_blend = df.copy()\n            else:\n                df_blend = (df_blend * i_blend  + df )/ (i_blend + 1)\n            i_blend += 1\n            #display(df_blend.head(2) )\n        \nif flag_add_blend_submits:\n    list_ids.append('BlendAll')  \n    list_df.append(df_blend )    \n    \nprint('All files loaded');\nprint()\n\n# %%time\nfrom sklearn.preprocessing import StandardScaler\nscaler = StandardScaler()\n\nlist_np = []\nif 1:\n    for k in range(len(list_df)): # dict_df:\n        df = list_df[k] \n        d = scaler.fit_transform(df)\n        list_np.append(d)\nelse:\n    for k in dict_df:\n        df = dict_df[k] \n        d = scaler.fit_transform(df)\n        list_np.append(d)\n    \nprint(len(list_np)) ","metadata":{"execution":{"iopub.status.busy":"2023-11-22T21:16:58.524488Z","iopub.execute_input":"2023-11-22T21:16:58.524997Z","iopub.status.idle":"2023-11-22T21:18:52.686117Z","shell.execute_reply.started":"2023-11-22T21:16:58.524955Z","shell.execute_reply":"2023-11-22T21:18:52.684707Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Compute correlations\n\nCorrelations can be compute as elementwise products averaged over columns (because the data has beed standard-scaled). Then average correlation for all genes. \n","metadata":{}},{"cell_type":"code","source":"%%time\n\nfig = plt.figure(figsize = (20,10))\ndf_stat = pd.DataFrame()\ndf_corr_averages_for_submits = pd.DataFrame()\ni_total = 0\nfor i0,d0 in  enumerate(list_np):\n    for i1,d1 in  enumerate(list_np):\n        i_total += 1\n        v = np.mean( d0*d1, axis = 0 )\n        if i0<i1:\n            plt.hist(v, label = str(i0)+ ' vs ' + str(i1))\n        dt = pd.Series(v).describe().to_frame()\n        dt.columns = [str(i0)+ ' vs  ' + str(i1)]\n        df_stat = pd.concat( [df_stat, dt] , axis = 1)\n        #print( i_total )\n        # c = np.mean( np.abs(v) )\n        c = np.mean( (v) )\n        if i_total<3:\n            print(c)\n        df_corr_averages_for_submits.loc[i0,i1] = c\n        \ndisplay(df_stat)        \n# plt.legend(fontsize = 12 )\nplt.grid()\nplt.title('Distribution of correlations (over genes) for each prediction pair ', fontsize = 20)\nplt.show()        \ndf_corr_averages_for_submits.columns = [t.replace('_',' ') for t in  list_ids]\ndf_corr_averages_for_submits.index = [t.replace('_',' ') for t in  list_ids]\n\nprint('Show 5x5 part of the correlation matrix:')\ndisplay(df_corr_averages_for_submits.iloc[:5,:5].round(2) )  ","metadata":{"execution":{"iopub.status.busy":"2023-11-22T21:18:52.687769Z","iopub.execute_input":"2023-11-22T21:18:52.688120Z","iopub.status.idle":"2023-11-22T21:19:02.837173Z","shell.execute_reply.started":"2023-11-22T21:18:52.688091Z","shell.execute_reply":"2023-11-22T21:19:02.836168Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_corr_averages_for_submits.to_csv('df_corr_averages_for_submits.csv')","metadata":{"execution":{"iopub.status.busy":"2023-11-22T21:19:02.840244Z","iopub.execute_input":"2023-11-22T21:19:02.840608Z","iopub.status.idle":"2023-11-22T21:19:02.849359Z","shell.execute_reply.started":"2023-11-22T21:19:02.840579Z","shell.execute_reply":"2023-11-22T21:19:02.848205Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Clustermap","metadata":{}},{"cell_type":"code","source":"cm = sns.clustermap( df_corr_averages_for_submits,figsize = (20,15),  annot=True, cmap='coolwarm')\nplt.show()\n\n\nif len(df_corr_averages_for_submits)>50:\n    sns.clustermap( df_corr_averages_for_submits,figsize = (20,12), annot=False, cmap='coolwarm')\n    plt.show()\n    \nlabels = cm.ax_heatmap.yaxis.get_majorticklabels()\nfor label in labels:\n    print(label.get_text())    ","metadata":{"execution":{"iopub.status.busy":"2023-11-22T21:19:02.850735Z","iopub.execute_input":"2023-11-22T21:19:02.851074Z","iopub.status.idle":"2023-11-22T21:19:04.580556Z","shell.execute_reply.started":"2023-11-22T21:19:02.851046Z","shell.execute_reply":"2023-11-22T21:19:04.579760Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if work_mode == 'submits_factory':\n    l = [i for i in range(len(list_ids)) if '608' in list_ids[i]]\n    # list_ids[l[0]], list_ids[l[1]]\n    print(l)\n    for i, k in enumerate(l):\n        if i == 0:\n            df_average = list_df[k]\n        else:\n            df_average += list_df[k]\n    df_average /= len(l)\n    for i, d in enumerate(list_df):\n        c = df_average.corrwith(d).mean()\n        print(i, c, list_ids[i] )\n    # l_corr = [ df_average.corrwith(d).mean() for d in list_df ]\n    # l_corr\n","metadata":{"execution":{"iopub.status.busy":"2023-11-22T21:19:04.581522Z","iopub.execute_input":"2023-11-22T21:19:04.581834Z","iopub.status.idle":"2023-11-22T21:19:04.591916Z","shell.execute_reply.started":"2023-11-22T21:19:04.581808Z","shell.execute_reply":"2023-11-22T21:19:04.590390Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Candidates to uplift the top\n\nLowest correlated with it ","metadata":{}},{"cell_type":"code","source":"df_corr_averages_for_submits.iloc[:,0].sort_values().head(30)","metadata":{"execution":{"iopub.status.busy":"2023-11-22T21:19:04.593152Z","iopub.execute_input":"2023-11-22T21:19:04.593527Z","iopub.status.idle":"2023-11-22T21:19:04.612947Z","shell.execute_reply.started":"2023-11-22T21:19:04.593483Z","shell.execute_reply":"2023-11-22T21:19:04.611900Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nplt.figure(figsize = (20,4) )\nv = df_corr_averages_for_submits.iloc[:,0].sort_values()\nplt.plot(v.values,'*-' )\nplt.xticks(range(len(v)), v.index)\nplt.xticks(rotation=90)\nplt.grid()\nplt.title('Least correlated - candidates to uplift the top', fontsize = 20 )\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-11-22T21:19:04.614899Z","iopub.execute_input":"2023-11-22T21:19:04.615418Z","iopub.status.idle":"2023-11-22T21:19:05.056519Z","shell.execute_reply.started":"2023-11-22T21:19:04.615373Z","shell.execute_reply":"2023-11-22T21:19:05.055222Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# KDE (s) ","metadata":{}},{"cell_type":"code","source":"genes_mitochondrial = [t for t in list_df[0].columns if t.startswith('MT')]\nprint(len(genes_mitochondrial), genes_mitochondrial)\n\n# Cell proliferation Tirosh genes:\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']\ngenes_Tirosh = G1S_genes_Tirosh + G2M_genes_Tirosh\ngenes_Tirosh = [t for t in list_df[0].columns if t in genes_Tirosh]\n","metadata":{"execution":{"iopub.status.busy":"2023-11-22T21:19:05.058151Z","iopub.execute_input":"2023-11-22T21:19:05.058558Z","iopub.status.idle":"2023-11-22T21:19:05.112189Z","shell.execute_reply.started":"2023-11-22T21:19:05.058525Z","shell.execute_reply":"2023-11-22T21:19:05.111114Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ncc = 0\nfor gene_groups in [ ['TP53'], ['SRC'], ['MYC'], ['CDK1'], ['CDK2'], ['CDK4'] ] + [ genes_Tirosh ] + [ genes_mitochondrial] + [ list(list_df[0].columns[:100])] :\n# for gene_name in ['TP53', 'SRC', 'MYC',  'CDK1', 'CDK2','CDK4']:\n    cc += 1\n    fig = plt.figure(figsize = (20,5))\n    for k in range(len(list_df)): # dict_df:\n        str_inf = list_ids[k]\n        df = list_df[k] \n        sns.kdeplot(df[gene_groups].values.ravel(), label = str_inf)\n        #sns.kdeplot(df.values.ravel(), label = str_inf)\n        #print(k, str_inf)\n\n    plt.title(str(gene_groups)[:30]+' Genes predictions', fontsize = 20 )    \n    if cc == 1:\n        plt.legend()# fontsize = 20)\n    plt.grid()\n    plt.show()    \n","metadata":{"execution":{"iopub.status.busy":"2023-11-22T21:19:05.113753Z","iopub.execute_input":"2023-11-22T21:19:05.114079Z","iopub.status.idle":"2023-11-22T21:19:22.627271Z","shell.execute_reply.started":"2023-11-22T21:19:05.114051Z","shell.execute_reply":"2023-11-22T21:19:22.626139Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Final timing","metadata":{}},{"cell_type":"code","source":"print('%.1f seconds passed total '%(time.time()-t0start) )\nprint('%.1f minutes passed total '%( (time.time()-t0start)/60)  )\nprint('%.2f hours passed total '%( (time.time()-t0start)/3600)  )","metadata":{"execution":{"iopub.status.busy":"2023-11-22T21:19:22.629092Z","iopub.execute_input":"2023-11-22T21:19:22.629770Z","iopub.status.idle":"2023-11-22T21:19:22.637455Z","shell.execute_reply.started":"2023-11-22T21:19:22.629714Z","shell.execute_reply":"2023-11-22T21:19:22.636265Z"},"trusted":true},"execution_count":null,"outputs":[]}]}