{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":51294,"databundleVersionId":6923401,"sourceType":"competition"}],"dockerImageVersionId":30587,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\n# for dirname, _, filenames in os.walk('/kaggle/input'):\n#     for filename in filenames:\n#         print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-11-27T05:55:59.673011Z","iopub.execute_input":"2023-11-27T05:55:59.673986Z","iopub.status.idle":"2023-11-27T05:56:00.176019Z","shell.execute_reply.started":"2023-11-27T05:55:59.673941Z","shell.execute_reply":"2023-11-27T05:56:00.174783Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To predict the reactivity of an RNA sequence to chemical modifiers :\n- RNA Secondary Structure Prediction: RNAfold from the ViennaRNA package or RNAstructure to predict the secondary structure of your RNA sequence. The secondary structure can influence the reactivity of specific regions.\n- Accessibility Prediction: Predict the accessibility or solvent accessibility of different regions in the secondary structure. More accessible regions may be more prone to chemical modification. Tools like DMS-MaPseq (Dimethyl Sulfate-Mutational Profiling with Sequencing) or SHAPE-MaP (Selective 2'-Hydroxyl Acylation analyzed by Primer Extension-MaP) for experimental accessibility data.\n- Biophysical Modeling:biophysical models that incorporate known chemical modifiers and their preferences for different RNA structures. These models might consider factors like base pairing, stacking, and loop structures.","metadata":{}},{"cell_type":"code","source":"from Bio.Seq import Seq \nx = Seq(\"ATGGCGGGAAAATGA\") \nx","metadata":{"execution":{"iopub.status.busy":"2023-11-27T05:06:15.649473Z","iopub.execute_input":"2023-11-27T05:06:15.649909Z","iopub.status.idle":"2023-11-27T05:06:15.657483Z","shell.execute_reply.started":"2023-11-27T05:06:15.649875Z","shell.execute_reply":"2023-11-27T05:06:15.656148Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#GC content of a sequence is the percentage of the nucleotides in the sequence that are either G or C\nfrom Bio.SeqUtils import GC \nfrom Bio.Seq import Seq \nx = Seq(\"ATGGCGGGAAAATGA\") \nGC(x) ","metadata":{"execution":{"iopub.status.busy":"2023-11-27T05:11:40.254895Z","iopub.execute_input":"2023-11-27T05:11:40.255851Z","iopub.status.idle":"2023-11-27T05:11:40.264760Z","shell.execute_reply.started":"2023-11-27T05:11:40.255792Z","shell.execute_reply":"2023-11-27T05:11:40.262971Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# with Biopython we can create objects corresponding to DNA, RNA, and Protein sequences\nfrom Bio.Seq import Seq \nDNA = Seq(\"ATGCTGGGATATTGA\") \nmRNA = DNA.transcribe() \nprint(mRNA) \nmRNA \nprotein = mRNA.translate() \nprint(protein) \nprotein","metadata":{"execution":{"iopub.status.busy":"2023-11-27T05:15:48.200997Z","iopub.execute_input":"2023-11-27T05:15:48.202096Z","iopub.status.idle":"2023-11-27T05:15:48.211552Z","shell.execute_reply.started":"2023-11-27T05:15:48.202048Z","shell.execute_reply":"2023-11-27T05:15:48.210437Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#search for an exact string or pattern within a larger sequence\n#when using a restriction enzyme to cut a larger sequence at sequence-specific sites that match a particular pattern,\n#we would like know where this pattern occurs in the larger sequence. \n#This task can be called string matching. \nfrom Bio.Seq import Seq \nfrom Bio import SeqUtils \npattern = Seq(\"ACG\") \nsequence = Seq(\"ATGCGCGACGGCGTGATCAGCTTATAGCCGTACGACTGCTGCAACGTGACTGAT\") \nresults = SeqUtils.nt_search(str(sequence),pattern) \nprint(results) ","metadata":{"execution":{"iopub.status.busy":"2023-11-27T05:16:37.889042Z","iopub.execute_input":"2023-11-27T05:16:37.889475Z","iopub.status.idle":"2023-11-27T05:16:37.896220Z","shell.execute_reply.started":"2023-11-27T05:16:37.889444Z","shell.execute_reply":"2023-11-27T05:16:37.894934Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from Bio.Seq import Seq\nfrom Bio.SeqUtils import GC\n\n# Your RNA sequence\nrna_sequence = \"AUGUAGCGUAGCGUACGUACGUAGCUAGCUAGCUAGCUACGUA\"\n\n# Create a Biopython Seq object\nrna_seq = Seq(rna_sequence)\n\n# Calculate GC content\ngc_content = GC(rna_seq)\n\nprint(f\"GC Content: {gc_content}%\")\n","metadata":{"execution":{"iopub.status.busy":"2023-11-27T05:21:46.534433Z","iopub.execute_input":"2023-11-27T05:21:46.534850Z","iopub.status.idle":"2023-11-27T05:21:46.541863Z","shell.execute_reply.started":"2023-11-27T05:21:46.534819Z","shell.execute_reply":"2023-11-27T05:21:46.540374Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install ViennaRNA","metadata":{"execution":{"iopub.status.busy":"2023-11-27T05:22:38.048772Z","iopub.execute_input":"2023-11-27T05:22:38.049203Z","iopub.status.idle":"2023-11-27T05:22:52.095026Z","shell.execute_reply.started":"2023-11-27T05:22:38.049168Z","shell.execute_reply":"2023-11-27T05:22:52.093535Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.metrics import mean_squared_error\nimport matplotlib.pyplot as plt\n\nnp.random.seed(42)\nn_samples = 1000\n\n# RNA sequence features\nrna_sequences = np.random.rand(n_samples, 10)\n\n# Chemical modifier features\nmodifier1 = np.random.rand(n_samples, 1)\nmodifier2 = np.random.rand(n_samples, 1)\n\n# Combine features\nX = np.hstack((rna_sequences, modifier1, modifier2))\n\n# Simulate reactivity values (target variable)\n# replace with actual dataset\ntrue_coefficients = np.array([0.5, -0.2, 0.3, 0.1, -0.4, 0.2, 0.8, 0.6, -0.3, -0.7, 0.4, 0.2])\nnoise = 0.1 * np.random.randn(n_samples)\ny = np.dot(X, true_coefficients) + noise\n\n# Split the dataset into training and testing sets\nX_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)\n\nmodel = LinearRegression()\nmodel.fit(X_train, y_train)\n\ny_pred = model.predict(X_test)\n\n# Evaluate the model\nmse = mean_squared_error(y_test, y_pred)\nprint(f\"Mean Squared Error: {mse}\")\n\n# Visualize the predicted vs. actual values\nplt.scatter(y_test, y_pred)\nplt.xlabel(\"True Reactivity\")\nplt.ylabel(\"Predicted Reactivity\")\nplt.title(\"Predicted vs. Actual Reactivity\")\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-11-27T05:30:07.249704Z","iopub.execute_input":"2023-11-27T05:30:07.250221Z","iopub.status.idle":"2023-11-27T05:30:07.562436Z","shell.execute_reply.started":"2023-11-27T05:30:07.250170Z","shell.execute_reply":"2023-11-27T05:30:07.559467Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import random\nimport RNA\n\n# Generate a random RNA sequence\nsequence_length = 100\nsequence = ''.join(random.choice('ACGU') for _ in range(sequence_length))\n\n# Predict the secondary structure of the RNA\nstructure = RNA.fold(sequence)\n\n# Calculate the reactivity of each nucleotide to DMS and CMCT\ndms_reactivity = []\ncmct_reactivity = []\n\n# Ensure that the length of the sequence matches the length of the structure\nsequence_length = min(sequence_length, len(structure))\n\nfor i in range(sequence_length):\n    # Check if the nucleotide is unpaired in the secondary structure\n    if structure[i] == '.':\n        # If the nucleotide is unpaired, it is more likely to be reactive to DMS\n        dms_reactivity.append(1)\n        # If the nucleotide is unpaired, it is also more likely to be reactive to CMCT\n        cmct_reactivity.append(1)\n    else:\n        # If the nucleotide is paired, it is less likely to be reactive to DMS\n        dms_reactivity.append(0)\n        # If the nucleotide is paired, it is also less likely to be reactive to CMCT\n        cmct_reactivity.append(0)\n\n# Print the results\nprint(\"DMS reactivity:\", dms_reactivity)\nprint(\"CMCT reactivity:\", cmct_reactivity)\n","metadata":{"execution":{"iopub.status.busy":"2023-11-27T05:36:15.989424Z","iopub.execute_input":"2023-11-27T05:36:15.989842Z","iopub.status.idle":"2023-11-27T05:36:16.006939Z","shell.execute_reply.started":"2023-11-27T05:36:15.989810Z","shell.execute_reply":"2023-11-27T05:36:16.005802Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import random\n\ndef generate_random_rna_sequence(length):\n    return ''.join(random.choice('ACGU') for _ in range(length))\n\nrna_sequence = generate_random_rna_sequence(50) \nprint(\"Random RNA Sequence:\", rna_sequence)\n","metadata":{"execution":{"iopub.status.busy":"2023-11-27T05:38:42.959099Z","iopub.execute_input":"2023-11-27T05:38:42.959543Z","iopub.status.idle":"2023-11-27T05:38:42.966604Z","shell.execute_reply.started":"2023-11-27T05:38:42.959512Z","shell.execute_reply":"2023-11-27T05:38:42.965294Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.metrics import mean_squared_error\nimport matplotlib.pyplot as plt\n\nnp.random.seed(42)\ndata_size = 100\nrna_sequences = np.random.rand(data_size, 10)  # 10 features as an example\nreactivity_dms = 3 * rna_sequences[:, 0] + 2 * np.random.randn(data_size)\nreactivity_2a3 = 2 * rna_sequences[:, 1] + 3 * np.random.randn(data_size)\n\n# DataFrame from the generated data\ndf = pd.DataFrame({\n    'RNA_Sequence': [''.join(map(str, seq)) for seq in rna_sequences],\n    'Reactivity_DMS': reactivity_dms,\n    'Reactivity_2A3': reactivity_2a3\n})\n\n# Split the dataset into training and testing sets\nX = rna_sequences  # Features\ny_dms = reactivity_dms  # Target for DMS\ny_2a3 = reactivity_2a3  # Target for 2A3\n\nX_train, X_test, y_dms_train, y_dms_test, y_2a3_train, y_2a3_test = train_test_split(\n    X, y_dms, y_2a3, test_size=0.2, random_state=42\n)\n\n# Train linear regression models\nmodel_dms = LinearRegression()\nmodel_2a3 = LinearRegression()\n\nmodel_dms.fit(X_train, y_dms_train)\nmodel_2a3.fit(X_train, y_2a3_train)\n\n# predictions on the test set\ny_dms_pred = model_dms.predict(X_test)\ny_2a3_pred = model_2a3.predict(X_test)\n\n# Evaluate the models\nmse_dms = mean_squared_error(y_dms_test, y_dms_pred)\nmse_2a3 = mean_squared_error(y_2a3_test, y_2a3_pred)\n\nprint(f'Mean Squared Error (DMS): {mse_dms}')\nprint(f'Mean Squared Error (2A3): {mse_2a3}')\n\nplt.scatter(y_dms_test, y_dms_pred, label='DMS')\nplt.scatter(y_2a3_test, y_2a3_pred, label='2A3')\nplt.xlabel('True Values')\nplt.ylabel('Predictions')\nplt.legend()\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-11-27T05:41:17.589969Z","iopub.execute_input":"2023-11-27T05:41:17.590393Z","iopub.status.idle":"2023-11-27T05:41:17.867429Z","shell.execute_reply.started":"2023-11-27T05:41:17.590362Z","shell.execute_reply":"2023-11-27T05:41:17.866028Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"file = \"/kaggle/input/stanford-ribonanza-rna-folding/train_data.csv\"\ndf = pd.read_csv(file,nrows=1000)\ndf.head()","metadata":{"execution":{"iopub.status.busy":"2023-11-27T05:56:05.582619Z","iopub.execute_input":"2023-11-27T05:56:05.583207Z","iopub.status.idle":"2023-11-27T05:56:05.769476Z","shell.execute_reply.started":"2023-11-27T05:56:05.583168Z","shell.execute_reply":"2023-11-27T05:56:05.768346Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# # train a Multilayer Perceptron (MLP) regressor model\n# from sklearn.model_selection import train_test_split\n# from sklearn.preprocessing import StandardScaler\n# from sklearn.neural_network import MLPRegressor\n\n# def load_data(filename):\n#     data = pd.read_csv(filename)\n#     X = data[['sequence']]\n#     y = data[['DMS_reactivity', '2A3_reactivity']]\n#     return X, y\n\n# def preprocess_data(X, y):\n#     # One-hot encode RNA sequences\n#     sequences = X['sequence'].str.split('')\n#     one_hot_encoded_sequences = []\n#     for sequence in sequences:\n#         encoded_sequence = np.zeros((4, len(sequence)))\n#         for i, nucleotide in enumerate(sequence):\n#             if nucleotide == 'A':\n#                 encoded_sequence[0, i] = 1\n#             elif nucleotide == 'C':\n#                 encoded_sequence[1, i] = 1\n#             elif nucleotide == 'G':\n#                 encoded_sequence[2, i] = 1\n#             elif nucleotide == 'U':\n#                 encoded_sequence[3, i] = 1\n#         one_hot_encoded_sequences.append(encoded_sequence)\n#     X = np.array(one_hot_encoded_sequences)\n\n#     # Standardize reactivity values\n#     scaler = StandardScaler()\n#     scaler.fit(y)\n#     y = scaler.transform(y)\n\n#     return X, y\n\n# def train_model(X, y):\n#     # Split data into training and testing sets\n#     X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)\n\n#     # Create and train the MLPRegressor model\n#     model = MLPRegressor(hidden_layer_sizes=(128, 64), max_iter=1000)\n#     model.fit(X_train, y_train)\n\n#     return model\n\n# def evaluate_model(model, X_test, y_test):\n#     # Evaluate the model on the test set\n#     predictions = model.predict(X_test)\n#     MSE = np.mean((predictions - y_test) ** 2)\n#     print('Mean Squared Error:', MSE)\n\n# def main():\n#     # Load data\n#     X, y = load_data('/kaggle/input/stanford-ribonanza-rna-folding/train_data.csv')\n\n#     # Preprocess data\n#     X, y = preprocess_data(X, y)\n\n#     model = train_model(X, y)\n\n#     evaluate_model(model, X, y)\n\n# if __name__ == '__main__':\n#     main()\n","metadata":{"execution":{"iopub.status.busy":"2023-11-27T05:56:43.454120Z","iopub.execute_input":"2023-11-27T05:56:43.454584Z","iopub.status.idle":"2023-11-27T05:56:43.462388Z","shell.execute_reply.started":"2023-11-27T05:56:43.454535Z","shell.execute_reply":"2023-11-27T05:56:43.461114Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}