Computational Analysis of Sequencing and Proteomic Data to Understand Sub-Cellular Protein Localizations, Functions, and Clinical Value
Advisor: Dr. Sharma
Internship Location: UCSF Mission Bay Campus
Onsite Mentor: Biao Wang
Abstract: In this era of big data, vast biomedical datasets, including genomic, transcriptomic, and clinical, are being generated daily. Thus, there is an increasing need to develop robust computational tools to analyze these datasets to extract biological value, whether understanding the spread of disease in populations or modeling complex cellular processes. At the pure-biology Wang Lab at the University of California, San Francisco, I will work with Principal Investigator Dr. Biao Wang to analyze lab-produced and public datasets to understand phosphorylation regulatory mechanisms, protein shuttling, and metabolic homeostasis. I will work to understand how specific protein sequences, called nuclear localization sequences (NLSs), are phosphorylated by kinases to regulate protein nucleocytoplasmic shuttling, i.e., the movement of proteins between the nucleus and cytoplasm. My project involves conducting an extensive literature review to understand pathways of nuclear import and export. It involves handling numerous classes of data, including protein sequence, phosphoproteomic, and annotation data, as well as corresponding genomic data. Understanding phosphorylation regulatory mechanisms can help us understand how protein sub-cellular localizations are impaired in disease.
Project Blog
Post 1: Decoding ADAR1: Unraveling RNA Editing in Fat Cells and Obesity
Hello! My name is Aashna. Welcome to my introductory blog post! Here, I’d like to talk to you about what excites me about my research project, how I got involved, and provide a bit of background on the biology.
My journey into the UCSF Wang Lab has been a long-winded one – from volunteering at the Sikh Temple medical clinic to connecting with a UCSF post-doc at a medical conference and shadowing under her to conducting a data analysis project under the Principal Investigator this past summer! My overarching goal at this pure-biology lab is to analyze their vast fat cell datasets to give them significance in the context of metabolic research.
One of the lab’s key projects is to understand the mechanism of action of the enzyme ADAR1. ADAR1 is the star of the show in this research project, so make sure to remember that it is an RNA-editing enzyme (operates on double-stranded RNAs (dsRNAs), replacing select “A” nucleotides with “G” nucleotides to mark them as “self”). Its role has been well-explored in viral infections and cancer, but not so much in fat cell biology. That’s where we come in!
This past summer, at the UCSF Wang Lab, I had a lot of spare time as I waited for computationally intensive code segments to run on the lab’s mouse fat data to identify ADAR1’s RNA-editing sites. What else to do with this time than look for public datasets that might complement the lab’s data? So, I found a Finnish human metabolic dataset and analyzed it in different ways. What I definitely was not expecting to find was that the inflammatory pathway, interferon-alpha signaling, regulated by the enzyme ADAR1, is enriched in obesity! Eureka! Doors of possibilities opened for me and the lab about investigating connections between inflammation and obesity!
Armed with the knowledge I gained through last summer’s trial and error exploration of the lab’s data, I now return to the task of statistically analyzing the lab’s data to understand the mechanism of action of ADAR1 in mouse brown and white fat cells by identifying its RNA-editing sites. I hope such results can give us a glimpse into ADAR1’s complex cellular function and how perturbations to its function can contribute to obesity.
Week 1: Literature Review & Coding!
Hello everyone! Welcome to the first week of my Senior Project.
This week, I was busy reading research papers, drawing maps of signaling pathways, writing (and debugging) code, and learning how to use the Wynton HPC Compute Cluster (!!) (as well as vibing to music on the BART to UCSF…)
My project has taken me into the vibrant, complex world of the RNA-editing enzyme ADAR1. By reading around 8 original research and review papers this week, I learned that ADAR1 contains two primary isoforms (ADAR1p110 and ADAR1p150), each with a distinct cellular localization, range of functions, and disease implications. I explored recent studies investigating how ADAR1 correlates with obesity (its overexpression inhibits fat cell differentiation, which helps inhibit the progression of obesity) and its role in autoimmune diseases, viral infections, and cancer. I have also learned how ADAR1 interacts with other cellular proteins, notably MDA5, PKR, and DICER.
The papers I read this week illuminate ADAR1’s interactions with different, sometimes seemingly unrelated, proteins, which made it difficult for me to really understand the full picture. So, I embarked on the daunting task of synthesizing these papers to paint a clear (albeit complicated) web of pathway interactions. After a lot of trial and error (many pencil lines were erased and redrawn…) and revision by my mentor, I was able to draw up a map of ADAR1’s complex and multifaceted interactions, which will help me interpret my computational results later on!
My computational project aims to understand ADAR1’s mechanism of action by identifying its RNA-editing sites and RNA-editing-site clusters. Identification of RNA-editing sites from FASTQ files is a multi-step, lengthy process–one look at a never-ending FASTQ file and I knew I was in for a ride. This week, I surveyed multiple statistical packages that can help me do a fraction of this task, including REDItools, SPRINT, and RNAEditor, and I identified their pros and cons. I then explored a more recent tool developed by Torkler et al., called local differential editing index (LoDEI), that improves upon these existing packages in multiple ways–for one, it is a window-based detection approach, which accounts for the fact that ADAR1 differential editing between two experimental conditions may not be at the exact same site. Also, since it directly identifies differential editing sites, manual removal of SNPs is not needed! Wonderful! So, I will endeavor to use this tool to help me with RNA-editing site detection, and then do further analysis on its output files.
What I did this week in the way of coding was mostly data preprocessing and learning to use the compute cluster. Attending UCSF Wynton HPC Cluster office hours on Tuesday, I learned how to set up VSCode in the virtual environment, enabling easy editing/running of code! I generated genome indices for the lab’s mouse reference genome on UNIX–essentially, this means I annotated the genome with its functional characteristics–coding sequences, noncoding sequences, introns, exons, 5’ UTRs, 3’ UTRs...I am in the middle of preprocessing the lab’s knockout and wildtype ADAR1 mouse brown and white fat cell data (128 Gigabytes…it’s taking a while) – what this entails is aligning the raw FASTQ reads to the mouse reference genome to produce BAM files and then sorting these BAM files. Once this code finishes running, I can begin implementing LoDEI.
My PI also introduced me to a new project the lab is doing: understanding significant protein translocations between the nucleus and cytoplasm in response to environmental and chemical stimuli. This connects to my research on ADAR1, as the cellular localization of ADAR1 isoforms is tightly regulated: the isoform I am studying, ADAR1p150, has both a nuclear localization sequence (NLS) and a nuclear export sequence (NES). So, this is another project on my radar–I will work to identify proteins with NLS and NES sequences from a Nature mass spectrometry databank since this can potentially deliver insight into ADAR1’s cellular localization and activity. For now, this project is on the back burner as I work on the ADAR1 RNA-editing sites project, but I’ve started reading mass spectrometry papers to gain background knowledge on the topic.
Until next week!
Week 2: Detecting RNA Editing Sites and Understanding Protein Shuttling
Hello, and welcome back to my Senior Project! This week was full of reading papers, writing Bash and Python scripts, and brainstorming algorithmic workflows (as well as discovering new music genres on BART!).
Let’s dive right in.
For my ADAR1 RNA-editing site detection project, I realized that I generated some files incorrectly last week. The reason for this is that the GFF file (a descriptive file with genome annotations) that I was using was for the wrong mouse genome build. So, I first downloaded the correct GFF file from the Ensembl database. Then, I re-generated the reference genome indices and re-created the sequencing BAM files for the knockout and wildtype ADAR1 brown fat cells. These scripts took a few days to run. From there, I was finally ready to use the LoDEI (local differential editing index) tool to identify differential A to G RNA-editing sites between knockout and wildtype conditions! This process involved multiple iterations–for example, I learned that LoDEI only handles single-read RNA-seq data, so I had to modify its code to handle the lab’s paired-end data. I also had to process the mouse reference genome’s GFF file (removing the sex chromosome annotations) so that it was compatible with the LoDEI code.
I obtained preliminary results using LoDEI, showing differential A to T editing between knockout and wildtype ADAR1 brown fat cells. However, these results are not fully in line with expectations of ADAR1 editing–the A to G differential editing sites (the type of editing ADAR1 does) that I identified have high false discovery rates, which is unexpected. After meeting with my mentor, Dr. Wang, we decided that I will test LoDEI on the white fat cells to see how they compare to the brown fat results. Also, Dr. Wang shared details regarding how the mouse fat samples were prepared for RNA-sequencing. Learning about these, I realized there are some parameter adjustments that I need to make to LoDEI to run the code correctly (importantly, the RNA-seq data is unstranded, not reverse-stranded as I had assumed). So, while I have shown a basic proof-of-concept working of LoDEI, the results will likely be revised over the next week.
Switching gears to the protein translocation project I introduced last week! I read many lengthy (50+ page) papers on the topic, which helped me build background knowledge on pathways of nuclear import and export of proteins and how phosphorylation regulates these shuttling events. My mentor shared a Nature mass spectrometry dataset consisting of phosphosites across the entire proteome. After a lot of brainstorming, I devised a multi-step methodology to analyze this dataset to identify phosphorylation-regulated shuttling proteins, specifically focusing on shuttling events from the cytoplasm to the nucleus. I will first identify nuclear localization signals (NLSs) in protein amino acid sequences and then find phosphorylated serine/threonine residues near or within these sequences. As part of this workflow, I will use a bioinformatics tool called cNLS mapper (which predicts NLSs in an input protein sequence), developed originally for budding yeast cells! Luckily, the nuclear import pathway we are most interested in (importin α/β) is highly conserved in eukaryotes, so it will still make high accuracy predictions on mouse or human protein data. I tested cNLS mapper on well-known phosphorylation-regulated shuttling proteins (HDAC4 and FOXO1), verifying that it works. I also wrote code to parse through the Nature mass spectrometry dataset and filter it down from 20,000+ to 10,000 functionally-significant phosphosites to analyze. I conducted filtering based on p-value and q-value calculations.
I expect my code to continue running over the weekend, as the Wynton cluster is currently down (my terminal is full of “queue, waiting” messages). I’m excited to return next week to share updates on this project!
Week 3: Automating Data Collection and An Interesting Research Paper!
Thanks for joining me in week 3! Let’s get into some interesting science.
Thankfully, the Wynton HPC Cluster was working again on Monday! No more endless “queue, waiting” messages in the terminal. Hooray!
Last week, my mentor and I decided that I would try to run LoDEI on the white fat data to see how it compares to the brown fat results that I obtained. So, I carried out the same process to unpack and preprocess the white fat data as I did for the brown fat data – unzipping FASTQ files, aligning them to the reference genome, and generating BAM files. The process went smoothly this time, as I had hashed out all the errors with the brown fat data. Then, I re-ran LoDEI on both the brown and white fat datasets, this time with the revised parameters (the BAM files are unstranded, not stranded). The results I obtained are still inconclusive–LoDEI shows q-value significant differential editing from A to T and G to T nucleotides in the brown fat, and no significant differential editing in the white fat. Meeting with my mentor, we decided that I will try to find a public dataset consisting of 2 distinct ADAR1-related conditions in mice (ex: ADAR1 wildtype and ADAR1 deficient) and see if LoDEI produces anticipated A to G differential editing results. If it does, this might indicate an error in how the RNA-seq data I’m working with was collected in the first place. Furthermore, LoDEI was designed originally for human data, and it uses q-values as a ranking measure. It is possible that these design aspects are introducing incompatibilities with the lab's data, which is mouse fat. For now, these are the results I have produced, and it is possible we may need to move on and focus more on the next (below!) project.
Onto exciting developments! This week, I made significant progress on the protein shuttling project. While last week, I developed a framework for how I plan to go about this project, the challenge I was confronted with at the beginning of this week was: I am using tools (like cNLS Mapper and UniProt protein sequence retriever) that have web interfaces. To test them out on a few proteins, I can simply search up uniprot IDs in the UniProt website, copy-paste them into cNLS Mapper, obtain HTML results, and visually inspect the Nature phosphoproteome database to identify phosphosites located in/next to the predicted NLSs identified by cNLS Mapper. How do I automate the process of finding significant phosphosites in the massive (10,000+ proteins!) Nature databank that I’m working with?
The answer: APIs and a lot of web-scraping! I used the UniProt API to generate fasta files for the 10,000+ proteins in the Nature dataset. From there, inputting and retrieving data from cNLS Mapper was a challenge because it doesn’t have an API. So, I wrote code to mimic a user submitting a form (with the fasta files I generated) on the website for all of the proteins. Next, I saved the resulting HTML file results (I used parallel processing to speed this up). Then, I used web-scraping to extract the necessary data from these HTML files: specifically, the predicted NLS sequences, their locations, their confidence measures, and their class (monopartite or bipartite). Finally, I wrote a Python function to identify phosphorylated serine/threonine sites in the Nature dataset that lie within or near these NLSs. I obtained a total of 9,136 significant phosphosites across 2,746 unique proteins from this analysis. After sharing these results with my mentor, we decided to now check the sub-cellular localizations of the proteins that contain these significant phosphosites using the Human Protein Atlas (the proteins must be localized to the nucleus, cytosol, or both to be relevant), and the next exciting project will be to conduct sequence motif analysis (finding consensus phosphorylation sites in these proteins!) By doing this, we can better understand if these are certain conserved sequences related to nuclear import that protein kinases operate on to regulate protein shuttling. As I produce more filtered data, new research questions open for how we can use this data to investigate biological phenomena.
And a serendipitous happening: I stumbled across an insightful research article from the Koch Institute (published just 4 days ago!). Researchers developed a new machine-learning tool, ProtGPS (https://github.com/pgmikhael/protgps), to predict protein localization in the cell from its amino acid sequence. It was a very interesting read, and I’m starting to think about ways to apply this tool to the Nature dataset I’m analyzing. In this way, my project can not only identify NLS sequences in proteins but relate these sequences to proteins' predicted intra-cellular localizations.
Until next week!
Week 4: Exploring New Tools, Writing Code that Opens Google Chrome, and More
Welcome back to my Senior Project! This week, I learned about sequence motif analysis and debugged my LoDEI workflow for RNA-editing-site detection.
Let’s begin with the protein translocation project!
Last week, I identified a list of around 800 phosphosites, located in or near protein nuclear localization sequences (NLSs), that are potentially crucial in regulating nuclear import of proteins. My task, from there, was to identify which kinases phosphorylate these sites, and how these kinases differ between classes of phosphosites. To do this, I pursued sequence motif analysis.
What is sequence motif analysis? Sequence motifs are amino acid sequence patterns widespread and conserved across proteins. Highly conserved motifs are likely key to biological function and can elucidate insights into evolution. As sequence motif analysis was a completely new area for me, I began with literature review. One notable paper that my mentor introduced me to was “An Atlas of Substrate Specificities for the Human Serine/Threonine Kinome,” in which the authors developed a web tool that takes as input an amino acid sequence and outputs kinases (with scores indicating likelihood/match) most likely to phosphorylate it. I also explored tools like WebLogo, which generates visualizations (called “sequence logos”) of sequence motifs, and Clustal Omega, which uses Hidden Markov Models to align sequences to each other to identify patterns.
I began by applying Clustal Omega and WebLogo to my filtered list of 800 phosphosites. This required me to first parse through 700 or so protein fasta files to identify flanking sequences of these phosphosites, and then try to align the sequences to each other. Clustal Omega showed poor alignment, and Web Logo showed some conserved amino acids in the centers of the sequences (mostly serine and threonine). After talking to my research mentor, I learned that I will actually need to apply these tools at the end of the sequence motif analysis process – it is unlikely that there is a single consensus sequence in my list of 800 phosphosites. Instead, we first need to group the phosphosites into distinct classes. I sought to do this by identifying the kinases that bind to them.
To predict the kinases most likely to bind to each phosphorylation sequence, I used a tool called the Kinase Library. This was challenging – it only has a web interface and no API, and the website crashes whenever I upload a text file with greater than 15 rows. So, I had to be a little creative – it is impossible (or at least, highly unappealing) to generate 800/15 ~ 53 text files and manually upload each one of them to the Kinase Library to output TSV files of most likely kinases. I discovered a Python package called Selenium that helps me automate this task! Using this package, I can write code that mimics web browser submissions. After a lot of trial and error and debugging, I was able to write this script. It essentially loops through my phosphorylation sequence text files. For each one, it opens a Google Chrome browser, uploads the text file, displays the output, and scrolls down the webpage to click the "download TSV file results." This all happens in the background, allowing me to use my time to do more important things, like read papers and interpret the results of my code in a biological context. The code is still running, and I will hopefully provide updates on the output next week!
For the ADAR1 RNA-editing sites project, I found a few datasets I can use to establish a positive control. The first is a 2021 PubMed RNA-seq dataset on Parkinson’s Disease in healthy and diseased individuals, which was used to identify differential A-to-I editing between these 2 conditions. The second dataset is a 2020 RNA-seq dataset on patients with and without Huntington’s Disease. Both are human datasets; I am still searching for mouse datasets. Also, upon closer re-reading of the LoDEI paper, I realized there are some quality filtering steps I missed, like FASTQC trimming, used when there is noise in the data. So, I will return to the lab data and try additional quality filtering on the raw FASTQ reads to see if this improves the results.
Finally, I’m working to configure the ProtGPS tool I introduced last week–it is difficult to implement on the cluster because it requires GPU resources, but I’m getting help from the Wynton HPC support team to figure this out. Using ProtGPS would add greater context to my findings for the protein translocation project.
Thanks, and see you again next week!
Week 5: Steady Progress!
Welcome back to my Senior Project!
Last week, I started running code to collect kinase data for my phosphorylation sites using The Kinase Library website tool. My focus at the beginning of this week was to finish collecting all of this kinase data for further analysis. I encountered some challenges with the Kinase Library tool, as it would occasionally crash because of the volume of the data I was inputting, despite my parallelization attempts. However, at the end, I was able to produce around 90 kinase text files for my 800+ phosphosites. I then put all the relevant data from these text files into a format that is easier to read and analyze: a single CSV file, with columns for kinase names and rows for phosphosites. I realized that the Kinase Library tool missed characterizing around 30 phosphosites because their phosphorylation sequences were asymmetric (the phosphosites were located either at the beginning or the end of a protein sequence) or because the sequence was located on a non-canonical protein isoform. So, I had to go back and manually find these phosphosites’ phosphorylation sequences, retrieve data from the Kinase Library, and add them to my master CSV file. This manual addition of data took a significant amount of time.
The CSV file served as the main subject of my downstream analysis! For one, I calculated the average percentile value of each kinase for all 800+ phosphosites. This analysis led to a dead end–no clear enrichment for a specific kinase. Next, I tried analyzing kinases with scores only in the top 90th percentile. This path led to some interesting results. I found kinases highly enriched across the phosphosite dataset (such as RSK4 and CRIK), as well as kinases key in regulating specific phosphosites (for example, CDK7 for phosphosite located at residue 224 in the protein SNW1). At the end of this analysis, I had produced 2 key text files in tree-like format for clear and structured viewing.
The next task was to figure out–how can we best visualize these results? After talking to my mentor, I learned about kinase tree visualization tools, like KinMap and Kinome Render. I spent a few days reading their corresponding papers and tool documentation, as well as reading studies that have used these tools. I was able to generate a preliminary interactive visualization of highly enriched kinases from my dataset using KinMap. The kinase tree is presented in Figure 1. Essentially, this tree groups kinases by classes, where the size of the circle representing a kinase is proportional to its abundance in my dataset. Using this tool was not without some hitches–the tool did not recognize 100 of the kinase names I inputted, so I had to manually input alternative gene names recognizable by the tool. My KinMap plot shows significant enrichment of the kinase classes AGC and CAMK. I also noticed that some of these kinases correlate with kinases enriched in metabolic diseases (these are the kinases highlighted in green, while the ones in my dataset are red).
Next week, I will try to generate a kinase tree with even more stringent filtration metrics (only including kinases present in the top 95th percentile across the phosphosites). Also, I am currently working to put together every step of my protein translocation analysis algorithm into a single pipeline–I’m first tracing out the key steps on paper. My hope is to bring together all of my analysis steps into a single customizable “package” that can be useful to pure-biology researchers. More to come on this in the following weeks.
Until next week!
Figure 1: KinMap Static Plot
Week 6: Opening Doors of Possibilities!
Hi everyone! Welcome back to my project. This week, I worked on debugging and rerunning the ADAR1 RNA-editing-site detection pipeline. I also made significant progress in the protein translocation project through bioinformatics approaches and algorithm brainstorming. Let’s dive right in!
While last week I took a hiatus from the ADAR1 project to focus on the protein translocation project, this week I made another attempt to troubleshoot and bolster the RNA-editing-site detection method. Upon close rereading of the original LoDEI paper, I learned about additional methods of data preprocessing that I should apply to my raw FASTQ RNA-Seq data: namely, FASTQC trimming and manual removal of variants that are likely to be single nucleotide polymorphisms, or SNPs, by cross-comparing with established SNP databases like gnomAD and dbSNP. This removal of SNPs helps differentiate true editing events from commonly occurring variations to the genome that exist naturally in the human population. I spent half of the week working on these changes and re-running the detection pipeline, which involved writing, running, and debugging code. However, I obtained the same results that I obtained the first time I ran the pipeline (likely because the RNA-Seq company had already implemented quality control practices on the mouse fat data), showing that data preprocessing was not the point of error in the pipeline. Nevertheless, I still finished up my project workflow by grouping RNA-editing sites into clusters, or RNA-editing-site islands, using K-means clustering and DBSCAN, and presented these findings. Alas, this marks the end of my current efforts on the ADAR1 project’s end.
Nevertheless, doors of possibilities opened up for understanding protein shuttling! Last week, I generated a KinMap plot for the protein kinases that have greater than 90th percentile abundance in my phosphoproteome dataset. This week, I made the conditions more stringent, changing the threshold to the 95th percentile, and built a new plot (See Figure 1). Comparing this plot to the one from last week, I saw the same relative distribution of kinases across the different kinase classes (no significant change in the clustering of kinases on different branches). This suggests that increasing the threshold does not offer greater delineation or meaningful filtration of relevant kinases. So, I needed to figure out a new way to narrow down and classify functionally-significant phosphosites. I describe this method below!
After meeting with my external advisor on Wednesday, I learned about a very exciting direction we can explore: the clinical relevance of protein phosphorylation sequences. Single nucleotide polymorphisms (SNPs) are, as the name suggests, are genomic variations at single base positions in the DNA sequence. Some SNPs do not have any effect on the protein produced due to the redundancy in the genetic code (the same amino acid is encoded), while others can have significant effects on protein function. My mentor suggested that I use public SNP databases (like gnomAD) to search my filtered phosphosite dataset for SNPs that occur within these phosphorylation sequences, which themselves are located near nuclear localization sequences (NLSs). If such SNPs result in a modified amino acid, they can be of great clinical significance by affecting a protein's sub-cellular localization, and this function! This week was mostly about better understanding this area – reading review papers, exploring the gnomAD database, and investigating APIs I can use to collect SNP data for 800+ genes in an automated fashion using the Python Selenium package. I was able to accomplish the first step – collecting SNP data in the form of 800+ CSV files. I also explored API tools from the Ensembl genome database, and I planned out the remainder of the workflow for identifying variants of clinical significance. I hope to begin implementing it and interpreting results in the following weeks.
See you next week!
Figure 1
Week 7: Investigating Clinical Significance
Hello! Welcome back to my Senior Project. This week, I made significant progress in an exciting, clinically relevant area: identifying and characterizing SNPs that are located within phosphorylation sequences!
First, let’s recap what my protein translocation project has been addressing so far so that you can see where it leads. I began with a Nature phosphoproteome dataset, which aggregates thousands of mass spectrometry experiments to comprehensively characterize phosphosites across the proteome. Using the tool cNLS Mapper, I predicted the monopartite and bipartite nuclear localization sequences (NLSs) of these proteins, which are sequences involved in regulating the protein’s transport from the cytoplasm to the nucleus. I also filtered this dataset to only include proteins that have been characterized by the literature to have dual cellular localization. From there, I identified phosphosites (phosphorylated protein residues) located within or near these NLSs, constructed their sequence motifs, and predicted the kinases that act on these sequences. Finally, I built kinase tree plots to visualize these results.
With these filtered datasets and results, the next question is: how do we connect this phosphosite information to actual disease phenotypes? As I explained last week, one key method of doing this is by identifying single nucleotide polymorphisms (SNPs) near these phosphosites—these single-base variations can impair the kinase’s ability to recognize and phosphorylate that NLS-implicated sequence, which can impair the ability of the protein to move between the cytoplasm and nucleus to carry out its normal biological function.
This week, I spent a significant amount of time writing code to retrieve and map SNP data from the gnomAD database to my phosphorylation sequence dataset. But before I could even do this, a very careful consideration needed to be made: my phosphorylation sequence dataset is in terms of protein residue coordinates, not genomic coordinates, whereas SNP data is always in genomic coordinates! Therefore, I had to carry out a few extra steps to convert my protein coordinates into genomic coordinates—I first had to retrieve the protein Ensembl IDs and use the Ensembl API to retrieve the corresponding genomic coordinates of my protein phosphorylation sequences. This process involved a lot of trial and error—bugs in the API calls, unavailable data for a few proteins, and validating the results on a few test cases. After that, I used the Selenium Python package to automate the retrieval of SNP data for 800+ genes (the genes that my phosphosites are located within), leaving me with 800+ CSV files of characterized SNPs and their protein consequences. Then, I scanned these CSVs for SNPs that occur within phosphorylation sequences in my phosphosite dataset, and returned the results as a new CSV file with extensive information on each SNP, including its protein consequence, rsID, clinical classification (benign, pathogenic, unknown, etc.), and CADD (predicted deleteriousness level). My resulting file had around 60,000 SNPs!
After meeting with my mentor, we decided on filters we would use to narrow down this list further—we are looking for a small number of SNPs that actually have clinical value and that we should test further in the lab. So, I removed SNPs that did not change the resulting amino acid (silent mutations) as well as SNPs that cause early termination of the protein (nonsense mutations). I also removed frameshift mutations, as we are specifically looking for missense mutations that might be clinically significant. At the end, I produced a list of around 25,000 SNPs. Next week, I will work with my mentor to understand other methods of filtration we can use. We will narrow down this list to focus mostly on pathogenic or likely pathogenic SNPs, cross-check them with known data sources, and maybe develop a plan to test them experimentally.
Thanks for reading!
Week 8: Revising and Re-running the Entire Pipeline
Hi everyone! Thanks for tuning in to my Senior Project. This week, I brought together my analysis steps into a coherent pipeline and re-ran my phosphoproteome data analysis code with revised parameters.
Before diving into code this week, I wanted to visualize my algorithmic approach in the form of a flowchart. This flowchart summarizes my methodology with relevant intermediate statistics (see Figure 1). Building this chart was very useful – it helped me keep the analysis steps clear in my head and also communicate my research work, as the goal of my project is not only to produce results for the lab’s specific research area but to create a generalizable method that researchers can customize to their needs (when it comes to analyzing phosphosproteome data for clinical relevance). I will include this flowchart (or a more sophisticated version) in my final research poster.
After meeting with my mentor with the flowchart in hand, he made some comments on how I can revise the workflow. While in the past weeks, I implemented this workflow as a proof-of-concept with preliminary parameters, he suggested that I make certain modifications: for example, including tyrosine phosphosite residues, lowering the threshold for best localization probability of phosphosites, and decreasing the cut-off score for NLSs to 4.0. So, this week, after modularizing and organizing my Python functions, as well as making them more interpretable with commenting (I often overlook this important aspect of programming…), I re-ran the workflow from the beginning: from filtering the Nature dataset to running cNLS mapper to parsing HTML results to finding SNPs. As my dataset now has 20,000+ phosphosites (located within 4,500+ genes), it is taking a long time to retrieve genomic and SNPs data (the code is still running). I tried parallelizing the code – sending multiple web requests at once using threads – but it is still slow, and I am experiencing some challenges with this, as the server keeps crashing. Hopefully, I will be able to share the results early next week.
Furthermore, I decided to break up the SNP analysis into three branches, which is something I did not do last time. This time, I am checking for SNPs in three categories, which will help us rank them by importance: SNPs that correspond with phosphosite locations, SNPs located in phosphosite flanking sequences, and SNPs located in NLSs. By the end of the week, I obtained some results (SNPs in phosphosite flanking sequences), which I shared with my mentor. Specifically, I found 7 SNPs (missense and 1 deletion) located in these regions that are likely pathogenic (classified using established bioinformatic tools). The genes these phosphosites are located within are H1-4, SASH1, PNPO, YARS1, ATXN7, and SPAST. This narrowed-down list of SNPs is important because computationally, I have demonstrated that they are clinically relevant. As the next step, we can create a plan to test them in vitro to confirm our hypotheses. As I finish running the aforementioned revised workflow (with less stringent parameters) and obtain SNPs results, I will similarly interpret the results biologically.
Thanks!
Figure 1
Week 9: Coding & Working on the Poster
Hi everyone! Welcome back to my Senior Project.
First off this week, I added some key information to the methodological flowchart that I made last week. Namely, I added intermediate statistics (number of proteins/phosphosites at each step) so that readers can clearly see how the filtration method is progressing. I also added relevant statistics for the two SNPs finding methods–SNPs that correspond with phosphosites, and SNPs located in phosphosites’ flanking sequences. Overall, this work helped make the flowchart more informative and useful.
This week also had some challenges with regard to the coding aspect. One of the things I have learned over the course of this project is how much iteration is necessary in programming and research in general. After obtaining SNPs results, the next key step was to compare them to positive control data. The positive control we are using is the gene HDAC4, known to be involved in nuclear import and to contain SNPs near NLS regions. For one, this week, I spent a lot of time researching HDAC4 – a key paper I read and annotated was “Missense substitutions at a conserved 14-3-3 binding site in HDAC4 cause a novel intellectual disability syndrome” (Wakeling et al.). I also looked at HDAC4 mass spectrometry data, which gave me a sense of the functional scores/importance of different amino acid sites in the protein. My mentor also shared a list of SNPs in HDAC4 from the ClinVar database that have known pathogenic consequences.
Comparing this provided list to the SNPs produced from my computational analysis, I found none of them to be present in my results! This signaled that something was wrong, either in my code or in the databases I was using. After a lot of troubleshooting and investigation, I found a few things–the API I was using to extract and download SNPs data was connected to the gnomAD database. This database provides external links to SNPs information, but unfortunately is not a comprehensive database, meaning that it did not contain a full list of SNPs for every gene. And, turns out, it was missing the pathogenic SNPs of HDAC4 that my mentor shared from the ClinVar database. Furthermore, going back to the code I wrote, I realized that not all of the SNPs data for the 4,000+ proteins had been downloaded as I intended due to the overwhelming number of requests to the server.
These observations made me go back to my workflow to see how I could improve it. For one, I decided to re-download SNPs data from the ClinVar database (instead of gnomAD), which has a comprehensive list of SNPs and their functional/clinical consequences (it is where my mentor found the positive control data for HDAC4). This code took a few days to run, as I increased the timeout between one gene and the next to ensure all the SNPs data was downloaded. Next, I decided to simplify the workflow I was using to find SNPs located in phosphorylation sequences, seeing as my current approach was missing out on some key known SNPs. Before, I was converting the phosphosite data (which is in terms of amino acid coordinates) into genomic coordinates, and then finding SNPs whose genomic coordinates fall within these regions. However, I realized that I could directly find the SNPs protein consequence data, which includes the amino acid residue they impact, and thus compare the phosphosite data and SNPs data directly without going through the laborious process of converting into genomic coordinates! This made the workflow a lot simpler and more intuitive.
In addition to running this code, I worked on interpreting my results biologically. What this means is, after getting the SNPs results, I investigated the functions of the genes they are located within to identify potentially interesting SNPs for further investigation. One interesting result I found was a SNP (rs672601307) located in the gene USP8, which is involved in membrane trafficking and cargo sorting–as my entire investigation centers around protein import and export, this gene could play a pivotal role. Additionally, I worked on planning my poster and assembling the necessary materials for it–the literature review I did at the beginning, the flowchart, and some of the results I have mentioned here. Currently, I am making a Google Docs outline of my poster, which I will later convert into the poster format.
Thanks for reading!
Week 10: The Final Post!
Hi everyone! Welcome to my last Senior Project blog post. I can’t believe these 10 weeks have come to an end–time really does fly when you’re having fun! I hope this post will be a memorable end to my project.
This week, I focused on creating my final product (a research poster) and presentation. After obtaining so many results and insights over the course of these 10 weeks, the task became to figure out how best to convey them. The research process was super nonlinear, iterative, and at times confusing (as is clear from my past blog posts…), and I needed to distill it to a clear, comprehensible format for my poster! Organizing and communicating my findings was one of the most challenging parts of the entire project.
I spent a lot of time making and revising my poster. I revised my methodology flowchart from last week to highlight key statistics and include some new analysis I did this week (filtering of pathogenic SNPs, inclusion of data on the number of unique phosphorylated NLSs, and researching SNPs’ clinical significance). I wrote the introduction for my poster by synthesizing the literature reviews I did at the beginning of the Senior Project. I created and found visualizations to help communicate the biological phenomena under investigation, both for my poster and my presentation. For example, I found figures on HDAC4 (the project’s positive control gene) that show its NLSs and key phosphosites. I put some images of my algorithm's step-by-step outputs alongside the methodological workflow so that readers can easily follow along. Additionally, I included visualizations of the kinase trees from KinMap.
I also did some further literature review this week to put my project results in a biological context. An interesting result that I identified from my SNPs discovery pipeline was 2 high functional score SNPs in the gene USP8, which is responsible for deubiquitinating proteins, protecting cell-surface receptors like EGFR and CXCR4 from degradation, as well as regulating nucleocytoplasmic shuttling. I focused on this result and the validation of the positive control (HDAC4 gene) in my discussion section. I also identified areas for further improvement to my methodology, such as testing different NLS prediction algorithms to increase the confidence of the results.
My poster went through multiple iterations as I regularly met with my external advisor and implemented his feedback. There were some important project results that I could not directly include, as they would not fit in the space–for example, the database of SNPs I built for further bioinformatic analysis. As such, I used QR codes to link to these external files.
After finishing my poster, I worked on my presentation, focusing on the research question, methodology, results, and most importantly, my learnings from the research process. I rehearsed my presentation multiple times to get it within the 8-minute limit. It was definitely challenging to condense my work into a poster and an 8-minute presentation!
Thank you all for following along with my project. I’m excited to present on May 17th!