(c) 2019 Justin Bois. With the exception of pasted graphics, where the source is noted, this work is licensed under a Creative Commons Attribution License CC-BY 4.0. All code contained herein is licensed under an MIT license.
This lesson was generated from an Jupyter notebook. You can download the notebook here.
a) The trick here is to do what we did in Lesson 7, except use
[::-1] indexing instead of the
def complement_base(base): """Returns the Watson-Crick complement of a base.""" if base == 'A' or base == 'a': return 'T' elif base == 'T' or base == 't': return 'A' elif base == 'G' or base == 'g': return 'C' else: return 'G' def reverse_complement(seq): """Compute reverse complement of a sequence.""" # Initialize reverse complement rev_seq = '' # Loop through and populate list with reverse complement for base in seq: rev_seq += complement_base(base) return rev_seq[::-1]
And we'll do a quick test with the same sequence as in lesson 7.
b) We can eliminate the
for loop by using the
replace() method of strings.
def reverse_complement(seq): """Compute reverse complement of a sequence.""" # Initialize rev_seq to a lowercase seq rev_seq = seq.lower() # Substitute bases rev_seq = rev_seq.replace('t', 'A') rev_seq = rev_seq.replace('a', 'T') rev_seq = rev_seq.replace('g', 'C') rev_seq = rev_seq.replace('c', 'G') return rev_seq[::-1]
And let's give it a test!
Note: We haven't learned about it yet, but some Googling would allow you to use the
maketrans() string methods.
maketrans() makes a translation table for characters in a string, and then the
translate() functions uses it to mutate the characters in the list.
def reverse_complement(seq): """Compute reverse complement of a sequence.""" return seq.translate(str.maketrans('ATGCatgc', 'TACGTACG'))[::-1] reverse_complement('GCAGTTGCA')
So, we were able to do it in one line!
Write a function that takes two sequences and returns the longest common substring. A substring is a contiguous portion of a string. For example:
TGCA T TAT
Not substrings of
AGCA # Skipped T CCATA # Added another C Hello, world. # Has nothing to do with the input sequence
There may be more than one longest common substring; you only need to return one of them.
The call signature of the function should be
Here are some return values you should get.
def longest_common_substring(s1, s2): """Return one of the longest common substrings""" # Make sure s1 is the shorter if len(s1) > len(s2): s1, s2 = s2, s1 # Start with the entire sequence and shorten substr_len = len(s1) while substr_len > 0: # Try all substrings for i in range(len(s1) - substr_len + 1): if s1[i:i+substr_len] in s2: return s1[i:i+substr_len] substr_len -= 1 # If we haven't returned, there is no common substring return ''
Let's try our function out.
In this problem, we will write a function that takes an RNA sequence and an RNA secondary structure and decides if the secondary structure is possible given the sequence. Remember, single stranded RNA can fold back on itself and form base pairs. An RNA secondary structure is simply the list of base pairs that are present. We will represent the base pairs in dot-parentheses notation. For example, a sequence/secondary structure pair would be
0123456789 GCAUCUAUGC (((....)))
For convenience of discussion, I have labeled the indices of the bases on the top row. In this case, base
G, pairs with base
1 pairs with base
8, and base
2 pairs with base
6 are unpaired. (This structure is aptly called a "hairpin.")
I hope the dot-parentheses notation is clear. An open parenthesis is paired with the parenthesis that closes it. Dots are unpaired.
So, the goal of our function is to check all base pairs present in a secondary structure and see if they are with
A-U, or (optionally)
a) Write a function to make sure that the number of closed parentheses is equal to the number of open parentheses, a requirement for a valid secondary structure. It should return
True if the parentheses are valid and
b) Write a function that converts the dot-parens notation to a tuple of 2-tuples representing the base pairs. We'll call this function
dotparen_to_bp(). An example input/output of this function would be:
dotparen_to_bp('(((....)))') ((0, 9), (1, 8), (2, 7))
Hint: You should look at methods that are available for lists. You might find the
pop() methods useful.
c) Because of sterics, the minimal length of a hairpin loop is three bases. A hairpin loop is a series of unpaired bases that are closed by a base pair. For example, the secondary structure
(.(....).) has a single hairpin loop of length 4. So, the structure
((((..)))) is not valid because it has a hairpin loop of only two bases.
Write a function that verifies that a list of base pairs (as outputted by
dotparen_to_bp()) satisfies the minimal hairpin length requirement.
d) Now write your validator function. The function definition should look like this:
def rna_ss_validator(seq, sec_struc, wobble=True):
It should return
True if the sequence is commensurate with a valid secondary structure and
False otherwise. The
wobble keyword argument is
True if we allow wobble pairs (
G paired with
U). Here are some expected results:
rna_ss_validator('GCAUCUAUGC', '(((....)))') rna_ss_validator('GCAUCUAUGU', '(((....)))') rna_ss_validator('GCAUCUAUGU', '(.(....).)')
rna_ss_validator('GCAUCUACGC', '(((....)))') rna_ss_validator('GCAUCUAUGU', '(((....)))', wobble=False) rna_ss_validator('GCAUCUAUGU', '(.(....)).') rna_ss_validator('GCCCUUGGCA', '(.((..))).')
a) This part of the validation is simple. We just need to make sure the number of open and closed parentheses are equal.
def parens_count(struc): """ Ensures there are equal number of open and closed parentheses in structure. """ return struc.count('(') == struc.count(')')
Let's give it a try.
b) As we scan a dot-parens structure from left to right, we can keep a list of the positions of open parentheses. Whenever we encounter a closed one, we have closed the last open one we added. So, we can just scan through the dot-parens string and pop out base pairs. If this procedure fails, we know that there was an error in the input structure (i.e., a closed parenthesis appeared without a corresponding open one before it).
def dot_parens_to_bp(struc): """ Convert a dot-parens structure to a list of base pairs. Return False if the structure is invalid. """ if not parens_count(struc): print('Error in input structure.') return False # Initialize list of open parens and list of base pairs open_parens =  bps =  # Scan through string for i, x in enumerate(struc): if x == '(': open_parens.append(i) elif x == ')': if len(open_parens) > 0: bps.append((open_parens.pop(), i)) else: print('Error in input structure.') return False # Return the result as a tuple return tuple(sorted(bps))
Let's try it on some legitimate sequences and on some bad ones.
# Good structure dot_parens_to_bp('..(((...)))(((((....)))).)..')
((2, 10), (3, 9), (4, 8), (11, 25), (12, 23), (13, 22), (14, 21), (15, 20))
# Good structure dot_parens_to_bp('(((..(((...)).))))')
((0, 17), (1, 16), (2, 15), (5, 14), (6, 12), (7, 11))
# Bad structure dot_parens_to_bp('((....)))(')
Error in input structure.
Error in input structure.
c) It is quite easy to detect short hairpins once we have a list of base pairs. We just need to make sure the difference in index of any pair of paired bases is not less than three.
def hairpin_check(bps): """Check to make sure no hairpins are too short.""" for bp in bps: if bp - bp < 4: print('A hairpin is too short.') return False # Everything checks out return True
d) Most everything is in place. We just need to check the sequence.
def rna_ss_validator(seq, sec_struc, wobble=True): """Validate and RNA structure""" # Convert structure to base pairs bps = dot_parens_to_bp(sec_struc) # If this failed, the structure was invalid if not bps: return False # Do the hairpin check if not hairpin_check(bps): return False # Possible base pairs if wobble: ok_bps = ('gc', 'cg', 'au', 'ua', 'gu', 'ug') else: ok_bps = ('gc', 'cg', 'au', 'ua') # Check complementarity for bp in bps: bp_str = (seq[bp] + seq[bp]).lower() if bp_str not in ok_bps: print('Invalid base pair.') return False # Everything passed return True
Let's test it on the test cases from the problem statement.
print('Should be True:') print(rna_ss_validator('GCAUCUAUGC', '(((....)))')) print(rna_ss_validator('GCAUCUAUGU', '(((....)))')) print(rna_ss_validator('GCAUCUAUGU', '(.(....).)')) print(rna_ss_validator('AUUGAUGCACGUGCAUCCCCAGCGGGUCCCGCGAGCUCACCCCCUUCCAAAAGCACCACGUGCCAGGCCUCGCCCCCGGAAGUAUACCUGUGAGCCAGA', '...(((((....)))))....((((...))))..((((((...(((((....((((...))))..(((...)))...))))).......))))))....')) print('\nShould be False:') print(rna_ss_validator('GCAUCUACGC', '(((....)))'), '\n') print(rna_ss_validator('GCAUCUAUGU', '(((....)))', wobble=False), '\n') print(rna_ss_validator('GCAUCUAUGU', '(.(....)).'), '\n') print(rna_ss_validator('GCCCUUGGCA', '(.((..))).'),'\n') print(rna_ss_validator('AUUGAUGCACGUGCAUCCCCAGCGGGUCCCGCGAGCCCACCCCCUUCCAAAAGCACCACGUGCCAGGCCUCGCCCCCGGAAGUAUACCUGUGAGCCAGA', '...(((((....)))))....((((...))))..((((((...(((((....((((...))))..(((...)))...))))).......))))))....'))
Should be True: True True True True Should be False: Invalid base pair. False Invalid base pair. False Invalid base pair. False A hairpin is too short. False Invalid base pair. False
In this problem, you will play around with the command line on your machine and get more familiar with it.
a) Let's play around with some options for the
ls command. First
cd into a directory that has some interesting files in it (like
~git/bootcamp/command_line_tutorial). Try the following if you are using
ls -F ls -G # Might not be as cool with Git Bash on Windows ls -l ls -lh ls -lS ls -FGLh
You should be able to infer what these different options do, but you should talk with the TAs as well.
Normally, files that begin with a dot (
.) are omitted when listing things. They are also generally omitted when you use your OS's GUI-based file handling system (like Finder on Macs). To see them, use
ls -a. So,
cd into your home directory (you remember how to do that, right?), and then do
b) The nuclear option to delete everything in a directory is
rm -rf. The
r means to delete recursively, and the
f means to "force" deletion. I was going to give you an exercise that uses the nuclear option, but I'm not going to do that. So, just forget I said anything. For this part of the problem, I want you to discuss with your neighbor when the nuclear option might be used, and what needs to be in place before exercising it.
c) Try doing this if you are using macOS or Linux:
cd-ing there and seeing what's in there. Do not delete anything!
.bashrc file allows you to configure your
bash shell how you like.
a) Create a
.bashrc file in your home directory. If you already have a
.bashrc file, open it up for editing using Jupyter's text editor.
b) It is often useful to
alias functions to other functions. For example, I am always worried I will accidentally delete things by accident. I therefore have the following line in my
alias rm="rm -i"
You should create aliases for commands like
ls based on the flags you like to always use. Do the same for
mv (I use the
-i flag with these). To figure out what flags are available, you can look at the
man pages. Asking Google will usually give you the information you need on flags.
If you like, you can use my
.bashrc file, available in
c) Depending on your operating system, your
~/.bashrc file may or may not be properly loaded upon opening a new bash shell. You may, e.g. for new macOS versions, need to explicitly source your
.bashrc file in your
~/.bash_profile file. Therefore, you should add the following to the bottom of your
if [ -f $HOME/.bashrc ]; then . $HOME/.bashrc fi
Again, this was mostly you messing around with the command line. The contents of my
.bashrc files are shown below.
# Give me a nice prompt that tells me my pwd export PS1="\[\e[1;32m\]\u\[\e[0m\]@\e[1;36m\]\h\[\e[0m\] [\w]\n% " # Keep me out of trouble! alias rm="rm -i" alias mv="mv -i" alias cp="cp -i" # customize list output alias ls="ls -FGh" export LSCOLORS="gxfxcxdxCxegedabagacad"
%load_ext watermark %watermark -v -p jupyterlab
CPython 3.7.3 IPython 7.1.1 jupyterlab 0.35.5