{"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":"# Introduction of C/C++ calculation and CITEseq visualization","metadata":{}},{"cell_type":"markdown","source":"# Summary\n\nIn my previous [notebook](https://www.kaggle.com/code/konomuabe/eda-multiome-targets), we found that zeros in inputs/targets should be considered missing. While I tried correlation of non-zero pairs, I found calculating this in python is very slow. Then, I tried to calculate non-zero pair correlations in C/C++ language <br>\n\nI also visualized the result of CITEseq correlations\n\n# resutlt\n\nWith C/C++ on Saturn Cloud, I could calculate CITEseq non-zero pair correlation in 30 seconds, and Multiome non-zero pair correlations in 2 to 3 hours, thanks to Satrun Cloud's 32-core CPU. We can also execute the C++ program in this notebook but takes 5 minutes for CITEseq calculation.\n \nI could find some gene expressions which have strong positive correlations with many proteins, but I also found proteins that have weak negative correlations with many gene expressions.\n \n# remaining issue\nVisualizing multiome correlation. It's huge.","metadata":{}},{"cell_type":"code","source":"import os, gc, pickle\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport matplotlib as mpl\nimport numpy as np\nfrom scipy.sparse import *\n\ntry:\n    os.environ['SATURN_IMAGE']\n    kernel = \"saturn_cloud\"\nexcept KeyError:\n    kernel = \"kaggle\"\n\nDATA_DIR = \"/kaggle/input/open-problems-multimodal/\"\nFP_CELL_METADATA = os.path.join(DATA_DIR,\"metadata.csv\")\n\nif kernel==\"saturn_cloud\":\n    SPARSE_DIR = \"./\"\n    META_DIR = \"./\"\nelse: # kaggle\n#    SPARSE_DIR = \"../input/opmsci-sparse-data/\"\n    META_DIR = \"../input/open-problems-multimodal-singlecell-metadata/\"\n    SPARSE_DIR = \"../input/open-problems-multimodal-singlecell-metadata/\"\n\nMETA_CITE_TRAIN = os.path.join(META_DIR,\"cite_train_meta.csv\")\nMETA_CITE_TEST = os.path.join(META_DIR,\"cite_test_meta.csv\")\nMETA_MULTI_TRAIN = os.path.join(META_DIR,\"multiome_train_meta.csv\")\nMETA_MULTI_TEST = os.path.join(META_DIR,\"multiome_test_meta.csv\")\n\nSPS_CITE_TRAIN_INPUTS = os.path.join(SPARSE_DIR,\"train_cite_inputs.npz\")\nSPS_CITE_TRAIN_TARGETS = os.path.join(SPARSE_DIR,\"train_cite_targets.npz\")\nSPS_CITE_TEST_INPUTS = os.path.join(SPARSE_DIR,\"test_cite_inputs.npz\")\n\nSPS_MULTIOME_TRAIN_INPUTS = os.path.join(SPARSE_DIR,\"train_multi_inputs.npz\")\nSPS_MULTIOME_TRAIN_TARGETS = os.path.join(SPARSE_DIR,\"train_multi_targets.npz\")\nSPS_MULTIOME_TEST_INPUTS = os.path.join(SPARSE_DIR,\"test_multi_inputs.npz\")\n\nFP_CITE_TRAIN_INPUTS = os.path.join(DATA_DIR,\"train_cite_inputs.h5\")\nFP_CITE_TRAIN_TARGETS = os.path.join(DATA_DIR,\"train_cite_targets.h5\")\nFP_CITE_TEST_INPUTS = os.path.join(DATA_DIR,\"test_cite_inputs.h5\")\n\nFP_MULTIOME_TRAIN_INPUTS = os.path.join(DATA_DIR,\"train_multi_inputs.h5\")\nFP_MULTIOME_TRAIN_TARGETS = os.path.join(DATA_DIR,\"train_multi_targets.h5\")\nFP_MULTIOME_TEST_INPUTS = os.path.join(DATA_DIR,\"test_multi_inputs.h5\")\n\nFP_SUBMISSION = os.path.join(DATA_DIR,\"sample_submission.csv\")\nFP_EVALUATION_IDS = os.path.join(DATA_DIR,\"evaluation_ids.csv\")\n\nif not os.path.exists('/opt/conda/lib/python3.7/site-packages/tables'):  ## AMBROSM's trick\n    !pip install --quiet tables","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"jupyter":{"source_hidden":true},"execution":{"iopub.status.busy":"2022-09-22T13:48:00.160171Z","iopub.execute_input":"2022-09-22T13:48:00.161438Z","iopub.status.idle":"2022-09-22T13:48:12.914797Z","shell.execute_reply.started":"2022-09-22T13:48:00.161243Z","shell.execute_reply":"2022-09-22T13:48:12.913587Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Export sparse matrix data in raw binary\nWe save the sparse matrix data as raw binary data so that we can read them easily in C/C++.","metadata":{}},{"cell_type":"code","source":"%%time\ncite_y = load_npz(SPS_CITE_TRAIN_TARGETS)\ncite_y = cite_y.tocsc()\ncite_y.indptr.tofile(\"cite_y_indptr.bin\")\ncite_y.indices.tofile(\"cite_y_indices.bin\")    \ncite_y.data.tofile(\"cite_y_data.bin\")\nprint(len(cite_y.indices))\nprint(len(cite_y.data))\n\ncite_x = load_npz(SPS_CITE_TRAIN_INPUTS)\ncite_x = cite_x.tocsc()\ncite_x.indptr.tofile(\"cite_x_indptr.bin\")\ncite_x.indices.tofile(\"cite_x_indices.bin\")\ncite_x.data.tofile(\"cite_x_data.bin\")\nprint(len(cite_x.indices))\nprint(len(cite_x.data))\n","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-09-22T13:48:12.918278Z","iopub.execute_input":"2022-09-22T13:48:12.918591Z","iopub.status.idle":"2022-09-22T13:48:37.485634Z","shell.execute_reply.started":"2022-09-22T13:48:12.918561Z","shell.execute_reply":"2022-09-22T13:48:37.484324Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Equivalent python code\nBelow is the equivalent python code of what I did in C++. Search the non-zero pairs and calculate correlations. This takes 2 or 3 minutes to process 1 target out of 140 targets. It also shows the first 10 correlations as a reference.","metadata":{}},{"cell_type":"code","source":"%%time \ncorr_matrix = np.zeros((cite_x.shape[1],cite_y.shape[1]),dtype = np.float32)\n\nnon_zero_y = np.zeros(cite_y.shape[0])\nnon_zero_x = np.zeros(cite_x.shape[0])\n\n#for y_col_idx in tqdm.tqdm(range(cite_y.shape[1])):\nfor y_col_idx in range(1):\n    for x_col_idx in range(cite_x.shape[1]):\n        nzx = cite_x[:,x_col_idx]!=0\n        nzy = cite_y[:,y_col_idx]!=0 \n        bothnz = (nzx.todense()&nzy.todense())  \n        if bothnz.sum()>=10:  \n            corr_matrix[x_col_idx,y_col_idx] = np.corrcoef(cite_x[:,x_col_idx].toarray()[bothnz],cite_y[:,y_col_idx].toarray()[bothnz])[0,1]\n        else:\n            corr_matrix[x_col_idx,y_col_idx] = -2\nfor x_col_idx in range(10):\n    print(corr_matrix[x_col_idx,0])\ndel cite_y,cite_x\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-09-22T13:48:37.486823Z","iopub.execute_input":"2022-09-22T13:48:37.487160Z","iopub.status.idle":"2022-09-22T13:50:52.223822Z","shell.execute_reply.started":"2022-09-22T13:48:37.487128Z","shell.execute_reply":"2022-09-22T13:50:52.222785Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# C++ program","metadata":{}},{"cell_type":"markdown","source":"I store the C++ program in python raw string. It's lengthy. If you are interested, please expand the below cell.","metadata":{}},{"cell_type":"code","source":"ccode = r'''\n// g++ -O3 full_corr.cc -lpthread\n\n#include <stdio.h>\n#include <sys/stat.h>\n#include <math.h>\n#include <pthread.h>\n#include <time.h>\n//#include <gsl/gsl_statistics.h>\n\n#define CITESEQ \n#define NUM_THREADS (30)\n#define CORR_MIN_DATA   (10)\n\n//#define OUTPUT_INT8\n//#define OUTPUT_TYPE signed char\n//#define OUTPUT_INT16\n//#define OUTPUT_TYPE int\n#define OUTPUT_FLOAT32\n#define OUTPUT_TYPE float\n\n#ifdef CITESEQ\nconst char Y_INDPTR[] = \"cite_y_indptr.bin\";\nconst char Y_INDICES[] = \"cite_y_indices.bin\";\nconst char Y_DATA[] = \"cite_y_data.bin\";\nconst char X_INDPTR[] = \"cite_x_indptr.bin\";\nconst char X_INDICES[] = \"cite_x_indices.bin\";\nconst char X_DATA[] = \"cite_x_data.bin\";\n\nunsigned int x_sz = 22050;\nunsigned int y_sz = 140;\nunsigned int row_sz = 70988;\n#else\nconst char Y_INDPTR[] = \"multi_y_indptr.bin\";\nconst char Y_INDICES[] = \"multi_y_indices.bin\";\nconst char Y_DATA[] = \"multi_y_data.bin\";\nconst char X_INDPTR[] = \"multi_x_indptr.bin\";\nconst char X_INDICES[] = \"multi_x_indices.bin\";\nconst char X_DATA[] = \"multi_x_data.bin\";\n\nunsigned int x_sz = 228942;\nunsigned int y_sz = 23418;\nunsigned int row_sz = 105942;\n#endif\n\nunsigned int y_indptr_len ;\nunsigned int y_indices_len ;\nunsigned int y_data_len ;\n\nunsigned int *y_indptr ;\nunsigned int *y_indices ;\nfloat *y_data ;\n\nunsigned int x_indptr_len ;\nunsigned int x_indices_len ;\nunsigned int x_data_len ;\n\nunsigned int *x_indptr ;\nunsigned int *x_indices ;\nfloat *x_data ;\n\nOUTPUT_TYPE **cor_matrix ;\n\npthread_mutex_t lock;\n\nlong GetFileSize(const char *file)\n{\n    struct stat statBuf;\n\n    if (stat(file, &statBuf) == 0)\n        return statBuf.st_size;\n    return -1L;\n}\n\nfloat correlation(float *x,float *y,int n)\n{\n    double sx = 0, sy=0, sxx = 0, syy = 0 , sxy = 0;\n    for(int i=0;i<n;i++){\n        sx += x[i];\n        sy += y[i];\n        sxx += ((double)x[i])*x[i];\n        syy += ((double)y[i])*y[i];\n        sxy += ((double)x[i])*y[i];\n    }\n    float Sxy = sxy - sx*sy/n;\n    float Sxx = sxx - (sx*sx)/n;\n    float Syy = syy - (sy*sy)/n;\n    return(Sxy/sqrt(Sxx*Syy));\n}\n\n\nint read_files(void)\n{\n//    unsigned int y_indptr_len = GetFileSize(\"y_indptr.bin\");\n//    printf(\"sizeof(unsigned int) %d \\n\",(int)sizeof(unsigned int));\n//    printf(\"%d \\n\",(int)sizeof(float));\n\n    y_indptr_len = GetFileSize(Y_INDPTR);\n    y_indices_len = GetFileSize(Y_INDICES);\n    y_data_len = GetFileSize(Y_DATA);\n    y_indptr_len /= 4;\n    y_indices_len /= 4;\n    y_data_len /= 4;\n\n    printf(\"y_indptr_len %d \\n\",y_indptr_len);\n    printf(\"y_indices_len %d \\n\",y_indices_len);\n    printf(\"y_data_len %d \\n\",y_data_len);\n    y_indptr  = new unsigned int[y_indptr_len];\n    y_indices  = new unsigned int[y_indices_len];\n    y_data  = new float[y_data_len];\n\n    FILE *fp = fopen(Y_INDPTR,\"rb\");\n    unsigned int n_read = fread(y_indptr,4,y_indptr_len,fp);\n    printf(\"%s read %d \\n\",Y_INDPTR,n_read);\n    fclose(fp);\n\n    fp = fopen(Y_INDICES,\"rb\");\n    n_read = fread(y_indices,4,y_indices_len,fp);\n    printf(\"%s read %d \\n\",Y_INDICES,n_read);\n    fclose(fp);\n\n    fp = fopen(Y_DATA,\"rb\");\n    n_read = fread(y_data,4,y_data_len,fp);\n    printf(\"%s read %d \\n\",Y_DATA,n_read);\n    fclose(fp);\n\n    x_indptr_len = GetFileSize(X_INDPTR);\n    x_indices_len = GetFileSize(X_INDICES);\n    x_data_len = GetFileSize(X_DATA);\n\n    x_indptr_len /= 4;\n    x_indices_len /= 4;\n    x_data_len /= 4;\n\n    printf(\"x_indptr_len %d \\n\",x_indptr_len);\n    printf(\"x_indices_len %d \\n\",x_indices_len);\n    printf(\"x_data_len %d \\n\",x_data_len);\n    x_indptr  = new unsigned int[x_indptr_len];\n    x_indices  = new unsigned int[x_indices_len];\n    x_data  = new float[x_data_len];\n\n    fp = fopen(X_INDPTR,\"rb\");\n    n_read = fread(x_indptr,4,x_indptr_len,fp);\n    printf(\"%s read %d \\n\",X_INDPTR,n_read);\n    fclose(fp);\n\n    fp = fopen(X_INDICES,\"rb\");\n    n_read = fread(x_indices,4,x_indices_len,fp);\n    printf(\"%s read %d \\n\",X_INDICES,n_read);\n    fclose(fp);\n\n    fp = fopen(X_DATA,\"rb\");\n    n_read = fread(x_data,4,x_data_len,fp);\n    printf(\"%s read %d \\n\",X_DATA,n_read);\n    fclose(fp);\n\n    return(0);\n}\n\nvoid *calc_y_correlation(void *threadid)\n{   \n\n    long y_col_idx = (long)threadid;\n    unsigned int x_col_idx;\n    float *non_zero_y = new float[row_sz];\n    float *non_zero_x = new float[row_sz];\n    double d;\n\n    for(x_col_idx=0;x_col_idx<x_sz;x_col_idx++){   \n        unsigned int x_idx = x_indptr[x_col_idx];\n        unsigned int x_idx_end = x_indptr[x_col_idx+1];\n        unsigned int y_idx = y_indptr[y_col_idx];\n        unsigned int y_idx_end = y_indptr[y_col_idx+1];\n        unsigned int non_zeros = 0;\n    \n        while(1){\n            if (x_indices[x_idx] == y_indices[y_idx]){\n                non_zero_x[non_zeros] = x_data[x_idx];\n                non_zero_y[non_zeros] = y_data[y_idx];\n                non_zeros += 1;\n                x_idx += 1;\n                y_idx += 1;\n            }\n            while ((x_indices[x_idx] < y_indices[y_idx])\n                    && (x_idx < x_idx_end))\n                x_idx += 1;\n            while ((x_indices[x_idx] > y_indices[y_idx])\n                    && (y_idx < y_idx_end))\n                y_idx += 1;\n            \n\n            if( x_idx >= x_idx_end || y_idx >= y_idx_end)\n                break;\n        }\n        if(non_zeros>=CORR_MIN_DATA){\n            d = correlation(non_zero_x,non_zero_y,non_zeros);\n        }else{\n            d = -2.0;\n        }\n        pthread_mutex_lock(&lock);\n#ifdef OUTPUT_INT8\n        if(d>-2.0)\n            cor_matrix[y_col_idx][x_col_idx] = d * 100;\n        else\n            cor_matrix[y_col_idx][x_col_idx] = 0x80;\n#endif\n#ifdef OUTPUT_INT16\n        if(d>-2.0)\n            cor_matrix[y_col_idx][x_col_idx] = d * 2500;\n        else\n            cor_matrix[y_col_idx][x_col_idx] = 0x8000;\n#endif\n#ifdef OUTPUT_FLOAT32\n            cor_matrix[y_col_idx][x_col_idx] = d;\n#endif\n        pthread_mutex_unlock(&lock);\n    }\n    delete [] non_zero_x;\n    delete [] non_zero_y;\n    pthread_exit(NULL);\n}\n\nint main(void)\n{\n    int rc;\n    time_t start_time, now; \n    read_files();\n    pthread_t threads[NUM_THREADS];\n    for(int i=0;i<NUM_THREADS;i++)threads[i] = 0;\n    if (pthread_mutex_init(&lock, NULL) != 0) {\n        printf(\"\\n mutex init has failed\\n\");\n        return 1;\n    }\n\n    cor_matrix = new OUTPUT_TYPE* [y_sz];\n    for(int i=0;i<y_sz;i++)\n        cor_matrix[i] = new OUTPUT_TYPE[x_sz];\n\n    long y_col_idx = 0;\n    int finished = 0;\n    start_time = time(NULL);\n    while(1){\n        for(int i=0;i<NUM_THREADS;i++){\n            if(threads[i]==0){\n                rc = pthread_create(&threads[i], NULL, calc_y_correlation, (void *)y_col_idx);\n//                rc = pthread_create(&threads[i], NULL, calc_y_correlation, (void *)(y_col_idx/10));\n//                printf(\"thread started %ld at %ld\\n\",y_col_idx,threads[i]);\n                if(rc==0) {\n                    y_col_idx++;\n                }else{\n                    printf(\"thread failed to start. %ld at %ld\\n\",y_col_idx,threads[i]);\n                }\n            }else{\n                rc = pthread_tryjoin_np(threads[i], NULL);\n                if(rc==0){\n//                     printf(\"thread %ld finished\\n\",threads[i]);\n                     threads[i] = 0;\n                     finished++;\n                     if(finished%1000==0){\n                        now = time(NULL);\n                        float sec_per_y = (now-start_time)/(float)finished;\n                        printf(\"%.1f hours %d finished. %.2f sec/y remaining %d %.1f hours\\n\",(now-start_time)/3600.0,finished,sec_per_y,y_sz-finished,(y_sz-finished)*sec_per_y/3600);\n                     }\n                }\n            }\n        }\n        if(y_col_idx>=y_sz)\n            break;\n    }\n\n    now = time(NULL);\n    printf(\"Calculation %d sec.\\n\",(int)(now-start_time));\n    for(int i=0;i<NUM_THREADS;i++){\n        if(threads[i]!=0){\n            rc = pthread_join(threads[i], NULL);\n            if(rc==0) threads[i] = 0;\n        }\n    }\n    pthread_mutex_destroy(&lock);\n\n    FILE *fp = fopen(\"correlation.dat\",\"wb\");\n    for(int y=0;y<y_sz;y++){\n            fwrite(cor_matrix[y],sizeof(OUTPUT_TYPE),x_sz,fp); \n    }\n    fclose(fp);\n\n    for(int x=0;x<10;x++){\n        float d;\n            d = cor_matrix[0][x];\n#ifdef OUTPUT_INT8\n            d /= 100.0;\n#endif\n#ifdef OUTPUT_INT16\n            d /= 2500.0;\n#endif\n            printf(\"%.9f \",d);\n    }\n    printf(\"\\n\");\n\n    for(int i=0;i<y_sz;i++)\n        delete [] cor_matrix[i];\n    delete[] cor_matrix;\n\n    delete [] y_indptr ;\n    delete [] y_indices ;\n    delete [] y_data ;\n    return(1);\n\n}\n\n'''","metadata":{"_kg_hide-input":true,"jupyter":{"source_hidden":true},"execution":{"iopub.status.busy":"2022-09-22T13:50:52.226475Z","iopub.execute_input":"2022-09-22T13:50:52.226745Z","iopub.status.idle":"2022-09-22T13:50:52.237674Z","shell.execute_reply.started":"2022-09-22T13:50:52.226721Z","shell.execute_reply":"2022-09-22T13:50:52.236653Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Make a c++ source file from a python string.","metadata":{}},{"cell_type":"code","source":"f = open('temp.cc',\"wt\")\nf.write(ccode)\nf.close()","metadata":{"execution":{"iopub.status.busy":"2022-09-22T13:50:52.239702Z","iopub.execute_input":"2022-09-22T13:50:52.240006Z","iopub.status.idle":"2022-09-22T13:50:52.255535Z","shell.execute_reply.started":"2022-09-22T13:50:52.239981Z","shell.execute_reply":"2022-09-22T13:50:52.254339Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Then compile it.","metadata":{}},{"cell_type":"code","source":"! g++ temp.cc -lpthread -O3","metadata":{"execution":{"iopub.status.busy":"2022-09-22T13:50:52.257089Z","iopub.execute_input":"2022-09-22T13:50:52.257428Z","iopub.status.idle":"2022-09-22T13:50:53.246524Z","shell.execute_reply.started":"2022-09-22T13:50:52.257400Z","shell.execute_reply":"2022-09-22T13:50:53.245210Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Then run. This takes about 5 to 6 minutes, although it takes just 30 seconds on Saturn Cloud(You need to modify NUM_THREADS to 30).<br>\nThe last 10 numbers are the first 10 correlations for comparison with the reference.","metadata":{}},{"cell_type":"code","source":"%%time\n! ./a.out","metadata":{"execution":{"iopub.status.busy":"2022-09-22T13:50:53.247664Z","iopub.execute_input":"2022-09-22T13:50:53.247997Z","iopub.status.idle":"2022-09-22T13:53:37.986976Z","shell.execute_reply.started":"2022-09-22T13:50:53.247963Z","shell.execute_reply":"2022-09-22T13:53:37.985756Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Reference from slow python code again,<br>\n-0.052215796\n-0.013762924\n0.3773076\n-0.061093338\n-0.043617174\n-0.0055685155\n-0.030543977\n0.005519003\n-0.034018744\n-0.025657881\n<br>\nThe result seems to be OK.","metadata":{}},{"cell_type":"markdown","source":"# Back import\nLet's back import the result. This is also in raw binary data.","metadata":{}},{"cell_type":"code","source":"#corr = np.fromfile(\"correlation.dat\",dtype=np.uint8)\ncorr = np.fromfile(\"correlation.dat\",dtype=np.float32)\ncorr = corr.reshape(140,22050)\ncorr[0,0:10]","metadata":{"execution":{"iopub.status.busy":"2022-09-22T13:53:37.988250Z","iopub.execute_input":"2022-09-22T13:53:37.988587Z","iopub.status.idle":"2022-09-22T13:53:38.001110Z","shell.execute_reply.started":"2022-09-22T13:53:37.988556Z","shell.execute_reply":"2022-09-22T13:53:38.000322Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In the C++ program, if non-zero pairs are less than 10, which are too few samples, I marked the correlation with -2. Replace these data with 0. (10 is just my rough guess.)","metadata":{}},{"cell_type":"code","source":"print((corr<-1.5).sum())\ncorr[corr<-1.5] = 0","metadata":{"execution":{"iopub.status.busy":"2022-09-22T13:53:38.002180Z","iopub.execute_input":"2022-09-22T13:53:38.002492Z","iopub.status.idle":"2022-09-22T13:53:38.015686Z","shell.execute_reply.started":"2022-09-22T13:53:38.002465Z","shell.execute_reply":"2022-09-22T13:53:38.014445Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Visualization\nFirst, I plot the correlations of the first 20 proteins and the first gene expressions.<br>\nWe can find some strong positive correlations(red) and some strong negative correlations(blue) ","metadata":{}},{"cell_type":"code","source":"fignum = 1\nfig,ax = plt.subplots(1,1,figsize=(12,5))\n#plt.subplots_adjust(hspace = 1)\nfig.suptitle(\"fig.%d:Correlation of top left corner\"%fignum)\nfignum+=1\npos = ax.imshow(corr[0:20,0:20],cmap=mpl.colormaps['seismic'],clim=(-1, 1))\nax.set_title('correlation of target[0:20] - gene expression[0:20]')\n_ = fig.colorbar(pos, ax=ax)\n","metadata":{"execution":{"iopub.status.busy":"2022-09-22T13:53:38.019094Z","iopub.execute_input":"2022-09-22T13:53:38.019810Z","iopub.status.idle":"2022-09-22T13:53:38.311954Z","shell.execute_reply.started":"2022-09-22T13:53:38.019775Z","shell.execute_reply":"2022-09-22T13:53:38.311292Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Next, let's plot the whole combination. Because the data is 22050 gene expressions x 140 proteins, the aspect ratio is very bad. So I split the data into a dozen fragments.","metadata":{}},{"cell_type":"code","source":"bksz = 2000\nfig,ax = plt.subplots(22050//bksz + 1,1,figsize=(40,25))\nplt.subplots_adjust(hspace=0.8)\nfig.suptitle(\"fig.%d:Correlation of all protein vs. all gene expressions\"%fignum)\nfignum+=1\n\nfor i in range(22050//bksz + 1):\n    pos = ax[i].imshow(corr[:,bksz*(i):bksz*(i+1)],cmap=mpl.colormaps['seismic'],clim=(-1, 1))\n    ax[i].set_aspect(1.0)\n    ax[i].set_title(\"gene expression %d - %d \"%(bksz*(i),bksz*(i+1) if i<22050//bksz else 22050))\n    ax[i].set_ylabel('protein column')\n","metadata":{"execution":{"iopub.status.busy":"2022-09-22T13:53:38.313193Z","iopub.execute_input":"2022-09-22T13:53:38.314544Z","iopub.status.idle":"2022-09-22T13:53:40.232954Z","shell.execute_reply.started":"2022-09-22T13:53:38.314510Z","shell.execute_reply":"2022-09-22T13:53:40.231801Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can find bluish horizontal stripes, which means some(or many) proteins have negative correlations with many gene expressions. On the other hand, positive correlations(red) tend to show vertical lines, which means some gene expressions have positive correlations with many proteins.<br>\nWe can also find areas where colors are thicker. I magnify a couple of these areas.","metadata":{}},{"cell_type":"code","source":"fig,ax = plt.subplots(1,3,figsize=(20,5))\nplt.subplots_adjust(wspace=1.0)\nfig.suptitle(\"fig.%d:Magnified Area with strong correlation\"%fignum)\nfignum+=1\npos = ax[0].imshow(corr[:,3500:3600],cmap=mpl.colormaps['seismic'],clim=(-1, 1))\nax[0].set_title('correlation of all protein vs.\\n - gene expression[3500:3600]')\nax[0].set_ylabel('protein column')\npos = ax[1].imshow(corr[:,6600:6700],cmap=mpl.colormaps['seismic'],clim=(-1, 1))\nax[1].set_title('correlation of all protein vs.\\n - gene expression[6600:6700]')\nax[1].set_ylabel('protein column')\npos = ax[2].imshow(corr[:,16800:16900],cmap=mpl.colormaps['seismic'],clim=(-1, 1))\nax[2].set_title('correlation of of all protein vs.\\n - gene expression[16800:16900]')\nax[2].set_ylabel('protein column')\n\n_ = fig.colorbar(pos, ax=ax)","metadata":{"execution":{"iopub.status.busy":"2022-09-22T13:53:40.233995Z","iopub.execute_input":"2022-09-22T13:53:40.234884Z","iopub.status.idle":"2022-09-22T13:53:40.774808Z","shell.execute_reply.started":"2022-09-22T13:53:40.234851Z","shell.execute_reply":"2022-09-22T13:53:40.773250Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The center graph in fig.3 has especially strong reddish lines. Let's check the gene expression names.","metadata":{}},{"cell_type":"code","source":"gene_columns = pd.read_hdf(FP_CITE_TRAIN_INPUTS,start=0,stop=1).columns.values\ngene_columns[6608:6618],gene_columns[6655:6665],gene_columns[6683:6693]","metadata":{"execution":{"iopub.status.busy":"2022-09-22T13:53:40.776603Z","iopub.execute_input":"2022-09-22T13:53:40.777103Z","iopub.status.idle":"2022-09-22T13:53:40.913670Z","shell.execute_reply.started":"2022-09-22T13:53:40.777061Z","shell.execute_reply":"2022-09-22T13:53:40.912754Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The right graph in fig.3 has a mixture of strong positive and negative correlations. Let's check the gene expression names, too.","metadata":{}},{"cell_type":"code","source":"gene_columns[16865:16880]","metadata":{"execution":{"iopub.status.busy":"2022-09-22T13:53:40.915490Z","iopub.execute_input":"2022-09-22T13:53:40.915839Z","iopub.status.idle":"2022-09-22T13:53:40.923802Z","shell.execute_reply.started":"2022-09-22T13:53:40.915808Z","shell.execute_reply":"2022-09-22T13:53:40.922485Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Finally, delete the work files.","metadata":{}},{"cell_type":"code","source":"! rm *.bin\n! rm temp.cc\n! rm a.out","metadata":{"execution":{"iopub.status.busy":"2022-09-22T13:53:40.925937Z","iopub.execute_input":"2022-09-22T13:53:40.926310Z","iopub.status.idle":"2022-09-22T13:53:42.027623Z","shell.execute_reply.started":"2022-09-22T13:53:40.926246Z","shell.execute_reply":"2022-09-22T13:53:42.026310Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}