I take things too seriously sometimes. As far as I'm aware, I've always been like that. When I was a kid, I knew, deep down, that walkthroughs were bad. Maybe they weren't quite morally wrong, but they were certainly pushing it. You can probably ask my brother how I felt about walkthroughs. He was often on the receiving end of my scorn when he used a website to figure out a particularly difficult puzzle. "That's not the way you're SUPPOSED to do it!" Children can be so vicious. Little did I know I would eventually be publishing my own walkthrough.
Hopefully I have done my fair share of maturing. Games, I realize almost two decades later than everyone else, are intended to be fun by definition. If utilizing a guide makes the game more entertaining for you, then by all means, use it. Enjoy the little things in life the way you like them, not the way that you're "supposed" to, or the way that others tell you you should. It's better to learn these things late than never, right?
I've had a tremendous amount of fun putting this game together. That being said, I am still taking things very seriously, and have also viewed this as a learning experience. I've never tried to market anything before. To put it lightly, I failed terribly. Most of the code that I wrote for this game is rather different than the style of code I write for work. Anyone with any web development experience would have a few things to say about my email script.
"Whatever. This is my first time, and my next one will be better." Tell yourself that and go try something new. Learn something difficult.
Needless to say, the following link contains spoilers. It also contains a pretty awesome story. While I fully encourage the use of this document, only use it if you know that's what you want to do.
[This link has been temporarily removed]
Until next time,
Cheers!
Code
Monday, July 13, 2015
Wednesday, July 8, 2015
Excision Winner Announcement
Congratulations to Ankur Sundara for being the first person to complete Excision! For his dedication and problem-solving abilities, he's won the first place price of $50! Ankur, a New Hampshire resident, is a rising senior in high school.
Ankur had this to say about the game: "Excision was an awesome competition! It was a fun mix of knowledge of biology, computer security, and programming skills. My prior participation in hacking/programming competitions definitely helped me out a lot. Some of the challenges had me stumped for a while, but they all had pretty elegant solutions. I hope to see a bigger, better sequel to this competition in the future!"
I want to thank everyone who participated in the game or even gave it a glance. Also, for those of you who still want to play, there's an awesome second place price: a sense of adventure and the feeling of a job well done! Until next time!
Saturday, February 21, 2015
Base Camp
Welcome to Excision.
Before we set off on our adventure, I think it's only fair to give you some provisions for your journey. They come mostly in the form of the basic philosophy of Excision. You're going to be playing the role of an Agent who runs into some unique circumstances which are going to require some equally unique talents. But I'm sure you're up for the job.
The basic objective of Excision is to complete a number of stages. There will also be an accompanying story, but I'll explain up front that the basic template for each stage is (with some variations):
1. Read the story. It provides important context for the stage, and it’s entertaining!
2. Work through the puzzle at that stage.
3. Uncover the password necessary to unlock the next stage/clue.
4. Find the location of the next stage/clue.
5. Repeat.
The game itself has a few main themes, which closely echo the themes of this blog. Computers and programming are really exciting subjects, not to mention increasingly ubiquitous in our everyday lives. Unfortunately, they can seem intimidating to learn. I'm hoping these stages can walk you through a few of the basics. Of course, since this is a bioinformatics blog, there's an emphasis on computing in biology, but the skills necessary to complete the stages are useful in any field. There is almost always more than one way to get to the next stage. Think creatively. Think logically. Talk to other people; sometimes people think about things in different ways.
Of course, what's a game without a prize? The first person to complete the final stage of the puzzle will be the winner of $50! The $50 can be sent directly to you (cash or check) or donated in your name to your favorite charity.
Additionally, here's a brief reminder that, while anyone can play Excision, to be eligible for the prize you have to like the new Facebook page for the blog: www.facebook.com/WhatIsBioinformatics
I'm trying to separate my personal Facebook account from the blog, and this is a good way of doing it. More importantly for the game, the Facebook page is the platform of the referral/reward system.
The Facebook page is how you can earn personalized hints in order to complete Excision more quickly. Games are always more entertaining when played with friends, so invite yours to like the page and help you solve the challenges. Below are the specifics of how to earn and spend hints:
• If someone mentions your name in a comment or a post on the Facebook page, you earn 1 hint.
• If you share the Facebook page on your Facebook account, you earn 2 hints.
• Special “Bonus Hints” times will give Agents unique ways to earn extra hints.
• You can spend hints by sending a private message to the Facebook page.
• In the message, describe how you earned your hint (Comment/Share/Bonus Hints).
• If you earned your hint via a comment, include the name of the person who mentioned you.
• Include whatever question you want to have answered.
• While responses may not be immediate, one of our team members will get back to you within 24 hours.
• Questions asked can be retracted, for whatever reason, at any point before a clue has been given by sending the message “Please ignore the previous question”.
• Retracted questions will not cost a hint.
• Hints can be stockpiled indefinitely and do not need to be spent any given time, even during Bonus Hints.
A copy of the official rules for the Excision contest may be found at http://tinyurl.com/oq3aw4r. Continuing to the first level from this blog post implies that you have read and accepted the rules and that you know what it takes to be eligible for the prize.
Note that this is Base Camp. A few of the stages are located here, though they're encrypted and password-protected. When you’re not sure where to go next, always try returning here.
A final tip: a good Agent, like a good scientist, always cHecks their source, not jusT the source of the inforMation he/she is currently looking at but also aLl the way back to the beginning.
Have fun and enjoy Excision!
Before we set off on our adventure, I think it's only fair to give you some provisions for your journey. They come mostly in the form of the basic philosophy of Excision. You're going to be playing the role of an Agent who runs into some unique circumstances which are going to require some equally unique talents. But I'm sure you're up for the job.
The basic objective of Excision is to complete a number of stages. There will also be an accompanying story, but I'll explain up front that the basic template for each stage is (with some variations):
1. Read the story. It provides important context for the stage, and it’s entertaining!
2. Work through the puzzle at that stage.
3. Uncover the password necessary to unlock the next stage/clue.
4. Find the location of the next stage/clue.
5. Repeat.
The game itself has a few main themes, which closely echo the themes of this blog. Computers and programming are really exciting subjects, not to mention increasingly ubiquitous in our everyday lives. Unfortunately, they can seem intimidating to learn. I'm hoping these stages can walk you through a few of the basics. Of course, since this is a bioinformatics blog, there's an emphasis on computing in biology, but the skills necessary to complete the stages are useful in any field. There is almost always more than one way to get to the next stage. Think creatively. Think logically. Talk to other people; sometimes people think about things in different ways.
Of course, what's a game without a prize? The first person to complete the final stage of the puzzle will be the winner of $50! The $50 can be sent directly to you (cash or check) or donated in your name to your favorite charity.
Additionally, here's a brief reminder that, while anyone can play Excision, to be eligible for the prize you have to like the new Facebook page for the blog: www.facebook.com/WhatIsBioinformatics
I'm trying to separate my personal Facebook account from the blog, and this is a good way of doing it. More importantly for the game, the Facebook page is the platform of the referral/reward system.
The Facebook page is how you can earn personalized hints in order to complete Excision more quickly. Games are always more entertaining when played with friends, so invite yours to like the page and help you solve the challenges. Below are the specifics of how to earn and spend hints:
• If someone mentions your name in a comment or a post on the Facebook page, you earn 1 hint.
• If you share the Facebook page on your Facebook account, you earn 2 hints.
• Special “Bonus Hints” times will give Agents unique ways to earn extra hints.
• You can spend hints by sending a private message to the Facebook page.
• In the message, describe how you earned your hint (Comment/Share/Bonus Hints).
• If you earned your hint via a comment, include the name of the person who mentioned you.
• Include whatever question you want to have answered.
• While responses may not be immediate, one of our team members will get back to you within 24 hours.
• Questions asked can be retracted, for whatever reason, at any point before a clue has been given by sending the message “Please ignore the previous question”.
• Retracted questions will not cost a hint.
• Hints can be stockpiled indefinitely and do not need to be spent any given time, even during Bonus Hints.
A copy of the official rules for the Excision contest may be found at http://tinyurl.com/oq3aw4r. Continuing to the first level from this blog post implies that you have read and accepted the rules and that you know what it takes to be eligible for the prize.
Note that this is Base Camp. A few of the stages are located here, though they're encrypted and password-protected. When you’re not sure where to go next, always try returning here.
A final tip: a good Agent, like a good scientist, always cHecks their source, not jusT the source of the inforMation he/she is currently looking at but also aLl the way back to the beginning.
Have fun and enjoy Excision!
Sunday, February 15, 2015
The Protein Folding Answer
I need to issue a disclaimer before anything else. The script below is most definitely NOT the answer to protein folding. In fact, I make no claims at all about the efficacy of this code beyond it's rudimentary ability to fold the beta-hairpin shown here. This code is just a homework assignment and learning experience for me. I believe the code is decently well commented for anyone who's working with PyRosetta themselves. MAYBE this can be used at a starting point.
I'm also not going to go through the code section by section as I did for the sequence similarity code, but here's a brief discussion of what the script is supposed to accomplish:
1. First the peptide has to be generated. It's what we are going to fold in this program.
2. Since we know we're making a beta-hairpin, we use knowledge based design to set the basic starting point angles of the backbone.
3. Then we go through a million or so cycles of what amounts to three steps:
i. Randomly shift the peptide by some amount.
ii. Calculate the energy score of the new position compared to the old one. (Since we know we're looking for hydrogen bonds, we add an extra bonus for moves that potentially form h-bonds)
iii. Decide whether to accept the new position of the peptide, or to reject it and keep the old position. This decision is based on the Metropolis Criteria.
4. As the folding simulation progresses we periodically output the energy score to see how things are going.
5. At the end of the run, we save the new peptide so we can view it later.
In a way, the steps above are deceptively simply. There are only a few processes and they're fairly straightforward. On the other hand, the theory underlying these steps are very powerful. More importantly, the amount of refinement that can be build on top of these steps is enormous.
At the beginning of the simulation, the peptide is just a straight chain:
I'm also not going to go through the code section by section as I did for the sequence similarity code, but here's a brief discussion of what the script is supposed to accomplish:
1. First the peptide has to be generated. It's what we are going to fold in this program.
2. Since we know we're making a beta-hairpin, we use knowledge based design to set the basic starting point angles of the backbone.
3. Then we go through a million or so cycles of what amounts to three steps:
i. Randomly shift the peptide by some amount.
ii. Calculate the energy score of the new position compared to the old one. (Since we know we're looking for hydrogen bonds, we add an extra bonus for moves that potentially form h-bonds)
iii. Decide whether to accept the new position of the peptide, or to reject it and keep the old position. This decision is based on the Metropolis Criteria.
4. As the folding simulation progresses we periodically output the energy score to see how things are going.
5. At the end of the run, we save the new peptide so we can view it later.
In a way, the steps above are deceptively simply. There are only a few processes and they're fairly straightforward. On the other hand, the theory underlying these steps are very powerful. More importantly, the amount of refinement that can be build on top of these steps is enormous.
At the beginning of the simulation, the peptide is just a straight chain:
By the end, the peptide has folded and formed hydrogen bond. The cartoon representation of the peptide shows that it does appear to be a beta-hairpin.
Below is dumped the code used to make this transformation. If you have questions about any of the specifics, I'd be happy to try to answer them. For those of you looking for the actual answer to the protein folding problem sorry for the disappointment.
Maybe next time. Cheers!!
"""Import necessary libraries and initialize rosetta"""
import random
import math
from rosetta import *
from rosetta.PyMolLink import *
init()
"""Generate the beta-hairpin pose and score it"""
strand = "AAAAA" #alanine end residues
turn = "GG" #glycine turn residues
hairpin = strand+turn+strand #combine to make hairpin
pose = pose_from_sequence(hairpin,"fa_standard") #Create Rosetta pose
fa = get_fa_scorefxn() #Initialize a scoring function
fa(pose) #Score the hairpin
pose.dump_pdb("start1.pdb") #Generate an image of the starting structure
pose.set_psi(5,90) #Set starting angle
pose.set_phi(6,60) #Set starting angle
pose.set_psi(6,-120) #Set starting angle
pose.set_phi(7,-80) #Set starting angle
pose.set_psi(7,0) #Set starting angle
pose.set_phi(8,-120) #Set starting angle
"""Initialize variables"""
old_pose = Pose() #Initialize empty poses for later
low_pose = Pose() #Initialize empty poses for later
max_temp = 100 #Initial temperature for metropolis and annealing
current_temp = max_temp #Current temperature will drop through cycles
max_cycles = 1000000 #How many moves you want to make
low_energy = 10000 #Will store the value of lowest energy conformation
constraint_weight = 0.5 #Parameter for how much h-bond distances matter
ideal_dist = 0.29 #The proper distance for an h-bond between the alanines
"""Begin iterating through moves"""
for i in range(0,max_cycles): #Loop through the all the cycles
"""Update accept flag, and the old pose/energy/score"""
old_pose.assign(pose) #Set the old_pose
fa(old_pose) #Score the old pose
old_energy = fa.score(old_pose) #Define variable for the energy of old pose
old_score = old_energy #Alias old_energy, to use old_score for hbond constraining
"""Move given angle or angles by a chosen amount."""
current_residue = random.randint(1,pose.total_residue()) #Choose the residue to move from a uniform distribution
curr_phi = pose.phi(current_residue) #Set the current phi angle to be moved
curr_psi = pose.psi(current_residue) #Set the current psi angle to be moved
angle_shift = random.gauss(0,25*(float(i)/float(max_cycles))) #As the folding progresses, move the angle(s) by a smaller amount
if(random.random()<(float(i)/float(max_cycles))): #As the folding progresses, make the probability of a shear move increase
pose.set_phi(current_residue,curr_phi+angle_shift) #Set the new phi angle (shear move)
pose.set_psi(current_residue,curr_psi-angle_shift) #Set the new psi angle (shear move)
elif(random.random()<0 data-blogger-escaped--0.01="" data-blogger-escaped--="" data-blogger-escaped-.5="" data-blogger-escaped-0:="" data-blogger-escaped-0="" data-blogger-escaped-1000="" data-blogger-escaped-3="" data-blogger-escaped-a="" data-blogger-escaped-above="" data-blogger-escaped-accept:="" data-blogger-escaped-accept="" data-blogger-escaped-accepted="" data-blogger-escaped-accepting="" data-blogger-escaped-actions="" data-blogger-escaped-addition="" data-blogger-escaped-after="" data-blogger-escaped-alculate="" data-blogger-escaped-and="" data-blogger-escaped-angle="" data-blogger-escaped-angle_shift="" data-blogger-escaped-appropriate="" data-blogger-escaped-are="" data-blogger-escaped-as="" data-blogger-escaped-atom="" data-blogger-escaped-atoms="" data-blogger-escaped-ave="" data-blogger-escaped-bad="" data-blogger-escaped-between="" data-blogger-escaped-bonding="" data-blogger-escaped-change="" data-blogger-escaped-chosen="" data-blogger-escaped-code="" data-blogger-escaped-constraining="" data-blogger-escaped-constraint_weight="" data-blogger-escaped-core="" data-blogger-escaped-criteria="" data-blogger-escaped-cterm="pose.residue(12-2*k)" data-blogger-escaped-cterm_pos="Cterm.xyz(" data-blogger-escaped-curr_phi="" data-blogger-escaped-curr_psi="" data-blogger-escaped-current="" data-blogger-escaped-current_residue="" data-blogger-escaped-current_temp="max_temp*(math.exp(-1*max_cycles**-1))" data-blogger-escaped-delta_e="" data-blogger-escaped-dimentions="" data-blogger-escaped-dist="dist" data-blogger-escaped-distance="" data-blogger-escaped-each="" data-blogger-escaped-ecide="" data-blogger-escaped-edo="" data-blogger-escaped-efine="" data-blogger-escaped-else:="" data-blogger-escaped-end="" data-blogger-escaped-energy:="" data-blogger-escaped-energy="" data-blogger-escaped-epeat="" data-blogger-escaped-eset="" data-blogger-escaped-estore="" data-blogger-escaped-et="" data-blogger-escaped-f="" data-blogger-escaped-fa="" data-blogger-escaped-file="" data-blogger-escaped-final="" data-blogger-escaped-find="" data-blogger-escaped-first="" data-blogger-escaped-folding="" data-blogger-escaped-for="" data-blogger-escaped-hbond="" data-blogger-escaped-hbonding="" data-blogger-escaped-heck="" data-blogger-escaped-here="" data-blogger-escaped-i="" data-blogger-escaped-ideal_dist-old_dist="" data-blogger-escaped-if="" data-blogger-escaped-in="" data-blogger-escaped-inal="" data-blogger-escaped-ind="" data-blogger-escaped-is="" data-blogger-escaped-isn="" data-blogger-escaped-iterations="" data-blogger-escaped-j="" data-blogger-escaped-just="" data-blogger-escaped-k="" data-blogger-escaped-lias="" data-blogger-escaped-low_energy="" data-blogger-escaped-low_pose.assign="" data-blogger-escaped-low_pose.dump_pdb="" data-blogger-escaped-lower="" data-blogger-escaped-lowest="" data-blogger-escaped-math.fabs="" data-blogger-escaped-math.pow="" data-blogger-escaped-max_cycles="" data-blogger-escaped-metropolis="" data-blogger-escaped-move="" data-blogger-escaped-moves="" data-blogger-escaped-new="" data-blogger-escaped-new_dist="math.sqrt(dist)" data-blogger-escaped-new_energy="" data-blogger-escaped-new_score="new_score+math.fabs(ideal_dist-new_dist)*constraint_weight" data-blogger-escaped-not="" data-blogger-escaped-nterm="pose.residue(1+2*k)" data-blogger-escaped-nterm_pos="Nterm.xyz(" data-blogger-escaped-of="" data-blogger-escaped-old="" data-blogger-escaped-old_dist="math.sqrt(dist)" data-blogger-escaped-old_energy="" data-blogger-escaped-old_pose="" data-blogger-escaped-old_score="" data-blogger-escaped-one="" data-blogger-escaped-onte_carlo_2-blog2.pdb="" data-blogger-escaped-or="" data-blogger-escaped-other="" data-blogger-escaped-outcome.="" data-blogger-escaped-pair="" data-blogger-escaped-pdate="" data-blogger-escaped-pdb="" data-blogger-escaped-phi="" data-blogger-escaped-pose.="" data-blogger-escaped-pose.assign="" data-blogger-escaped-pose.set_phi="" data-blogger-escaped-pose.set_psi="" data-blogger-escaped-pose="" data-blogger-escaped-poses="" data-blogger-escaped-pre="" data-blogger-escaped-print="" data-blogger-escaped-probability="" data-blogger-escaped-process="" data-blogger-escaped-progress="" data-blogger-escaped-psi="" data-blogger-escaped-r="" data-blogger-escaped-range="" data-blogger-escaped-residue="" data-blogger-escaped-residues="" data-blogger-escaped-rint="" data-blogger-escaped-runs="" data-blogger-escaped-save="" data-blogger-escaped-score.="" data-blogger-escaped-score="" data-blogger-escaped-script="" data-blogger-escaped-second="" data-blogger-escaped-set="" data-blogger-escaped-shear="" data-blogger-escaped-so="" data-blogger-escaped-spacial="" data-blogger-escaped-statements="" data-blogger-escaped-str="" data-blogger-escaped-structure="" data-blogger-escaped-t="" data-blogger-escaped-take="" data-blogger-escaped-temperature="" data-blogger-escaped-term_pos="" data-blogger-escaped-than="" data-blogger-escaped-the="" data-blogger-escaped-there="" data-blogger-escaped-to="" data-blogger-escaped-until="" data-blogger-escaped-update="" data-blogger-escaped-use="" data-blogger-escaped-utoff="" data-blogger-escaped-value="" data-blogger-escaped-values="" data-blogger-escaped-variable="" data-blogger-escaped-very="" data-blogger-escaped-whether="" data-blogger-escaped-with="">
Friday, February 6, 2015
Thursday, February 5, 2015
The Protein Folding Problem
First, if you've never had the opportunity to visit xkcd, I highly suggest you do it. It's a wonderful website. My personal favorite is the What If? feature. Read one of them; you'll like it. The comic relevant to this post is:
If you're unfamiliar with the protein folding problem, let me briefly set the stage. Imagine a protein made of 100 amino acids linked together in a chain. The protein has to figure out where to put each of these amino acids relative to each other, in order to work properly. (As an aside, proteins which fail to fold properly have implications to many horrible diseases.) For simplicity, let's also image that each of the amino acids can take one of 2 positions. This implies that the protein can take 2100 positions. This may or may not seem like a lot of choices to you. Try putting 2100 into WolframAlpha. It works out to be 1267650600228229401496703205376 different positions, and only one of them is right.
To think about this number a different way, consider having 100 computers that will each check a trillion of these different possibilities a second. That's pretty good, right? Actually, it's not even close. It's going to take you approximately ~400,000,000 years to be able to look at all possibilities. It's absolutely impossible to fold a protein this way. It's not how nature does it, and it isn't how we attempt to do it either. Instead, the general hypothesis is that as the protein begins to fold into a favorable/proper state, it lowers it's energy, which basically crosses certain possibilities off the list without even having to check. The usual analogy is a folding funnel:
So, while all spaces on the funnel are technically possible, the protein is almost always just going to fall down the funnel. When scientist try to fold proteins, they use this funnel analogy as a method for only checking some of the possible conformations. But they check intelligently, mostly looking at the positions that have a high chance of being the natural shape of the protein.
The specifics of how to implement such a search on a computer can get as technical as you like. But a teaching tool called PyRosetta has been made that allows students like me (and you?) to get a pretty good idea of whats going on and what tools we have available to us. In my next post I'll show you a very simple script written for class that allows for the folding of a short peptide.
Until next time, cheers.
If you're unfamiliar with the protein folding problem, let me briefly set the stage. Imagine a protein made of 100 amino acids linked together in a chain. The protein has to figure out where to put each of these amino acids relative to each other, in order to work properly. (As an aside, proteins which fail to fold properly have implications to many horrible diseases.) For simplicity, let's also image that each of the amino acids can take one of 2 positions. This implies that the protein can take 2100 positions. This may or may not seem like a lot of choices to you. Try putting 2100 into WolframAlpha. It works out to be 1267650600228229401496703205376 different positions, and only one of them is right.
To think about this number a different way, consider having 100 computers that will each check a trillion of these different possibilities a second. That's pretty good, right? Actually, it's not even close. It's going to take you approximately ~400,000,000 years to be able to look at all possibilities. It's absolutely impossible to fold a protein this way. It's not how nature does it, and it isn't how we attempt to do it either. Instead, the general hypothesis is that as the protein begins to fold into a favorable/proper state, it lowers it's energy, which basically crosses certain possibilities off the list without even having to check. The usual analogy is a folding funnel:
So, while all spaces on the funnel are technically possible, the protein is almost always just going to fall down the funnel. When scientist try to fold proteins, they use this funnel analogy as a method for only checking some of the possible conformations. But they check intelligently, mostly looking at the positions that have a high chance of being the natural shape of the protein.
The specifics of how to implement such a search on a computer can get as technical as you like. But a teaching tool called PyRosetta has been made that allows students like me (and you?) to get a pretty good idea of whats going on and what tools we have available to us. In my next post I'll show you a very simple script written for class that allows for the folding of a short peptide.
Until next time, cheers.
Thursday, January 15, 2015
Simple Sequence Similarity
I've mostly refrained from posting any of the code that I've been writing in the past months because the code seems scraped together in a "just get it done" method that wouldn't be particularly interesting or useful to anyone but me. The programs are simple scripts only applicable to the single task at hand. But I've finally gained enough confidence to show off a few lines that I've been using recently.
How similar is the phrase "Get it finished" to "Finished get it"? Or in protein terminology, how homologous is GETITFINISHED to FINISHEDGETIT? Many similarity calculations assume sequence linearity, and therefore would not detect much likeness between these two phrases/sequences. The script below is a very simple way of detecting similarity in these sort of "mixed up" strings.
Even with the comments, this is a lot to digest at once, so let's go through it section by section.
The first section is pretty straightforward. It simply give Python access to the necessary tools that I'm going to need later on in the script.
The next section defines a couple of functions. "k_mers" takes the name of a data file (currently this must be a FASTA file, but it would be easy to generalize the function to be able to read other strings as well), the name of the file to which you would like to store your data, and the size of the substring/kmer you want to analyze. It will then produce a file telling you how many times each kmer appeared in each of the the sequences you want to compare. Notice that the "occurrences" function is defined exclusively inside of the "k_mers" function. It is the function that actually does the counting. If you know much of anything about Python, you probably want to point out that Python has a built in count method. The count method Python provides, however, does not count overlapping sub-strings. For example if you used the default count method on "AGTC" with a size of 2, it would return "AG" and "TC". "occurrences", on the other hand, will return "AG", "GT", and "TC"'.
Next, we need to know how many transcripts we are going to analyze. So, we open the file and read through it with a counter that adds one each time a new transcript is encountered. I'm not going to get into generators here, but Python needs to reopen a file to be able to read through it again, so I have to recall the same line.
Now we can make an 2-by-2 array that's eventually going to be filled with the number of kmers that are present in each transcript and begin cycling through them. For each transcript, we call the "occurrences" function we defined above. We send "occurrences" the current transcript we are working on as well as the size of the kmers we are looking for. It will return with a dictionary of the kmers and their corresponding counts.
The next few lines might make more sense if you've looked ahead to what sort of statistical method we are eventually going to use for comparison of the sequences. But consider the trivial case where we want to compare "AGTA" to "CAAA" using a kmer size of 1. The first transcript is going to produce "A"=2, "G"=1, and "T"=1 while the second will be "A"=4 and "C"=1. But what is needed for "an apples-to-apples" comparison is the first "A"=2, "G"=1, "T"=1, "C"=0 while the second will now be "A"=3, "G"=0, "T"=0, and "C"=1. These next lines are producing those necessary placeholders.
The last two lines of this function are pretty self explanatory. It would have been a bit easier to keep the matrix in memory and continue on to the correlation calculations, but here is a natural breaking point which allows the user to do other things with the data if they wish. "k_mers_corr" is a fairly simple function, which takes the names of the file produced in the previous function, as well as another output file to store this new data. The new data will be a symmetric matrix where each transcript is compared to every other transcript in the file, and their similarity are saved. After getting things going by bringing in the data and creating a square matrix of the proper size, we start iterating through the rows of our data.
For each row we go through all of the rows (including the row we are currently on) and compare their similarity via the Pearson correlation. Once we have that number, we store it in the proper cell in the matrix we initialized earlier. Once we have filled all of the cells, we save our data.
And we're finished! The functions I've walked through here are actually stripped down versions of what I have been using recently (don't want to give away ALL my secrets!), so it would appear that I've made my first real Bioinformatics post. I hope I made it clear but if any clarifications are needed, feel free to ask a question in the comments below. Cheers!
How similar is the phrase "Get it finished" to "Finished get it"? Or in protein terminology, how homologous is GETITFINISHED to FINISHEDGETIT? Many similarity calculations assume sequence linearity, and therefore would not detect much likeness between these two phrases/sequences. The script below is a very simple way of detecting similarity in these sort of "mixed up" strings.
import numpy as np
from Bio import SeqIO
from itertools import product
from scipy.stats.stats import pearsonr
def k_mers(data, outfile,k):
def occurrences(s,k1):
dic = {} #Initialize dictionary
c = 0 #Initialize counter
while c+k1 <= len(s): #Iterate through the sequence until the end
k_mer = s[c:c+k1] #Kmer defined as a slice of sequence s and size k
if k_mer not in dic: #If the kmer hasn't been seen before
dic.update({k_mer:1}) #Add the kmer as a key to the dictionary.
elif k_mer in dic: #If the kmer has been seen before
dic[k_mer] += 1 #Increment the kmer's value by one.
c += 1 #Move over one letter in the sequence
return dic #Return the full dictionary
transcripts = SeqIO.parse(open(data, "r"),'fasta') #Open the fasta file
c = 0 #Initialize a counter to determine size of storage array
for i in transcripts: #For each transcript
c+=1 #Add 1 to the counter
transcripts = SeqIO.parse(open(data, "r"),'fasta') #Reopen the fasta file
c_mat = np.zeros([c,4**k]) #Initialize array to store the counts
for i, rna in enumerate(transcripts): #For each transcript
counts = occurrences(str(rna.seq),k) #Get the kmer counts for sequence
counts_ar = np.zeros([4**k]) #Initialize array to store all kmers
for j, sub in enumerate(product("AGTC", repeat=k)): #Generate each kmer
k_mer = "".join(sub) #Turn kmer list into a string
if k_mer in counts: #If the kmer is in the transcript
counts_ar[j] = counts[k_mer] #Add the proper kmer value
c_mat[i] = counts_ar #Add new row of counts to count matrix
np.savetxt(open(outfile, "w"),c_mat, delimiter=",") #Save the matrix
def k_mers_corr(c_mat, outfile):
c_mat = np.genfromtxt(open(c_mat,"r"), delimiter=",") #Open the count matrix
corr_mat = np.zeros([c_mat.shape[0],c_mat.shape[0]]) #Initialize new matrix
for i, row1 in enumerate(c_mat): #For each row in count matrix
for j, row2 in enumerate(c_mat): #For each row in count matrix
dist = pearsonr(row1,row2)[0] #Get the corr coeff between rows
corr_mat[j][i] = dist #Add coeff to matrix
np.savetxt(open(outfile, "w"),corr_mat, delimiter=",") #Save the matrix
p = "A:/Library2/Data/" #My directory containing all my working files
i = p+"gencode_lncRNA_overlap_m.txt" #My transcript file
o1 = p+"blogger1.csv" #Output file for the count data
o2 = p+"blogger2.csv" #Output file for the correlation data
k_mers(i,o1,4) #Call to first function
k_mers_corr(o1,o2) #Call to second function
Even with the comments, this is a lot to digest at once, so let's go through it section by section.
The first section is pretty straightforward. It simply give Python access to the necessary tools that I'm going to need later on in the script.
import numpy as np from Bio import SeqIO from itertools import product from scipy.stats.stats import pearsonr
The next section defines a couple of functions. "k_mers" takes the name of a data file (currently this must be a FASTA file, but it would be easy to generalize the function to be able to read other strings as well), the name of the file to which you would like to store your data, and the size of the substring/kmer you want to analyze. It will then produce a file telling you how many times each kmer appeared in each of the the sequences you want to compare. Notice that the "occurrences" function is defined exclusively inside of the "k_mers" function. It is the function that actually does the counting. If you know much of anything about Python, you probably want to point out that Python has a built in count method. The count method Python provides, however, does not count overlapping sub-strings. For example if you used the default count method on "AGTC" with a size of 2, it would return "AG" and "TC". "occurrences", on the other hand, will return "AG", "GT", and "TC"'.
def k_mers(data, outfile,k):
def occurrences(s,k1):
dic = {} #Initialize dictionary
c = 0 #Initialize counter
while c+k1 <= len(s): #Iterate through the sequence until the end
k_mer = s[c:c+k1] #Kmer defined as a slice of sequence s and size k
if k_mer not in dic: #If the kmer hasn't been seen before
dic.update({k_mer:1}) #Add the kmer as a key to the dictionary.
elif k_mer in dic: #If the kmer has been seen before
dic[k_mer] += 1 #Increment the kmer's value by one.
c += 1 #Move over one letter in the sequence
return dic #Return the full dictionary
Next, we need to know how many transcripts we are going to analyze. So, we open the file and read through it with a counter that adds one each time a new transcript is encountered. I'm not going to get into generators here, but Python needs to reopen a file to be able to read through it again, so I have to recall the same line.
transcripts = SeqIO.parse(open(data, "r"),'fasta') #Open the fasta file
c = 0 #Initialize a counter to determine size of storage array
for i in transcripts: #For each transcript
c+=1 #Add 1 to the counter
transcripts = SeqIO.parse(open(data, "r"),'fasta') #Reopen the fasta file
Now we can make an 2-by-2 array that's eventually going to be filled with the number of kmers that are present in each transcript and begin cycling through them. For each transcript, we call the "occurrences" function we defined above. We send "occurrences" the current transcript we are working on as well as the size of the kmers we are looking for. It will return with a dictionary of the kmers and their corresponding counts.
c_mat = np.zeros([c,4**k]) #Initialize array to store the counts
for i, rna in enumerate(transcripts): #For each transcript
counts = occurrences(str(rna.seq),k) #Get the kmer counts for sequence
The next few lines might make more sense if you've looked ahead to what sort of statistical method we are eventually going to use for comparison of the sequences. But consider the trivial case where we want to compare "AGTA" to "CAAA" using a kmer size of 1. The first transcript is going to produce "A"=2, "G"=1, and "T"=1 while the second will be "A"=4 and "C"=1. But what is needed for "an apples-to-apples" comparison is the first "A"=2, "G"=1, "T"=1, "C"=0 while the second will now be "A"=3, "G"=0, "T"=0, and "C"=1. These next lines are producing those necessary placeholders.
for j, sub in enumerate(product("AGTC", repeat=k)): #Generate each kmer
k_mer = "".join(sub) #Turn kmer list into a string
if k_mer in counts: #If the kmer is in the transcript
counts_ar[j] = counts[k_mer] #Add the proper kmer value
The last two lines of this function are pretty self explanatory. It would have been a bit easier to keep the matrix in memory and continue on to the correlation calculations, but here is a natural breaking point which allows the user to do other things with the data if they wish. "k_mers_corr" is a fairly simple function, which takes the names of the file produced in the previous function, as well as another output file to store this new data. The new data will be a symmetric matrix where each transcript is compared to every other transcript in the file, and their similarity are saved. After getting things going by bringing in the data and creating a square matrix of the proper size, we start iterating through the rows of our data.
def k_mers_corr(c_mat, outfile):
c_mat = np.genfromtxt(open(c_mat,"r"), delimiter=",") #Open the count matrix
corr_mat = np.zeros([c_mat.shape[0],c_mat.shape[0]]) #Initialize new matrix
for i, row1 in enumerate(c_mat): #For each row in count matrix
For each row we go through all of the rows (including the row we are currently on) and compare their similarity via the Pearson correlation. Once we have that number, we store it in the proper cell in the matrix we initialized earlier. Once we have filled all of the cells, we save our data.
for j, row2 in enumerate(c_mat): #For each row in count matrix
dist = pearsonr(row1,row2)[0] #Get the corr coeff between rows
corr_mat[j][i] = dist #Add coeff to matrix
np.savetxt(open(outfile, "w"),corr_mat, delimiter=",") #Save the matrix
And we're finished! The functions I've walked through here are actually stripped down versions of what I have been using recently (don't want to give away ALL my secrets!), so it would appear that I've made my first real Bioinformatics post. I hope I made it clear but if any clarifications are needed, feel free to ask a question in the comments below. Cheers!
Subscribe to:
Posts (Atom)