{"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":"# DRAFT under construction \n\n\n\n# What is about ?\n\nSome first look Ribonanza Kaggle challenge. \n\nNotebook partly forked from IAFOSS https://www.kaggle.com/code/iafoss/rna-starter-0-186-lb - please upvote !\n\n\n## Some EDA notes:\n\n- Train: 821840 samples, Test: 1343823 samples. \n\n- Lengths of train rna seq are mostly equal to 177 (95%) with small number of 170 (30000/2?), 115 (27290/2?), 155 (13038/2?) 206 (4998/2?), in test we have 1 million of sequences of length 207 (probably private), 330+K - 177, and few thousands 457, 307 \n\n- begining 26 symbols of ALL train sequnces are the same: \"GGGAACGACUCGAGUAGAGUCGAAAA\", up to position 34 there is dominance of the sequence \"GGGAACGACUCGAGUAGAGUCGAAAAGAUAUGGA\", from the position 35 it becomes more random - see below and discussion:   https://www.kaggle.com/competitions/stanford-ribonanza-rna-folding/discussion/442561 . However for the train sequences of length 206 - they is domination by 3 sequences till 26+90 position. \n\n- Train labels are mostly placed at positions 26-125, other positions are mostly NANs. While we need to predict full length lables in test ! \n\n- tail sequence mostly \"AAAAGAAACAACAACAACAAC\" -  21 symbols \n\n- First 26 target lablels are NAN (it correposnd to exactly the same 26 symbols at start for all sequences). We mostly have labels for the next  100 positions, and small amount of labels later positions. That does NOT correspond to the fact that only last 21 symbols are exactly the same. Also always NAN in train targets for reactivity_IX ,  IX = 166..206 and the  same for reactivity_error\n\n\n    There are some non-nans at exceptional positions: \n    3 positions: reactivity_0127,reactivity_0128,reactivity_0129 -  47791 - exceptional lengths 170,155,206\n    5 positions: reactivity_0131 ...reactivity_0135, count:   17813 - exceptional lengths 155,206\n    30 positions at reactivity_0136 ... reactivity_0165 count: 4998 - exactly those of length 206\n\n\n- For 2A3 median value seems to decrease with position (SN_filter = 1)\n\n- There are 2150418 files in the competition folder see https://www.kaggle.com/code/alexandervc/ribonanza-show-2-millions-files\n\n- Only 29 (2A3), 26 (DMS) percents of targets are within \"good\" range 0,1 \n\n- 45 (2A3), 33 (DMS) percents of targets are EXACTLY zero\n\n- there are some duplicated sequences/sequence_ids - 806573 - unique out of 821840 -  probably experiments were made with different conditions or just source databases are different\n\n\n- Each sequence comes with \"family\" of similar sequences which size varies from 3 thousands to about just 50 - seen by Levenshtein distances - these similar rna should probably grouped togather for groupwise CV - similar to: https://www.kaggle.com/code/alexandervc/cafa5-23-groups-and-folds-diamond-igraph \n\n\n\n- IAFOSS's transformer model take about 3h for one fold on Kaggle GPU P100 - 30 epochs - see below logs on error decrease . \n","metadata":{}},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport time\nt0start = time.time() \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\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\ncc = 0\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        cc += 1\n        if cc < 20:\n            print(os.path.join(dirname, filename))\n        else:\n            break\n    if cc < 20:\n        pass\n    else:\n        break\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":{"execution":{"iopub.status.busy":"2023-09-30T16:05:28.053809Z","iopub.execute_input":"2023-09-30T16:05:28.054151Z","iopub.status.idle":"2023-09-30T16:05:29.565050Z","shell.execute_reply.started":"2023-09-30T16:05:28.054126Z","shell.execute_reply":"2023-09-30T16:05:29.563897Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Train Data Load\n\n\nFrom IAFOSS: \n\nThe primary training data is provided in train_data.csv, which contains 821840 RNA sequences and the corresponding reactivity measurements with 2A3_MaP and DMS_MaP methods. The reactivity is reported in columns reactivity_0001 - reactivity_0206 and is set to NaN for the first 26 and the last 21 nucleotides as well as padding for sequences shorter than 206. For faster loading and effective RAM use, I converted the data into a float32 parquet file. \n\nEvaluation in this competition is performed only on samples with SN_filter = 1 for both measurement methods. In this example, I perform training only on samples wiht SN_filter = 1, which gives a noticeable CV boost but uses only 1/4 of the data (i.e. training on noisy SN_filter = 0 data degrades the performance). A proper consideration of all data as well as reactivity errors may boost the performance.\n\nIn this example, I use a simple CV Kfold split. However, given a mismatch in the RNA length between train/public LB vs. private LB data, **it may be important to verify the effect of the sequence length** to avoid a significant shakeup at the private LB.\n\nOne of the tricks, well known in NLP community, which I use here, is length matching batch sampling: composing batches of samples of approximately the same length to minimize the overhead caused by padding tokens.","metadata":{}},{"cell_type":"code","source":"%%time \nPATH = '/kaggle/input/stanford-ribonanza-rna-folding-converted/'\ndf = pd.read_parquet(os.path.join(PATH,'train_data.parquet'))\nprint( df.shape )\ndf\n","metadata":{"execution":{"iopub.status.busy":"2023-09-30T16:05:29.566856Z","iopub.execute_input":"2023-09-30T16:05:29.567910Z","iopub.status.idle":"2023-09-30T16:05:41.752874Z","shell.execute_reply.started":"2023-09-30T16:05:29.567864Z","shell.execute_reply":"2023-09-30T16:05:41.751780Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df['len'] = [len(t) for t in df['sequence']]\nprint( df['len'].value_counts() )\nprint()\nm = df['SN_filter'] == 1\nprint('SN_filter == 1:', m.sum())\nprint( df['len'][m].value_counts() )\n","metadata":{"execution":{"iopub.status.busy":"2023-09-30T16:06:01.913023Z","iopub.execute_input":"2023-09-30T16:06:01.913419Z","iopub.status.idle":"2023-09-30T16:06:02.567272Z","shell.execute_reply.started":"2023-09-30T16:06:01.913390Z","shell.execute_reply":"2023-09-30T16:06:02.566414Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print( df['experiment_type'].value_counts() )\nprint()\nprint(\"len(set(df['sequence'])):\", len(set(df['sequence'])))\nprint();\n\nm = df['SN_filter'] == 1\nprint('SN_filter == 1:', m.sum())\nprint( df['experiment_type'][m].value_counts() )\n","metadata":{"execution":{"iopub.status.busy":"2023-09-30T16:06:03.687148Z","iopub.execute_input":"2023-09-30T16:06:03.688284Z","iopub.status.idle":"2023-09-30T16:06:04.299077Z","shell.execute_reply.started":"2023-09-30T16:06:03.688244Z","shell.execute_reply":"2023-09-30T16:06:04.297979Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.columns[6], df.columns[7]","metadata":{"execution":{"iopub.status.busy":"2023-09-30T16:06:04.673084Z","iopub.execute_input":"2023-09-30T16:06:04.673460Z","iopub.status.idle":"2023-09-30T16:06:04.680743Z","shell.execute_reply.started":"2023-09-30T16:06:04.673429Z","shell.execute_reply":"2023-09-30T16:06:04.679502Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.head(1)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T14:09:44.120667Z","iopub.execute_input":"2023-09-24T14:09:44.121086Z","iopub.status.idle":"2023-09-24T14:09:44.151666Z","shell.execute_reply.started":"2023-09-24T14:09:44.121055Z","shell.execute_reply":"2023-09-24T14:09:44.150600Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# PCA/UMAP Onehot encoded sequences","metadata":{}},{"cell_type":"code","source":"%%time\nm1 = df['experiment_type'] == 'DMS_MaP'\nprint(df[m1]['len'].value_counts() )\nprint()\n\ndict_len_to_X = {}\nfor L in [155 ,170,177,206]:# 170\n\n    m1 = df['experiment_type'] == 'DMS_MaP'\n    m = m1 & ( df['len'] == L)\n    print(L, m1.sum(), m.sum() )\n#     print(m.sum() )\n\n    import time\n    t0 = time.time()\n    rna_onehot_dict = {'A': [1,0,0,0], 'U': [0,1,0,0], 'G':[0,0,1,0], 'C': [0,0,0,1],}\n\n    v = df['sequence'][m].values\n\n    N = m.sum()\n    X = np.zeros( (N, L*4) , np.int8)\n    for k in range(N):\n        seq = v[k]\n#         if k<=1: print( k, 'len(seq):', len(seq) )\n#         if (k % 400_000 == 0): print(k, '%.1f seconds passed'%(time.time()-t0))\n        l = [rna_onehot_dict[seq[II]] for II in range(L) ]\n        X[k,:] = np.array(l).ravel()\n        \n#     dict_len_to_X[L] = X\n    \n    \n    from sklearn.decomposition import PCA\n    str_inf = 'PCA'\n    reducer = PCA(n_components=3)\n    Xr = reducer.fit_transform(X)\n    for i,j in [[0,1]]: #  ,[0,2],[1,2]]:\n        plt.figure(figsize = (15,6))\n        sns.scatterplot(x= Xr[:,i], y = Xr[:,j] ) # df['reads'])\n        plt.xlabel(str_inf+str(i))\n        plt.ylabel(str_inf+str(j))\n        plt.title('Onehot encoded seqences reduced. Train. Len ' + str(L) + ' Count ' +str(Xr.shape[0] ) , fontsize = 20)\n\n        plt.show()\n\n    import umap\n    str_inf = 'umap' \n    reducer = umap.UMAP(n_components=2)\n    Xr = reducer.fit_transform(X[:150_000,:])\n    for i,j in [[0,1]]: #  ,[0,2],[1,2]]:\n        plt.figure(figsize = (15,6))\n        sns.scatterplot(x= Xr[:,i], y = Xr[:,j] ) # df['reads'])\n        plt.xlabel(str_inf+str(i))\n        plt.ylabel(str_inf+str(j))\n        plt.title('Onehot encoded seqences reduced. Train. Len ' + str(L) + ' Count ' +str(Xr.shape[0]) , fontsize = 20)\n        plt.show()\n    \n    ","metadata":{"execution":{"iopub.status.busy":"2023-09-30T16:44:35.727746Z","iopub.execute_input":"2023-09-30T16:44:35.728040Z","iopub.status.idle":"2023-09-30T16:45:20.189367Z","shell.execute_reply.started":"2023-09-30T16:44:35.728015Z","shell.execute_reply":"2023-09-30T16:45:20.188121Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Statistics for each nucleotide depending on position","metadata":{}},{"cell_type":"code","source":"m1 = df['experiment_type'] == 'DMS_MaP'\ndf[m1]['len'].value_counts()","metadata":{"execution":{"iopub.status.busy":"2023-09-30T16:38:07.595361Z","iopub.execute_input":"2023-09-30T16:38:07.595751Z","iopub.status.idle":"2023-09-30T16:38:08.439884Z","shell.execute_reply.started":"2023-09-30T16:38:07.595721Z","shell.execute_reply":"2023-09-30T16:38:08.438678Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndict_len_to_dfstat = {}\nfor L in [155,170,177,206]:# 170\n\n    m1 = df['experiment_type'] == 'DMS_MaP'\n    m = m1 & ( df['len'] == L)\n    print(L, m1.sum(), m.sum() )\n#     print(m.sum() )\n\n\n    import time\n    t0 = time.time()\n    rsna_dict = {'A': 0, 'U': 1, 'G':2, 'C': 3}\n\n    v = df['sequence'][m].values\n\n    N = m.sum()\n    res = np.zeros( (N, L) , np.int8)\n    for k in range(N):\n        seq = v[k]\n#         if k<=1: print( k, 'len(seq):', len(seq) )\n#         if (k % 400_000 == 0): print(k, '%.1f seconds passed'%(time.time()-t0))\n        l = [rsna_dict[seq[II]] for II in range(L) ]\n        res[k,:] = l\n    res.shape  \n\n    d = pd.DataFrame(index = [0,1,2,3])\n    for k in range(0,res.shape[1]):\n        v = pd.Series(res[:,k]).value_counts()\n        v.name = str(k)\n        d =d.join(v)\n#     print( )\n\n    dict_num2symbol = { rsna_dict[k]:k for k in rsna_dict}\n    d.index = [dict_num2symbol[k] for k in d.index]\n    dict_len_to_dfstat[L] = d\n\n    plt.figure(figsize = (20,6) )\n    for i in range(4):\n        pass\n        v = d.iloc[i,:].values\n        lb =  dict_num2symbol[i]\n        plt.plot(v,'.', label = lb )\n    #     plt.plot(d.iloc[i,:-21].values , label = rsna_dict[i]  )\n    plt.legend()\n    plt.xlabel('position on RNA',fontsize = 20)\n    plt.title('Numbers of nucleotides for different positions. Train Len ' + str(L) + ' Count ' +str(m.sum() ) , fontsize = 20)\n    plt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:29:50.710211Z","iopub.execute_input":"2023-09-30T10:29:50.710786Z","iopub.status.idle":"2023-09-30T10:30:21.700185Z","shell.execute_reply.started":"2023-09-30T10:29:50.710746Z","shell.execute_reply":"2023-09-30T10:30:21.699347Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dict_len_to_dfstat[206].iloc[:,26:46]","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:39:27.065464Z","iopub.execute_input":"2023-09-30T10:39:27.065800Z","iopub.status.idle":"2023-09-30T10:39:27.082617Z","shell.execute_reply.started":"2023-09-30T10:39:27.065775Z","shell.execute_reply":"2023-09-30T10:39:27.081157Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Sequences on length 206 in train are dominated by 3 sequences at begining till 26+90 positions","metadata":{}},{"cell_type":"code","source":"m = df['len'] == 206\nprint(m.sum())\nfor i in range(10):\n    seq = df[m]['sequence'].iat[i]\n    print(seq[26:56],i)","metadata":{"execution":{"iopub.status.busy":"2023-09-30T16:07:40.720118Z","iopub.execute_input":"2023-09-30T16:07:40.720594Z","iopub.status.idle":"2023-09-30T16:07:40.767928Z","shell.execute_reply.started":"2023-09-30T16:07:40.720546Z","shell.execute_reply":"2023-09-30T16:07:40.766788Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"m = df['len'] == 206\nprint(m.sum())\nfor N in [80,90,100,150]:\n    l = [t[26: N] for t in df['sequence'][m] ]\n    print('First 26+N ', N)\n    print( pd.Series(l).value_counts().head(5) )","metadata":{"execution":{"iopub.status.busy":"2023-09-30T16:15:47.112302Z","iopub.execute_input":"2023-09-30T16:15:47.113326Z","iopub.status.idle":"2023-09-30T16:15:47.139382Z","shell.execute_reply.started":"2023-09-30T16:15:47.113290Z","shell.execute_reply":"2023-09-30T16:15:47.138161Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dict_len_to_dfstat[170].iloc[:,126:136]","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:30:21.701703Z","iopub.execute_input":"2023-09-30T10:30:21.702153Z","iopub.status.idle":"2023-09-30T10:30:21.714780Z","shell.execute_reply.started":"2023-09-30T10:30:21.702125Z","shell.execute_reply":"2023-09-30T10:30:21.713855Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dict_len_to_dfstat[177].iloc[:,128:136]","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:31:27.848863Z","iopub.execute_input":"2023-09-30T10:31:27.849283Z","iopub.status.idle":"2023-09-30T10:31:27.862785Z","shell.execute_reply.started":"2023-09-30T10:31:27.849252Z","shell.execute_reply":"2023-09-30T10:31:27.861926Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# sanity check:\ndict_len_to_dfstat[177].iloc[:,126:136].sum()","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:27:23.195283Z","iopub.execute_input":"2023-09-30T10:27:23.196323Z","iopub.status.idle":"2023-09-30T10:27:23.206415Z","shell.execute_reply.started":"2023-09-30T10:27:23.196272Z","shell.execute_reply":"2023-09-30T10:27:23.204817Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# PCA/UMAP for sequences ordinally encoded","metadata":{}},{"cell_type":"code","source":"%%time\ndict_len_to_dfstat = {}\nfor L in [155,170,177,206]:# 170\n\n    m1 = df['experiment_type'] == 'DMS_MaP'\n    m = m1 & ( df['len'] == L)\n    print(L, m1.sum(), m.sum() )\n#     print(m.sum() )\n\n\n    import time\n    t0 = time.time()\n    rsna_dict = {'A': 0, 'U': 1, 'G':2, 'C': 3}\n\n    v = df['sequence'][m].values\n\n    N = m.sum()\n    res = np.zeros( (N, L) , np.int8)\n    for k in range(N):\n        seq = v[k]\n#         if k<=1: print( k, 'len(seq):', len(seq) )\n#         if (k % 400_000 == 0): print(k, '%.1f seconds passed'%(time.time()-t0))\n        l = [rsna_dict[seq[II]] for II in range(L) ]\n        res[k,:] = l\n    res.shape  \n\n    from sklearn.decomposition import PCA\n    X = res\n    str_inf = 'PCA'\n    reducer = PCA(n_components=3)\n    Xr = reducer.fit_transform(X)\n    for i,j in [[0,1]]: #  ,[0,2],[1,2]]:\n        plt.figure(figsize = (15,6))\n        sns.scatterplot(x= Xr[:,i], y = Xr[:,j] ) # df['reads'])\n        plt.xlabel(str_inf+str(i))\n        plt.ylabel(str_inf+str(j))\n        plt.title('Ordinal encoded seqences reduced. Train. Len ' + str(L) + ' Count ' +str(m.sum() ) , fontsize = 20)\n\n        plt.show()\n\n    import umap\n    str_inf = 'umap' \n    reducer = umap.UMAP(n_components=2)\n    X = res\n    Xr = reducer.fit_transform(X)\n    for i,j in [[0,1]]: #  ,[0,2],[1,2]]:\n        plt.figure(figsize = (15,6))\n        sns.scatterplot(x= Xr[:,i], y = Xr[:,j] ) # df['reads'])\n        plt.xlabel(str_inf+str(i))\n        plt.ylabel(str_inf+str(j))\n        plt.title('Ordinal encoded seqences reduced. Train. Len ' + str(L) + ' Count ' +str(m.sum() ) , fontsize = 20)\n        plt.show()\n    ","metadata":{"execution":{"iopub.status.busy":"2023-09-30T11:06:59.604970Z","iopub.execute_input":"2023-09-30T11:06:59.605481Z","iopub.status.idle":"2023-09-30T12:36:51.648831Z","shell.execute_reply.started":"2023-09-30T11:06:59.605445Z","shell.execute_reply":"2023-09-30T12:36:51.647306Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Dependence of medians on position \n\nFor 2A3 median seems to decrease with position ","metadata":{}},{"cell_type":"code","source":"m = df['SN_filter'] == 1\nm1 = df['experiment_type'] == '2A3_MaP'\nprint( m.sum(), (m&m1).sum() )\nm1 = df['experiment_type'] == 'DMS_MaP'\nprint( m.sum(), (m&m1).sum() )","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:52:42.542994Z","iopub.execute_input":"2023-09-30T10:52:42.543488Z","iopub.status.idle":"2023-09-30T10:52:42.832514Z","shell.execute_reply.started":"2023-09-30T10:52:42.543458Z","shell.execute_reply":"2023-09-30T10:52:42.830781Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\nm = df['SN_filter'] == 1\nm1 = df['experiment_type'] == '2A3_MaP'\nv  = df[m&m1].iloc[:,33:137].median(axis = 0)\nv = v.round(3)\nprint(list(np.round(v.values,3)))\nplt.figure(figsize = (20,4))\nplt.plot(v.values[:90],'*-')\nplt.title('2A3',fontsize = 20)\nplt.show()\n\nm1 = df['experiment_type'] == 'DMS_MaP'\nv  = df[m&m1].iloc[:,33:137].median(axis = 0)\nv = v.round(3)\nprint(list(np.round(v.values,3)))\nplt.figure(figsize = (20,4))\nplt.plot(v.values[:90],'*-')\nplt.title('DMS',fontsize = 20)\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:52:42.997718Z","iopub.execute_input":"2023-09-30T10:52:42.998083Z","iopub.status.idle":"2023-09-30T10:52:45.460326Z","shell.execute_reply.started":"2023-09-30T10:52:42.998056Z","shell.execute_reply":"2023-09-30T10:52:45.458847Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Reactivity_IX always NAN for IX = 1..26, IX = 166..206, same for error","metadata":{}},{"cell_type":"code","source":"%%time\nv = (~np.isnan(df.iloc[:,7:].values)).sum(axis = 0 )\ns = pd.Series(index = df.columns[7:], data = v)\nl = [t for t in list(s[s==0].index ) if 'error' not in t ]\nprint( l)","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:52:46.437889Z","iopub.execute_input":"2023-09-30T10:52:46.439059Z","iopub.status.idle":"2023-09-30T10:52:51.059553Z","shell.execute_reply.started":"2023-09-30T10:52:46.439021Z","shell.execute_reply":"2023-09-30T10:52:51.058051Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"l = [t for t in list(s[s==0].index ) if 'error' in t ]\nprint( l)","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:52:51.061914Z","iopub.execute_input":"2023-09-30T10:52:51.062370Z","iopub.status.idle":"2023-09-30T10:52:51.068920Z","shell.execute_reply.started":"2023-09-30T10:52:51.062331Z","shell.execute_reply":"2023-09-30T10:52:51.067511Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Look on each position interval more closely ","metadata":{}},{"cell_type":"code","source":"# 1643680 \nprint(list(v[26:126])) # ","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:52:53.767610Z","iopub.execute_input":"2023-09-30T10:52:53.768006Z","iopub.status.idle":"2023-09-30T10:52:53.773781Z","shell.execute_reply.started":"2023-09-30T10:52:53.767979Z","shell.execute_reply":"2023-09-30T10:52:53.772569Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(list(v[126:166]))","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:52:54.482729Z","iopub.execute_input":"2023-09-30T10:52:54.483293Z","iopub.status.idle":"2023-09-30T10:52:54.490726Z","shell.execute_reply.started":"2023-09-30T10:52:54.483243Z","shell.execute_reply":"2023-09-30T10:52:54.489741Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"s.iloc[126:130] # 3 positions: reactivity_0127,reactivity_0128,reactivity_0129 -  47791","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:52:54.808517Z","iopub.execute_input":"2023-09-30T10:52:54.808966Z","iopub.status.idle":"2023-09-30T10:52:54.817698Z","shell.execute_reply.started":"2023-09-30T10:52:54.808928Z","shell.execute_reply":"2023-09-30T10:52:54.816676Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nm = df['reactivity_0127'].notnull()\nl = [len(t) for t in df[m]['sequence'] ]\npd.Series(l).value_counts()","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:52:55.544373Z","iopub.execute_input":"2023-09-30T10:52:55.544721Z","iopub.status.idle":"2023-09-30T10:52:55.617551Z","shell.execute_reply.started":"2023-09-30T10:52:55.544690Z","shell.execute_reply":"2023-09-30T10:52:55.616330Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# print( s.iloc[130:132] )\nprint( s.iloc[130:136] ) # 5 positions: reactivity_0131 ...reactivity_0135, count:   17813\n","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:52:56.107317Z","iopub.execute_input":"2023-09-30T10:52:56.107683Z","iopub.status.idle":"2023-09-30T10:52:56.115062Z","shell.execute_reply.started":"2023-09-30T10:52:56.107653Z","shell.execute_reply":"2023-09-30T10:52:56.113642Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nm = df['reactivity_0131'].notnull()\nl = [len(t) for t in df[m]['sequence'] ]\npd.Series(l).value_counts()","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:52:56.919123Z","iopub.execute_input":"2023-09-30T10:52:56.919523Z","iopub.status.idle":"2023-09-30T10:52:56.953584Z","shell.execute_reply.started":"2023-09-30T10:52:56.919497Z","shell.execute_reply":"2023-09-30T10:52:56.952221Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print( (s.iloc[130:170] == 4998 ).sum() )# 30 positions at reactivity_0136 ... reactivity_0165 count: 4998\ns.iloc[164:167]","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:52:57.573485Z","iopub.execute_input":"2023-09-30T10:52:57.573836Z","iopub.status.idle":"2023-09-30T10:52:57.583941Z","shell.execute_reply.started":"2023-09-30T10:52:57.573810Z","shell.execute_reply":"2023-09-30T10:52:57.582471Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nl = [len(t)==206 for t in df['sequence']]\nnp.sum(l)","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:52:58.341047Z","iopub.execute_input":"2023-09-30T10:52:58.341549Z","iopub.status.idle":"2023-09-30T10:52:58.820136Z","shell.execute_reply.started":"2023-09-30T10:52:58.341518Z","shell.execute_reply":"2023-09-30T10:52:58.818954Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nl2 = [ ~np.isnan(t) for t in df['reactivity_0136'] ]\nnp.sum(l2)","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:52:58.822159Z","iopub.execute_input":"2023-09-30T10:52:58.822626Z","iopub.status.idle":"2023-09-30T10:53:01.594196Z","shell.execute_reply.started":"2023-09-30T10:52:58.822597Z","shell.execute_reply":"2023-09-30T10:53:01.593110Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"l2 == l","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:53:01.596289Z","iopub.execute_input":"2023-09-30T10:53:01.596919Z","iopub.status.idle":"2023-09-30T10:53:03.275684Z","shell.execute_reply.started":"2023-09-30T10:53:01.596881Z","shell.execute_reply":"2023-09-30T10:53:03.274400Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# First 26 symbols are the same, disbalance in 27-34 , GAUAUGGA\n\n\n    Same 26 start for all \n    GGGAACGACUCGAGUAGAGUCGAAAA\n    Top represented 34 with frequency 232092 (next is 31334)\n    GGGAACGACUCGAGUAGAGUCGAAAAGAUAUGGA    \n    \n    From length 35 disbalance become to disappear\n","metadata":{"execution":{"iopub.status.busy":"2023-09-22T20:03:42.050851Z","iopub.execute_input":"2023-09-22T20:03:42.051378Z","iopub.status.idle":"2023-09-22T20:03:42.060357Z","shell.execute_reply.started":"2023-09-22T20:03:42.051339Z","shell.execute_reply":"2023-09-22T20:03:42.058308Z"}}},{"cell_type":"code","source":"%%time\n# len('GGGAACGACUCGAGUAGAGUCGAAAA') = 26\n\n\nfor N in [26,27, 28,33,34,35,36]:\n    print( N)\n    l = [s[:N] for s in df['sequence'] ]\n    print( pd.Series(l).value_counts().head(10) )\n    print()\n","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:53:03.278156Z","iopub.execute_input":"2023-09-30T10:53:03.279409Z","iopub.status.idle":"2023-09-30T10:53:09.499930Z","shell.execute_reply.started":"2023-09-30T10:53:03.279368Z","shell.execute_reply":"2023-09-30T10:53:09.498319Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Tail AAAAGAAACAACAACAACAAC  - mostly these 21 symbol ","metadata":{}},{"cell_type":"code","source":"l = [len(s) for s in df['sequence'] ]\npd.Series(l).value_counts()","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:53:09.502107Z","iopub.execute_input":"2023-09-30T10:53:09.502882Z","iopub.status.idle":"2023-09-30T10:53:10.281426Z","shell.execute_reply.started":"2023-09-30T10:53:09.502838Z","shell.execute_reply":"2023-09-30T10:53:10.280318Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nprint(len('AAAAGAAACAACAACAACAAC'))\nprint(df.shape)\nfor N in [155,156,160]:\n    print( N)\n    l = [s[N:] for s in df['sequence'] ]\n    print( pd.Series(l).value_counts().head(10) )\n    print()","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:53:10.283523Z","iopub.execute_input":"2023-09-30T10:53:10.283841Z","iopub.status.idle":"2023-09-30T10:53:12.408201Z","shell.execute_reply.started":"2023-09-30T10:53:10.283813Z","shell.execute_reply":"2023-09-30T10:53:12.407096Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.columns[7]","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:53:12.409723Z","iopub.execute_input":"2023-09-30T10:53:12.410772Z","iopub.status.idle":"2023-09-30T10:53:12.418883Z","shell.execute_reply.started":"2023-09-30T10:53:12.410732Z","shell.execute_reply":"2023-09-30T10:53:12.417591Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Positions where number of NANs changes \n\nThe main change is of couse in position at position i = 25,26 - first 26 are totally NAN, from 27 - mostly not NAN ","metadata":{}},{"cell_type":"code","source":"%%time\nv = df.iloc[:,7:].isnull().sum(axis = 0 )\nfor k in range(len(v)-1):\n    if v[k]!=v[k+1]:print(k,v[k],v[k+1], v.index[k], v.index[k+1])","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:53:12.420496Z","iopub.execute_input":"2023-09-30T10:53:12.421381Z","iopub.status.idle":"2023-09-30T10:53:14.381004Z","shell.execute_reply.started":"2023-09-30T10:53:12.421350Z","shell.execute_reply":"2023-09-30T10:53:14.379633Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load test data","metadata":{}},{"cell_type":"code","source":"%%time\nfn = '/kaggle/input/stanford-ribonanza-rna-folding-converted/test_sequences.parquet'\ndft = pd.read_parquet(fn)\nprint(dft.shape)\ndft\n","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:53:14.382522Z","iopub.execute_input":"2023-09-30T10:53:14.382812Z","iopub.status.idle":"2023-09-30T10:53:18.187423Z","shell.execute_reply.started":"2023-09-30T10:53:14.382787Z","shell.execute_reply":"2023-09-30T10:53:18.186319Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndft['len'] = [len(t) for t in dft['sequence']]\ndft['len'].value_counts()","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:53:18.188784Z","iopub.execute_input":"2023-09-30T10:53:18.189728Z","iopub.status.idle":"2023-09-30T10:53:18.770117Z","shell.execute_reply.started":"2023-09-30T10:53:18.189690Z","shell.execute_reply":"2023-09-30T10:53:18.768696Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Test. PCA/UMAP onehot encoded sequences ","metadata":{}},{"cell_type":"code","source":"%%time\n\ndict_len_to_X = {}\nfor L in [177,207,307,457]:# 170\n\n#     m1 = dft['experiment_type'] == 'DMS_MaP'\n    m =  ( dft['len'] == L)\n    print(L, m1.sum(), m.sum() )\n#     print(m.sum() )\n\n    import time\n    t0 = time.time()\n    rna_onehot_dict = {'A': [1,0,0,0], 'U': [0,1,0,0], 'G':[0,0,1,0], 'C': [0,0,0,1],}\n\n    v = dft['sequence'][m].values\n\n    N = m.sum()\n    X = np.zeros( (N, L*4) , np.int8)\n    for k in range(N):\n        seq = v[k]\n#         if k<=1: print( k, 'len(seq):', len(seq) )\n#         if (k % 400_000 == 0): print(k, '%.1f seconds passed'%(time.time()-t0))\n        l = [rna_onehot_dict[seq[II]] for II in range(L) ]\n        X[k,:] = np.array(l).ravel()\n        \n#     dict_len_to_X[L] = X\n    \n    \n    from sklearn.decomposition import PCA\n    str_inf = 'PCA'\n    reducer = PCA(n_components=3)\n    Xr = reducer.fit_transform(X)\n    for i,j in [[0,1]]: #  ,[0,2],[1,2]]:\n        plt.figure(figsize = (15,6))\n        sns.scatterplot(x= Xr[:,i], y = Xr[:,j] ) # df['reads'])\n        plt.xlabel(str_inf+str(i))\n        plt.ylabel(str_inf+str(j))\n        plt.title('Onehot encoded seqences reduced. Test. Len ' + str(L) + ' Count ' +str(Xr.shape[0] ) , fontsize = 20)\n\n        plt.show()\n\n    import umap\n    str_inf = 'umap' \n    reducer = umap.UMAP(n_components=2)\n    Xr = reducer.fit_transform(X[:50_000,:])\n    for i,j in [[0,1]]: #  ,[0,2],[1,2]]:\n        plt.figure(figsize = (15,6))\n        sns.scatterplot(x= Xr[:,i], y = Xr[:,j] ) # df['reads'])\n        plt.xlabel(str_inf+str(i))\n        plt.ylabel(str_inf+str(j))\n        plt.title('Onehot encoded seqences reduced. Test. Len ' + str(L) + ' Count ' +str(Xr.shape[0]) , fontsize = 20)\n        plt.show()\n    \n    import gc\n    gc.collect()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Nucleotide count statistics for each position. Test","metadata":{}},{"cell_type":"code","source":"%%time\ndict_len_to_dfstat = {}\nfor L in [177,207,307,457]:# 170\n\n#     m1 = df['experiment_type'] == 'DMS_MaP'\n    m =( dft['len'] == L)\n    print(L,  m.sum() )\n#     print(m.sum() )\n\n\n    import time\n    t0 = time.time()\n    rsna_dict = {'A': 0, 'U': 1, 'G':2, 'C': 3}\n\n    v = dft['sequence'][m].values\n\n    N = m.sum()\n    res = np.zeros( (N, L) , np.int8)\n    for k in range(N):\n        seq = v[k]\n#         if k<=1: print( k, 'len(seq):', len(seq) )\n#         if (k % 400_000 == 0): print(k, '%.1f seconds passed'%(time.time()-t0))\n        l = [rsna_dict[seq[II]] for II in range(L) ]\n        res[k,:] = l\n    res.shape  \n\n    d = pd.DataFrame(index = [0,1,2,3])\n    for k in range(0,res.shape[1]):\n        v = pd.Series(res[:,k]).value_counts()\n        v.name = str(k)\n        d =d.join(v)\n#     print( )\n\n    dict_num2symbol = { rsna_dict[k]:k for k in rsna_dict}\n    d.index = [dict_num2symbol[k] for k in d.index]\n    dict_len_to_dfstat[L] = d\n\n    plt.figure(figsize = (20,6) )\n    for i in range(4):\n        pass\n        v = d.iloc[i,:].values\n        lb =  dict_num2symbol[i]\n        plt.plot(v,'.', label = lb )\n    #     plt.plot(d.iloc[i,:-21].values , label = rsna_dict[i]  )\n    plt.legend()\n    plt.xlabel('position on RNA',fontsize = 20)\n    plt.title('Numbers of nucleotides for different positions. Test. Len ' + str(L) + ' Count ' +str(m.sum() ) , fontsize = 20)\n    plt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:53:18.773595Z","iopub.execute_input":"2023-09-30T10:53:18.774208Z","iopub.status.idle":"2023-09-30T10:54:13.727647Z","shell.execute_reply.started":"2023-09-30T10:53:18.774154Z","shell.execute_reply":"2023-09-30T10:54:13.726713Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Test. Dim-reductions sequences ordinal encoded","metadata":{}},{"cell_type":"code","source":"%%time\ndict_len_to_dfstat = {}\nfor L in [177,207,307,457]:# 170\n\n#     m1 = df['experiment_type'] == 'DMS_MaP'\n    m =( dft['len'] == L)\n    print(L,  m.sum() )\n#     print(m.sum() )\n\n\n    import time\n    t0 = time.time()\n    rsna_dict = {'A': 0, 'U': 1, 'G':2, 'C': 3}\n\n    v = dft['sequence'][m].values\n\n    N = m.sum()\n    res = np.zeros( (N, L) , np.int8)\n    for k in range(N):\n        seq = v[k]\n#         if k<=1: print( k, 'len(seq):', len(seq) )\n#         if (k % 400_000 == 0): print(k, '%.1f seconds passed'%(time.time()-t0))\n        l = [rsna_dict[seq[II]] for II in range(L) ]\n        res[k,:] = l\n    res.shape  \n\n    from sklearn.decomposition import PCA\n    X = res\n    str_inf = 'PCA'\n    reducer = PCA(n_components=3)\n    Xr = reducer.fit_transform(X)\n    for i,j in [[0,1]]: #  ,[0,2],[1,2]]:\n        plt.figure(figsize = (15,6))\n        sns.scatterplot(x= Xr[:,i], y = Xr[:,j] ) # df['reads'])\n        plt.xlabel(str_inf+str(i))\n        plt.ylabel(str_inf+str(j))\n        plt.title('Ordinal encoded seqences reduced. Train. Len ' + str(L) + ' Count ' +str(m.sum() ) , fontsize = 20)\n\n        plt.show()\n\n    import umap\n    str_inf = 'umap' \n    reducer = umap.UMAP(n_components=2)\n    X = res[:50_000,:]\n    Xr = reducer.fit_transform(X)\n    for i,j in [[0,1]]: #  ,[0,2],[1,2]]:\n        plt.figure(figsize = (15,6))\n        sns.scatterplot(x= Xr[:,i], y = Xr[:,j] ) # df['reads'])\n        plt.xlabel(str_inf+str(i))\n        plt.ylabel(str_inf+str(j))\n        plt.title('Ordinal encoded seqences reduced. Train. Len ' + str(L) + ' Count ' +str(X.shape[0] ) , fontsize = 20)\n        plt.show()\n    ","metadata":{"execution":{"iopub.status.busy":"2023-09-30T12:36:57.519742Z","iopub.execute_input":"2023-09-30T12:36:57.520086Z","iopub.status.idle":"2023-09-30T12:39:55.539399Z","shell.execute_reply.started":"2023-09-30T12:36:57.520058Z","shell.execute_reply":"2023-09-30T12:39:55.537790Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Check we need to predict all positions , i.e. id_max = id_min + len(sequence)-1","metadata":{}},{"cell_type":"code","source":"%%time\nl =  [ (dft['id_max'].iat[k] - dft['id_min'].iat[k]+1) == len(dft['sequence'].iat[k]) for k in range(len(dft)) ]\nnp.sum(l)","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:54:13.729167Z","iopub.execute_input":"2023-09-30T10:54:13.729734Z","iopub.status.idle":"2023-09-30T10:54:34.417194Z","shell.execute_reply.started":"2023-09-30T10:54:13.729702Z","shell.execute_reply":"2023-09-30T10:54:34.416070Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len(dft['sequence'].iat[0]), len(dft['sequence'].iat[0])+len(dft['sequence'].iat[1]), len(dft['sequence'].iat[0])+len(dft['sequence'].iat[1]) + len(dft['sequence'].iat[2]) ","metadata":{"execution":{"iopub.status.busy":"2023-09-24T19:46:39.198900Z","iopub.execute_input":"2023-09-24T19:46:39.199747Z","iopub.status.idle":"2023-09-24T19:46:39.208974Z","shell.execute_reply.started":"2023-09-24T19:46:39.199702Z","shell.execute_reply":"2023-09-24T19:46:39.208067Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# future - (Boolean) sequences whose data will be collected after the start of the competition (but before final scoring) are labeled as 1.\ndft['future'].value_counts()","metadata":{"execution":{"iopub.status.busy":"2023-09-24T19:46:52.335834Z","iopub.execute_input":"2023-09-24T19:46:52.336652Z","iopub.status.idle":"2023-09-24T19:46:52.359528Z","shell.execute_reply.started":"2023-09-24T19:46:52.336605Z","shell.execute_reply":"2023-09-24T19:46:52.358127Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Lengths in test - mostly 207 (1million)  ","metadata":{}},{"cell_type":"code","source":"%%time\nl = [len(t) for t in dft['sequence']] \npd.Series(l).value_counts()","metadata":{"execution":{"iopub.status.busy":"2023-09-24T19:46:56.373770Z","iopub.execute_input":"2023-09-24T19:46:56.374177Z","iopub.status.idle":"2023-09-24T19:46:57.351562Z","shell.execute_reply.started":"2023-09-24T19:46:56.374144Z","shell.execute_reply":"2023-09-24T19:46:57.350482Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Check how sequences are grouped by length - it will be helpful to prepare submit by X.ravel() - four different subgroups","metadata":{}},{"cell_type":"code","source":"%%time\nl = [len(t) for t in dft['sequence']] \n# pd.Series(l).value_counts()\n\nfor k in range(len(l)-1):\n    if l[k]!=l[k+1]:\n        print(k,l[k],l[k+1])","metadata":{"execution":{"iopub.status.busy":"2023-09-24T20:12:53.132675Z","iopub.execute_input":"2023-09-24T20:12:53.133078Z","iopub.status.idle":"2023-09-24T20:12:53.829009Z","shell.execute_reply.started":"2023-09-24T20:12:53.133033Z","shell.execute_reply":"2023-09-24T20:12:53.828011Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Test 26 start is the same, but no dominance till 34","metadata":{}},{"cell_type":"code","source":"%%time\n# len('GGGAACGACUCGAGUAGAGUCGAAAA') = 26\n\nprint( 'GGGAACGACUCGAGUAGAGUCGAAAA' == 'GGGAACGACUCGAGUAGAGUCGAAAA') # First 26 are exactly the same in train and test \nfor N in [26,27, 28,33,34,35,36]:\n    print( N)\n    l = [s[:N] for s in dft['sequence'] ]\n    print( pd.Series(l).value_counts().head(10) )\n    print()\n","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:54:52.395649Z","iopub.execute_input":"2023-09-30T10:54:52.396049Z","iopub.status.idle":"2023-09-30T10:54:58.140419Z","shell.execute_reply.started":"2023-09-30T10:54:52.396019Z","shell.execute_reply":"2023-09-30T10:54:58.139243Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Test Tail 21 position is typically same to train 'AAAAGAAACAACAACAACAAC'\n\nSeveral repeats of AAC at the end , and before AAAAG ","metadata":{}},{"cell_type":"code","source":"%%time\nprint(len('AAAAGAAACAACAACAACAAC'))\nprint( 'AAAAGAAACAACAACAACAAC' == 'AAAAGAAACAACAACAACAAC', ' test and train top tail are the same  ')\nprint(df.shape)\nfor N in [(207-22), (207-21),(207-20)]:\n    print('N = ', N)\n    l = [s[N:] for s in dft['sequence'] ]\n    print( pd.Series(l).value_counts().head(3) )\n    print()","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:54:58.144350Z","iopub.execute_input":"2023-09-30T10:54:58.144657Z","iopub.status.idle":"2023-09-30T10:54:59.866876Z","shell.execute_reply.started":"2023-09-30T10:54:58.144632Z","shell.execute_reply":"2023-09-30T10:54:59.865711Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# SN_filter seems almost randomly distributed by index, but may be not 100%","metadata":{}},{"cell_type":"code","source":"d = df[['SN_filter']].copy()\nd['Index'] = range(len(df))\nd.corr()","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:54:59.868301Z","iopub.execute_input":"2023-09-30T10:54:59.868606Z","iopub.status.idle":"2023-09-30T10:54:59.926470Z","shell.execute_reply.started":"2023-09-30T10:54:59.868581Z","shell.execute_reply":"2023-09-30T10:54:59.925317Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":" df[ 'SN_filter'].value_counts()","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:54:59.928361Z","iopub.execute_input":"2023-09-30T10:54:59.928655Z","iopub.status.idle":"2023-09-30T10:54:59.952901Z","shell.execute_reply.started":"2023-09-30T10:54:59.928631Z","shell.execute_reply":"2023-09-30T10:54:59.951303Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nplt.figure(figsize = (20,5))\nplt.plot( df[ 'SN_filter'].values)\nplt.show()\nplt.figure(figsize = (20,5))\nplt.plot( df[ 'SN_filter'].values[:5000])\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:54:59.954648Z","iopub.execute_input":"2023-09-30T10:54:59.955060Z","iopub.status.idle":"2023-09-30T10:55:03.996437Z","shell.execute_reply.started":"2023-09-30T10:54:59.955020Z","shell.execute_reply":"2023-09-30T10:55:03.995314Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Look on Targets","metadata":{}},{"cell_type":"code","source":"df.iloc[1:10,33:137]","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:55:03.999529Z","iopub.execute_input":"2023-09-30T10:55:04.000651Z","iopub.status.idle":"2023-09-30T10:55:04.033814Z","shell.execute_reply.started":"2023-09-30T10:55:04.000598Z","shell.execute_reply":"2023-09-30T10:55:04.032536Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nplt.figure(figsize = (20,5))\nfor i in range(33,40):\n    plt.plot(df.iloc[:,i].values)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:55:04.035162Z","iopub.execute_input":"2023-09-30T10:55:04.035494Z","iopub.status.idle":"2023-09-30T10:55:07.504558Z","shell.execute_reply.started":"2023-09-30T10:55:04.035459Z","shell.execute_reply":"2023-09-30T10:55:07.503405Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nplt.figure(figsize = (20,5))\nfor i in range(33,40):\n    plt.plot(np.clip( df.iloc[:,i].values,0,1)  )\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:55:07.506072Z","iopub.execute_input":"2023-09-30T10:55:07.507047Z","iopub.status.idle":"2023-09-30T10:55:51.034769Z","shell.execute_reply.started":"2023-09-30T10:55:07.507007Z","shell.execute_reply":"2023-09-30T10:55:51.033407Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfor i in range(33,50):\n    plt.figure(figsize = (20,5))\n    plt.plot(df.iloc[:,i].values, label = str(i))\n    plt.grid()\n    plt.title(df.columns[i],fontsize = 20 )\n    plt.xlabel('sample index',fontsize =20 )\n    plt.show()\n    ","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:55:51.038595Z","iopub.execute_input":"2023-09-30T10:55:51.039087Z","iopub.status.idle":"2023-09-30T10:56:03.899526Z","shell.execute_reply.started":"2023-09-30T10:55:51.039043Z","shell.execute_reply":"2023-09-30T10:56:03.898159Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nsns.pairplot(np.clip(df.iloc[:10_000,33:38],0,1)  )\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:56:03.901025Z","iopub.execute_input":"2023-09-30T10:56:03.901345Z","iopub.status.idle":"2023-09-30T10:56:11.463569Z","shell.execute_reply.started":"2023-09-30T10:56:03.901319Z","shell.execute_reply":"2023-09-30T10:56:11.462322Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# PCA for targets ","metadata":{}},{"cell_type":"code","source":"%%time\nN0,N1 = 0, 800_000\nfrom sklearn.decomposition import PCA\nX = df.iloc[N0:N1,33:137].fillna(0)\nstr_inf = 'PCA'\nreducer = PCA(n_components=3)\nXr = reducer.fit_transform(X)\nfor i,j in [[0,1],[0,2],[1,2]]:\n    sns.scatterplot(x= Xr[:,i], y = Xr[:,j], hue =  df[ 'SN_filter'].iloc[N0:N1] ) # df['reads'])\n    plt.xlabel(str_inf+str(i))\n    plt.ylabel(str_inf+str(j))\n    \n    plt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-09-30T10:56:11.465046Z","iopub.execute_input":"2023-09-30T10:56:11.465382Z","iopub.status.idle":"2023-09-30T10:57:26.447017Z","shell.execute_reply.started":"2023-09-30T10:56:11.465356Z","shell.execute_reply":"2023-09-30T10:57:26.445599Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":" df[ 'SN_filter'].value_counts()","metadata":{"execution":{"iopub.status.busy":"2023-09-23T18:41:48.367830Z","iopub.execute_input":"2023-09-23T18:41:48.368764Z","iopub.status.idle":"2023-09-23T18:41:48.395991Z","shell.execute_reply.started":"2023-09-23T18:41:48.368730Z","shell.execute_reply":"2023-09-23T18:41:48.393545Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfrom sklearn.decomposition import PCA\n\nfor N0,N1,str_inf1 in [ [0, 821840,'2A4'],[821840, 1643680, 'DMS']]:\n    X = np.clip(df.iloc[N0:N1,33:137].fillna(0),0, 1)\n    str_inf = 'PCA' \n    reducer = PCA(n_components=3)\n    Xr = reducer.fit_transform(X)\n    for i,j in [[0,1],[0,2],[1,2]]:\n        plt.figure(figsize = (20,5))\n        plt.subplot(1,2,1)\n        sns.scatterplot(x= Xr[:,i], y = Xr[:,j], hue =  df[ 'SN_filter'].iloc[N0:N1] ) # df['reads'])\n        plt.xlabel(str_inf+str(i))\n        plt.ylabel(str_inf+str(j))\n        plt.title(str_inf1 + ' SN_filter', fontsize = 20 )\n\n        plt.subplot(1,2,2)\n        v = df[  'signal_to_noise' ].iloc[N0:N1].values\n        v = np.clip( v, np.percentile(v,5), np.percentile(v,95)  )\n        v = pd.Series(v); v.name = 'signal_to_noise' \n        sns.scatterplot(x= Xr[:,i], y = Xr[:,j], hue = v  ) # df['reads'])\n        plt.xlabel(str_inf+str(i))\n        plt.ylabel(str_inf+str(j))\n        plt.title(str_inf1+' signal_to_noise', fontsize = 20 )\n        plt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-09-23T18:41:48.398372Z","iopub.execute_input":"2023-09-23T18:41:48.398792Z","iopub.status.idle":"2023-09-23T18:47:24.289760Z","shell.execute_reply.started":"2023-09-23T18:41:48.398762Z","shell.execute_reply":"2023-09-23T18:47:24.287903Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfrom sklearn.decomposition import PCA\n\n# for N0,N1,str_inf1 in [ [0, 821840,'2A4'],[821840, 1643680, 'DMS']]:\nfor N0,N1,str_inf1 in [ [0, 10_000,'2A4'],[821840, 821840 + 10_000, 'DMS']]:\n    \n    X = np.clip(df.iloc[N0:N1,33:137].fillna(0),0, 1)\n    str_inf = 'PCA' \n    reducer = PCA(n_components=3)\n    Xr = reducer.fit_transform(X)\n    for i,j in [[0,1],[0,2],[1,2]]:\n        plt.figure(figsize = (20,5))\n        plt.subplot(1,2,1)\n        sns.scatterplot(x= Xr[:,i], y = Xr[:,j], hue =  df[ 'SN_filter'].iloc[N0:N1] ) # df['reads'])\n        plt.xlabel(str_inf+str(i))\n        plt.ylabel(str_inf+str(j))\n        plt.title(str_inf1 + ' SN_filter', fontsize = 20 )\n\n        plt.subplot(1,2,2)\n        v = df[  'signal_to_noise' ].iloc[N0:N1].values\n        v = np.clip( v, np.percentile(v,5), np.percentile(v,95)  )\n        v = pd.Series(v); v.name = 'signal_to_noise' \n        sns.scatterplot(x= Xr[:,i], y = Xr[:,j], hue = v  ) # df['reads'])\n        plt.xlabel(str_inf+str(i))\n        plt.ylabel(str_inf+str(j))\n        plt.title(str_inf1+' signal_to_noise', fontsize = 20 )\n        plt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-09-23T18:47:24.291553Z","iopub.execute_input":"2023-09-23T18:47:24.291896Z","iopub.status.idle":"2023-09-23T18:47:32.775451Z","shell.execute_reply.started":"2023-09-23T18:47:24.291871Z","shell.execute_reply":"2023-09-23T18:47:32.774131Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfrom sklearn.decomposition import PCA\nimport umap \n\nfor N0,N1,str_inf1 in [ [0, 10_000,'2A4'],[821840, 821840 + 10_000, 'DMS']]:\n    X = np.clip(df.iloc[N0:N1,33:137].fillna(0),0, 1)\n    str_inf = 'umap' \n#     reducer = PCA(n_components=3)\n    reducer = umap.UMAP(n_components=3)\n    Xr = reducer.fit_transform(X)\n    for i,j in [[0,1],[0,2],[1,2]]:\n        plt.figure(figsize = (20,5))\n        plt.subplot(1,2,1)\n        sns.scatterplot(x= Xr[:,i], y = Xr[:,j], hue =  df[ 'SN_filter'].iloc[N0:N1] ) # df['reads'])\n        plt.xlabel(str_inf+str(i))\n        plt.ylabel(str_inf+str(j))\n        plt.title(str_inf1 + ' SN_filter', fontsize = 20 )\n\n        plt.subplot(1,2,2)\n        v = df[  'signal_to_noise' ].iloc[N0:N1].values\n        v = np.clip( v, np.percentile(v,5), np.percentile(v,95)  )\n        v = pd.Series(v); v.name = 'signal_to_noise' \n        sns.scatterplot(x= Xr[:,i], y = Xr[:,j], hue = v  ) # df['reads'])\n        plt.xlabel(str_inf+str(i))\n        plt.ylabel(str_inf+str(j))\n        plt.title(str_inf1+' signal_to_noise', fontsize = 20 )\n        plt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-09-23T18:47:32.777220Z","iopub.execute_input":"2023-09-23T18:47:32.777601Z","iopub.status.idle":"2023-09-23T18:48:56.597734Z","shell.execute_reply.started":"2023-09-23T18:47:32.777573Z","shell.execute_reply":"2023-09-23T18:48:56.596200Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfrom sklearn.decomposition import PCA\nimport umap \n\nfor N0,N1,str_inf1 in [ [0, 50_000,'2A4'],[821840, 821840 + 50_000, 'DMS']]:\n    X = np.clip(df.iloc[N0:N1,33:137].fillna(0),0, 1)\n    str_inf = 'umap' \n#     reducer = PCA(n_components=3)\n    reducer = umap.UMAP(n_components=2)\n    Xr = reducer.fit_transform(X)\n    for i,j in [[0,1] ]: # ,[0,2],[1,2]]:\n        plt.figure(figsize = (20,5))\n        plt.subplot(1,2,1)\n        sns.scatterplot(x= Xr[:,i], y = Xr[:,j], hue =  df[ 'SN_filter'].iloc[N0:N1] ) # df['reads'])\n        plt.xlabel(str_inf+str(i))\n        plt.ylabel(str_inf+str(j))\n        plt.title(str_inf1 + ' SN_filter', fontsize = 20 )\n\n        plt.subplot(1,2,2)\n        v = df[  'signal_to_noise' ].iloc[N0:N1].values\n        v = np.clip( v, np.percentile(v,5), np.percentile(v,95)  )\n        v = pd.Series(v); v.name = 'signal_to_noise' \n        sns.scatterplot(x= Xr[:,i], y = Xr[:,j], hue = v  ) # df['reads'])\n        plt.xlabel(str_inf+str(i))\n        plt.ylabel(str_inf+str(j))\n        plt.title(str_inf1+' signal_to_noise', fontsize = 20 )\n        plt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-09-23T18:48:56.599465Z","iopub.execute_input":"2023-09-23T18:48:56.600245Z","iopub.status.idle":"2023-09-23T18:50:12.653476Z","shell.execute_reply.started":"2023-09-23T18:48:56.600210Z","shell.execute_reply":"2023-09-23T18:50:12.651900Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfrom sklearn.decomposition import PCA\nimport umap \n\nfor N0,N1,str_inf1 in [ [0, 500_000,'2A4'],[821840, 821840 + 500_000, 'DMS']]:\n    X = np.clip(df.iloc[N0:N1,33:137].fillna(0),0, 1)\n    str_inf = 'umap' \n#     reducer = PCA(n_components=3)\n    reducer = umap.UMAP(n_components=2)\n    Xr = reducer.fit_transform(X)\n    for i,j in [[0,1] ]: # ,[0,2],[1,2]]:\n        plt.figure(figsize = (20,5))\n        plt.suptitle(str(N0)+' - '  +str(N1), fontsize = 20 )\n        plt.subplot(1,2,1)\n        sns.scatterplot(x= Xr[:,i], y = Xr[:,j], hue =  df[ 'SN_filter'].iloc[N0:N1] ) # df['reads'])\n        plt.xlabel(str_inf+str(i))\n        plt.ylabel(str_inf+str(j))\n        plt.title(str_inf1 + ' SN_filter', fontsize = 20 )\n\n        plt.subplot(1,2,2)\n        v = df[  'signal_to_noise' ].iloc[N0:N1].values\n        v = np.clip( v, np.percentile(v,5), np.percentile(v,95)  )\n        v = pd.Series(v); v.name = 'signal_to_noise' \n        sns.scatterplot(x= Xr[:,i], y = Xr[:,j], hue = v  ) # df['reads'])\n        plt.xlabel(str_inf+str(i))\n        plt.ylabel(str_inf+str(j))\n        plt.title(str_inf1+' signal_to_noise '  , fontsize = 20 )\n        plt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-09-23T18:50:12.655348Z","iopub.execute_input":"2023-09-23T18:50:12.655712Z","iopub.status.idle":"2023-09-23T18:50:28.695547Z","shell.execute_reply.started":"2023-09-23T18:50:12.655683Z","shell.execute_reply":"2023-09-23T18:50:28.693886Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Sequence, experiment_type, dataset_name , etc","metadata":{}},{"cell_type":"code","source":"len( set(df['sequence']) ), len( set(df['sequence_id']) ), ","metadata":{"execution":{"iopub.status.busy":"2023-09-23T18:50:28.702934Z","iopub.execute_input":"2023-09-23T18:50:28.703235Z","iopub.status.idle":"2023-09-23T18:50:29.533021Z","shell.execute_reply.started":"2023-09-23T18:50:28.703211Z","shell.execute_reply":"2023-09-23T18:50:29.531937Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df['experiment_type'].value_counts()","metadata":{"execution":{"iopub.status.busy":"2023-09-23T18:50:29.534226Z","iopub.execute_input":"2023-09-23T18:50:29.534534Z","iopub.status.idle":"2023-09-23T18:50:29.607635Z","shell.execute_reply.started":"2023-09-23T18:50:29.534509Z","shell.execute_reply":"2023-09-23T18:50:29.606428Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df['dataset_name'].value_counts().head(5)","metadata":{"execution":{"iopub.status.busy":"2023-09-23T18:50:29.609153Z","iopub.execute_input":"2023-09-23T18:50:29.609526Z","iopub.status.idle":"2023-09-23T18:50:29.679440Z","shell.execute_reply.started":"2023-09-23T18:50:29.609495Z","shell.execute_reply":"2023-09-23T18:50:29.678291Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"m = df['sequence_id'] == '1728e0d67d7c'\nprint(m.sum() )\nprint( df[m]['sequence'].value_counts() )\ndf[m]\n","metadata":{"execution":{"iopub.status.busy":"2023-09-23T18:50:29.681341Z","iopub.execute_input":"2023-09-23T18:50:29.681747Z","iopub.status.idle":"2023-09-23T18:50:29.817614Z","shell.execute_reply.started":"2023-09-23T18:50:29.681713Z","shell.execute_reply":"2023-09-23T18:50:29.816186Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"d = df[['reads','signal_to_noise','SN_filter']].copy()\nd['mean targets'] = df.iloc[:,33:137].mean(axis = 1 )\nd['median targets'] = df.iloc[:,33:137].median(axis = 1 )\nd['std targets'] = df.iloc[:,33:137].std(axis = 1 )\nd['var targets'] = df.iloc[:,33:137].var(axis = 1 )\n\n\n\nd.corr().round(2)","metadata":{"execution":{"iopub.status.busy":"2023-09-23T18:50:29.819125Z","iopub.execute_input":"2023-09-23T18:50:29.819491Z","iopub.status.idle":"2023-09-23T18:50:43.091617Z","shell.execute_reply.started":"2023-09-23T18:50:29.819461Z","shell.execute_reply":"2023-09-23T18:50:43.089289Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df['sequence_id'].value_counts()","metadata":{"execution":{"iopub.status.busy":"2023-09-23T18:50:43.093626Z","iopub.execute_input":"2023-09-23T18:50:43.095626Z","iopub.status.idle":"2023-09-23T18:50:44.572568Z","shell.execute_reply.started":"2023-09-23T18:50:43.095572Z","shell.execute_reply":"2023-09-23T18:50:44.571597Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df['sequence_id'].value_counts().value_counts()","metadata":{"execution":{"iopub.status.busy":"2023-09-23T18:50:44.574640Z","iopub.execute_input":"2023-09-23T18:50:44.575055Z","iopub.status.idle":"2023-09-23T18:50:45.862662Z","shell.execute_reply.started":"2023-09-23T18:50:44.575025Z","shell.execute_reply":"2023-09-23T18:50:45.861181Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len( set( df['sequence_id'] ) ), df.shape, len(df) - len( set( df['sequence_id'] ) )","metadata":{"execution":{"iopub.status.busy":"2023-09-23T18:50:45.864292Z","iopub.execute_input":"2023-09-23T18:50:45.865085Z","iopub.status.idle":"2023-09-23T18:50:46.576095Z","shell.execute_reply.started":"2023-09-23T18:50:45.865049Z","shell.execute_reply":"2023-09-23T18:50:46.574583Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Targets EDA","metadata":{}},{"cell_type":"code","source":"%%time\nfor N0,N1, str_info in [ [0,821840,'2A3'], [821840,  1643680,'DMS']] :\n    #N1 = 821840\n    \n    print( str_info )\n    v = df.iloc[N0:N1,33:136].values.ravel()\n    print('Percent Nan:', np.isnan(v).sum() / len(v)*100 )\n    print('Percent == 0 :', (v==0).sum() / len(v)*100 )\n\n    print('Mean:', np.mean(v[~np.isnan(v)]))\n    print('Percent < 0 (bad)', (v <0).sum()/len(v)*100 )\n    print('Percent < 1 (very bad)', (v <-1).sum()/len(v)*100 )\n    print('Percent - \"good\" - within 0,1:',  ((v >0)&(v<1)).sum() /len(v) ) \n    print( 'Percent >1 (bad)', (v >1).sum()/len(v) *100 )\n    print( 'Percent >2 (very bad)', (v >2).sum()/len(v) *100 )\n    plt.figure(figsize = (20,5))\n    plt.subplot(1,3,1)\n    plt.hist(v.ravel(),bins = 100)\n    plt.title(str_info, fontsize = 20 )\n    plt.subplot(1,3,2)\n    plt.hist(v[ (v>0)&(v<3)] ,bins = 100)\n    plt.title(str_info + ' cut not >0 , <3', fontsize = 20 )\n    plt.subplot(1,3,3)\n    plt.hist(v[ (v>-1)&(v<2)] ,bins = 100)\n    plt.title(str_info+ ' cut including 0', fontsize = 20 )\n    \n    plt.show()\n    print()","metadata":{"execution":{"iopub.status.busy":"2023-09-23T18:50:46.577782Z","iopub.execute_input":"2023-09-23T18:50:46.578139Z","iopub.status.idle":"2023-09-23T18:50:54.586678Z","shell.execute_reply.started":"2023-09-23T18:50:46.578111Z","shell.execute_reply":"2023-09-23T18:50:54.584793Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nplt.figure(figsize = (20,4) )\nfor N0,N1, str_info in [ [0,821840,'2A3'], [821840,  1643680,'DMS']] :\n    #N1 = 821840\n    \n    print( str_info )\n\n    d = df.iloc[N0:N1,33:137].mean()#.values\n    d.plot(label = str_info)\n    plt.legend()\n    plt.grid()\n    plt.title(str_info + ' targets mean for positions', fontsize = 20  )\nplt.show()\n\nplt.figure(figsize = (20,4) )\nfor N0,N1, str_info in [ [0,821840,'2A3'], [821840,  1643680,'DMS']] :\n    #N1 = 821840\n    \n    print( str_info )\n\n    d = df.iloc[N0:N1,33:137].std()#.values\n    d.plot(label = str_info)\n    plt.legend()\n    plt.grid()\n    plt.title(str_info + ' targets std for positions', fontsize = 20  )\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-09-23T18:50:54.588164Z","iopub.execute_input":"2023-09-23T18:50:54.588554Z","iopub.status.idle":"2023-09-23T18:50:56.743830Z","shell.execute_reply.started":"2023-09-23T18:50:54.588501Z","shell.execute_reply":"2023-09-23T18:50:56.742857Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Correlations for Targets\n\nSeems the first (by index) RNAs in the list - are NOT random, but quite similar\nand we see certain correlations for neigbour-placed targets.\n\nIf we take similar for say indexed 40_000 - 50_000 - that disappears \n","metadata":{}},{"cell_type":"code","source":"%%time\ncm = df.iloc[:1_000,33:137].corr()\nsns.heatmap(cm.abs()[cm.abs()>0.2])\nplt.title('First 1000 - 2A3', fontsize = 20)\nplt.show()\ndisplay( cm.iloc[:10,:10] )\n\ncm = df.iloc[:10_000,33:137].corr()\nsns.heatmap(cm.abs()[cm.abs()>0.2])\nplt.title('First 10_000 - 2A3', fontsize = 20)\nplt.show()\ndisplay( cm.iloc[:10,:10] )\n\ncm = df.iloc[40_000:50_000,33:137].corr()\nsns.heatmap(cm.abs()[cm.abs()>0.2])\nplt.title('40_000:50_000 - 2A3', fontsize = 20)\nplt.show()\ndisplay( cm.iloc[:10,:10] )\n\nN0 = 821840\ncm = df.iloc[N0:N0+10_000,33:137].corr()\nsns.heatmap(cm.abs()[cm.abs()>0.2])\nplt.title('Fist 10_000  - DMS', fontsize = 20)\n\nplt.show()\ndisplay( cm.iloc[:10,:10] )","metadata":{"execution":{"iopub.status.busy":"2023-09-23T18:50:56.745204Z","iopub.execute_input":"2023-09-23T18:50:56.745692Z","iopub.status.idle":"2023-09-23T18:50:59.460979Z","shell.execute_reply.started":"2023-09-23T18:50:56.745662Z","shell.execute_reply":"2023-09-23T18:50:59.459887Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Lengths","metadata":{"execution":{"iopub.status.busy":"2023-09-12T09:58:21.281807Z","iopub.execute_input":"2023-09-12T09:58:21.282197Z","iopub.status.idle":"2023-09-12T09:58:22.007053Z","shell.execute_reply.started":"2023-09-12T09:58:21.282166Z","shell.execute_reply":"2023-09-12T09:58:22.005942Z"}}},{"cell_type":"code","source":"%%time\nl = [len(t) for t in df['sequence']]\nprint( pd.Series(l).describe() )\n\nplt.hist(l, bins = 100)\nplt.show()\npd.Series(l).value_counts()","metadata":{"execution":{"iopub.status.busy":"2023-09-23T18:50:59.462423Z","iopub.execute_input":"2023-09-23T18:50:59.462824Z","iopub.status.idle":"2023-09-23T18:51:04.993346Z","shell.execute_reply.started":"2023-09-23T18:50:59.462790Z","shell.execute_reply":"2023-09-23T18:51:04.992090Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"1568354/len(df)*100","metadata":{"execution":{"iopub.status.busy":"2023-09-23T18:51:04.994853Z","iopub.execute_input":"2023-09-23T18:51:04.995181Z","iopub.status.idle":"2023-09-23T18:51:05.003322Z","shell.execute_reply.started":"2023-09-23T18:51:04.995152Z","shell.execute_reply":"2023-09-23T18:51:05.001743Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Levenshtein distances ","metadata":{}},{"cell_type":"code","source":"from Levenshtein import distance\nedit_dist = distance(\"ah\", \"aho\")\nedit_dist","metadata":{"execution":{"iopub.status.busy":"2023-09-23T18:51:05.005411Z","iopub.execute_input":"2023-09-23T18:51:05.005863Z","iopub.status.idle":"2023-09-23T18:51:05.114163Z","shell.execute_reply.started":"2023-09-23T18:51:05.005824Z","shell.execute_reply":"2023-09-23T18:51:05.112536Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ni = 0\ns = df['sequence'].iat[i]\nl = [distance(s,t) for t in df['sequence']]\nprint( pd.Series(l).describe() )\n\nplt.hist(l, bins = 100)\nplt.show()\npd.Series(l).value_counts()","metadata":{"execution":{"iopub.status.busy":"2023-09-23T18:51:05.115787Z","iopub.execute_input":"2023-09-23T18:51:05.116156Z","iopub.status.idle":"2023-09-23T18:51:15.370025Z","shell.execute_reply.started":"2023-09-23T18:51:05.116127Z","shell.execute_reply":"2023-09-23T18:51:15.368305Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"m = pd.Series(l) < 40\nprint(m.sum() )\nplt.plot(df[m].index)\ndf[m]","metadata":{"execution":{"iopub.status.busy":"2023-09-23T18:51:15.372631Z","iopub.execute_input":"2023-09-23T18:51:15.373092Z","iopub.status.idle":"2023-09-23T18:51:15.867475Z","shell.execute_reply.started":"2023-09-23T18:51:15.373060Z","shell.execute_reply":"2023-09-23T18:51:15.864972Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfor j in range(10):\n    i = 10_000*j\n    s = df['sequence'].iat[i]\n    l = [distance(s,t) for t in df['sequence']]\n    sr = pd.Series(l)\n    print()\n    print(i, 'Count, percent dist < 40 : ' , (sr<40).sum() , (sr<40).sum() / len(sr)*100 ) \n    print()\n    print( sr.describe() )\n    print(sr.value_counts().head(3) )\n    \n    plt.hist(l, bins = 100)\n    plt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-09-23T18:51:15.869619Z","iopub.execute_input":"2023-09-23T18:51:15.869994Z","iopub.status.idle":"2023-09-23T18:52:58.345104Z","shell.execute_reply.started":"2023-09-23T18:51:15.869963Z","shell.execute_reply":"2023-09-23T18:52:58.343515Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n#for j in range(10):\ni = 10_000*1\ns = df['sequence'].iat[i]\nl = [distance(s,t) for t in df['sequence']]\nsr = pd.Series(l)\nprint()\nprint(i, 'Count, percent dist < 40 : ' , (sr<40).sum() , (sr<40).sum() / len(sr)*100 ) \nprint()\nprint( sr.describe() )\nprint(sr.value_counts().head(3) )\n\nplt.hist(l, bins = 100)\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-09-23T18:52:58.346891Z","iopub.execute_input":"2023-09-23T18:52:58.347893Z","iopub.status.idle":"2023-09-23T18:53:08.557040Z","shell.execute_reply.started":"2023-09-23T18:52:58.347860Z","shell.execute_reply":"2023-09-23T18:53:08.555448Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nm = sr<40\nprint( m.sum() )\nv = df[m].iloc[:,33:137].std()\n# v = np.clip(df[m].iloc[:,33:137],0,1).std()\nplt.plot(v)\n\n# v = np.clip(df.iloc[:,33:137],0,1).std()\nv = df.iloc[:,33:137].std()\nplt.plot(v)\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-09-23T18:53:08.559510Z","iopub.execute_input":"2023-09-23T18:53:08.559969Z","iopub.status.idle":"2023-09-23T18:53:10.498710Z","shell.execute_reply.started":"2023-09-23T18:53:08.559936Z","shell.execute_reply":"2023-09-23T18:53:10.497144Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nm = sr<40\nprint( m.sum() )\n# v = df[m].iloc[:,33:137].std()\nv = np.clip(df[m].iloc[:,33:137],0,1).std()\nplt.plot(v)\n\nv = np.clip(df.iloc[:,33:137],0,1).std()\nplt.plot(v)\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-09-23T18:53:10.500062Z","iopub.execute_input":"2023-09-23T18:53:10.500413Z","iopub.status.idle":"2023-09-23T18:53:15.063296Z","shell.execute_reply.started":"2023-09-23T18:53:10.500368Z","shell.execute_reply":"2023-09-23T18:53:15.062249Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np.sort(l)[2:60]","metadata":{"execution":{"iopub.status.busy":"2023-09-23T18:53:15.064619Z","iopub.execute_input":"2023-09-23T18:53:15.064915Z","iopub.status.idle":"2023-09-23T18:53:15.209941Z","shell.execute_reply.started":"2023-09-23T18:53:15.064890Z","shell.execute_reply":"2023-09-23T18:53:15.208776Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Description\n\nWelcome to Stanford Ribonanza RNA Folding challenge. The task of this competition is predicting the chemical reactivity at each position of an RNA molecule. These data are extremely sensitive to the structure that each RNA forms, and an algorithm that could perfectly predict these chemical reactivities would need to have an implicit ‘understanding’ of RNA structure. Such an oracle could be then utilized to predictively model structures of novel RNA molecules. A better understanding of how to manipulate RNA could help usher in an age of programmable medicine, including first cures for pancreatic cancer and Alzheimer’s disease as well as much-needed antibiotics and new biotechnology approaches for climate change. \n\nThis notebook provides a simple baseline that may be used as a starting point for further experiments. Improvment of the baseline may include: \n1. Use of proper loss function to incorporate SN_filter = 0 samples as well as reactivity errors into training\n2. Model improvement and use of additional data, e.g. Ribonanza_bpp_files\n\nFinally, working on this competition keep in mind that train/public LB have different sequence length distribution from private LB, i.e. 115-206 vs. 207-457. Therefore, to avoid a strong shakeup one may need to look into performance vs. sequence length end ensure generalizability.","metadata":{}},{"cell_type":"markdown","source":"# Data\n\nThe primary training data is provided in train_data.csv, which contains 821840 RNA sequences and the corresponding reactivity measurements with 2A3_MaP and DMS_MaP methods. The reactivity is reported in columns reactivity_0001 - reactivity_0206 and is set to NaN for the first 26 and the last 21 nucleotides as well as padding for sequences shorter than 206. For faster loading and effective RAM use, I converted the data into a float32 parquet file. \n\nEvaluation in this competition is performed only on samples with SN_filter = 1 for both measurement methods. In this example, I perform training only on samples wiht SN_filter = 1, which gives a noticeable CV boost but uses only 1/4 of the data (i.e. training on noisy SN_filter = 0 data degrades the performance). A proper consideration of all data as well as reactivity errors may boost the performance.\n\nIn this example, I use a simple CV Kfold split. However, given a mismatch in the RNA length between train/public LB vs. private LB data, **it may be important to verify the effect of the sequence length** to avoid a significant shakeup at the private LB.\n\nOne of the tricks, well known in NLP community, which I use here, is length matching batch sampling: composing batches of samples of approximately the same length to minimize the overhead caused by padding tokens.","metadata":{}},{"cell_type":"markdown","source":"    About 3hours with GPU\n    \n    epoch\ttrain_loss\tvalid_loss\tmae\ttime\n    0\t0.259645\t0.242505\t0.242824\t05:32\n    1\t0.232985\t0.239556\t0.240022\t05:17\n    2\t0.224673\t0.230146\t0.230736\t05:17\n    3\t0.217918\t0.224602\t0.225077\t05:17\n    4\t0.211552\t0.214193\t0.214585\t05:17\n    5\t0.203101\t0.208112\t0.208401\t05:16\n    6\t0.197911\t0.199414\t0.199593\t05:17\n    7\t0.193171\t0.193157\t0.193358\t05:17\n    8\t0.189614\t0.188340\t0.188529\t05:17\n    9\t0.186561\t0.187671\t0.187903\t05:17\n    10\t0.183876\t0.183966\t0.184206\t05:16\n    11\t0.181590\t0.182045\t0.182316\t05:17\n    12\t0.180300\t0.180280\t0.180582\t05:17\n    13\t0.177671\t0.177043\t0.177349\t05:17\n    14\t0.175484\t0.174906\t0.175268\t05:17\n    15\t0.174364\t0.172992\t0.173380\t05:17\n    16\t0.173053\t0.171768\t0.172208\t05:17\n    17\t0.171538\t0.169716\t0.170161\t05:17\n    18\t0.169739\t0.168699\t0.169170\t05:17\n    19\t0.169160\t0.168604\t0.169072\t05:17\n    20\t0.167819\t0.166619\t0.167131\t05:17\n    21\t0.167098\t0.165825\t0.166352\t05:17\n    22\t0.165852\t0.164539\t0.165069\t05:17\n    23\t0.165208\t0.164545\t0.165073\t05:17\n    24\t0.164834\t0.163786\t0.164341\t05:17\n    25\t0.164269\t0.163232\t0.163790\t05:17\n    26\t0.163931\t0.163240\t0.163799\t05:17\n    27\t0.163718\t0.162777\t0.163343\t05:18\n    28\t0.164019\t0.162552\t0.163110\t05:17\n    29\t0.163497\t0.162560\t0.163127\t05:17\n    30\t0.163550\t0.162594\t0.163162\t05:17\n    31\t0.163162\t0.162560\t0.163128\t05:17\n","metadata":{}},{"cell_type":"code","source":"# Broadcast: \n# A = np.zeros( (2,3) )\n# b = np.ones(3)\n# b\n# A+b[np.newaxis,:]","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Submission prepare by column mean\n\n","metadata":{}},{"cell_type":"code","source":"%%time\nl = [len(t) for t in dft['sequence']] \nprint( pd.Series(l).value_counts() )\n\nfor k in range(len(l)-1):\n    if l[k]!=l[k+1]:\n        print(k,l[k],l[k+1])","metadata":{"execution":{"iopub.status.busy":"2023-09-24T20:21:08.145226Z","iopub.execute_input":"2023-09-24T20:21:08.145691Z","iopub.status.idle":"2023-09-24T20:21:09.536895Z","shell.execute_reply.started":"2023-09-24T20:21:08.145655Z","shell.execute_reply":"2023-09-24T20:21:09.535832Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\nll = ['reactivity_'+'%04d'%i for i in range(1,178)]\nprint(ll[:2])\n\nm = df['SN_filter'] == 1\nm1 = df['experiment_type'] == '2A3_MaP'\nv  = df[m&m1][ll].median(axis = 0).fillna(0)\nv[v<0] = 0\nv[v>0.35] = 0.35\nv = v.round(3)\nprint(list(np.round(v.values,3)))\nplt.figure(figsize = (20,4))\nplt.plot(v.values,'*-')\nplt.title('2A3',fontsize = 20)\nplt.show()\n\nv = np.round(v.values, 3)\nY_2A3_pred = np.zeros( (335823, 177), dtype = np.float16 )\nY_2A3_pred += v[np.newaxis,:]\nY_2A3_pred[:3,25:28]","metadata":{"execution":{"iopub.status.busy":"2023-09-24T20:26:07.316935Z","iopub.execute_input":"2023-09-24T20:26:07.317388Z","iopub.status.idle":"2023-09-24T20:26:10.622143Z","shell.execute_reply.started":"2023-09-24T20:26:07.317339Z","shell.execute_reply":"2023-09-24T20:26:10.620957Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\nll = ['reactivity_'+'%04d'%i for i in range(1,178)]\nprint(ll[:2])\n\nm = df['SN_filter'] == 1\nm1 = df['experiment_type'] == 'DMS_MaP'\nv  = df[m&m1][ll].median(axis = 0).fillna(0)\nv[v<0] = 0\nv[v>0.35] = 0.35\nv = v.round(3)\nprint(list(np.round(v.values,3)))\nplt.figure(figsize = (20,4))\nplt.plot(v.values,'*-')\nplt.title('2A3',fontsize = 20)\nplt.show()\n\nv = np.round(v.values, 3)\nY_DMS_pred = np.zeros( (335823, 177) , dtype = np.float16)\nY_DMS_pred += v[np.newaxis,:]\nY_DMS_pred[:3,25:28]","metadata":{"execution":{"iopub.status.busy":"2023-09-24T20:26:25.931311Z","iopub.execute_input":"2023-09-24T20:26:25.931762Z","iopub.status.idle":"2023-09-24T20:26:28.603912Z","shell.execute_reply.started":"2023-09-24T20:26:25.931727Z","shell.execute_reply":"2023-09-24T20:26:28.602612Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nN = 269796671\ndf_submit = pd.DataFrame(index = range(N) )\n","metadata":{"execution":{"iopub.status.busy":"2023-09-24T20:29:00.876000Z","iopub.execute_input":"2023-09-24T20:29:00.876519Z","iopub.status.idle":"2023-09-24T20:29:00.885432Z","shell.execute_reply.started":"2023-09-24T20:29:00.876483Z","shell.execute_reply":"2023-09-24T20:29:00.884395Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndf_submit.index.name = 'id'\ndf_submit['reactivity_DMS_MaP'] = np.zeros( N, dtype = np.float16 )\ndf_submit['reactivity_2A3_MaP'] = np.zeros( N, dtype = np.float16 )\nK = len(Y_DMS_pred.ravel() )\ndf_submit['reactivity_DMS_MaP'].iloc[:K] = Y_DMS_pred.ravel()\ndf_submit['reactivity_2A3_MaP'].iloc[:K] = Y_2A3_pred.ravel()\n","metadata":{"execution":{"iopub.status.busy":"2023-09-24T20:31:54.224857Z","iopub.execute_input":"2023-09-24T20:31:54.225830Z","iopub.status.idle":"2023-09-24T20:31:55.260919Z","shell.execute_reply.started":"2023-09-24T20:31:54.225789Z","shell.execute_reply":"2023-09-24T20:31:55.259830Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_submit.iloc[30:35,:]","metadata":{"execution":{"iopub.status.busy":"2023-09-24T20:33:50.379589Z","iopub.execute_input":"2023-09-24T20:33:50.380037Z","iopub.status.idle":"2023-09-24T20:33:50.392251Z","shell.execute_reply.started":"2023-09-24T20:33:50.379997Z","shell.execute_reply":"2023-09-24T20:33:50.391124Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndf_submit.to_csv('submission.csv',float_format='%.3f')","metadata":{"execution":{"iopub.status.busy":"2023-09-24T20:35:24.683339Z","iopub.execute_input":"2023-09-24T20:35:24.683743Z","iopub.status.idle":"2023-09-24T21:14:46.655082Z","shell.execute_reply.started":"2023-09-24T20:35:24.683712Z","shell.execute_reply":"2023-09-24T21:14:46.651783Z"},"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":[]},{"cell_type":"code","source":"print('%.1f seconds passed total '%(time.time()-t0start) )\nprint('%.1f minutes passed total '%( (time.time()-t0start)/60)  )\nprint('%.2f hours passed total '%( (time.time()-t0start)/3600)  )","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}