{"metadata":{"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"name":"R","codemirror_mode":"r","pygments_lexer":"r","mimetype":"text/x-r-source","file_extension":".r","version":"4.0.5"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"library(tidyverse) # metapackage of all tidyverse packages\nlibrary(seqinr) # for reading .fasta files\nlibrary(ggthemes) # for ggplot themes\nlist.files(path = \"../input/stanford-ribonanza-rna-folding\")\noptions(warn=-1)","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","execution":{"iopub.status.busy":"2023-09-11T10:59:00.675981Z","iopub.execute_input":"2023-09-11T10:59:00.709012Z","iopub.status.idle":"2023-09-11T10:59:00.736194Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Read the train data -- 2.37 GB!\ntrain_data <- read.csv(\"../input/stanford-ribonanza-rna-folding/train_data.csv\")\n\ntrain_data %>% head()","metadata":{"execution":{"iopub.status.busy":"2023-09-11T10:59:04.885657Z","iopub.execute_input":"2023-09-11T10:59:04.887523Z","iopub.status.idle":"2023-09-11T11:03:38.994060Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# This is what we have from the competition hosts regarding the train_data.csv:\n\n* sequence_id - (string) An arbitrary identifier like 8cdfeef00 for each sequence.\n* sequence - (string) Describes the RNA sequence, a combination of A, G, U, and C for each sample. Should be 115 to 457 characters long).\n* experiment_type - (string) Either DMS_MaP or 2A3_MaP to describe the type of chemical mapping experiment that was used to generate each profile. References: DMS, 2A3.\n* dataset_name - (string) name of high throughput sequencing dataset from which the reactivity profile was extracted.\n* reads - (integer) Number of reads in the high throughput sequencing experiment that were assigned to the RNA sequence, and whose mutations were tabulated to compile the reactivity profile.\n* signal_to_noise - (float) Signal/noise value for the profile, defined as mean( measurement value over probed nts )/mean( statistical error in measurement value over probed nts). \n* SN_filter - (Boolean) 0 or 1 depending on whether the profile has signal_to_noise>1.0 and reads>100. For evaluation, only sequences whose DMS_MaP and 2A3_MaP profiles both pass this filter will be used to score submissions.\n* reactivity_0001, reactivity_0002,… - (float) An array of floating point numbers of the train data, should have the same length as the RNA sequence, which defines the reactivity profile for the RNA. For sequences shorter than the maximum RNA length, positions that go beyond the sequence length have null. Several positions early and late in the sequence also cannot be probed due to technical reasons, and their reactivity values are null.\n* reactivity_error_0001, reactivity_error_0002,… - (float) An array of floating point numbers, should have the same length as the corresponding reactivity_* columns, calculated errors in experimental values obtained in reactivity derived from counting statistics in the high-throughput sequencing experiment.\n* reactivity_DMS_MaP, reactivity_2A3_MaP - (float) sample submission values.\n* future - (Boolean) sequences whose data will be collected after the start of the competition (but before final scoring) are labeled as 1.\n","metadata":{}},{"cell_type":"markdown","source":"### When we first look at the dataset, we see there are **1,643,680 rows** of data and **806,573 unique sequence IDs**.  We know from the competition material there are two experiment types, do each sequence ID have both experiements?","metadata":{}},{"cell_type":"code","source":"train_data %>% nrow()\ntrain_data %>% select(sequence_id) %>% unique() %>% nrow()","metadata":{"execution":{"iopub.status.busy":"2023-09-11T11:06:44.347512Z","iopub.execute_input":"2023-09-11T11:06:44.349631Z","iopub.status.idle":"2023-09-11T11:06:44.835259Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 806,573 * 2 is **1,613,146** so some sequence_ids are showing up more than twice for the two experiments.","metadata":{}},{"cell_type":"code","source":"# If we look at the following sequence_IDs, we see that some have the same two experiemnt types, just more of them:\ntrain_data %>% filter(sequence_id == '00005a0b365f') %>% select(sequence_id, experiment_type,reads, dataset_name)\n\ntrain_data %>% filter(sequence_id == '000d960a19ca') %>% select(sequence_id, experiment_type,reads, dataset_name)\n\n# Upon closer examination, it looks like the dataset_name is what brings back more rows.  So, some sequence IDs will have more profiling than others.","metadata":{"execution":{"iopub.status.busy":"2023-09-11T11:06:47.452107Z","iopub.execute_input":"2023-09-11T11:06:47.453827Z","iopub.status.idle":"2023-09-11T11:06:47.599670Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### So how do the dataset_names feature in the dataset?","metadata":{}},{"cell_type":"code","source":"# Adjust plot size\noptions(repr.plot.width = 20, repr.plot.height =12)\n\ntrain_data %>% group_by(dataset_name) %>% summarize(total = n()) %>% select(dataset_name, total) %>% \nggplot(aes(x=reorder(dataset_name, total))) + geom_col(aes(y=total), fill = 'steelblue') + coord_flip() + theme_clean()","metadata":{"execution":{"iopub.status.busy":"2023-09-11T11:06:50.660242Z","iopub.execute_input":"2023-09-11T11:06:50.662084Z","iopub.status.idle":"2023-09-11T11:06:52.214119Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Do the experiment types feature evenly in the dataset?","metadata":{}},{"cell_type":"code","source":"# Yes! Looks like there are even amounts of 2A3_MaP & DMS_MaP experiment types\ntrain_data %>% group_by(experiment_type) %>% summarize(total = n())","metadata":{"execution":{"iopub.status.busy":"2023-09-11T11:06:57.999001Z","iopub.execute_input":"2023-09-11T11:06:58.000738Z","iopub.status.idle":"2023-09-11T11:06:58.075900Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Let's take a look at the reads column -- this may need some more in-depth analysis to determine if 0 is an appropriate value.  From the data, it looks like read == 0 features more prominently in 2A3_MaP experiment types. ","metadata":{}},{"cell_type":"code","source":"# Looking at reads column\ntrain_data %>% select(reads) %>% summary()\n\n# There are 37,867 reads == 0.  If there were no reads assigned to the RNA sequence, did that experiment fail, or did it complete successfully and yield no positive read result?\ntrain_data %>% filter(reads == 0) %>% nrow()\n\n# What do reads == 0 look like from the experiment_type perspective?\ntrain_data %>% filter(reads == 0) %>% group_by(experiment_type) %>% summarize(total = n())","metadata":{"execution":{"iopub.status.busy":"2023-09-11T11:07:01.321262Z","iopub.execute_input":"2023-09-11T11:07:01.322967Z","iopub.status.idle":"2023-09-11T11:07:01.976280Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Summary of signal_to_noise\ntrain_data %>% select(signal_to_noise) %>% summary()","metadata":{"execution":{"iopub.status.busy":"2023-09-11T11:07:05.655905Z","iopub.execute_input":"2023-09-11T11:07:05.657818Z","iopub.status.idle":"2023-09-11T11:07:05.779852Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Let's look at the SN_filter.  From the material, sequence profiles both need to pass this filter to be used to score submissions.  This variable is **TRUE (or 1)** when **signal_to_noise > 1.0 and reads > 100**.  There are approximately 2.75 times more sequence profiles that failed to pass this filter than passed.","metadata":{}},{"cell_type":"code","source":"train_data %>% group_by(SN_filter) %>% summarize(total = n()) %>% select(SN_filter, total)\n\n# Adjust plot size\noptions(repr.plot.width = 8, repr.plot.height =4)\ntrain_data %>% group_by(SN_filter) %>% summarize(total = n()) %>% select(SN_filter, total) %>% ggplot(aes(x=SN_filter)) + geom_col(aes(y=total), fill = 'steelblue') + theme_clean()","metadata":{"execution":{"iopub.status.busy":"2023-09-11T11:07:09.640899Z","iopub.execute_input":"2023-09-11T11:07:09.642934Z","iopub.status.idle":"2023-09-11T11:07:10.139888Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Let's take a look at the reactivity_ and reactivity_error columns: there are 206 of each.","metadata":{}},{"cell_type":"code","source":"# If we want to see all the variables & types, we can use the list.len = ncol(train_data) argument.  \n# This is where I initially saw the different variable types within both reactivity_ & reactivity_error variable sets\n# str(train_data, list.len = ncol(train_data))\n\n# Grab only the 'reactivity_' columns and bring back unique class types\nsapply(train_data[,which(str_detect(colnames(train_data), 'reactivity_0'))], class) %>% unique()\n\n# Grab only the 'reactivity_error' columns and bring back unique class types\nsapply(train_data[,which(str_detect(colnames(train_data), 'reactivity_error'))], class) %>% unique()\n\n# Display the two types for variable sets:\ncolumns <- c('reactivity_0001', 'reactivity_0027', 'reactivity_error_0001', 'reactivity_error_0027')\ntrain_data[,columns] %>% head()\n\n# Upon further scrutiny, even though the competition material states those columns are (float), some of the columns appear to be exclusively NA, thus those are being assigned 'logical' types in R.  We will most likely need to account for those in the future.","metadata":{"execution":{"iopub.status.busy":"2023-09-11T11:07:20.995399Z","iopub.execute_input":"2023-09-11T11:07:20.997127Z","iopub.status.idle":"2023-09-11T11:07:21.047474Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# This is what we have from the competition hosts regarding the additional data folders:\n*     sequence_libraries - all 2,150,401 sequences, with titles, in FASTA format, grouped into the actual collections that were synthesized together. Note: some files list DNA not RNA sequences (T instead of U).\n*     supplementary_silico_predictions - CSV files of predicted RNA secondary structures in dot-parentheses notations. Note: not all packages were run for all sequences due to computational expense of prediction packages that are able to predict pseudoknots (non-nested pairings).\n*     eterna_openknot_metadata - TSV files of id, name (author), title (author-provided design title), body (description), and sequence for Eterna OpenKnot sequences. Note: some files have additional tabs in title or description that may affect readin.\n*     Ribonanza_bpp_files - TXT file of base pair probabilities from LinearPartition-EternaFold package for train and test sequences. Indexed by sequence_id. Note: this package simulates RNA secondary structure ensembles without pseudoknots or other tertiary structure features.","metadata":{}},{"cell_type":"markdown","source":"### Initial thoughts on the additional data folders: reading in the .fasta files does take some time.  Will require some additional analysis / research to determine a time-efficient method of reading in the additional files.  Some draft approaches are below:","metadata":{}},{"cell_type":"code","source":"# Read the sequence_libraries...there are 9. Reading the list of files doesn't take too long, but reading the actual files took quite some time....so be prepared!\n#fasta_list <- list.files(path = \"../input/stanford-ribonanza-rna-folding/sequence_libraries\",full.names = TRUE, recursive = TRUE)\n\n# Placeholder for reading the 9 fasta files","metadata":{"execution":{"iopub.status.busy":"2023-09-11T10:38:18.674657Z","iopub.execute_input":"2023-09-11T10:38:18.676068Z","iopub.status.idle":"2023-09-11T10:38:18.687879Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Read the eterna_openknot_metadata tsv flies...there are 6.  These don't take too long.\ntsv_list <- list.files(path = \"../input/stanford-ribonanza-rna-folding/eterna_openknot_metadata\",full.names = TRUE, recursive = TRUE)\n\n# Set up the tsv_df\ncol_names <- c('id','title','name','body','sequence')\ntsv_df <- data.frame(matrix(nrow=0,ncol=5))\ncolnames(tsv_df) <- col_names\n\nfor(i in 1:length(tsv_list)){\n    tsv_file <- read_tsv(tsv_list[4], col_names = TRUE,show_col_types = FALSE) %>% as.data.frame()\n    tsv_df <- rbind(tsv_df, tsv_file)\n}","metadata":{"execution":{"iopub.status.busy":"2023-09-11T11:07:40.673873Z","iopub.execute_input":"2023-09-11T11:07:40.675527Z","iopub.status.idle":"2023-09-11T11:07:41.260307Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Read the Ribonanza_bpp_files.  They are .txt, and there are 2150396!!\n#ribo_list <- list.files(path = \"../input/stanford-ribonanza-rna-folding/Ribonanza_bpp_files\",full.names = TRUE, recursive = TRUE)\n\n# Place holder for reading the 2150396 .txt files","metadata":{"execution":{"iopub.status.busy":"2023-09-11T10:38:18.926716Z","iopub.execute_input":"2023-09-11T10:38:18.928269Z","iopub.status.idle":"2023-09-11T10:38:18.940651Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Read the test_sequences.csv\n#test_sequences <- read.csv(\"../input/stanford-ribonanza-rna-folding/test_sequences.csv\")","metadata":{"execution":{"iopub.status.busy":"2023-09-11T10:38:18.944588Z","iopub.execute_input":"2023-09-11T10:38:18.946061Z","iopub.status.idle":"2023-09-11T10:38:18.958304Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Read the sample submission -- 3.67 GB!\n#sample_submission <-read.csv(\"../input/stanford-ribonanza-rna-folding/sample_submission.csv\")","metadata":{"execution":{"iopub.status.busy":"2023-09-11T11:07:51.088172Z","iopub.execute_input":"2023-09-11T11:07:51.090295Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Write the sample submission\n#write_csv(sample_submission, 'submission.csv')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Hope this was helpful!  If so, please ***upvote*** or leave a ***comment*** to make it better!  Will continue to add more as time permits!","metadata":{}}],"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"name":"R","codemirror_mode":"r","pygments_lexer":"r","mimetype":"text/x-r-source","file_extension":".r","version":"4.0.5"}}