{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# What is about ?\n\n\nAnalyse correlations between CD proteins and RNA. Choose top N correlated and save to csv.\n\nFor the CITE-seq single cell RNA seq + CD-proteins data from the Kaggle / NIPS22 competition: \nhttps://www.kaggle.com/competitions/open-problems-multimodal\n\n\n**Conclusions:** For CD36 protein:\n\n    0.73 Pearson Correlation (0.5 Spearman) with CD36 RNA\n    Interesection of the top100 - quite big - 79 \n\n    RNA for CD36 is practically zero on HSC, NeuP, MasP cell types - so their consideration should probably be avoided \n    HSC\t97.8\n    NeuP\t96.5\n    EryP\t32.2\n    MasP\t93.5\n    MoP \t35.5\n\n    Percent Zero RNA\tCorr Pearson CD36 Protein vs RNA\t\n                            Count\tMean Protein\tMean RNA\n    HSC\t97.8\t0.046302\t29879.0\t0.814702\t0.083006\n    NeuP\t96.5\t0.032865\t12493.0\t0.728760\t0.132491\n    EryP\t32.2\t0.622676\t14241.0\t11.629602\t3.278012\n    MasP\t93.5\t0.411630\t8242.0\t1.053324\t0.260617\n    MoP\t35.5\t0.664600\t591.0\t13.004719\t3.341950\n\n            Counts:\n    HSC     29879\n    EryP    14241\n    NeuP    12493\n    MasP     8242\n    MkP      5382\n    MoP       591\n    BP        160\n    \n    For random vectors of such dimension the std of correlations is about 0.0037, thus observed correlations for most of the rna are above 10*sigma - explanation for that:\n    When we fix cell type, donor, day - that effect disappears - we see std of correlations between random vectors and biological vectors is approximately the same. The most important is fixing the cell type. \n    That gives an opportunity to estimate threshold such that correlations above it are NOT random. \n    For gaussian probability to be beyond 4 sigma is 6.334E-05  , thus with Bonferony correction for 1000 samples we get 0.063\n    For gaussian probability to be beyond 3.5 sigma is 4.653E-04  , thus with Bonferony correction for 100 samples we get 0.046\n    \n    So we can take about top100 correlations 3.5 sigma, or 1000 correlations beyond 4 sigma - as statistically meaningful. Possible even more is Okay. \n    \n**Outcome::** For the cell type EryP we have taking top100 correlated  would be stastically meaningful by all possible estimates of std.\n\n**Lists Intersection analysis:** probably the list of meaningfully correlated genes is quite big - 2000+-  the linear approximation of the intersection becomes worse than quadratic around topK K= 2000. Intersection is taken for Pearson and Spearman lists.  The slope is about 3/4. \n\n**Enrichment:** enriching with transcription facrtors lists and even with 0.01 p-value, and even interserting lists - we get about 100 tfs - too many ? \n\n**V9,13: LGB:** correlations with LGB from each feature - non-linear dependences for individual features - it is very slow \n\n**Correlation RNA-Protein CD36 (on EryP):** decreases with the day, but in different way for different donors - see same named section \"So correlation protein-rna decreases with the day, but in different manner for different donors - reasons are unclear\" \n    \n    Technical: Pearson correlations by the direct loop 23 seconds for calculations, Spearman - 6.5 minutes\n\n","metadata":{}},{"cell_type":"markdown","source":"# Key Params","metadata":{"execution":{"iopub.status.busy":"2022-12-18T21:18:08.728359Z","iopub.execute_input":"2022-12-18T21:18:08.729709Z","iopub.status.idle":"2022-12-18T21:18:08.735176Z","shell.execute_reply.started":"2022-12-18T21:18:08.729662Z","shell.execute_reply":"2022-12-18T21:18:08.733772Z"}}},{"cell_type":"code","source":"target_name = 'CD36'\n\nn_top = 100\n\nslow_calculations_on = 1 # calculations which take hours e.g. LGB for each feauture and corrlation  based on it - it might take 6 hours ","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-01-06T14:14:57.675199Z","iopub.execute_input":"2023-01-06T14:14:57.675604Z","iopub.status.idle":"2023-01-06T14:14:57.681918Z","shell.execute_reply.started":"2023-01-06T14:14:57.675571Z","shell.execute_reply":"2023-01-06T14:14:57.680684Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Preparations","metadata":{}},{"cell_type":"code","source":"import pandas as pd\n\n# Main results will be stored here: \ndf_results = pd.DataFrame()","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:14:57.687252Z","iopub.execute_input":"2023-01-06T14:14:57.688206Z","iopub.status.idle":"2023-01-06T14:14:57.694928Z","shell.execute_reply.started":"2023-01-06T14:14:57.688171Z","shell.execute_reply":"2023-01-06T14:14:57.693817Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pd.set_option('display.max_columns', 500)\n# pd.set_option('display.max_rows', 500)\n# pd.set_option('display.width', 1000)","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:14:57.696845Z","iopub.execute_input":"2023-01-06T14:14:57.697141Z","iopub.status.idle":"2023-01-06T14:14:57.705396Z","shell.execute_reply.started":"2023-01-06T14:14:57.697113Z","shell.execute_reply":"2023-01-06T14:14:57.704443Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nimport time\nt0start = time.time()\n\nimport os","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:14:57.709584Z","iopub.execute_input":"2023-01-06T14:14:57.710142Z","iopub.status.idle":"2023-01-06T14:14:57.717056Z","shell.execute_reply.started":"2023-01-06T14:14:57.710100Z","shell.execute_reply":"2023-01-06T14:14:57.715879Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load data for Proteins","metadata":{}},{"cell_type":"code","source":"%%time\ndf_y = pd.read_hdf('/kaggle/input/open-problems-multimodal/train_cite_targets.h5')\ndf_y","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:14:57.740335Z","iopub.execute_input":"2023-01-06T14:14:57.740961Z","iopub.status.idle":"2023-01-06T14:14:58.798818Z","shell.execute_reply.started":"2023-01-06T14:14:57.740925Z","shell.execute_reply":"2023-01-06T14:14:58.797697Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#  Load RNA expression data \n","metadata":{}},{"cell_type":"code","source":"%%time\ndf_rna = pd.read_hdf('/kaggle/input/open-problems-multimodal/train_cite_inputs.h5')\ndisplay(df_rna) \n\n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:14:58.801030Z","iopub.execute_input":"2023-01-06T14:14:58.801571Z","iopub.status.idle":"2023-01-06T14:16:07.157713Z","shell.execute_reply.started":"2023-01-06T14:14:58.801526Z","shell.execute_reply":"2023-01-06T14:16:07.156888Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Correlations (Pearson) CD vs RNA","metadata":{}},{"cell_type":"code","source":"%%time\nl = []\ny = df_y[target_name]\nfor col in df_rna.columns:\n    c = np.corrcoef(y,df_rna[col])[0,1]\n    l.append(c)\n    if (len(l)%5000) == 1:\n        print(len(l))\n        \ncol = 'Corr Pearson NIPS22'\ndf_stat = pd.DataFrame(index = df_rna.columns, data = l , columns = [col] )\ndf_stat['Abs '+col] = df_stat[col].abs()\ndf_stat.sort_values('Abs '+col, ascending = False).head(n_top)\n        ","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:16:07.158984Z","iopub.execute_input":"2023-01-06T14:16:07.159535Z","iopub.status.idle":"2023-01-06T14:16:29.037639Z","shell.execute_reply.started":"2023-01-06T14:16:07.159502Z","shell.execute_reply":"2023-01-06T14:16:29.036456Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"col = 'Corr Pearson NIPS22'\nv = df_stat.sort_values('Abs '+col, ascending = False)['Abs '+col].head(n_top*2)\nplt.figure(figsize = (20,5))\nplt.plot(v.values,'*-')\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:16:29.039966Z","iopub.execute_input":"2023-01-06T14:16:29.040269Z","iopub.status.idle":"2023-01-06T14:16:29.265599Z","shell.execute_reply.started":"2023-01-06T14:16:29.040241Z","shell.execute_reply":"2023-01-06T14:16:29.264496Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_results['Corr Pearson NIPS22'] = df_stat.sort_values('Abs '+col, ascending = False).index[:n_top]\ndf_results","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:16:29.266711Z","iopub.execute_input":"2023-01-06T14:16:29.266991Z","iopub.status.idle":"2023-01-06T14:16:29.286183Z","shell.execute_reply.started":"2023-01-06T14:16:29.266964Z","shell.execute_reply":"2023-01-06T14:16:29.285160Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"display( df_stat.describe() )\nimport matplotlib.pyplot as plt\nplt.hist(df_stat[col], bins = 100 )\nplt.show()\n\nfor t in [0.01, 0.1, 0.2,0.3,0.4,0.5,0.6,0.7]:\n    m = df_stat['Abs '+col] > t\n    print(t, m.sum(), np.round( m.sum()/len(df_stat) * 100,4)  )","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:16:29.287610Z","iopub.execute_input":"2023-01-06T14:16:29.288185Z","iopub.status.idle":"2023-01-06T14:16:29.647513Z","shell.execute_reply.started":"2023-01-06T14:16:29.288141Z","shell.execute_reply":"2023-01-06T14:16:29.646510Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Spearman Correlations","metadata":{}},{"cell_type":"code","source":"%%time\nfrom scipy import stats\nl = []\ny = df_y[target_name]\nfor col in df_rna.columns:\n    c = stats.spearmanr(y,df_rna[col])[0]\n    l.append(c)\n    if (len(l)%5000) == 1:\n        print(len(l))\n        \ncol = 'Corr Spearman NIPS22'\n#df_stat = pd.DataFrame(index = df_rna.columns, data = l , columns = [col] )\ndf_stat[col] = l\ndf_stat['Abs '+col] = df_stat[col].abs()\ndf_stat.sort_values('Abs '+col, ascending = False).head(n_top)\n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:16:29.648716Z","iopub.execute_input":"2023-01-06T14:16:29.649025Z","iopub.status.idle":"2023-01-06T14:22:21.560860Z","shell.execute_reply.started":"2023-01-06T14:16:29.648985Z","shell.execute_reply":"2023-01-06T14:22:21.559510Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize = (20,5))\nplt.plot(df_stat.sort_values('Abs '+col, ascending = False).head(n_top)['Abs '+col].values, '*-' )\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:22:21.562601Z","iopub.execute_input":"2023-01-06T14:22:21.563393Z","iopub.status.idle":"2023-01-06T14:22:21.787532Z","shell.execute_reply.started":"2023-01-06T14:22:21.563352Z","shell.execute_reply":"2023-01-06T14:22:21.786390Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_stat","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:22:21.789311Z","iopub.execute_input":"2023-01-06T14:22:21.789756Z","iopub.status.idle":"2023-01-06T14:22:21.807190Z","shell.execute_reply.started":"2023-01-06T14:22:21.789714Z","shell.execute_reply":"2023-01-06T14:22:21.805864Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_results['Corr Spearman NIPS22'] = df_stat.sort_values('Abs '+col, ascending = False).index[:n_top]\ndf_results        ","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:22:21.813386Z","iopub.execute_input":"2023-01-06T14:22:21.814676Z","iopub.status.idle":"2023-01-06T14:22:21.834835Z","shell.execute_reply.started":"2023-01-06T14:22:21.814622Z","shell.execute_reply":"2023-01-06T14:22:21.833747Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"display( df_stat.describe() )\nimport matplotlib.pyplot as plt\nplt.hist(df_stat[col], bins = 100 )\nplt.show()\n\nfor t in [0.1, 0.2,0.3,0.4,0.5,0.6,0.7]:\n    m = df_stat['Abs '+col] > t\n    print(t, m.sum(), np.round( m.sum()/len(df_stat) * 100,4)  )","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:22:21.836209Z","iopub.execute_input":"2023-01-06T14:22:21.836520Z","iopub.status.idle":"2023-01-06T14:22:22.214557Z","shell.execute_reply.started":"2023-01-06T14:22:21.836493Z","shell.execute_reply":"2023-01-06T14:22:22.213233Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"s = set(df_results.iloc[:,0]) & set(df_results.iloc[:,1])\n\nprint(len(s))\nprint( (s))","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:22:22.216199Z","iopub.execute_input":"2023-01-06T14:22:22.216656Z","iopub.status.idle":"2023-01-06T14:22:22.223791Z","shell.execute_reply.started":"2023-01-06T14:22:22.216613Z","shell.execute_reply":"2023-01-06T14:22:22.222525Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_results.head(n_top)","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:22:22.225877Z","iopub.execute_input":"2023-01-06T14:22:22.226669Z","iopub.status.idle":"2023-01-06T14:22:22.242633Z","shell.execute_reply.started":"2023-01-06T14:22:22.226624Z","shell.execute_reply":"2023-01-06T14:22:22.241694Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_results.to_csv(target_name + '_Top'+str(n_top)+'Correlated.csv' ) ","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:22:22.243916Z","iopub.execute_input":"2023-01-06T14:22:22.244299Z","iopub.status.idle":"2023-01-06T14:22:22.251398Z","shell.execute_reply.started":"2023-01-06T14:22:22.244270Z","shell.execute_reply":"2023-01-06T14:22:22.250497Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#df_results.to_csv(target_name + '_Top'+str(n_top)+'Correlated.csv' ) \n\ndf_stat.to_csv(target_name + '_all_rna_correlations.csv' ) ","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:22:22.253030Z","iopub.execute_input":"2023-01-06T14:22:22.253447Z","iopub.status.idle":"2023-01-06T14:22:22.422452Z","shell.execute_reply.started":"2023-01-06T14:22:22.253406Z","shell.execute_reply":"2023-01-06T14:22:22.421660Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Analysis of intersection lenghts for Pearson and Spearman \n\nThe quadratic (parabolic) growth of intersection lenght - indication of randomness.\n\nWhile linear growth - non randomness.\n\n\nhttps://mathoverflow.net/q/437569/10446\n\nhttps://mathoverflow.net/q/437545/10446","metadata":{"execution":{"iopub.status.busy":"2023-01-06T11:17:43.512356Z","iopub.execute_input":"2023-01-06T11:17:43.512702Z","iopub.status.idle":"2023-01-06T11:17:43.518483Z","shell.execute_reply.started":"2023-01-06T11:17:43.512671Z","shell.execute_reply":"2023-01-06T11:17:43.517653Z"}}},{"cell_type":"code","source":"%%time\nlist_interesection_lengths = []\nfor k in range(len(df_stat)):\n    for i,col in enumerate(['Abs Corr Pearson NIPS22', 'Abs Corr Spearman NIPS22']):\n        if i == 0:\n            s = set( df_stat[col].sort_values( ascending = False).index[:k] )\n        else:\n            s = s & set(  df_stat[col].sort_values( ascending = False).index[:k] )\n    list_interesection_lengths.append(len(s))\n    ","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:22:22.424114Z","iopub.execute_input":"2023-01-06T14:22:22.424837Z","iopub.status.idle":"2023-01-06T14:26:13.920873Z","shell.execute_reply.started":"2023-01-06T14:22:22.424794Z","shell.execute_reply":"2023-01-06T14:26:13.919854Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ll = list_interesection_lengths\nplt.figure(figsize = (20,6))\nplt.plot(ll, label = 'real data')# , label = 'intersection lengths ' )\n\ny = np.array(ll)\nx = np.arange(len(y))\np = np.polyfit(x,y,2)\nprint('Polynom coefficients:', p)\nplt.plot(np.polyval(p,x), label = 'Parabolic approximation' )\n    \nplt.title('intersection lengths  for topX Pearson and Spearman  correlated', fontsize = 20)\nplt.grid()\nplt.legend(fontsize = 20)\nplt.show()\n\n\nplt.figure(figsize = (20,8))\nplt.plot(ll, label = 'real data ' )\n\ny = np.array(ll)\nx = np.arange(len(y))\np = np.polyfit(x,y,2)\nprint('Polynom coefficients:', p)\nplt.plot(np.polyval(p,x), label = 'Parabolic approximation' )\n\ny = np.array(ll)\nx = np.arange(len(y))\np = np.polyfit(x[:100],y[:100],1)\nprint('Polynom coefficients:', p)\nplt.plot(np.polyval(p,x), label = 'Linear approximation' )\n\nplt.ylim([0,2500])\nplt.xlim([0,2500])\n\nplt.title('intersection lengths  for topX Pearson and Spearman  correlated', fontsize = 20)\nplt.grid()\nplt.legend(fontsize = 20)\nplt.show()\n\n\n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:26:13.922454Z","iopub.execute_input":"2023-01-06T14:26:13.922808Z","iopub.status.idle":"2023-01-06T14:26:14.575300Z","shell.execute_reply.started":"2023-01-06T14:26:13.922776Z","shell.execute_reply":"2023-01-06T14:26:14.574521Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Calculate Correlation with LGB from each feature","metadata":{}},{"cell_type":"code","source":"%%time\n\nif slow_calculations_on:\n    import lightgbm as lgbm\n    from sklearn.model_selection import cross_val_predict\n    from sklearn.model_selection import KFold\n    import time\n    t0 = time.time()\n\n    n_splits_for_cross_valdition = 2\n\n    model = lgbm.LGBMRegressor(random_state = 0) # Ridge(alpha = alpha_selected )\n\n    for random_state_cross_validation in [0,1]:\n        kf = KFold(n_splits=n_splits_for_cross_valdition,  shuffle=True, random_state= random_state_cross_validation )\n\n        l = []\n        y = df_y[target_name]\n        for col in df_rna.columns:\n            y_pred = cross_val_predict(model, df_rna[[col]], y, cv=kf)\n\n            c = np.corrcoef(y,y_pred)[0,1]\n            l.append(c)\n            if (len(l)%500) == 1:\n                print(len(l), '%.1f secs passed'%(time.time()-t0))\n\n        col = 'LGB Corr Pearson NIPS22 rs'+str(random_state_cross_validation )\n        df_stat[col] = l\n        df_stat.sort_values(col, ascending = False).head(10) # n_top)\n\n    df_results['LGB Corr Pearson NIPS22 rs0'] = df_stat.sort_values('LGB Corr Pearson NIPS22 rs0', ascending = False).index[:n_top]\n    df_results['LGB Corr Pearson NIPS22 rs1'] = df_stat.sort_values('LGB Corr Pearson NIPS22 rs1', ascending = False).index[:n_top]\n    df_results        \n\n    df_results.to_csv(target_name + '_Top'+str(n_top)+'Correlated.csv' ) \n\n    df_stat.to_csv(target_name + '_all_rna_correlations.csv' )         ","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:26:14.576600Z","iopub.execute_input":"2023-01-06T14:26:14.577094Z","iopub.status.idle":"2023-01-06T14:26:14.587203Z","shell.execute_reply.started":"2023-01-06T14:26:14.577060Z","shell.execute_reply":"2023-01-06T14:26:14.586115Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_stat.describe()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\nif slow_calculations_on:\n\n    list_interesection_lengths = []\n    for k in range(len(df_stat)):\n        for i,col in enumerate(['Abs Corr Pearson NIPS22', 'LGB Corr Pearson NIPS22 rs0']):\n\n            if i == 0:\n                s = set( df_stat[col].sort_values( ascending = False).index[:k] )\n            else:\n                s = s & set(  df_stat[col].sort_values( ascending = False).index[:k] )\n        list_interesection_lengths.append(len(s))\n\n\n    ll = list_interesection_lengths\n    plt.figure(figsize = (20,6))\n    plt.plot(ll, label = 'real data')# , label = 'intersection lengths ' )\n\n    y = np.array(ll)\n    x = np.arange(len(y))\n    p = np.polyfit(x,y,2)\n    print('Polynom coefficients:', p)\n    plt.plot(np.polyval(p,x), label = 'Parabolic approximation' )\n\n    plt.title('intersection lengths  for topX Pearson and Pearson-LGB  correlated', fontsize = 20)\n    plt.grid()\n    plt.legend(fontsize = 20)\n    plt.show()\n\n\n    plt.figure(figsize = (20,8))\n    plt.plot(ll, label = 'real data ' )\n\n    y = np.array(ll)\n    x = np.arange(len(y))\n    p = np.polyfit(x,y,2)\n    print('Polynom coefficients:', p)\n    plt.plot(np.polyval(p,x), label = 'Parabolic approximation' )\n\n    y = np.array(ll)\n    x = np.arange(len(y))\n    p = np.polyfit(x[:100],y[:100],1)\n    print('Polynom coefficients:', p)\n    plt.plot(np.polyval(p,x), label = 'Linear approximation' )\n\n    plt.ylim([0,2500])\n    plt.xlim([0,2500])\n\n    plt.title('intersection lengths  for topX Pearson and Pearson-LGB  correlated', fontsize = 20)\n    plt.grid()\n    plt.legend(fontsize = 20)\n    plt.show()\n\n\n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:26:14.588816Z","iopub.execute_input":"2023-01-06T14:26:14.589497Z","iopub.status.idle":"2023-01-06T14:26:14.605511Z","shell.execute_reply.started":"2023-01-06T14:26:14.589427Z","shell.execute_reply":"2023-01-06T14:26:14.604225Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\nif slow_calculations_on:\n\n    list_interesection_lengths = []\n    for k in range(len(df_stat)):\n        for i,col in enumerate(['Abs Corr Spearman NIPS22', 'LGB Corr Pearson NIPS22 rs0']):\n\n            if i == 0:\n                s = set( df_stat[col].sort_values( ascending = False).index[:k] )\n            else:\n                s = s & set(  df_stat[col].sort_values( ascending = False).index[:k] )\n        list_interesection_lengths.append(len(s))\n\n\n    ll = list_interesection_lengths\n    plt.figure(figsize = (20,6))\n    plt.plot(ll, label = 'real data')# , label = 'intersection lengths ' )\n\n    y = np.array(ll)\n    x = np.arange(len(y))\n    p = np.polyfit(x,y,2)\n    print('Polynom coefficients:', p)\n    plt.plot(np.polyval(p,x), label = 'Parabolic approximation' )\n\n    plt.title('intersection lengths  for topX Spearman and Pearson-LGB  correlated', fontsize = 20)\n    plt.grid()\n    plt.legend(fontsize = 20)\n    plt.show()\n\n\n    plt.figure(figsize = (20,8))\n    plt.plot(ll, label = 'real data ' )\n\n    y = np.array(ll)\n    x = np.arange(len(y))\n    p = np.polyfit(x,y,2)\n    print('Polynom coefficients:', p)\n    plt.plot(np.polyval(p,x), label = 'Parabolic approximation' )\n\n    y = np.array(ll)\n    x = np.arange(len(y))\n    p = np.polyfit(x[:100],y[:100],1)\n    print('Polynom coefficients:', p)\n    plt.plot(np.polyval(p,x), label = 'Linear approximation' )\n\n    plt.ylim([0,2500])\n    plt.xlim([0,2500])\n\n    plt.title('intersection lengths  for topX Spearman and Pearson-LGB  correlated', fontsize = 20)\n    plt.grid()\n    plt.legend(fontsize = 20)\n    plt.show()\n\n\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\nif slow_calculations_on:\n\n    list_interesection_lengths = []\n    for k in range(len(df_stat)):\n        for i,col in enumerate(['LGB Corr Pearson NIPS22 rs0', 'LGB Corr Pearson NIPS22 rs1']):\n\n            if i == 0:\n                s = set( df_stat[col].sort_values( ascending = False).index[:k] )\n            else:\n                s = s & set(  df_stat[col].sort_values( ascending = False).index[:k] )\n        list_interesection_lengths.append(len(s))\n\n\n    ll = list_interesection_lengths\n    plt.figure(figsize = (20,6))\n    plt.plot(ll, label = 'real data')# , label = 'intersection lengths ' )\n\n    y = np.array(ll)\n    x = np.arange(len(y))\n    p = np.polyfit(x,y,2)\n    print('Polynom coefficients:', p)\n    plt.plot(np.polyval(p,x), label = 'Parabolic approximation' )\n\n    plt.title('intersection lengths  for topX LGB rs0 and LGB rs1  correlated', fontsize = 20)\n    plt.grid()\n    plt.legend(fontsize = 20)\n    plt.show()\n\n\n    plt.figure(figsize = (20,8))\n    plt.plot(ll, label = 'real data ' )\n\n    y = np.array(ll)\n    x = np.arange(len(y))\n    p = np.polyfit(x,y,2)\n    print('Polynom coefficients:', p)\n    plt.plot(np.polyval(p,x), label = 'Parabolic approximation' )\n\n    y = np.array(ll)\n    x = np.arange(len(y))\n    p = np.polyfit(x[:100],y[:100],1)\n    print('Polynom coefficients:', p)\n    plt.plot(np.polyval(p,x), label = 'Linear approximation' )\n\n    plt.ylim([0,2500])\n    plt.xlim([0,2500])\n\n    plt.title('intersection lengths  for topX LGB rs0 and LGB rs1 correlated', fontsize = 20)\n    plt.grid()\n    plt.legend(fontsize = 20)\n    plt.show()\n\n\n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:26:14.607660Z","iopub.execute_input":"2023-01-06T14:26:14.608098Z","iopub.status.idle":"2023-01-06T14:26:14.625973Z","shell.execute_reply.started":"2023-01-06T14:26:14.608067Z","shell.execute_reply":"2023-01-06T14:26:14.624680Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Enrichment Analysis","metadata":{}},{"cell_type":"code","source":"!pip install gseapy\nimport gseapy as gp","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:26:14.627613Z","iopub.execute_input":"2023-01-06T14:26:14.628140Z","iopub.status.idle":"2023-01-06T14:26:29.505427Z","shell.execute_reply.started":"2023-01-06T14:26:14.628106Z","shell.execute_reply":"2023-01-06T14:26:29.502778Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_stat.sort_values('Corr Pearson NIPS22', ascending = True).head(20)","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:26:29.509730Z","iopub.execute_input":"2023-01-06T14:26:29.510450Z","iopub.status.idle":"2023-01-06T14:26:29.552231Z","shell.execute_reply.started":"2023-01-06T14:26:29.510362Z","shell.execute_reply":"2023-01-06T14:26:29.550365Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## KEGG Enrichment","metadata":{}},{"cell_type":"code","source":"%%time\n#list_genes  = df_results['Corr Pearson NIPS22'].tolist()\nlist_genes  = df_results['Corr Spearman NIPS22'].tolist()\n\nlist_genes  = list( df_stat.sort_values('Corr Pearson NIPS22', ascending = False).index[:300] )\n\nlist_genes = [t.split('_')[1] for t in list_genes]\nprint(list_genes[:100])\n\nenr = gp.enrichr(\n    gene_list=list_genes,\n    gene_sets=['KEGG_2016','KEGG_2021_Human'],\n    organism='human',\n    outdir=None,\n)\ndisplay( enr.results.sort_values('P-value').head(20) )\n\n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:26:29.553615Z","iopub.execute_input":"2023-01-06T14:26:29.554058Z","iopub.status.idle":"2023-01-06T14:26:34.848857Z","shell.execute_reply.started":"2023-01-06T14:26:29.554022Z","shell.execute_reply":"2023-01-06T14:26:34.847419Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Transciption factor related enrichement","metadata":{}},{"cell_type":"markdown","source":"### Lists directly related to transcription factors ","metadata":{}},{"cell_type":"code","source":"%%time\n#list_genes  = df_results['Corr Pearson NIPS22'].tolist()\n#list_genes  = df_results['Corr Spearman NIPS22'].tolist()\n\nlist_genes  = list( df_stat.sort_values('Corr Pearson NIPS22', ascending = False).index[:300] )\n\nlist_genes = [t.split('_')[1] for t in list_genes]\nprint(list_genes[:100])\n\nenr = gp.enrichr(\n    gene_list=list_genes,\n    gene_sets= ['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    organism='human',\n    outdir=None,\n)\ndisplay( enr.results.sort_values('P-value').head(20) )\n\ncol = 'Adjusted P-value'\nm = enr.results[col] < 0.01\nprint(col, '<0.01 count:', m.sum() )\nv = enr.results['Term'][m]\nl = []\nfor k in v:\n    c = k.split(' ')[0].upper()\n    l.append(c)\ns = set(l)\nprint(len(s), s)\nsave_tf1 =  s\n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:26:34.851238Z","iopub.execute_input":"2023-01-06T14:26:34.851649Z","iopub.status.idle":"2023-01-06T14:26:49.927719Z","shell.execute_reply.started":"2023-01-06T14:26:34.851615Z","shell.execute_reply":"2023-01-06T14:26:49.926535Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Lists non-directly related to transcription factors  (coexpression, ppi)","metadata":{}},{"cell_type":"code","source":"%%time\n#list_genes  = df_results['Corr Pearson NIPS22'].tolist()\n#list_genes  = df_results['Corr Spearman NIPS22'].tolist()\n\nlist_genes  = list( df_stat.sort_values('Corr Pearson NIPS22', ascending = False).index[:300] )\n\nlist_genes = [t.split('_')[1] for t in list_genes]\nprint(list_genes[:100])\n\nenr = gp.enrichr(\n    gene_list=list_genes,\n    gene_sets= [    'ARCHS4_TFs_Coexp', 'Enrichr_Submissions_TF-Gene_Coocurrence',\n    'Transcription_Factor_PPIs'  ],\n    organism='human',\n    outdir=None,\n)\ndisplay( enr.results.sort_values('P-value').head(20) )\n\n\ncol = 'Adjusted P-value'\nm = enr.results[col] < 0.01\nprint(col, '<0.01 count:', m.sum() )\nv = enr.results['Term'][m]\nl = []\nfor k in v:\n    c = k.split(' ')[0].upper()\n    l.append(c)\ns = set(l)\nprint(len(s), s)\nsave_tf2 =  s\n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:26:49.929209Z","iopub.execute_input":"2023-01-06T14:26:49.929567Z","iopub.status.idle":"2023-01-06T14:26:57.061914Z","shell.execute_reply.started":"2023-01-06T14:26:49.929534Z","shell.execute_reply":"2023-01-06T14:26:57.060246Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"s = save_tf1 & save_tf2\nprint('Intersection of tfs lists:', len(s), s )","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:26:57.064004Z","iopub.execute_input":"2023-01-06T14:26:57.064531Z","iopub.status.idle":"2023-01-06T14:26:57.072074Z","shell.execute_reply.started":"2023-01-06T14:26:57.064481Z","shell.execute_reply":"2023-01-06T14:26:57.070794Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Look on non-zero RNA by cell types  ","metadata":{}},{"cell_type":"markdown","source":"## Preparations","metadata":{}},{"cell_type":"code","source":"%%time\n\nrna_name = ''\nfor t in df_rna.columns:\n    if target_name in t: rna_name = t\nprint(rna_name)\n\n\nfn = '/kaggle/input/open-problems-multimodal/metadata.csv'\ndf_meta = pd.read_csv(fn, index_col = 0 )\ndf_meta\n# Cut only train cite-seq part: \nd = pd.DataFrame(index = df_y.index)\nprint(d.shape)\ndf_meta = d.join(df_meta, how = 'left')\ndf_meta","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:26:57.073425Z","iopub.execute_input":"2023-01-06T14:26:57.073754Z","iopub.status.idle":"2023-01-06T14:26:57.711533Z","shell.execute_reply.started":"2023-01-06T14:26:57.073725Z","shell.execute_reply":"2023-01-06T14:26:57.710283Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Results","metadata":{}},{"cell_type":"code","source":"%%time\nfrom scipy import stats\n\nd = pd.DataFrame(  ); \nfor cell_type in ['HSC', 'NeuP', 'EryP', 'MasP','MoP']:\n    m = df_meta['cell_type'] == cell_type # 'NeuP'\n    v1 = df_y[m][target_name];     \n    v2 = df_rna[m][rna_name];     \n    d.loc[cell_type,'Corr Pearson '+target_name + ' Protein vs RNA'] = np.corrcoef(v1,v2)[0,1] \n    d.loc[cell_type,'Corr Spearman '+target_name + ' Protein vs RNA'] = stats.spearmanr(v1,v2)[0]\n    d.loc[cell_type,'Count'] = m.sum()#  np.corrcoef(v1,v2)[0,1] \n    d.loc[cell_type,'Mean Protein'] = v1.mean() #  np.corrcoef(v1,v2)[0,1] \n    d.loc[cell_type,'Mean RNA'] = v2.mean() #  np.corrcoef(v1,v2)[0,1] \n    d.loc[cell_type,'Count Zero RNA'] = (v2 == 0).sum() #  np.corrcoef(v1,v2)[0,1] \n    d.loc[cell_type,'Percent Zero RNA'] = np.round( (v2 == 0).sum() / len(v2) * 100 , 1) #  np.corrcoef(v1,v2)[0,1] \n    \n    \n#display( d.describe() )\ndisplay(d)\ndisplay(d[ ['Percent Zero RNA','Corr Pearson CD36 Protein vs RNA', 'Count', 'Mean Protein', 'Mean RNA']])","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:26:57.721561Z","iopub.execute_input":"2023-01-06T14:26:57.722894Z","iopub.status.idle":"2023-01-06T14:27:05.831859Z","shell.execute_reply.started":"2023-01-06T14:26:57.722851Z","shell.execute_reply":"2023-01-06T14:27:05.830809Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfrom scipy import stats\nprint('Results where RNA is NON zero ' )\nd = pd.DataFrame(  ); \nfor cell_type in ['HSC', 'NeuP', 'EryP', 'MasP','MoP']:\n    m = df_meta['cell_type'] == cell_type # 'NeuP'\n    m = m & (df_rna[rna_name] != 0  )\n    v1 = df_y[m][target_name];     \n    v2 = df_rna[m][rna_name];     \n    d.loc[cell_type,'Corr Pearson '+target_name + ' Protein vs RNA'] = np.corrcoef(v1,v2)[0,1] \n    d.loc[cell_type,'Corr Spearman '+target_name + ' Protein vs RNA'] = stats.spearmanr(v1,v2)[0]\n    d.loc[cell_type,'Count'] = m.sum()#  np.corrcoef(v1,v2)[0,1] \n    d.loc[cell_type,'Mean Protein'] = v1.mean() #  np.corrcoef(v1,v2)[0,1] \n    d.loc[cell_type,'Mean RNA'] = v2.mean() #  np.corrcoef(v1,v2)[0,1] \n    d.loc[cell_type,'Count Zero RNA'] = (v2 == 0).sum() #  np.corrcoef(v1,v2)[0,1] \n    d.loc[cell_type,'Percent Zero RNA'] = np.round( (v2 == 0).sum() / len(v2) * 100 , 1) #  np.corrcoef(v1,v2)[0,1] \n    \n    \n#display( d.describe() )\ndisplay(d)\ndisplay(d[ ['Percent Zero RNA','Corr Pearson CD36 Protein vs RNA', 'Count', 'Mean Protein', 'Mean RNA']])","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:27:05.833260Z","iopub.execute_input":"2023-01-06T14:27:05.833578Z","iopub.status.idle":"2023-01-06T14:27:07.925561Z","shell.execute_reply.started":"2023-01-06T14:27:05.833549Z","shell.execute_reply":"2023-01-06T14:27:07.924325Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Correlations with random vector - to estimate p-values","metadata":{}},{"cell_type":"markdown","source":"## Compare correlations of random vectors with protein vector\n\nWe see std is order of magnitude higher for biological version - that requires explanation.\nExplanation - we need to fix cell types, days, donors  - then correlations will have approximately same std. \n\nHere we take random vectors just by gaussian N(0,1) and other option - as random permutation of the protein vector - we get the same results for both - std is around 0.0038, while for biological vector we get std - 0.05 , which is order of magnitude higher. \n\nSimilar results for RNA - we get similar result. ","metadata":{}},{"cell_type":"code","source":"rna_name = ''\nfor t in df_rna.columns:\n    if target_name in t: rna_name = t\nprint(rna_name)","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:27:07.927125Z","iopub.execute_input":"2023-01-06T14:27:07.928163Z","iopub.status.idle":"2023-01-06T14:27:07.939679Z","shell.execute_reply.started":"2023-01-06T14:27:07.928113Z","shell.execute_reply":"2023-01-06T14:27:07.938547Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ny_bio = df_y[target_name].values\ny = df_y[target_name] \nIX = np.random.permutation(len(y))\ny = y.values[IX]\ny2 = np.random.randn( len(df_y) ) #  df_y[target_name]\nl = []; l2 = []; l_bio = []; l_rna = []; l_rna_perm = []\nprint('y.shape', y.shape )\nfor col in df_rna.columns:\n    c = np.corrcoef(y,df_rna[col])[0,1]\n    l.append(c)\n    c = np.corrcoef(y2,df_rna[col])[0,1]\n    l2.append(c)\n    c = np.corrcoef(y_bio,df_rna[col])[0,1]\n    l_bio.append(c)\n    c = np.corrcoef(df_rna[rna_name] ,df_rna[col])[0,1]\n    l_rna.append(c)\n    c = np.corrcoef(df_rna[rna_name].values[IX] ,df_rna[col])[0,1]\n    l_rna_perm.append(c)\n    \n    if (len(l)%5000) == 1:\n        print(len(l))\nd = pd.DataFrame()\nd['Corr Protein ' + target_name ] = l_bio\nd['Corr Random Permutation ' + target_name ] = l2\nd['Corr Random N(0,1) ' ] = l\nd['Corr RNA ' +rna_name ] = l_rna\nd['Corr RNA Permutation ' +rna_name ] = l_rna_perm\n\ndisplay( d.describe() )\nplt.hist(l, bins = 100 )\nplt.hist(l2, bins = 100 )\nplt.hist(l_bio, bins = 100 )\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:27:07.941302Z","iopub.execute_input":"2023-01-06T14:27:07.941647Z","iopub.status.idle":"2023-01-06T14:28:43.928482Z","shell.execute_reply.started":"2023-01-06T14:27:07.941616Z","shell.execute_reply":"2023-01-06T14:28:43.927268Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Even with restriction rna != 0  we get bio-std 5-6 times higher than random - it is too much ","metadata":{}},{"cell_type":"code","source":"%%time\ny_bio = df_y[target_name].values\ny = df_y[target_name] \nIX = np.random.permutation(len(y))\ny = y.values[IX]\ny2 = np.random.randn( len(df_y) ) #  df_y[target_name]\nl = []; l2 = []; l_bio = []; l_rna = []; l_rna_perm = []\nprint('y.shape', y.shape )\nfor col in df_rna.columns:\n    m = df_rna[col] != 0\n    if m.sum() <= 1000: continue\n    v2 = df_rna[col][m]    \n    c = np.corrcoef(y[m],v2)[0,1]\n    l.append(c)\n    c = np.corrcoef(y2[m],v2)[0,1]\n    l2.append(c)\n    c = np.corrcoef(y_bio[m],v2)[0,1]\n    l_bio.append(c)\n    c = np.corrcoef(df_rna[rna_name][m] ,v2)[0,1]\n    l_rna.append(c)\n    c = np.corrcoef(df_rna[rna_name].values[IX][m] ,v2)[0,1]\n    l_rna_perm.append(c)\n    \n    if (len(l)%5000) == 1:\n        print(len(l))\nd = pd.DataFrame()\nd['Corr Protein ' + target_name ] = l_bio\nd['Corr Random Permutation ' + target_name ] = l2\nd['Corr Random N(0,1) ' ] = l\nd['Corr RNA ' +rna_name ] = l_rna\nd['Corr RNA Permutation ' +rna_name ] = l_rna_perm\n\ndisplay( d.describe() )\nplt.hist(l, bins = 100 )\nplt.hist(l2, bins = 100 )\nplt.hist(l_bio, bins = 100 )\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:28:43.929875Z","iopub.execute_input":"2023-01-06T14:28:43.930198Z","iopub.status.idle":"2023-01-06T14:30:26.046484Z","shell.execute_reply.started":"2023-01-06T14:28:43.930169Z","shell.execute_reply":"2023-01-06T14:30:26.044888Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Just correlate two completely random gaussian N(0,1) vectors \n\nGet the same std = 0.0037+- ","metadata":{}},{"cell_type":"code","source":"%%time\nl = []\nfor i in range(10000): # col in df_rna.columns:\n    v1 = np.random.randn( len(y) ) #  df_y[target_name]\n    v2 = np.random.randn( len(y) ) #  df_y[target_name]\n    c = np.corrcoef(v1,v2)[0,1]\n    l.append(c)\n    if (len(l)%5000) == 1:\n        print(len(l))\n        \ndisplay( pd.Series(l).describe() )\nplt.hist(l, bins = 100 )\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:30:26.048041Z","iopub.execute_input":"2023-01-06T14:30:26.048379Z","iopub.status.idle":"2023-01-06T14:31:15.836902Z","shell.execute_reply.started":"2023-01-06T14:30:26.048347Z","shell.execute_reply":"2023-01-06T14:31:15.835597Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Fix real protein expression and correlate with random vectors\n","metadata":{}},{"cell_type":"code","source":"%%time\ny = df_y[target_name] \nl = []\nfor i in range(10000): # col in df_rna.columns:\n    v = np.random.randn( len(y) ) #  df_y[target_name]\n    c = np.corrcoef(y,v)[0,1]\n    l.append(c)\n    if (len(l)%5000) == 1:\n        print(len(l))\n        \ndisplay( pd.Series(l).describe() )\nplt.hist(l, bins = 100 )\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:31:15.838610Z","iopub.execute_input":"2023-01-06T14:31:15.839489Z","iopub.status.idle":"2023-01-06T14:31:45.120009Z","shell.execute_reply.started":"2023-01-06T14:31:15.839422Z","shell.execute_reply":"2023-01-06T14:31:45.118894Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Correlate random vector with ALL rna data\n\nStd of the result is similar to what we get above - so probably correspond to random vectors \n","metadata":{}},{"cell_type":"code","source":"%%time\nl = []\ny = np.random.randn( len(df_y) ) #  df_y[target_name]\nfor col in df_rna.columns:\n    c = np.corrcoef(y,df_rna[col])[0,1]\n    l.append(c)\n    if (len(l)%5000) == 1:\n        print(len(l))\n        \ndisplay( pd.Series(l).describe() )\nplt.hist(l, bins = 100 )\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:31:45.121113Z","iopub.execute_input":"2023-01-06T14:31:45.121401Z","iopub.status.idle":"2023-01-06T14:32:03.644607Z","shell.execute_reply.started":"2023-01-06T14:31:45.121374Z","shell.execute_reply":"2023-01-06T14:32:03.643813Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"l2 = np.abs(l)\n\n\nfor t in [0.005,0.01, 0.015, ]:\n    m = l2 > t\n    print(t, m.sum(), np.round( m.sum()/len(df_stat) * 100,4)  )","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:32:03.646119Z","iopub.execute_input":"2023-01-06T14:32:03.647312Z","iopub.status.idle":"2023-01-06T14:32:03.659106Z","shell.execute_reply.started":"2023-01-06T14:32:03.647268Z","shell.execute_reply":"2023-01-06T14:32:03.657927Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Somewhat a bit faster way to calculate correlations - by normalization and multiplications\n\nNo a big improvement","metadata":{}},{"cell_type":"code","source":"%%time\nfrom sklearn.preprocessing import StandardScaler\nscaler = StandardScaler()\nX = df_rna\nX = scaler.fit_transform(X)\nX = X/np.sqrt(len(X))\nX\n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:32:03.660548Z","iopub.execute_input":"2023-01-06T14:32:03.661090Z","iopub.status.idle":"2023-01-06T14:32:30.650373Z","shell.execute_reply.started":"2023-01-06T14:32:03.661041Z","shell.execute_reply":"2023-01-06T14:32:30.649210Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nl = []\nfor i in range(10):\n    v = np.random.randn(len(df_y))\n    v = (v-np.mean(v)) /np.std(v) /np.sqrt(len(v))\n    c = np.matmul(v, X ) \n    l += (list(np.ravel(c)) )\nprint(len(l))\ndisplay( pd.Series(l).describe() )\nplt.hist(l, bins = 100 )\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:32:30.651849Z","iopub.execute_input":"2023-01-06T14:32:30.652182Z","iopub.status.idle":"2023-01-06T14:34:34.096750Z","shell.execute_reply.started":"2023-01-06T14:32:30.652152Z","shell.execute_reply":"2023-01-06T14:34:34.095650Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"l2 = np.abs(l)\nfor t in [0.005,0.01, 0.015, 0.02]:\n    m = l2 > t\n    print(t, m.sum(), np.round( m.sum()/len(l2) * 100,4)  )","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:34:34.098195Z","iopub.execute_input":"2023-01-06T14:34:34.098544Z","iopub.status.idle":"2023-01-06T14:34:34.124296Z","shell.execute_reply.started":"2023-01-06T14:34:34.098509Z","shell.execute_reply":"2023-01-06T14:34:34.123223Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# So std = 0.0037\n# 3*std = 0.01 \n# so 10* (3*std) = \n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:34:34.126021Z","iopub.execute_input":"2023-01-06T14:34:34.126450Z","iopub.status.idle":"2023-01-06T14:34:34.130537Z","shell.execute_reply.started":"2023-01-06T14:34:34.126415Z","shell.execute_reply":"2023-01-06T14:34:34.129506Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:34:34.132266Z","iopub.execute_input":"2023-01-06T14:34:34.132659Z","iopub.status.idle":"2023-01-06T14:34:34.167817Z","shell.execute_reply.started":"2023-01-06T14:34:34.132626Z","shell.execute_reply":"2023-01-06T14:34:34.166789Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Analysis on segments - cell types, days , donors ","metadata":{}},{"cell_type":"code","source":"%%time\nfn = '/kaggle/input/open-problems-multimodal/metadata.csv'\ndf_meta = pd.read_csv(fn, index_col = 0 )\ndf_meta\n# Cut only train cite-seq part: \nd = pd.DataFrame(index = df_y.index)\nprint(d.shape)\ndf_meta = d.join(df_meta, how = 'left')\ndf_meta","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:34:34.169179Z","iopub.execute_input":"2023-01-06T14:34:34.169523Z","iopub.status.idle":"2023-01-06T14:34:34.580671Z","shell.execute_reply.started":"2023-01-06T14:34:34.169490Z","shell.execute_reply":"2023-01-06T14:34:34.579439Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Cut only train cite-seq part: \nd = pd.DataFrame(index = df_y.index)\nprint(d.shape)\ndf_meta = d.join(df_meta, how = 'left')\ndf_meta","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:34:34.582712Z","iopub.execute_input":"2023-01-06T14:34:34.583536Z","iopub.status.idle":"2023-01-06T14:34:34.601958Z","shell.execute_reply.started":"2023-01-06T14:34:34.583492Z","shell.execute_reply":"2023-01-06T14:34:34.600755Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## HSC","metadata":{}},{"cell_type":"code","source":"m = df_meta['cell_type'] == 'HSC'\nm = m & ( df_meta['day'] == 2 )\nm = m & ( df_meta['donor'] == 32606 )\nprint(m.sum())","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:34:34.603663Z","iopub.execute_input":"2023-01-06T14:34:34.604603Z","iopub.status.idle":"2023-01-06T14:34:34.619175Z","shell.execute_reply.started":"2023-01-06T14:34:34.604556Z","shell.execute_reply":"2023-01-06T14:34:34.617959Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y = df_y[m][target_name]\nX = df_rna[m]\nprint(X.shape, y.shape)","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:34:34.620458Z","iopub.execute_input":"2023-01-06T14:34:34.620774Z","iopub.status.idle":"2023-01-06T14:34:34.963219Z","shell.execute_reply.started":"2023-01-06T14:34:34.620738Z","shell.execute_reply":"2023-01-06T14:34:34.962446Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"l = []\nfor col in X.columns:\n    c = np.corrcoef(y,X[col])[0,1]\n    l.append(c)\n    if (len(l)%5000) == 1:\n        print(len(l))\n        \ndisplay( pd.Series(l).describe() )\nplt.hist(l, bins = 100 )\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:34:34.964345Z","iopub.execute_input":"2023-01-06T14:34:34.964840Z","iopub.status.idle":"2023-01-06T14:34:41.281672Z","shell.execute_reply.started":"2023-01-06T14:34:34.964784Z","shell.execute_reply":"2023-01-06T14:34:41.280512Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nl = []\nfor i in range(10000): # col in df_rna.columns:\n    v1 = np.random.randn( len(y) ) #  df_y[target_name]\n    v2 = np.random.randn( len(y) ) #  df_y[target_name]\n    c = np.corrcoef(v1,v2)[0,1]\n    l.append(c)\n    if (len(l)%5000) == 1:\n        print(len(l))\n        \ndisplay( pd.Series(l).describe() )\nplt.hist(l, bins = 100 )\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:34:41.283125Z","iopub.execute_input":"2023-01-06T14:34:41.283505Z","iopub.status.idle":"2023-01-06T14:34:45.227145Z","shell.execute_reply.started":"2023-01-06T14:34:41.283447Z","shell.execute_reply":"2023-01-06T14:34:45.225589Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_meta['cell_type'].value_counts(), df_meta['day'].value_counts(), ","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:34:45.229039Z","iopub.execute_input":"2023-01-06T14:34:45.229538Z","iopub.status.idle":"2023-01-06T14:34:45.245804Z","shell.execute_reply.started":"2023-01-06T14:34:45.229490Z","shell.execute_reply":"2023-01-06T14:34:45.244559Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## All cell types","metadata":{}},{"cell_type":"code","source":"%%time\nfor day in [2,3,4]:\n    d = pd.DataFrame(  ); \n    for cell_type in ['HSC', 'NeuP', 'EryP', 'MasP','MoP']:\n        m = df_meta['cell_type'] == cell_type # 'NeuP'\n        m = m & ( df_meta['day'] == day ); m = m & ( df_meta['donor'] == 32606 )\n\n        y = df_y[m][target_name];     X = df_rna[m]\n        #print(m.sum(), X.shape, y.shape)\n\n        l_bio = [];    l_random = []\n        for col in X.columns:\n            l_bio.append(  np.corrcoef(y,X[col])[0,1] )\n            l_random.append(  np.corrcoef(np.random.randn( len(y) ), np.random.randn( len(y) )  )[0,1] )\n            #if (len(l)%5000) == 1:         print(len(l))\n        d['Corr Bio ' +cell_type] = l_bio; d['Corr Random '+cell_type] = l_random   \n    print('day:',day)\n    display( d.describe() )\n\n\n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:34:45.247113Z","iopub.execute_input":"2023-01-06T14:34:45.247487Z","iopub.status.idle":"2023-01-06T14:37:32.310644Z","shell.execute_reply.started":"2023-01-06T14:34:45.247439Z","shell.execute_reply":"2023-01-06T14:37:32.309459Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# plt.hist(d, bins = 100 )\n# #plt.hist(l_random, bins = 100 )\n# plt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:37:32.312400Z","iopub.execute_input":"2023-01-06T14:37:32.313616Z","iopub.status.idle":"2023-01-06T14:37:32.318701Z","shell.execute_reply.started":"2023-01-06T14:37:32.313567Z","shell.execute_reply":"2023-01-06T14:37:32.317350Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Vary only cell type ","metadata":{}},{"cell_type":"code","source":"%%time\nd = pd.DataFrame(  ); \nfor cell_type in ['HSC', 'NeuP', 'EryP', 'MasP','MoP']:\n#for day in [2,3,4]:\n    m = pd.Series(index = df_meta.index, data = True )\n    m = df_meta['cell_type'] == cell_type # 'NeuP'\n    #m = m & ( df_meta['day'] == day ); \n    #m = m & ( df_meta['donor'] == donor )\n\n    y = df_y[m][target_name];     X = df_rna[m]\n    #print(m.sum(), X.shape, y.shape)\n\n    l_bio = [];    l_random = []\n    for col in X.columns:\n        l_bio.append(  np.corrcoef(y,X[col])[0,1] )\n        l_random.append(  np.corrcoef(np.random.randn( len(y) ), np.random.randn( len(y) )  )[0,1] )\n        #if (len(l)%5000) == 1:         print(len(l))\n    postfix = str(cell_type)\n    d['Corr Bio ' +postfix] = l_bio; d['Corr Random '+postfix] = l_random   \n#print('day:',day)\ndisplay( d.describe() )\n\n\n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:37:32.319994Z","iopub.execute_input":"2023-01-06T14:37:32.320322Z","iopub.status.idle":"2023-01-06T14:40:22.261922Z","shell.execute_reply.started":"2023-01-06T14:37:32.320292Z","shell.execute_reply":"2023-01-06T14:40:22.260765Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Vary only day","metadata":{}},{"cell_type":"code","source":"%%time\nd = pd.DataFrame(  ); \n#for cell_type in ['HSC', 'NeuP', 'EryP', 'MasP','MoP']:\nfor day in [2,3,4]:\n    m = pd.Series(index = df_meta.index, data = True )\n    #m = df_meta['cell_type'] == cell_type # 'NeuP'\n    m = m & ( df_meta['day'] == day ); \n    #m = m & ( df_meta['donor'] == donor )\n\n    y = df_y[m][target_name];     X = df_rna[m]\n    #print(m.sum(), X.shape, y.shape)\n\n    l_bio = [];    l_random = []\n    for col in X.columns:\n        l_bio.append(  np.corrcoef(y,X[col])[0,1] )\n        l_random.append(  np.corrcoef(np.random.randn( len(y) ), np.random.randn( len(y) )  )[0,1] )\n        #if (len(l)%5000) == 1:         print(len(l))\n    postfix = str(day)\n    d['Corr Bio ' +postfix] = l_bio; d['Corr Random '+postfix] = l_random   \n#print('day:',day)\ndisplay( d.describe() )\n\n\n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:40:22.263839Z","iopub.execute_input":"2023-01-06T14:40:22.264449Z","iopub.status.idle":"2023-01-06T14:43:01.160168Z","shell.execute_reply.started":"2023-01-06T14:40:22.264414Z","shell.execute_reply":"2023-01-06T14:43:01.158887Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Vary only donor ","metadata":{}},{"cell_type":"code","source":"df_meta['donor'].value_counts()\ndf_meta['donor'].unique()\n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:43:01.161635Z","iopub.execute_input":"2023-01-06T14:43:01.162519Z","iopub.status.idle":"2023-01-06T14:43:01.173743Z","shell.execute_reply.started":"2023-01-06T14:43:01.162480Z","shell.execute_reply":"2023-01-06T14:43:01.172221Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n#for day in [2,3,4]:\nd = pd.DataFrame(  ); \n# for cell_type in ['HSC', 'NeuP', 'EryP', 'MasP','MoP']:\n# for day in [2,3,4]:\nfor donor in [32606, 13176, 31800]:    \n    m = pd.Series(index = df_meta.index, data = True )\n    \n    #m = df_meta['cell_type'] == cell_type # 'NeuP'\n    #m = m & ( df_meta['day'] == day ); \n    m = m & ( df_meta['donor'] == donor )\n\n    y = df_y[m][target_name];     X = df_rna[m]\n    #print(m.sum(), X.shape, y.shape)\n\n    l_bio = [];    l_random = []\n    for col in X.columns:\n        l_bio.append(  np.corrcoef(y,X[col])[0,1] )\n        l_random.append(  np.corrcoef(np.random.randn( len(y) ), np.random.randn( len(y) )  )[0,1] )\n        #if (len(l)%5000) == 1:         print(len(l))\n    postfix = str(donor)        \n    d['Corr Bio ' +postfix] = l_bio; d['Corr Random '+postfix] = l_random   \n#print('day:',day)\ndisplay( d.describe() )\n\n\n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:43:01.175793Z","iopub.execute_input":"2023-01-06T14:43:01.177180Z","iopub.status.idle":"2023-01-06T14:45:40.108214Z","shell.execute_reply.started":"2023-01-06T14:43:01.177143Z","shell.execute_reply":"2023-01-06T14:45:40.107002Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Estimate how many correlations are not random ","metadata":{"execution":{"iopub.status.busy":"2022-12-29T21:40:30.690231Z","iopub.execute_input":"2022-12-29T21:40:30.691365Z","iopub.status.idle":"2022-12-29T21:40:30.696521Z","shell.execute_reply.started":"2022-12-29T21:40:30.691319Z","shell.execute_reply":"2022-12-29T21:40:30.694884Z"}}},{"cell_type":"markdown","source":"## Cell types for fixed day , donor ","metadata":{}},{"cell_type":"code","source":"%%time\nfor day in [2,3,4]:\n    d = pd.DataFrame(  ); \n    for cell_type in ['HSC', 'NeuP', 'EryP', 'MasP','MoP']:\n        m = df_meta['cell_type'] == cell_type # 'NeuP'\n        m = m & ( df_meta['day'] == day ); \n        m = m & ( df_meta['donor'] == 32606 )\n\n        y = df_y[m][target_name];     X = df_rna[m]\n        #print(m.sum(), X.shape, y.shape)\n\n        l_bio = [];    l_random = []; l_random_perm = []\n        for col in X.columns:\n            l_bio.append(  np.corrcoef(y,X[col])[0,1] )\n            l_random.append(  np.corrcoef(np.random.randn( len(y) ), np.random.randn( len(y) )  )[0,1] )\n            l_random_perm.append(  np.corrcoef(y.values[np.random.permutation(len(y)) ] ,X[col])[0,1] )\n            #if (len(l)%5000) == 1:         print(len(l))\n        d['Corr Bio ' +cell_type] = l_bio; d['Corr Random '+cell_type] = l_random ; d['Corr Random Perm '+cell_type] = l_random_perm ; \n    print('day:',day)\n    display( d.describe() )\n\n    for i in range(int( d.shape[1]/3) ) :\n        col1 = d.columns[3*i]\n        col2 = d.columns[3*i+1]\n        col3 = d.columns[3*i+2]\n        s1 = d[col1].std(); s2 = d[col2].std(); s3 = d[col3].std(); \n        #print(col1)\n        m = d[col1]>3.5*s2\n        print(col1, 'Beyond Random (3.5 and 4)*std:', (d[col1]>3.5*s2).sum() , (d[col1]>4*s2).sum(), 'std1 = %.4f , std2 = %.4f, std3 = %.4f '%(s1,s2,s3) )\n        print(col1, 'Beyond Random Perm (3.5 and 4)*std   :', (d[col1]>3.5*s3).sum() , (d[col1]>4*s3).sum())# , 'std1 = %.4f , std2 = %.4f '%(s1,s2) )\n        print(col1, 'Beyond Bio (3.5 and 4)*std   :', (d[col1]>3.5*s1).sum() , (d[col1]>4*s1).sum())# , 'std1 =  '%(s1,s2) )\n\n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:45:40.110049Z","iopub.execute_input":"2023-01-06T14:45:40.110378Z","iopub.status.idle":"2023-01-06T14:49:41.422267Z","shell.execute_reply.started":"2023-01-06T14:45:40.110348Z","shell.execute_reply":"2023-01-06T14:49:41.421097Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.hist(d, bins = 100 )\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:49:41.423572Z","iopub.execute_input":"2023-01-06T14:49:41.423892Z","iopub.status.idle":"2023-01-06T14:49:43.952544Z","shell.execute_reply.started":"2023-01-06T14:49:41.423855Z","shell.execute_reply":"2023-01-06T14:49:43.951327Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Correlation Protein/RNA statistics by cell type ","metadata":{}},{"cell_type":"code","source":"rna_name = ''\nfor t in df_rna.columns:\n    if target_name in t: rna_name = t\nprint(rna_name)","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:49:43.953710Z","iopub.execute_input":"2023-01-06T14:49:43.954010Z","iopub.status.idle":"2023-01-06T14:49:43.965167Z","shell.execute_reply.started":"2023-01-06T14:49:43.953983Z","shell.execute_reply":"2023-01-06T14:49:43.963894Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfrom scipy import stats\n\nd = pd.DataFrame(  ); \nfor cell_type in ['HSC', 'NeuP', 'EryP', 'MasP','MoP']:\n    m = df_meta['cell_type'] == cell_type # 'NeuP'\n    v1 = df_y[m][target_name];     \n    v2 = df_rna[m][rna_name];     \n    d.loc[cell_type,'Corr Pearson '+target_name + ' Protein vs RNA'] = np.corrcoef(v1,v2)[0,1] \n    d.loc[cell_type,'Corr Spearman '+target_name + ' Protein vs RNA'] = stats.spearmanr(v1,v2)[0]\n    d.loc[cell_type,'Count'] = m.sum()#  np.corrcoef(v1,v2)[0,1] \n    d.loc[cell_type,'Mean Protein'] = v1.mean() #  np.corrcoef(v1,v2)[0,1] \n    d.loc[cell_type,'Mean RNA'] = v2.mean() #  np.corrcoef(v1,v2)[0,1] \n    d.loc[cell_type,'Count Zero RNA'] = (v2 == 0).sum() #  np.corrcoef(v1,v2)[0,1] \n    d.loc[cell_type,'Percent Zero RNA'] = np.round( (v2 == 0).sum() / len(v2) * 100 , 1) #  np.corrcoef(v1,v2)[0,1] \n    \n#display( d.describe() )\ndisplay(d)","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:49:43.967159Z","iopub.execute_input":"2023-01-06T14:49:43.967673Z","iopub.status.idle":"2023-01-06T14:49:51.877086Z","shell.execute_reply.started":"2023-01-06T14:49:43.967630Z","shell.execute_reply":"2023-01-06T14:49:51.876320Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_meta['donor'].unique()","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:49:51.878642Z","iopub.execute_input":"2023-01-06T14:49:51.878949Z","iopub.status.idle":"2023-01-06T14:49:51.887148Z","shell.execute_reply.started":"2023-01-06T14:49:51.878921Z","shell.execute_reply":"2023-01-06T14:49:51.886034Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nd = pd.DataFrame(  ); \nd2 = pd.DataFrame(  );\nfor donor in [32606, 13176, 31800]:\n    for day in [2,3,4]:\n        for cell_type in ['EryP', ]: #  'HSC', 'NeuP', 'EryP', 'MasP','MoP']:\n            m = df_meta['cell_type'] == cell_type # 'NeuP'\n            m = m & ( df_meta['day'] == day ); \n            m = m & ( df_meta['donor'] == donor )\n\n            v1 = df_y[m][target_name];     \n            v2 = df_rna[m][rna_name];     \n            d.loc[str(donor) + ' ' + str(day),cell_type + ' Corr Pearson '+target_name + ' Protein vs RNA'] = np.corrcoef(v1,v2)[0,1] \n            d.loc[str(donor) + ' ' + str(day),cell_type + ' Count'] = m.sum()\n            d.loc[str(donor) + ' ' + str(day),cell_type + ' Mean Protein'] = v1.mean()\n            d.loc[str(donor) + ' ' + str(day),cell_type + ' Mean RNA'] = v2.mean()\n            \n            d2.loc[donor, day] = np.corrcoef(v1,v2)[0,1] \ndisplay( d )\ndisplay(d2)\nd.describe() ","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:49:51.889730Z","iopub.execute_input":"2023-01-06T14:49:51.890451Z","iopub.status.idle":"2023-01-06T14:49:53.701564Z","shell.execute_reply.started":"2023-01-06T14:49:51.890407Z","shell.execute_reply":"2023-01-06T14:49:53.700136Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# So correlation protein-rna decreases with the day (on EryP), but in different manner for different donors  - reasons are unclear ","metadata":{}},{"cell_type":"code","source":"print('Correlations RNA - Protein for ', target_name , 'for days and donors: ')\ndisplay(d2)\n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T14:58:41.361546Z","iopub.execute_input":"2023-01-06T14:58:41.362200Z","iopub.status.idle":"2023-01-06T14:58:41.375701Z","shell.execute_reply.started":"2023-01-06T14:58:41.362156Z","shell.execute_reply":"2023-01-06T14:58:41.374296Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}