{"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\n### Version 1,2: \n\n    CD36: calculate average correlations for segements related to donor and day for fixed CELL TYPE = EryP (since that is the cell type where CD36 is mainly expressed). \n    \n    Order genes by such averaged correlation.\n    \n    Make enricments for top100,top200 ordered by Abs/positive/negative correlations.\n    \n    Correlation CD36 with rnas decreases with the day - systematically for all donors. How that can be explained ? \n    CD36 is getting. CD36 is increasing with days for EryP - is is unclear why, and it is more strange that correlation decreases.  Also it is systematiclly higher for EryP of one donor 32606 - also unclear why \n\n\n### Attempt to understand how many genes are correlated in non random way\n    \n    See the section: \"# Attempts to estimations number of non-random correlations by Bonferoni correction \n\"\n## Analysis of interesetion of All EryP segments, and additionaly intersection with All cell , and MoP\n\n    See section \"# Add analysis of correlation over all cells and over monocytes progenitors \"\n\n    All EryP intersection of topK varies 14-20 percent\n    Intersection with MoP for top First - is only about 10% , but later grows to 20%\n    Interesection EryP1 with topK for all cells is about 45% \n    \n    On MoP the sign of some correlations changes comparing to EryP \n    \n    The \"All cells\" are more or less consistent with EryP - which is quite good for us   - taking in mind generalization to other datasets possibly without split to cell types. \n    \n    We observe phenomena of stabilization of the slopes of the curve - that might give a choice for selection of genes. \n    We can take stabilized linear curve - and take that many genes where that linear approximation is good - the deviation point from \n    the linear approximation is stopping criteria. \n    \n    Saved results with top100 genes which are in the intersection of topN for large N of the EryP for 9 segments: days & donors.\n    The ordering of that list is obtained by median correlation for each segment.  \n    See subsection \"## Save results \" and https://www.kaggle.com/code/alexandervc/mmscel-correlations-cd-vs-rna-v2?scriptVersionId=115287304&cellId=95\n    \n## Analysis of correlations on random subsamples\n\n    See section '# Correlations on random subsamples '\n    \n    For random subsamples we see quite great concordance between correlations - much more than for splits of EryP by days & donors \n    \n    Even if take 10% size of subsamples - the average correlation would be 0.86 within results on different subsamples\n\n    When we make 50 trials we see that slope of the curve is stabilized with 0.48. \n    But with 100 trials it goes to 0.43. \n    The point of deviation from linear trend also changes - 400, 300, 200. \n    Thus the stabilization seems not quite exist. \n    Version 13 - contains 100 trials \n    \n  \n    Results with top100 saved:\n    https://www.kaggle.com/code/alexandervc/mmscel-correlations-cd-vs-rna-v2/data?scriptVersionId=115305377\n    Another version with intersection over 50 random sabsamples\n    https://www.kaggle.com/code/alexandervc/mmscel-correlations-cd-vs-rna-v2?scriptVersionId=115310325&cellId=109    \n    \n    V15 - run with 250 trials\n    V16 - run with 500 trials\n    V17 - run with 300 trials\n    v18 - run with 10 trials again \n    \n","metadata":{"execution":{"iopub.status.busy":"2022-12-30T14:32:51.115380Z","iopub.execute_input":"2022-12-30T14:32:51.115913Z","iopub.status.idle":"2022-12-30T14:32:51.124048Z","shell.execute_reply.started":"2022-12-30T14:32:51.115845Z","shell.execute_reply":"2022-12-30T14:32:51.122397Z"}}},{"cell_type":"markdown","source":"# Key params","metadata":{}},{"cell_type":"code","source":"target_name = 'CD36'\n\n#n_top = 100 ","metadata":{"execution":{"iopub.status.busy":"2023-01-06T22:23:03.610320Z","iopub.execute_input":"2023-01-06T22:23:03.610824Z","iopub.status.idle":"2023-01-06T22:23:03.615553Z","shell.execute_reply.started":"2023-01-06T22:23:03.610786Z","shell.execute_reply":"2023-01-06T22:23:03.614748Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Preparations","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-06T22:23:05.192292Z","iopub.execute_input":"2023-01-06T22:23:05.193062Z","iopub.status.idle":"2023-01-06T22:23:05.223064Z","shell.execute_reply.started":"2023-01-06T22:23:05.193022Z","shell.execute_reply":"2023-01-06T22:23:05.222311Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport seaborn as sns","metadata":{"execution":{"iopub.status.busy":"2023-01-06T22:23:05.713148Z","iopub.execute_input":"2023-01-06T22:23:05.713961Z","iopub.status.idle":"2023-01-06T22:23:06.418302Z","shell.execute_reply.started":"2023-01-06T22:23:05.713918Z","shell.execute_reply":"2023-01-06T22:23:06.417219Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Explore some idea\n\nConsider bootstrap over subsamples - the mean appears to be the same as original correlation,  thus we cannot use bootstrap over random subsamples to better estimate the correlation. \n\nWhat will be done in the present notebook - non random subsamples related to donor and day ","metadata":{}},{"cell_type":"code","source":"%%time\nN = 10_000\nv1 = np.random.randn(N)\nv2 = np.random.randn(N)\nc = np.corrcoef(v1,v2)[0,1]\nprint(c)\n\nN2 = int(N*0.75)\nl = []\nfor k in range(1000):\n    IX = np.random.permutation(N)[:N2]\n    c = np.corrcoef(v1[IX], v2[IX])[0,1]\n    l.append(c)\ndisplay( pd.Series(l).describe()  )\nplt.hist(l)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-06T22:23:06.867607Z","iopub.execute_input":"2023-01-06T22:23:06.868029Z","iopub.status.idle":"2023-01-06T22:23:07.624029Z","shell.execute_reply.started":"2023-01-06T22:23:06.867993Z","shell.execute_reply":"2023-01-06T22:23:07.622356Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load data","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#%%time\ndf_y = pd.read_hdf('/kaggle/input/open-problems-multimodal/train_cite_targets.h5')\ndisplay(df_y)\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-06T22:23:10.141015Z","iopub.execute_input":"2023-01-06T22:23:10.141568Z","iopub.status.idle":"2023-01-06T22:24:06.932851Z","shell.execute_reply.started":"2023-01-06T22:23:10.141518Z","shell.execute_reply":"2023-01-06T22:24:06.931626Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_meta['donor'].unique()","metadata":{"execution":{"iopub.status.busy":"2023-01-06T22:32:15.968517Z","iopub.execute_input":"2023-01-06T22:32:15.969622Z","iopub.status.idle":"2023-01-06T22:32:15.979340Z","shell.execute_reply.started":"2023-01-06T22:32:15.969579Z","shell.execute_reply":"2023-01-06T22:32:15.978419Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Look on mean over cell types for protein levels  and count non-zero for rna \n\nThe main types where CD36 is actvive are EryP and MoP, but MoP has only 591 cell. \n\nFor the other cells mean expression is low and well rna has huge percent of zero elements ","metadata":{}},{"cell_type":"code","source":"d = df_y.copy()\nd = d.join(df_meta) \nprint('Counts:')\ndisplay(d['cell_type'].value_counts() )\nprint()\n\n#m = d['cell_type'] == 'EryP'\nprint('Mean protein:')\ndisplay( d.groupby('cell_type')[target_name].mean().sort_values() )\n#display( d[m].groupby('donor')[target_name].mean() )\n\nrna_name = ''\nfor t in df_rna.columns:\n    if target_name in t: rna_name = t\nprint(rna_name)\nprint()\nd['rna is non zero'] = (df_rna[rna_name] != 0)*100\nprint('Percent of non-zeros RNA:')\ndisplay( d.groupby('cell_type')['rna is non zero'].mean().sort_values() )\n\nprint()\nprint('Mean rna:')\nd['rna'] = df_rna[rna_name]# != 0\ndisplay( d.groupby('cell_type')['rna'].mean().sort_values() )\n\n\n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T22:32:19.973292Z","iopub.execute_input":"2023-01-06T22:32:19.973779Z","iopub.status.idle":"2023-01-06T22:32:20.103826Z","shell.execute_reply.started":"2023-01-06T22:32:19.973736Z","shell.execute_reply":"2023-01-06T22:32:20.102564Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## CD36 is increasing with days for EryP, and is systematiclly higher for EryP of one donor 32606 - both unclear why , however when look by pairs day x donor - the trends are not that much clear \n","metadata":{}},{"cell_type":"code","source":"d = df_y.copy()\nd = d.join(df_meta) \nm = d['cell_type'] == 'EryP'\ndisplay( d[m].groupby('day')[target_name].mean() )\ndisplay( d[m].groupby('donor')[target_name].mean() )\ndisplay( d[m].groupby(['donor', 'day'])[target_name].mean() )\n\n\n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T22:32:28.970656Z","iopub.execute_input":"2023-01-06T22:32:28.971066Z","iopub.status.idle":"2023-01-06T22:32:29.078293Z","shell.execute_reply.started":"2023-01-06T22:32:28.971033Z","shell.execute_reply":"2023-01-06T22:32:29.077081Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Calculate correlations for several segements separately","metadata":{"execution":{"iopub.status.busy":"2022-12-30T13:53:33.535624Z","iopub.execute_input":"2022-12-30T13:53:33.535991Z","iopub.status.idle":"2022-12-30T13:53:33.541258Z","shell.execute_reply.started":"2022-12-30T13:53:33.535960Z","shell.execute_reply":"2022-12-30T13:53:33.540275Z"}}},{"cell_type":"code","source":"%%time\ndict_corrs = {}\nd_bio = pd.DataFrame(index = df_rna.columns)\nd_rand_perm = pd.DataFrame(index = df_rna.columns)\n\nfor cell_type in ['EryP']:#  ['HSC', 'NeuP', 'EryP', 'MasP','MoP']:\n    for donor in [32606, 13176, 31800]: \n        for day in [2,3,4]:\n            #d = pd.DataFrame(  ); \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];     \n            X = df_rna[m]\n            #print(m.sum(), X.shape, y.shape)\n            l_bio = []; l_random_perm = []\n            for col in X.columns:\n                l_bio.append(  np.corrcoef(y,X[col])[0,1] )\n                l_random_perm.append(  np.corrcoef(y.values[np.random.permutation(len(y)) ] ,X[col])[0,1] )\n                if (len(l_bio)%5000) == 1:         print(len(l_bio))\n                \n            d_bio[str(day)+' '+str(donor) + ' ' + str(cell_type) ] = l_bio\n            d_rand_perm[str(day)+' '+str(donor) + ' ' + str(cell_type) ] = l_random_perm\n            \n            \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-06T22:33:13.030828Z","iopub.execute_input":"2023-01-06T22:33:13.031234Z","iopub.status.idle":"2023-01-06T22:34:57.032025Z","shell.execute_reply.started":"2023-01-06T22:33:13.031202Z","shell.execute_reply":"2023-01-06T22:34:57.030551Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(); print('d_bio')\ndisplay(d_bio.sort_values(d_bio.columns[0], ascending = False).head(20) )\n\nprint(); print('d_bio')\ndisplay(d_bio)\nprint(); print('d_rand_perm')\ndisplay(d_rand_perm)\nprint(); print('d_bio')\ndisplay(d_bio.describe())\ndisplay(d_bio.describe().mean(axis = 1 ))\nprint(); print('d_rand_perm')\ndisplay(d_rand_perm.describe())\ndisplay(d_rand_perm.describe().mean(axis = 1 ))\n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T22:41:12.901468Z","iopub.execute_input":"2023-01-06T22:41:12.902660Z","iopub.status.idle":"2023-01-06T22:41:13.154563Z","shell.execute_reply.started":"2023-01-06T22:41:12.902598Z","shell.execute_reply":"2023-01-06T22:41:13.153389Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"(-d_bio).rank().sort_values(d_bio.columns[0], ascending = True ).head(20)","metadata":{"execution":{"iopub.status.busy":"2023-01-06T22:34:57.202450Z","iopub.execute_input":"2023-01-06T22:34:57.202928Z","iopub.status.idle":"2023-01-06T22:34:57.266725Z","shell.execute_reply.started":"2023-01-06T22:34:57.202882Z","shell.execute_reply":"2023-01-06T22:34:57.265571Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Most stably related to CD36 genes ","metadata":{}},{"cell_type":"code","source":"d3 = (-d_bio).rank().sort_values(d_bio.columns[0], ascending = True )\nd4 = (-d_bio).rank().sort_values(d_bio.columns[0], ascending = True )\nd4['Max Rank'] = d3.max(axis = 1)\nd4 = d4.sort_values('Max Rank')\nd4.head(40)","metadata":{"execution":{"iopub.status.busy":"2023-01-02T19:51:42.961971Z","iopub.execute_input":"2023-01-02T19:51:42.962496Z","iopub.status.idle":"2023-01-02T19:51:43.088879Z","shell.execute_reply.started":"2023-01-02T19:51:42.962447Z","shell.execute_reply":"2023-01-02T19:51:43.087634Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install gseapy\nimport gseapy as gp","metadata":{"execution":{"iopub.status.busy":"2023-01-02T19:51:43.090488Z","iopub.execute_input":"2023-01-02T19:51:43.091004Z","iopub.status.idle":"2023-01-02T19:51:57.719742Z","shell.execute_reply.started":"2023-01-02T19:51:43.090955Z","shell.execute_reply":"2023-01-02T19:51:57.718109Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"N = 400\nlist_genes = d4.index[:N]\nlist_genes = [t.split('_')[1] for t in list_genes]\nprint(len(list_genes ))\n\nprint(list_genes[:100])\n\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) )\ndisplay( enr.results.sort_values('Adjusted P-value').head(20) )","metadata":{"execution":{"iopub.status.busy":"2023-01-02T19:51:57.728386Z","iopub.execute_input":"2023-01-02T19:51:57.729161Z","iopub.status.idle":"2023-01-02T19:51:59.798784Z","shell.execute_reply.started":"2023-01-02T19:51:57.729099Z","shell.execute_reply":"2023-01-02T19:51:59.797585Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Attempts to estimations number of non-random correlations by Bonferoni correction \n\nWe have about  20 000 genes. We need Bonferoni correction. \n\n\nRoughly speaking each 0.5 sigma gives order of magnitude growth in p-value\n(see https://en.wikipedia.org/wiki/68%E2%80%9395%E2%80%9399.7_rule#Table_of_numerical_values )\n\n    μ ± 3.5σ 4.653E-04\n    μ ± 4σ   6.334E-05\n    μ ± 4.5σ 6.795E-06 \n    μ ± 5σ   5.733E-07 \n    μ ± 5.5σ 3.798E-08\n    μ ± 6σ   1.973E-09 \n    μ ± 6.5σ 8.032E-11\n    \nSo values beyond μ ± 4.5σ slighly not enough for 20 000 correction,\nAnd values beyond μ ± 5σ MORE THAN ENOGUH for 20 000 correction \n\nThe question is what SIGMA to consider - if take sigma from purely random vectors (or with random permutation),\nsuch sigma is very far from the sigma obtained from real correlations of CD36 to all RNA, that is somewhat puzzling, but probably is indication that there are indeed many genes which correspond to powerful differentiation process of EryP cells from stem like state to more differentiated state. CD36 might be in line with these process. \n\nHowever for the moment it is not that much clear, and we will give number of genes which goes K * sigma for both sigma from random vectors and sigma from actual distribution of correlations. \n\n\nIf we conservatively  μ ± 5σ  with uncorrected p-value 5.733E-07 ,  for sigma from random sample we will get\n266 - 566 genes, while taking sigma from actual sample we get 44 - 78 genes\n\n\n","metadata":{}},{"cell_type":"code","source":"dt1 = d_rand_perm.describe()\ndt2 = d_bio.describe()\n\nd = pd.DataFrame()\nd.index.name = 'K * sgima, K:'\ni = 0\ns1 = dt1.iloc[:,i]['std']\ns2 = dt2.iloc[:,i]['std']\nfor t in [4.5,5,5.5,6,6.5,7]:\n    m = (d_bio.iloc[:,i].abs() > t*s1 )\n    d.loc[str(t), 'Count beyond K * sigma-random ' ] = m.sum()\n    print(t, m.sum())\n    m = (d_bio.iloc[:,i].abs() > t*s2 )\n    d.loc[str(t), 'Count beyond K * sigma-bio ' ] = m.sum()\n    print(t, m.sum())\nd","metadata":{"execution":{"iopub.status.busy":"2023-01-02T19:51:59.800684Z","iopub.execute_input":"2023-01-02T19:51:59.801289Z","iopub.status.idle":"2023-01-02T19:51:59.908871Z","shell.execute_reply.started":"2023-01-02T19:51:59.801252Z","shell.execute_reply":"2023-01-02T19:51:59.907714Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dt1 = d_rand_perm.describe()\ndt2 = d_bio.describe()\n\nd = pd.DataFrame()\nd.index.name = 'K * sgima, K:'\nfor i in range(d_bio.shape[1]): # i = 0\n    s1 = dt1.iloc[:,i]['std']\n    s2 = dt2.iloc[:,i]['std']\n    for t in [4.5,5,5.5,6,6.5,7]:\n        m = (d_bio.iloc[:,i].abs() > t*s1 )\n        d.loc[str(t), 'Count beyond K * sigma-random ' + d_bio.columns[i] ] = m.sum()\n        #print(t, m.sum())\n        m = (d_bio.iloc[:,i].abs() > t*s2 )\n        d.loc[str(t), 'Count beyond K * sigma-bio ' + d_bio.columns[i] ] = m.sum()\n        #print(t, m.sum())\nd","metadata":{"execution":{"iopub.status.busy":"2023-01-02T19:51:59.910669Z","iopub.execute_input":"2023-01-02T19:51:59.911021Z","iopub.status.idle":"2023-01-02T19:52:00.117184Z","shell.execute_reply.started":"2023-01-02T19:51:59.910990Z","shell.execute_reply":"2023-01-02T19:52:00.116040Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np.min([d.loc['5'].iat[2*i+1] for i in range(9)  ] ), np.max([d.loc['5'].iat[2*i+1] for i in range(9)  ] )","metadata":{"execution":{"iopub.status.busy":"2023-01-02T19:52:00.118840Z","iopub.execute_input":"2023-01-02T19:52:00.119750Z","iopub.status.idle":"2023-01-02T19:52:00.132385Z","shell.execute_reply.started":"2023-01-02T19:52:00.119698Z","shell.execute_reply":"2023-01-02T19:52:00.130999Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"d.loc['5',::2].min(), d.loc['5',::2].max()\n","metadata":{"execution":{"iopub.status.busy":"2023-01-02T19:52:00.134547Z","iopub.execute_input":"2023-01-02T19:52:00.135260Z","iopub.status.idle":"2023-01-02T19:52:00.147884Z","shell.execute_reply.started":"2023-01-02T19:52:00.135209Z","shell.execute_reply":"2023-01-02T19:52:00.146597Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Analysis of intersections ","metadata":{}},{"cell_type":"code","source":"d = d_bio.std(axis = 1)\nd = d.to_frame()\nd.columns = ['Std']\nd['Corr'] = d_bio.mean(axis = 1)\nd['Abs'] = d_bio.mean(axis = 1).abs()\ndisplay(d)\nplt.title('Std of correlations over day x donor', fontsize = 20 )\nplt.plot( d.sort_values('Abs', ascending = False)['Std'].values[:2000] )\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-01-02T19:52:00.149572Z","iopub.execute_input":"2023-01-02T19:52:00.150477Z","iopub.status.idle":"2023-01-02T19:52:00.435376Z","shell.execute_reply.started":"2023-01-02T19:52:00.150391Z","shell.execute_reply":"2023-01-02T19:52:00.433843Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"d_bio.corr()","metadata":{"execution":{"iopub.status.busy":"2023-01-02T19:52:00.436785Z","iopub.execute_input":"2023-01-02T19:52:00.437152Z","iopub.status.idle":"2023-01-02T19:52:00.465432Z","shell.execute_reply.started":"2023-01-02T19:52:00.437121Z","shell.execute_reply":"2023-01-02T19:52:00.464500Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"d2 = pd.DataFrame(index = d_bio.index)\nd2['Index1'] = range(len(d2))\nd2['Index2'] = range(len(d2))\nd2 = d2.join(d_bio)\nd2 = d2.sort_values(d_bio.columns[0], ascending = False )\nd2['Index1'] = range(len(d2))\n\nd2 = d2.sort_values(d_bio.columns[1], ascending = False )\nd2['Index2'] = range(len(d2))\ndisplay(d2)\nv = (d2['Index1'] - d2['Index2'])\nplt.figure(figsize = (20,5) )\nplt.plot( v.values[:100] )\nplt.grid()\nplt.show()\n\nv.head(30)\n","metadata":{"execution":{"iopub.status.busy":"2023-01-02T19:52:00.466794Z","iopub.execute_input":"2023-01-02T19:52:00.467723Z","iopub.status.idle":"2023-01-02T19:52:00.752511Z","shell.execute_reply.started":"2023-01-02T19:52:00.467682Z","shell.execute_reply":"2023-01-02T19:52:00.751514Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pd.set_option('display.max_rows', 100)","metadata":{"execution":{"iopub.status.busy":"2023-01-02T19:52:00.753959Z","iopub.execute_input":"2023-01-02T19:52:00.754665Z","iopub.status.idle":"2023-01-02T19:52:00.760724Z","shell.execute_reply.started":"2023-01-02T19:52:00.754621Z","shell.execute_reply":"2023-01-02T19:52:00.759217Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"v2 = v.to_frame().join(d2)","metadata":{"execution":{"iopub.status.busy":"2023-01-02T19:52:00.762776Z","iopub.execute_input":"2023-01-02T19:52:00.763331Z","iopub.status.idle":"2023-01-02T19:52:00.782340Z","shell.execute_reply.started":"2023-01-02T19:52:00.763292Z","shell.execute_reply":"2023-01-02T19:52:00.780826Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Show most \"unstable\" genes - those whose order changes from day 2 to day 3 ","metadata":{}},{"cell_type":"code","source":"v2[v2[0]>100].head(100)\n","metadata":{"execution":{"iopub.status.busy":"2023-01-02T19:52:00.784088Z","iopub.execute_input":"2023-01-02T19:52:00.784609Z","iopub.status.idle":"2023-01-02T19:52:00.846022Z","shell.execute_reply.started":"2023-01-02T19:52:00.784559Z","shell.execute_reply":"2023-01-02T19:52:00.844755Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Linear growth for initial period\n\nThat confirms that permutation of genes ordered by correlations varying days and donors is not random.\nBecause for random permutations the growth is x^2/N, see discussions on Mathoverflow:\n\nhttps://mathoverflow.net/q/437569/10446\n\nhttps://mathoverflow.net/q/437545/10446\n","metadata":{}},{"cell_type":"code","source":"%%time \nl1 = []; l2 = []; l3 = []; l4 = []; l5 = []\n\nfor k in range(10, 1000):\n    s0 = set(  d_bio.iloc[:,0].abs().sort_values(ascending = False ).index[:k] )\n    for i in [1,2,3,4,5]:\n        s = s0 & set(  d_bio.iloc[:,i].abs().sort_values(ascending = False ).index[:k] )\n        if i==1:\n            l1.append(len(s))       \n        if i==2:\n            l2.append(len(s))       \n        if i==3:\n            l3.append(len(s))       \n        if i==4:\n            l4.append(len(s))       \n        if i==5:\n            l5.append(len(s))       \n    \nplt.figure(figsize = (20,6))\nplt.plot(l1, label = 'len intersection 1 & 2: day: ' + d_bio.columns[1])\nplt.plot(l2, label = 'len intersection 1 & 3: day: ' + d_bio.columns[2])\nplt.plot(l3, label = 'len intersection 1 & 4: day: ' + d_bio.columns[3])\nplt.plot(l4, label = 'len intersection 1 & 5: day: ' + d_bio.columns[4])\nplt.plot(l5, label = 'len intersection 1 & 6: day: ' + d_bio.columns[5])\nplt.legend(fontsize = 20 )\nplt.grid()\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-01-02T19:52:00.847627Z","iopub.execute_input":"2023-01-02T19:52:00.848039Z","iopub.status.idle":"2023-01-02T19:52:26.683241Z","shell.execute_reply.started":"2023-01-02T19:52:00.848004Z","shell.execute_reply":"2023-01-02T19:52:26.681821Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time \nl2 = []; l3 = []; l6 = []\nfor k in range(10, 1000):\n    for i in [0,1,2,3,4,5]:\n        if i == 0:\n            s = set(  d_bio.iloc[:,i].abs().sort_values(ascending = False ).index[:k] )\n        else:\n            s = s & set(  d_bio.iloc[:,i].abs().sort_values(ascending = False ).index[:k] )\n        if i==1:\n            l2.append(len(s))       \n        if i==2:\n            l3.append(len(s))       \n        if i==5:\n            l6.append(len(s))       \n    \nplt.figure(figsize = (20,6))\nplt.plot(l2, label = 'len 2 lists intersection')\nplt.plot(l3, label = 'len 3 lists intersection')\nplt.plot(l6, label = 'len 6 lists intersection')\nplt.legend(fontsize = 20 )\nplt.grid()\nplt.show()\nprint( l2[-1],  l3[-1], l6[-1] )\n","metadata":{"execution":{"iopub.status.busy":"2023-01-02T19:52:26.685655Z","iopub.execute_input":"2023-01-02T19:52:26.686578Z","iopub.status.idle":"2023-01-02T19:52:52.147281Z","shell.execute_reply.started":"2023-01-02T19:52:26.686511Z","shell.execute_reply":"2023-01-02T19:52:52.145977Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time \nl2 = []; l3 = []; l6 = []\nfor k in range(10, len(d_bio)):\n    for i in [0,1,2,3,4,5]:\n        if i == 0:\n            s = set(  d_bio.iloc[:,i].abs().sort_values(ascending = False ).index[:k] )\n        else:\n            s = s & set(  d_bio.iloc[:,i].abs().sort_values(ascending = False ).index[:k] )\n        if i==1:\n            l2.append(len(s))       \n        if i==2:\n            l3.append(len(s))       \n        if i==5:\n            l6.append(len(s))       \n    \nplt.figure(figsize = (20,6))\nplt.plot(l2, label = 'len 2 lists intersection')\nplt.plot(l3, label = 'len 3 lists intersection')\nplt.plot(l6, label = 'len 6 lists intersection')\nplt.legend(fontsize = 20 )\nplt.grid()\nplt.show()\nprint( l2[-1],  l3[-1], l6[-1] )\n","metadata":{"execution":{"iopub.status.busy":"2023-01-02T19:52:52.149141Z","iopub.execute_input":"2023-01-02T19:52:52.149570Z","iopub.status.idle":"2023-01-02T20:10:30.245691Z","shell.execute_reply.started":"2023-01-02T19:52:52.149533Z","shell.execute_reply":"2023-01-02T20:10:30.243918Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Compare to random shape \n\nFor random permutations the growth is x^2/N for intersection of the two lists, see discussions on Mathoverflow:\n\nhttps://mathoverflow.net/q/437569/10446\n\nhttps://mathoverflow.net/q/437545/10446\n","metadata":{}},{"cell_type":"code","source":"%%time\nN = 1000\nfor i in range(5):\n    v = np.arange(N)\n    p = np.random.permutation(N)\n    l = []\n    for k in range(N):\n        s = set(v[:k]) & set(v[p][:k])\n        l.append(len(s))\n    plt.plot(l)\nplt.show()\nplt.plot(l[:int(N/2)])\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:10:30.247283Z","iopub.execute_input":"2023-01-02T20:10:30.247689Z","iopub.status.idle":"2023-01-02T20:10:31.169256Z","shell.execute_reply.started":"2023-01-02T20:10:30.247656Z","shell.execute_reply":"2023-01-02T20:10:31.167935Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time \n\nd2  = d_bio.mean(axis = 1).to_frame()\nd2.columns = ['Corr']\nd2['Abs'] = d2.iloc[:,0].abs()\nl0 = d2.sort_values('Abs',ascending = False).index\n\nl2 = [];l2c = []; l3 = []; l6 = []; l6b = []; l6c = []\nfor k in range(10, 1000):\n    for i in [0,1,2,3,4,5]:\n        if i == 0:\n            s = set(  d_bio.iloc[:,i].abs().sort_values(ascending = False ).index[:k] )\n        else:\n            s = s & set(  d_bio.iloc[:,i].abs().sort_values(ascending = False ).index[:k] )\n        if i==1:\n            l2.append(len(s & set(l0[:k])         ))       \n            l2c.append(len(s & set(l0[:len(s)])         ))       \n        if i==2:\n            l3.append(len(s & set(l0[:k])         )) \n        if i==5:\n            l6.append(len(s & set(l0[:k])         ))       \n            l6b.append(len(s & set(l0[:100])      ))\n            l6c.append(len(s & set(l0[:len(s)])      ))\n            #print(k, len(s), len(s & set(l0[:len(s)])  ))\n    \nplt.figure(figsize = (20,6))\nplt.plot(l2, label = 'len 2 lists intersection')\nplt.plot(l2c, label = 'len 2 lists intersection with len(s) sorted way 1')\nplt.plot(l3, label = 'len 3 lists intersection')\nplt.plot(l6, label = 'len 6 lists intersection')\nplt.plot(l6b, label = 'len 6 lists intersection with 100 sorted way 1')\nplt.plot(l6b, label = 'len 6 lists intersection with len(s) sorted way 1')\nplt.legend(fontsize = 20 )\nplt.grid()\nplt.show()\nprint( l2[-1],  l3[-1], l6[-1], l6b[-1], l6c[-1] )\nprint()\n","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:10:31.170683Z","iopub.execute_input":"2023-01-02T20:10:31.171199Z","iopub.status.idle":"2023-01-02T20:10:57.607128Z","shell.execute_reply.started":"2023-01-02T20:10:31.171157Z","shell.execute_reply":"2023-01-02T20:10:57.605612Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Amazing stability and coincidence\n\nUp to  50 top K genes sorted by the average correlation (average taken over 9 day x donor segments) \n\nwill be the same as top K from the list obtained by  the intersection of top N lists. \n","metadata":{}},{"cell_type":"code","source":"%%time \n\nd2  = d_bio.mean(axis = 1).to_frame()\nd2.columns = ['Corr']\nd2['Abs'] = d2.iloc[:,0].abs()\nl0 = d2.sort_values('Abs',ascending = False).index\n\nl2 = [];l2c = []; l3 = []; l6 = []; l6b = []; l6c = []\nfor k in range(10, 500):\n    for i in [0,1,2,3,4,5]:\n        if i == 0:\n            s = set(  d_bio.iloc[:,i].abs().sort_values(ascending = False ).index[:k] )\n        else:\n            s = s & set(  d_bio.iloc[:,i].abs().sort_values(ascending = False ).index[:k] )\n        if i==1:\n            l2.append(len(s & set(l0[:k])         ))       \n            l2c.append(len(s & set(l0[:len(s)])         ))       \n        if i==2:\n            l3.append(len(s & set(l0[:k])         )) \n        if i==5:\n            l6.append(len(s & set(l0[:k])         ))       \n            l6b.append(len(s & set(l0[:100])      ))\n            l6c.append(len(s & set(l0[:len(s)])      ))\n            #print(k, len(s), len(s & set(l0[:len(s)])  ))\n    \nplt.figure(figsize = (20,6))\n# plt.plot(l2, label = 'len 2 lists intersection')\n# plt.plot(l2c, label = 'len 2 lists intersection with len(s) sorted way 1')\n# plt.plot(l3, label = 'len 3 lists intersection')\nplt.plot(l6, label = 'len 6 lists intersection')\n#plt.plot(l6b, label = 'len 6 lists intersection with 100 sorted way 1')\nplt.plot(l6b, label = 'len 6 lists intersection with len(s) sorted way 1')\nplt.legend(fontsize = 20 )\nplt.grid()\nplt.show()\n#print( l2[-1],  l3[-1], l6[-1], l6b[-1], l6c[-1] )\nprint()\n","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:10:57.608944Z","iopub.execute_input":"2023-01-02T20:10:57.610388Z","iopub.status.idle":"2023-01-02T20:11:10.807707Z","shell.execute_reply.started":"2023-01-02T20:10:57.610329Z","shell.execute_reply":"2023-01-02T20:11:10.806116Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"display(d_bio.describe())\ndisplay(d_rand_perm.describe())\n\nt = pd.DataFrame()\nt['Bio'] =d_bio.mean(axis = 1)\nt['Rand Perm'] =d_rand_perm.mean(axis = 1)\n\ndisplay( t.describe() )","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:11:10.809819Z","iopub.execute_input":"2023-01-02T20:11:10.810204Z","iopub.status.idle":"2023-01-02T20:11:10.954419Z","shell.execute_reply.started":"2023-01-02T20:11:10.810171Z","shell.execute_reply":"2023-01-02T20:11:10.953001Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"s1 = t['Bio'].std()\ns2 = t['Rand Perm'].std()\nprint(s1,s2)\nm = t['Bio'].abs() > 3.5*s1\nprint(m.sum() )\nm = t['Bio'].abs() > 4*s1\nprint(m.sum() )\nm = t['Bio'].abs() > 3.5*s2\nprint(m.sum() )\nm = t['Bio'].abs() > 4*s2\nprint(m.sum() )\n","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:11:10.958665Z","iopub.execute_input":"2023-01-02T20:11:10.959088Z","iopub.status.idle":"2023-01-02T20:11:10.981229Z","shell.execute_reply.started":"2023-01-02T20:11:10.959053Z","shell.execute_reply":"2023-01-02T20:11:10.980055Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"v = d_bio.mean(axis = 1)\nplt.hist(v,bins = 100)\nplt.hist(d_bio.iloc[:,2],bins = 100)\n\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:11:10.983482Z","iopub.execute_input":"2023-01-02T20:11:10.984590Z","iopub.status.idle":"2023-01-02T20:11:11.993801Z","shell.execute_reply.started":"2023-01-02T20:11:10.984535Z","shell.execute_reply":"2023-01-02T20:11:11.992338Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"d_bio.index = df_rna.columns\nd_bio","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:11:11.995609Z","iopub.execute_input":"2023-01-02T20:11:11.996037Z","iopub.status.idle":"2023-01-02T20:11:12.020874Z","shell.execute_reply.started":"2023-01-02T20:11:11.995984Z","shell.execute_reply":"2023-01-02T20:11:12.019461Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Correlations protein vs RNAs are systematically lower for days 3,4 comparing to day 2 - for all donors the same","metadata":{}},{"cell_type":"code","source":"d2 = d_bio.sort_values(d_bio.columns[0], ascending = False)\nfig = plt.figure(figsize = (20,4))\nplt.plot(d2.iloc[:200,0].values, label = d_bio.columns[0])\nplt.plot(d2.iloc[:200,1].values, label = d_bio.columns[1])\nplt.plot(d2.iloc[:200,3].values,  label = d_bio.columns[2])\nplt.legend(fontsize = 20)\nplt.show()\n\nd2 = d_bio.sort_values(d_bio.columns[3], ascending = False)\nfig = plt.figure(figsize = (20,4))\nfor i in [3,4,5]:\n    plt.plot(d2.iloc[:200,i].values, label = d_bio.columns[i])\nplt.legend(fontsize = 20)\nplt.show()\n\nd2 = d_bio.sort_values(d_bio.columns[6], ascending = False)\nfig = plt.figure(figsize = (20,4))\nfor i in [6,7,8]:\n    plt.plot(d2.iloc[:200,i].values, label = d_bio.columns[i])\nplt.legend(fontsize = 20)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:11:12.022845Z","iopub.execute_input":"2023-01-02T20:11:12.023321Z","iopub.status.idle":"2023-01-02T20:11:12.886604Z","shell.execute_reply.started":"2023-01-02T20:11:12.023268Z","shell.execute_reply":"2023-01-02T20:11:12.885273Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Correlations seems does not change much by donor, (in contrast to day considered above) ","metadata":{}},{"cell_type":"code","source":"for j in [0,1,2]:\n    d2 = d_bio.sort_values(d_bio.columns[j], ascending = False)\n    fig = plt.figure(figsize = (20,5))\n    for i in [0,3,6]:\n        plt.plot(d2.iloc[:200,i+j].values, label = d_bio.columns[i+j])\n    plt.legend(fontsize = 20)\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:11:12.888612Z","iopub.execute_input":"2023-01-02T20:11:12.889488Z","iopub.status.idle":"2023-01-02T20:11:13.774282Z","shell.execute_reply.started":"2023-01-02T20:11:12.889434Z","shell.execute_reply":"2023-01-02T20:11:13.773150Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_corrs = d_bio.mean(axis = 1) \ndf_corrs = df_corrs.to_frame()\ndf_corrs.columns = [ 'Corr']\ndf_corrs['Abs'] = df_corrs.iloc[:,0].abs()\ndf_corrs.sort_values('Abs',ascending = False)","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:11:13.775866Z","iopub.execute_input":"2023-01-02T20:11:13.777071Z","iopub.status.idle":"2023-01-02T20:11:13.806752Z","shell.execute_reply.started":"2023-01-02T20:11:13.777025Z","shell.execute_reply":"2023-01-02T20:11:13.805122Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#pd.set_option('display.max_columns', 500)\npd.set_option('display.max_rows', 100)\n# pd.set_option('display.width', 1000)","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:11:13.808382Z","iopub.execute_input":"2023-01-02T20:11:13.808789Z","iopub.status.idle":"2023-01-02T20:11:13.815064Z","shell.execute_reply.started":"2023-01-02T20:11:13.808754Z","shell.execute_reply":"2023-01-02T20:11:13.813611Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_corrs.sort_values('Abs',ascending = False).head(100)","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:11:13.816641Z","iopub.execute_input":"2023-01-02T20:11:13.817049Z","iopub.status.idle":"2023-01-02T20:11:13.854811Z","shell.execute_reply.started":"2023-01-02T20:11:13.816990Z","shell.execute_reply":"2023-01-02T20:11:13.853210Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"m = df_corrs['Abs'].notnull()\n\ndf_corrs[m].sort_values('Corr',ascending = True).head(100)","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:11:13.856806Z","iopub.execute_input":"2023-01-02T20:11:13.858230Z","iopub.status.idle":"2023-01-02T20:11:13.892994Z","shell.execute_reply.started":"2023-01-02T20:11:13.858167Z","shell.execute_reply":"2023-01-02T20:11:13.891613Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Enrichment","metadata":{}},{"cell_type":"code","source":"!pip install gseapy\nimport gseapy as gp","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:11:13.905571Z","iopub.execute_input":"2023-01-02T20:11:13.906482Z","iopub.status.idle":"2023-01-02T20:11:26.097232Z","shell.execute_reply.started":"2023-01-02T20:11:13.906438Z","shell.execute_reply":"2023-01-02T20:11:26.095869Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Top correlated by absolute value","metadata":{}},{"cell_type":"code","source":"%%time\nN = 100\n#list_genes  = df_results['Corr Pearson NIPS22'].tolist()\nlist_genes  = df_corrs['Abs'].sort_values(ascending = False).index[:N].tolist()\n\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\nN = 200\n#list_genes  = df_results['Corr Pearson NIPS22'].tolist()\nlist_genes  = df_corrs['Abs'].sort_values(ascending = False).index[:N].tolist()\n\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-02T20:11:26.098742Z","iopub.execute_input":"2023-01-02T20:11:26.099175Z","iopub.status.idle":"2023-01-02T20:11:28.699791Z","shell.execute_reply.started":"2023-01-02T20:11:26.099136Z","shell.execute_reply":"2023-01-02T20:11:28.698583Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Top positive correlated ","metadata":{}},{"cell_type":"code","source":"m = df_corrs['Abs'].notnull()\ndf_corrs = df_corrs[m]\n","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:11:28.701629Z","iopub.execute_input":"2023-01-02T20:11:28.702112Z","iopub.status.idle":"2023-01-02T20:11:28.711001Z","shell.execute_reply.started":"2023-01-02T20:11:28.702063Z","shell.execute_reply":"2023-01-02T20:11:28.709705Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nN = 100\n#list_genes  = df_results['Corr Pearson NIPS22'].tolist()\nlist_genes  = df_corrs['Corr'].sort_values(ascending = False).index[:N].tolist()\n\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\nN = 200\n#list_genes  = df_results['Corr Pearson NIPS22'].tolist()\nlist_genes  = df_corrs['Corr'].sort_values(ascending = False).index[:N].tolist()\n\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-02T20:11:28.712803Z","iopub.execute_input":"2023-01-02T20:11:28.713591Z","iopub.status.idle":"2023-01-02T20:11:31.180648Z","shell.execute_reply.started":"2023-01-02T20:11:28.713539Z","shell.execute_reply":"2023-01-02T20:11:31.178990Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Top ANTI correlated","metadata":{}},{"cell_type":"code","source":"df_corrs['Corr'].sort_values(ascending = True).head(10)","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:11:31.183295Z","iopub.execute_input":"2023-01-02T20:11:31.183725Z","iopub.status.idle":"2023-01-02T20:11:31.197176Z","shell.execute_reply.started":"2023-01-02T20:11:31.183686Z","shell.execute_reply":"2023-01-02T20:11:31.195831Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nN = 100\n#list_genes  = df_results['Corr Pearson NIPS22'].tolist()\nlist_genes  = df_corrs['Corr'].sort_values(ascending = True).index[:N].tolist()\n\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\nN = 200\n#list_genes  = df_results['Corr Pearson NIPS22'].tolist()\nlist_genes  = df_corrs['Corr'].sort_values(ascending = True).index[:N].tolist()\n\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-02T20:11:31.199088Z","iopub.execute_input":"2023-01-02T20:11:31.201083Z","iopub.status.idle":"2023-01-02T20:11:33.702316Z","shell.execute_reply.started":"2023-01-02T20:11:31.201035Z","shell.execute_reply":"2023-01-02T20:11:33.700857Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Add analysis of correlation over all cells and over monocytes progenitors ","metadata":{}},{"cell_type":"code","source":"d_bio2 = d_bio.copy()\nd_bio2","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:11:33.704703Z","iopub.execute_input":"2023-01-02T20:11:33.705222Z","iopub.status.idle":"2023-01-02T20:11:33.730899Z","shell.execute_reply.started":"2023-01-02T20:11:33.705174Z","shell.execute_reply":"2023-01-02T20:11:33.729606Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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)\nd_bio2['All cells'] = l\nd_bio2.head(2)    ","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:11:33.733162Z","iopub.execute_input":"2023-01-02T20:11:33.733566Z","iopub.status.idle":"2023-01-02T20:11:55.385952Z","shell.execute_reply.started":"2023-01-02T20:11:33.733514Z","shell.execute_reply":"2023-01-02T20:11:55.384594Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ny = df_y[target_name]\n\nfor cell_type in ['MoP']:#  ['HSC', 'NeuP', 'EryP', 'MasP','MoP']:\n    l = []\n    m = df_meta['cell_type'] == cell_type\n    X = df_rna[m]\n    y2 = y[m]\n    for col in df_rna.columns:\n        c = np.corrcoef(y2,X[col])[0,1]\n        l.append(c)\nd_bio2['MoP'] = l\nd_bio2.sort_values(d_bio2.columns[0], ascending = False).head(30)    ","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:11:55.387918Z","iopub.execute_input":"2023-01-02T20:11:55.388282Z","iopub.status.idle":"2023-01-02T20:12:01.340603Z","shell.execute_reply.started":"2023-01-02T20:11:55.388251Z","shell.execute_reply":"2023-01-02T20:12:01.339277Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"c = d_bio2.columns[0]\nplt.hist(d_bio2[c], bins = 100)\nplt.show()\nplt.hist(d_bio2['All cells'], bins = 100)\nplt.show()\nplt.hist(d_bio2['MoP'], bins = 100)\nplt.show()\n\nd_bio2.describe()","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:12:01.342532Z","iopub.execute_input":"2023-01-02T20:12:01.342942Z","iopub.status.idle":"2023-01-02T20:12:02.560576Z","shell.execute_reply.started":"2023-01-02T20:12:01.342903Z","shell.execute_reply":"2023-01-02T20:12:02.559152Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"d_bio2 = d_bio2.sort_values(d_bio2.columns[0], ascending = False)","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:12:02.562045Z","iopub.execute_input":"2023-01-02T20:12:02.563128Z","iopub.status.idle":"2023-01-02T20:12:02.575296Z","shell.execute_reply.started":"2023-01-02T20:12:02.563084Z","shell.execute_reply":"2023-01-02T20:12:02.573975Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize = (20,7)  )\nfor col in ['2 32606 EryP', 'All cells', 'MoP']:\n    v = d_bio2[col]\n    plt.plot(v.values[:200],'*-', label = col)\nplt.legend(fontsize = 20 )\nplt.grid()    \nplt.title('Correlation with ' + target_name , fontsize = 20)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:12:02.576902Z","iopub.execute_input":"2023-01-02T20:12:02.577918Z","iopub.status.idle":"2023-01-02T20:12:02.896181Z","shell.execute_reply.started":"2023-01-02T20:12:02.577863Z","shell.execute_reply":"2023-01-02T20:12:02.894880Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"r = (-d_bio2).rank().sort_values(d_bio2.columns[0])\nr","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:12:02.898042Z","iopub.execute_input":"2023-01-02T20:12:02.898577Z","iopub.status.idle":"2023-01-02T20:12:02.963108Z","shell.execute_reply.started":"2023-01-02T20:12:02.898499Z","shell.execute_reply":"2023-01-02T20:12:02.961855Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"v = r['All cells']\nfig = plt.figure(figsize = (20,5)  )\nplt.plot(v.values[:200],'*-')\nplt.show()\n\nfig = plt.figure(figsize = (20,5)  )\nplt.plot(v.values[:200],'*-')\nplt.ylim([0,500])\nplt.show()\n\ncol = 'MoP'\nv = r[col]\nfig = plt.figure(figsize = (20,5)  )\nplt.plot(v.values[:200],'*-')\nplt.title(col)\nplt.show()\n\nfig = plt.figure(figsize = (20,5)  )\nplt.plot(v.values[:200],'*-')\nplt.ylim([0,500])\nplt.title(col)\nplt.grid()\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:12:02.964728Z","iopub.execute_input":"2023-01-02T20:12:02.966088Z","iopub.status.idle":"2023-01-02T20:12:03.980942Z","shell.execute_reply.started":"2023-01-02T20:12:02.966037Z","shell.execute_reply":"2023-01-02T20:12:03.979590Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"d_bio2.columns","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:12:03.982339Z","iopub.execute_input":"2023-01-02T20:12:03.982707Z","iopub.status.idle":"2023-01-02T20:12:03.990541Z","shell.execute_reply.started":"2023-01-02T20:12:03.982675Z","shell.execute_reply":"2023-01-02T20:12:03.989485Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time \ndict_counts = {}\nfor k in range(10, 1000):\n    col = '2 32606 EryP'\n    s0 = set(  d_bio[col].abs().sort_values(ascending = False ).index[:k] )\n    for col in ['3 32606 EryP','4 32606 EryP', 'All cells', 'MoP' ]:\n        s = s0 & set(  d_bio2[col].abs().sort_values(ascending = False ).index[:k] )\n        if col not in dict_counts.keys(): \n            dict_counts[col] = []\n        dict_counts[col].append(len(s))\n    \nplt.figure(figsize = (20,6))\nfor col in dict_counts:\n    plt.plot(dict_counts[col], label = 'len intersection 2 32606 EryP and ' + col )\nplt.legend(fontsize = 20 )\nplt.grid()\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:12:03.992034Z","iopub.execute_input":"2023-01-02T20:12:03.992379Z","iopub.status.idle":"2023-01-02T20:12:26.223686Z","shell.execute_reply.started":"2023-01-02T20:12:03.992349Z","shell.execute_reply":"2023-01-02T20:12:26.222172Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for col in dict_counts:\n    plt.figure(figsize = (20,6))\n    y  = np.array(dict_counts[col])\n    plt.plot(y, label = 'len intersection 2 32606 EryP and ' + col )\n    x = np.arange(len(y))\n    p = np.polyfit(x,y,1)\n    print(p)\n    plt.plot(np.polyval(p, x)  )\n    plt.legend(fontsize = 20 )\n    plt.grid()\n    plt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:12:26.225463Z","iopub.execute_input":"2023-01-02T20:12:26.225964Z","iopub.status.idle":"2023-01-02T20:12:27.071153Z","shell.execute_reply.started":"2023-01-02T20:12:26.225918Z","shell.execute_reply":"2023-01-02T20:12:27.067933Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time \nl2 = []; l3 = []; l6 = []\nfor k in range(10, 1000):\n    for i in [0,1,2,3,4,5]:\n        if i == 0:\n            s = set(  d_bio.iloc[:,i].abs().sort_values(ascending = False ).index[:k] )\n        else:\n            s = s & set(  d_bio.iloc[:,i].abs().sort_values(ascending = False ).index[:k] )\n        if i==1:\n            l2.append(len(s))       \n        if i==2:\n            l3.append(len(s))       \n        if i==5:\n            l6.append(len(s))       \n    \nplt.figure(figsize = (20,6))\nplt.plot(l2, label = 'len 2 lists intersection')\nplt.plot(l3, label = 'len 3 lists intersection')\nplt.plot(l6, label = 'len 6 lists intersection')\nplt.legend(fontsize = 20 )\nplt.grid()\nplt.show()\nprint( l2[-1],  l3[-1], l6[-1] )\n","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:12:27.073893Z","iopub.execute_input":"2023-01-02T20:12:27.074583Z","iopub.status.idle":"2023-01-02T20:12:53.238529Z","shell.execute_reply.started":"2023-01-02T20:12:27.074510Z","shell.execute_reply":"2023-01-02T20:12:53.236923Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize = (20,6))\nplt.plot(l6, label = 'len 6 lists intersection')\ny = np.array(l6)\nx = np.arange(len(y))\np = np.polyfit(x,y,1)\nprint(p)\nplt.plot(np.polyval(p,x) )\nplt.legend(fontsize = 20 )\nplt.grid()\nplt.show()\n\nprint('Six list of EryP interessection leaves about 20%  from original genes ')","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:12:53.240448Z","iopub.execute_input":"2023-01-02T20:12:53.240981Z","iopub.status.idle":"2023-01-02T20:12:53.529686Z","shell.execute_reply.started":"2023-01-02T20:12:53.240931Z","shell.execute_reply":"2023-01-02T20:12:53.528068Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time \nl_all = [];\nfor k in range(10, 1000):\n    for i in range(d_bio.shape[1]):\n        if i == 0:\n            s = set(  d_bio.iloc[:,i].abs().sort_values(ascending = False ).index[:k] )\n        else:\n            s = s & set(  d_bio.iloc[:,i].abs().sort_values(ascending = False ).index[:k] )\n    l_all.append(len(s))       \n    \nplt.figure(figsize = (20,6))\nplt.plot(l_all, label = 'all EryP intersection')\nplt.legend(fontsize = 20 )\nplt.grid()\nplt.show()\n\nplt.figure(figsize = (20,6))\nplt.plot(l_all, label = 'all EryP intersection')\ny = np.array(l_all)\nx = np.arange(len(y))\np = np.polyfit(x,y,1)\nprint(p)\nplt.plot(np.polyval(p,x) )\nplt.legend(fontsize = 20 )\nplt.grid()\nplt.show()\n\nprint('All EryP intersection leaves - about 14-20%   ')","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:12:53.531634Z","iopub.execute_input":"2023-01-02T20:12:53.532313Z","iopub.status.idle":"2023-01-02T20:13:33.151441Z","shell.execute_reply.started":"2023-01-02T20:12:53.532270Z","shell.execute_reply":"2023-01-02T20:13:33.149975Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time \ndict_counts = {}\nfor k in range(0, d_bio.shape[0]):\n    col = '2 32606 EryP'\n    s0 = set(  d_bio[col].abs().sort_values(ascending = False ).index[:k] )\n    for col in ['3 32606 EryP']: # ,'4 32606 EryP', 'All cells', 'MoP' ]:\n        s = s0 & set(  d_bio[col].abs().sort_values(ascending = False ).index[:k] )\n        if col not in dict_counts.keys(): \n            dict_counts[col] = []\n        dict_counts[col].append(len(s))\n    \nplt.figure(figsize = (20,6))\nfor col in dict_counts:\n    plt.plot(dict_counts[col], label = 'len intersection 2 32606 EryP and ' + col )\nplt.legend(fontsize = 20 )\nplt.grid()\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:13:33.153179Z","iopub.execute_input":"2023-01-02T20:13:33.153731Z","iopub.status.idle":"2023-01-02T20:19:25.627378Z","shell.execute_reply.started":"2023-01-02T20:13:33.153688Z","shell.execute_reply":"2023-01-02T20:19:25.625894Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Quadratic approximation is quite good on big range - indication of randomness, linear approximation is good at the begining - indication of non-randomness \n","metadata":{}},{"cell_type":"code","source":"col = '3 32606 EryP'\nll = dict_counts[col]\nplt.figure(figsize = (20,6))\nplt.plot(ll, label = 'all EryP intersection')\ny = np.array(ll)\nx = np.arange(len(y))\np = np.polyfit(x,y,2)\nprint(p)\nplt.plot(np.polyval(p,x) )\nplt.legend(fontsize = 20 )\nplt.grid()\nplt.show()\n\n\nN = 2000\nplt.figure(figsize = (20,6))\nplt.plot(ll[:N], label = 'all EryP intersection')\nplt.plot(np.polyval(p,x[:N]) )\n\nM = 100\np2 = np.polyfit(x[:M],y[:M],1)\nprint(p2)\nplt.plot(np.polyval(p2,x[:N]) )\n\nplt.legend(fontsize = 20 )\nplt.grid()\nplt.show()\n\n","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:19:25.628787Z","iopub.execute_input":"2023-01-02T20:19:25.629150Z","iopub.status.idle":"2023-01-02T20:19:26.229537Z","shell.execute_reply.started":"2023-01-02T20:19:25.629118Z","shell.execute_reply":"2023-01-02T20:19:26.228087Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"d_bio.columns","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:19:26.231217Z","iopub.execute_input":"2023-01-02T20:19:26.231636Z","iopub.status.idle":"2023-01-02T20:19:26.240224Z","shell.execute_reply.started":"2023-01-02T20:19:26.231594Z","shell.execute_reply":"2023-01-02T20:19:26.238774Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dict_sorted = {} # pd.DataFrame()\nfor col in ['2 32606 EryP', '3 32606 EryP', '4 32606 EryP', '2 13176 EryP',\n   '3 13176 EryP', '4 13176 EryP', '2 31800 EryP', '3 31800 EryP',\n   '4 31800 EryP']:\n    dict_sorted[col] = d_bio[col].abs().sort_values(ascending = False )\ndict_sorted[col]","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:19:26.242224Z","iopub.execute_input":"2023-01-02T20:19:26.242627Z","iopub.status.idle":"2023-01-02T20:19:26.297939Z","shell.execute_reply.started":"2023-01-02T20:19:26.242591Z","shell.execute_reply":"2023-01-02T20:19:26.296796Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time \ndict_counts = {}\n\ndict_sorted = {} # pd.DataFrame()\nfor col in ['2 32606 EryP', '3 32606 EryP', '4 32606 EryP', '2 13176 EryP',\n   '3 13176 EryP', '4 13176 EryP', '2 31800 EryP', '3 31800 EryP',\n   '4 31800 EryP']:\n    dict_sorted[col] = d_bio[col].abs().sort_values(ascending = False )\n\nfor k in range(0, d_bio.shape[0]):\n    if (k%5000 == 1): print(k)\n    col = '2 32606 EryP'\n    s = set(  dict_sorted[col].index[:k] )\n    for col in ['2 32606 EryP', '3 32606 EryP', '4 32606 EryP', '2 13176 EryP',\n       '3 13176 EryP', '4 13176 EryP', '2 31800 EryP', '3 31800 EryP',\n       '4 31800 EryP']:\n        \n        s = s & set(dict_sorted[col].index[:k] )\n        if col not in dict_counts.keys(): \n            dict_counts[col] = []\n        dict_counts[col].append(len(s))\n\nplt.figure(figsize = (20,6))\nfor col in dict_counts:\n    plt.plot(dict_counts[col], label = 'intersection till: ' + col )\nplt.legend(fontsize = 20 )\nplt.grid()\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:19:26.299573Z","iopub.execute_input":"2023-01-02T20:19:26.299936Z","iopub.status.idle":"2023-01-02T20:31:52.708889Z","shell.execute_reply.started":"2023-01-02T20:19:26.299904Z","shell.execute_reply":"2023-01-02T20:31:52.707329Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Stabilization of intersection - we see it - might be criteria to choose meaningfully correlated genes \n\nThe slope of the line seems to  stabilize with the number of segments growing - so those part which follow linear pattern - might be considered as stable one \n\nWe see that more than 100 genes are in that region even on ","metadata":{}},{"cell_type":"code","source":"for col in  dict_counts: # = '3 32606 EryP'\n    ll = dict_counts[col]\n    plt.figure(figsize = (20,6))\n    plt.plot(ll, label = 'intersection till: ' + col)\n    y = np.array(ll)\n    x = np.arange(len(y))\n    p = np.polyfit(x,y,2)\n    print(p)\n    plt.plot(np.polyval(p,x) )\n    plt.legend(fontsize = 20 )\n    plt.grid()\n    plt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:31:52.710808Z","iopub.execute_input":"2023-01-02T20:31:52.711287Z","iopub.status.idle":"2023-01-02T20:31:55.578542Z","shell.execute_reply.started":"2023-01-02T20:31:52.711250Z","shell.execute_reply":"2023-01-02T20:31:55.577100Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i, col in  enumerate(dict_counts): # = '3 32606 EryP'\n    ll = dict_counts[col]\n    plt.figure(figsize = (20,6))\n    plt.plot(ll, label = 'intersection till: ' + col)\n    y = np.array(ll)\n    x = np.arange(len(y))\n    p = np.polyfit(x,y,2)\n    print(p)\n    plt.plot(np.polyval(p,x) )\n\n    M = 500\n    p = np.polyfit(x[:M],y[:M],1)\n    print(p)\n    plt.plot(np.polyval(p,x) )\n    \n    if i < 3:\n        plt.ylim([0,1000])\n    else:\n        plt.ylim([0,600])\n    plt.xlim([0,2000])\n    plt.legend(fontsize = 20 )\n    plt.grid()\n    plt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:31:55.580401Z","iopub.execute_input":"2023-01-02T20:31:55.580802Z","iopub.status.idle":"2023-01-02T20:31:58.540235Z","shell.execute_reply.started":"2023-01-02T20:31:55.580766Z","shell.execute_reply":"2023-01-02T20:31:58.538537Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndict_sorted = {} # pd.DataFrame()\nfor col in ['2 32606 EryP', '3 32606 EryP', '4 32606 EryP', '2 13176 EryP',\n   '3 13176 EryP', '4 13176 EryP', '2 31800 EryP', '3 31800 EryP',\n   '4 31800 EryP']:\n    dict_sorted[col] = d_bio[col].abs().sort_values(ascending = False )\n\nfor k in range(0, d_bio.shape[0]):\n    if (k%5000 == 1): print(k)\n    col = '2 32606 EryP'\n    s = set(  dict_sorted[col].index[:k] )\n    for col in ['2 32606 EryP', '3 32606 EryP', '4 32606 EryP', '2 13176 EryP',\n       '3 13176 EryP', '4 13176 EryP', '2 31800 EryP', '3 31800 EryP',\n       '4 31800 EryP']:\n        \n        s = s & set(dict_sorted[col].index[:k] )\n    if len(s) >=100: break\nprint(s)\n","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:31:58.542779Z","iopub.execute_input":"2023-01-02T20:31:58.544171Z","iopub.status.idle":"2023-01-02T20:31:58.878749Z","shell.execute_reply.started":"2023-01-02T20:31:58.544109Z","shell.execute_reply":"2023-01-02T20:31:58.877349Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"m = d_bio.index.isin(s)\nd = d_bio[m].sort_values(d_bio.columns[0],ascending = False)\n\nd['Median Corr'] = d.median(axis=1)\nd['Abs Median Corr'] = d['Median Corr'].abs()\n\nd = d.sort_values('Abs Median Corr', ascending = False)\nd.head(50)\n","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:31:58.880694Z","iopub.execute_input":"2023-01-02T20:31:58.881405Z","iopub.status.idle":"2023-01-02T20:31:58.946264Z","shell.execute_reply.started":"2023-01-02T20:31:58.881357Z","shell.execute_reply":"2023-01-02T20:31:58.944859Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Save results ","metadata":{}},{"cell_type":"code","source":"d.to_csv('Top100 Correlated Stable on EryP.csv')","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:31:58.948038Z","iopub.execute_input":"2023-01-02T20:31:58.948551Z","iopub.status.idle":"2023-01-02T20:31:58.962645Z","shell.execute_reply.started":"2023-01-02T20:31:58.948483Z","shell.execute_reply":"2023-01-02T20:31:58.961287Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Correlations on random subsamples ","metadata":{}},{"cell_type":"code","source":"n_trials = 10 # \n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T22:45:00.865932Z","iopub.execute_input":"2023-01-06T22:45:00.866564Z","iopub.status.idle":"2023-01-06T22:45:00.873752Z","shell.execute_reply.started":"2023-01-06T22:45:00.866518Z","shell.execute_reply":"2023-01-06T22:45:00.872007Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nN = int( len(y)*0.1 )\ny = df_y[target_name]\nIX_tmp = np.arange( len(y))\np = np.random.permutation(len(y))\nIX = IX_tmp[p][:N]\nprint( len(IX), len(IX_tmp) )\n\ndf_corr = pd.DataFrame(index = df_rna.columns )\nfor trial in range(n_trials):\n    p = np.random.permutation(len(y))\n    print(p[:10])\n    IX = IX_tmp[p][:N]\n    y2 = y.values[IX]\n    l = []\n    for i,col in enumerate(df_rna.columns):\n        x2 = df_rna[col].values[IX]\n        c = np.corrcoef(y2, x2 )[0,1]\n        l.append(c)\n        if (i%5000==1): print(i, trial);\n    df_corr['Trial '+str(trial)] = l\n    \ndf_corr.sort_values(df_corr.columns[0], ascending = False).head(20)    ","metadata":{"execution":{"iopub.status.busy":"2023-01-06T22:45:04.868825Z","iopub.execute_input":"2023-01-06T22:45:04.869291Z","iopub.status.idle":"2023-01-06T22:45:21.660012Z","shell.execute_reply.started":"2023-01-06T22:45:04.869228Z","shell.execute_reply":"2023-01-06T22:45:21.658879Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_corr.corr()","metadata":{"execution":{"iopub.status.busy":"2023-01-06T22:45:24.631774Z","iopub.execute_input":"2023-01-06T22:45:24.632227Z","iopub.status.idle":"2023-01-06T22:45:24.652239Z","shell.execute_reply.started":"2023-01-06T22:45:24.632193Z","shell.execute_reply":"2023-01-06T22:45:24.650945Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nd2 = df_corr.copy().sort_values(df_corr.columns[0], ascending= False)\nd3 = df_corr.copy().sort_values(df_corr.columns[0], ascending= False)\nfor col in d2.columns:\n    d2[col] = d2[col] - df_corr[df_corr.columns[0]]\n    d3[col] = d2[col].abs() + df_corr[df_corr.columns[0]].abs()\nd2 = d2/d3\nd2\n\nif d2.shape[1]<=100:\n    fig = plt.figure(figsize = (20,5) )\n    for col in d2.columns:\n        plt.plot(d2[col].values, label = col )\n    plt.grid()\n    plt.legend()\n    plt.show()    \n    fig = plt.figure(figsize = (20,5) )\n    for col in d2.columns:\n        plt.plot(d2[col].values[:1000], label = col )\n    plt.grid()\n    plt.legend()\n    plt.show()    \n    fig = plt.figure(figsize = (20,5) )\n    for col in d2.columns:\n        plt.plot(d2[col].values[:100], label = col )\n    plt.grid()\n    plt.legend()\n    plt.show()    \n    fig = plt.figure(figsize = (20,5) )\n    for col in d2.columns:\n        plt.plot(d2[col].values[:500], label = col )\n    plt.grid()\n    plt.legend()\n    plt.show()    ","metadata":{"execution":{"iopub.status.busy":"2023-01-06T22:45:28.718218Z","iopub.execute_input":"2023-01-06T22:45:28.718884Z","iopub.status.idle":"2023-01-06T22:45:31.307287Z","shell.execute_reply.started":"2023-01-06T22:45:28.718849Z","shell.execute_reply":"2023-01-06T22:45:31.305942Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"d = df_corr.sort_values(df_corr.columns[0], ascending = False)\nif d.shape[1]<=100:\n    for col in d.columns:\n        d[col].plot()\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-06T22:45:31.309460Z","iopub.execute_input":"2023-01-06T22:45:31.310374Z","iopub.status.idle":"2023-01-06T22:45:31.651655Z","shell.execute_reply.started":"2023-01-06T22:45:31.310333Z","shell.execute_reply":"2023-01-06T22:45:31.650529Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if d.shape[1]<=100:\n    for col in d.columns:\n        d[col].iloc[:200].plot()\n    plt.show()\n    for col in d.columns:\n        d[col].iloc[5000:5200].plot()\n    plt.show()\n    for col in d.columns:\n        d[col].iloc[1000:1200].plot()\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-06T22:45:33.199095Z","iopub.execute_input":"2023-01-06T22:45:33.199557Z","iopub.status.idle":"2023-01-06T22:45:33.780406Z","shell.execute_reply.started":"2023-01-06T22:45:33.199519Z","shell.execute_reply":"2023-01-06T22:45:33.779252Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_corr.head(1)","metadata":{"execution":{"iopub.status.busy":"2023-01-06T22:45:35.737574Z","iopub.execute_input":"2023-01-06T22:45:35.738381Z","iopub.status.idle":"2023-01-06T22:45:35.753470Z","shell.execute_reply.started":"2023-01-06T22:45:35.738336Z","shell.execute_reply":"2023-01-06T22:45:35.752018Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time \ndict_counts = {}\n\ndict_sorted = {} # pd.DataFrame()\nfor col in df_corr.columns:\n    dict_sorted[col] = df_corr[col].abs().sort_values(ascending = False )\n\nfor k in range(0, d_bio.shape[0]):\n    if (k%5000 == 1): print(k)\n    col = df_corr.columns[0]\n    s = set(  dict_sorted[col].index[:k] )\n    for col in df_corr.columns:\n        \n        s = s & set(dict_sorted[col].index[:k] )\n        if col not in dict_counts.keys(): \n            dict_counts[col] = []\n        dict_counts[col].append(len(s))\n\nplt.figure(figsize = (20,10))\nfor i,col in enumerate(dict_counts):\n    if len(dict_counts ) >= 50:\n        if (i%10) != 0: continue\n    plt.plot(dict_counts[col], label = 'intersection till: ' + col )\nplt.legend(fontsize = 20 )\nplt.grid()\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T22:45:37.248499Z","iopub.execute_input":"2023-01-06T22:45:37.249275Z","iopub.status.idle":"2023-01-06T22:50:57.858725Z","shell.execute_reply.started":"2023-01-06T22:45:37.249204Z","shell.execute_reply":"2023-01-06T22:50:57.857757Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i,col in  enumerate(dict_counts): # = '3 32606 EryP'\n    \n    if len(dict_counts ) >= 50:\n        if (i%10) != 0: continue \n\n    ll = dict_counts[col]\n    plt.figure(figsize = (20,6))\n    plt.plot(ll, label = 'intersection till: ' + col)\n    y = np.array(ll)\n    x = np.arange(len(y))\n    p = np.polyfit(x,y,2)\n    print(p)\n    plt.plot(np.polyval(p,x) )\n    plt.legend(fontsize = 20 )\n    plt.grid()\n    plt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T22:51:14.702233Z","iopub.execute_input":"2023-01-06T22:51:14.702672Z","iopub.status.idle":"2023-01-06T22:51:14.939527Z","shell.execute_reply.started":"2023-01-06T22:51:14.702638Z","shell.execute_reply":"2023-01-06T22:51:14.938354Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"list_slopes = []\nfor i, col in  enumerate(dict_counts): # = '3 32606 EryP'\n    ll = dict_counts[col]\n    \n    y = np.array(ll)\n    x = np.arange(len(y))\n    \n    M = 500\n    p = np.polyfit(x[:M],y[:M],1)\n    print(p)\n    list_slopes.append(p[0])\n    \n    if len(dict_counts ) >= 50:\n        if (i%10) != 0: continue \n        \n    plt.figure(figsize = (20,6))\n    plt.plot(ll, label = 'intersection till: ' + col)\n    y = np.array(ll)\n    x = np.arange(len(y))\n    p = np.polyfit(x,y,2)\n    print(p)\n    plt.plot(np.polyval(p,x) )\n\n    M = 500\n    p = np.polyfit(x[:M],y[:M],1)\n    print(p)\n    plt.plot(np.polyval(p,x) )\n\n    \n    if i < 3:\n        plt.ylim([0,5000])\n    else:\n        plt.ylim([0,600])\n    plt.xlim([0,5000])\n    plt.legend(fontsize = 20 )\n    plt.grid()\n    plt.show()\n    \n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T22:51:15.879834Z","iopub.execute_input":"2023-01-06T22:51:15.880526Z","iopub.status.idle":"2023-01-06T22:51:16.574495Z","shell.execute_reply.started":"2023-01-06T22:51:15.880487Z","shell.execute_reply":"2023-01-06T22:51:16.573308Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Slopes for linear approximation ","metadata":{}},{"cell_type":"code","source":"    \nplt.figure(figsize = (20,10))\nplt.plot( list_slopes     )\nplt.title('Slopes for linear approximation ', fontsize = 20)\nplt.grid()\nplt.show()\n\nplt.figure(figsize = (20,10))\nx = np.arange( 100,len(list_slopes) )\nplt.plot(x,  np.array(list_slopes)[x]     )\nplt.title('Slopes for linear approximation ' , fontsize = 20)\nplt.grid()\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-01-06T22:53:55.623815Z","iopub.execute_input":"2023-01-06T22:53:55.624346Z","iopub.status.idle":"2023-01-06T22:53:56.149914Z","shell.execute_reply.started":"2023-01-06T22:53:55.624303Z","shell.execute_reply":"2023-01-06T22:53:56.148989Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"v = df_corr.sort_values(df_corr.columns[0], ascending = False).std(axis = 1)\nplt.plot(v.values)\nplt.show()\nplt.plot(v.values[:100])\nplt.show()\nplt.plot(v.values[:1000])\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:46:15.972628Z","iopub.execute_input":"2023-01-02T20:46:15.973499Z","iopub.status.idle":"2023-01-02T20:46:16.635638Z","shell.execute_reply.started":"2023-01-02T20:46:15.973445Z","shell.execute_reply":"2023-01-02T20:46:16.633941Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndict_sorted = {} # pd.DataFrame()\nfor col in df_corr.columns:\n    dict_sorted[col] = df_corr[col].abs().sort_values(ascending = False )\n\nfor k in range(0, d_bio.shape[0]):\n    if (k%5000 == 1): print(k)\n    col = df_corr.columns[0]\n    s = set(  dict_sorted[col].index[:k] )\n    for col in df_corr.columns:\n        s = s & set(dict_sorted[col].index[:k] )\n    if len(s) >=100: break\nprint(s)\n\nm = df_corr.index.isin(s)\nd = df_corr[m].sort_values(df_corr.columns[0],ascending = False)\n\nd['Median Corr'] = d.median(axis=1)\nd['Abs Median Corr'] = d['Median Corr'].abs()\n\nd = d.sort_values('Abs Median Corr', ascending = False)\nd.head(50)\n","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:46:16.638300Z","iopub.execute_input":"2023-01-02T20:46:16.639264Z","iopub.status.idle":"2023-01-02T20:46:16.767629Z","shell.execute_reply.started":"2023-01-02T20:46:16.639212Z","shell.execute_reply":"2023-01-02T20:46:16.766350Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Save ordered results for many random subsamples ","metadata":{}},{"cell_type":"code","source":"d.to_csv('Top100_CD36_correlated_allcells_'+str(n_trials) + 'random_subsamples_intersection_ordered.csv')\nd.head(2)","metadata":{"execution":{"iopub.status.busy":"2023-01-02T20:46:16.769541Z","iopub.execute_input":"2023-01-02T20:46:16.770046Z","iopub.status.idle":"2023-01-02T20:46:16.797386Z","shell.execute_reply.started":"2023-01-02T20:46:16.769994Z","shell.execute_reply":"2023-01-02T20:46:16.796216Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}