Applied Bioinformatics

of Nucleic Acid Sequences

David A. Hendrix

October 3, 2019

Applied Bioinformatics by David A. Hendrix is licensed under a Creative Commons Attribution 4.0 International License, except where otherwise noted.

This book was produced with Pressbooks.

1

For the students and learners of the world.

Contents

- Dedication

- For the students and learners of the world.

- Foreword

- Foreword

- Chapter 1: Introduction to Biological Sequences, Biopython, and GNU/Linux

- 1.1 Nucleic Acid Bioinformatics

- 1.2 Sequences, Strings, and the Genetic Code

- 1.3 Sequences File Formats

- 1.4 Lab 1: Introduction to GNU/Linux and Fasta files

- 1.5 Biological Sequence Database

- 1.6 Lab 2: FASTQ and Quality Scores

- Chapter 2: Sequence Motifs

- 2.1 Introduction to Motifs

- 2.2 String Matching

- 2.3 Consensus Sequences

- 2.4 Motif Finding

- 2.5 Promoters

- 2.6 De novo Motif Finding

- 2.7 Lab 3: Introduction to Motifs

- Chapter 3: Sequence Alignments

- 3.1 Alignment Algorithms and Dynamic Programming

- 3.2 Alignment Software

- 3.3 Alignment Statistics

- 3.4 Short Read Mapping

- 3.5 Lab 4: Using BLAST on the command line

- Chapter 4: Multiple Sequence Alignments, Molecular Evolution, and Phylogenetics

- 4.1 Multiple Sequence Alignment

- 4.2 Phylogenetic Trees

- 4.3 Models of mutations

- 4.4 Lab 5: Phylogenetics

- Chapter 5: Genomics

- 5.1 The Three Fundamental “Gotchas” of Genomics

- 5.2 Genomic Data and File Formats

- 5.3 Genome Browsers

- 5.4 Lab 6: Genome Annotation Data

- Chapter 6: Transcriptomics

- 6.1 High-throughout Sequencing (HTS)

- 6.2 RNA Deep Sequencing

- 6.3 Small RNA sequencing

- 6.4 Long RNA sequencing

- 6.5 Single-Cell Transcriptomics

- 6.6 Transcription Initiation

- 6.7 Transcription

- 6.8 Elongation

- 6.9 Lab 7: RNA-seq

- Chapter 7: Noncoding RNAs

- 7.1 Small Noncoding RNAs (srcRNAs)

- 7.2. Long Noncoding RNAs

- 7.3 RNA Structure Prediction

- 7.4 Destabilizing energies

- 7.5 Lab 8: RNA Structure

- Chapter 8: Proteins

- 8.1 Protein Alignment

- 8.2 Functional Annotation of Proteins

- 8.3 Secondary Structure prediction

- 8.4 Gene Ontology

- 8.5 Lab 9: Proteins

- Chapter 9: Gene Regulation

- 9.1 Transcription Factors and ChIP-seq

- 9.2 MicroRNA regulation and Small RNA-seq

- 9.3 Regulatory Networks

- 9.4 Lab 10: ChIP-seq

- Appendix A: Mathematical Preliminaries

- Appendix B: Probability

- Bibliography

- Creative Commons License

- Recommended Citations

- Versioning

2

Foreword

Please note that this is the first edition of the text. Applied Bioinformatics is frequently updated, and the current edition may be downloaded from the Applied Bioinformatics site.

Chapter 1: Introduction to Biological Sequences, Biopython, and GNU/Linux

1.1 Nucleic Acid Bioinformatics

The tremendous growth of bioinformatics and computational biology in the late 20th and early 21st centuries has had an associated growth in software and algorithms for studying biological sequences. This book is designed to introduce students of the life sciences with little to no programming experience to the concepts and methodologies of bioinformatics, and contemporary software applications. The majority of this course will deal with the analysis of nucleic acid sequence data. Many of these applications have web interfaces, but where possible command line software, and strategies for dealing with sequence data with the GNU/Linux command line interface (the terminal), will be used. We won’t be content to simply run commands blindly, but rather we will also learn a great deal of the equations and theory behind these methods. In addition, these notes will provide an introduction to Biopython as a tool for bioinformatics and as a framework to connect all the concepts presented.

1.1.1 GNU/Linux and the command line

Bioinformatics is for the most part done best using the GNU/Linux command line interface, and therefore we will need a brief introduction. The GNU/Linux command line interface is well suited for working with the kinds of files commonly used in bioinformatics. After this introduction, we will continue to learn new GNU/Linux commands while we broaden our repertoire of software tools, Biopython commands, file formats, and databases.

The GNU/Linux command line interface is a text-based interface to files and commands. Such a terminal can be accessed from many terminal applications in a GNU/Linux distribution, or with the Terminal program (and others such as iTerminal) for Mac OS X, and other software available for Windows, such as Putty and MobaXterm. Advanced users of the terminal can accomplish difficult tasks with amazing efficiency, such as running complex commands on thousands of files in one fell-swoop while barely breaking a sweat. The combination of a strong knowledge of the terminal and shell scripting, along with knowledge of a scripting language like python or perl can be very powerful and worth learning for any student of bioinformatics. Scripting languages are named after the fact that they are used to write “scripts”, which are programs written to automate a sequence of commands and typically use an interpreter rather than a compiler.

Once you’ve opened your terminal application, you should be given a cursor in which to type text-based commands. For example, say you want to print a list of all files and directories in your current location (directory), you use the command . Directories are like folders for keeping files, but they can also be a location. In other words, you can be working in a directory, and any files you create will then be kept in that directory. A file, such as a text file for storing a sequence, is the unit in which data is stored. Although most people know what a file is from their experiences with computers, there is a lot that could be said about what a file actually is beyond the scope of these class notes. For our purposes right now, this intuitive concept of a file that most people have will suffice.

We can create a directory with the command . For example, let’s begin by creating two directories called “data1” and “data2”, followed by the command to see the contents of the current directory:

$ mkdir scripts
$ mkdir data
$ ls
data    scripts

Please note that the dollar sign at the beginning of each line is not to be typed, but rather signifies the line where commands are typed, hence the name “command line”. We will keep this notation throughout the book. As with most GNU/Linux commands, the command can take many options; for example, (that’s an lower case L) presents a detailed list, and (that’s the number one) puts each file on one line. The current directory, or “working directory} is the directory that we are in, and upon which all commands we run will act on. We can then see what directory we are currently in with the command , which is short for “print working directory”:

$ pwd
/home/dhendrix

We can then change directories with the command cd, which is short for “change directory”. For example, let’s use this command to go into the scripts directory:

$ cd scripts
$ pwd
/home/dhendrix/scripts

So with these examples, it is clear that directories can contain files and other directories, but they also represent an abstract location in which you are working.

We can create a file in countless ways, Let’s start with the creation of a basic text file. There are many programs called text editors that can be used to create a text file, such as , and . For , we can create a file by simply typing , where is the file we are creating. Since we are in the scripts directory, let’s create a simple python script called “helloworld.py”, which as the name suggests is a “hello world script”. This is the first script that people typically learn in a programming language because it is the simplest script that actually does something. To do this, while still in the scripts directory, let’s type , and type the following code:

print("hello world")

Then to save, we can type (or the Control key and X at the same time, followed by the Control key and the S at the same time). After it is saved, we can type to close the program. We can then run our new script by typing the command , which should produce the following output:

$ python helloworld.py
hello world

If you’ve never programmed before, you’ve just created your first script! The preceding was by no means a complete introduction to GNU/Linux. If anything, you should now have a rudimentary starting point. Throughout these class notes, we will learn new commands as we go, and new python code as we go. The intent here is to learn the essential pieces to do something simple, and keep learning as we go. In the end, we’ll have a decent set of commands to accomplish a lot.

1.2 Sequences, Strings, and the Genetic Code

1.2.1 Introduction to Sequences and Biopython

Central to the subject of bioinformatics is the study of sequences, and sequence analysis. So we begin with the notion of the biological sequence. Consider a character string x such as . Obviously, this sequence of characters is very different from a biological polymer such as DNA or protein, but this sequence is useful for the purposes of studying the bioinformatic properties of the molecule it represents. In the language of computer science, these character sequences are often called “strings”. We can represent a particular position i of the sequence as x[i]. For example, in this example, is . That is, when using a -based position, where the first position is numbered , the second position is , and so on. Although this may seem intuitive, most programming languages use as the start position.

Figure 1.1: A simple representation of a character string. Each position of the string can be accessed by its index number, which goes from 0 to the length minus one.

Let’s consider how sequences are represented in python. First, from the GNU/Linux terminal or shell, we can enter a “python terminal” by simply typing on the command line. This python terminal is like a sandbox in which we can create sequences and perform different sequence analyses. The same code we type in this python terminal should work in a python script, if we just save the same series of commands in a text file. That said, there might be some minor indentation/formatting changes needed. The python terminal is different from the GNU/Linux terminal, and the commands that work in one will typically not work in the other. The python terminal is indicated by three greater-than signs >>> , and we will use this notation in this book. We can define a sequence and access a position as follows:

>>> x = "ATGGCGGGAAAATGA"
>>> x
'ATGGCGGGAAAATGA'
>>> x[5]
'G'

Therefore, in python strings are -based and begin with position , then the second position is , the third is and so on. Defining the sequence in Biopython isn’t that much different:

>>> from Bio.Seq import Seq
>>> x = Seq("ATGGCGGGAAAATGA")
>>> x
Seq('ATGGCGGGAAAATGA', Alphabet())
>>> x[5]
'G'

Perhaps the difference is that Biopython takes the concept of a sequence, and makes it an “object”. Object oriented programming is a paradigm in computing where the central concept is the “object”, and there are defined “methods” that can be applied to these objects. Objects are defined as you might expect, as discrete units that fit some set of parameters or specifications as defined by a “class”. For example, the sequence object defined above has some similar features that a physical DNA molecule might have. We can transcribe our DNA sequence by simply typing:

>>> x.transcribe()
Seq('AUGGCGGGAAAAUGA', RNAAlphabet())

In this example, is the object, and is the method that is called on the object. Much of python is object oriented, and Biopython is very object oriented. We can see that this particular object is called a object. Note that in order to use Biopython we need to import it. We generally need to import a library or module with such a statement in order to use the object defined by them.

The advantage of importing libraries like this is that we get to use built in functions, methods, and objects. For example, a common quantity that one might be interested in regarding a nucleic acid sequence is its GC content. The GC content of a sequence is the percentage of the nucleotides in the sequence that are either G or C. Probably the easiest way to compute it is to import a function from another Biopython module:

>>> from Bio.SeqUtils import GC
>>> from Bio.Seq import Seq
>>> x = Seq("ATGGCGGGAAAATGA")
>>> GC(x)
46.666666666666664

1.2.2 The Central Dogma

Let’s now consider the central dogma of molecular biology. The central dogma states that biological information generally flows form DNA to RNA to proteins [1][2]. Similarly, with Biopython we can create objects corresponding to DNA, RNA, and Protein sequences. We can use methods on these objects to transcribe the DNA, and to translate the RNA to a protein using the and functions, respectively. Here’s how it’s done:

>>> from Bio.Seq import Seq
>>> DNA = Seq("ATGCTGGGATATTGA")
>>> print(DNA)
ATGCUGGGAUAUTGA
>>> DNA
Seq('ATGCTGGGATATTGA', Alphabet())
>>> mRNA = DNA.transcribe()
>>> print(mRNA)
AUGCUGGGAUAUUGA
>>> mRNA
Seq('AUGCUGGGAUAUUGA', RNAAlphabet())
>>> protein = mRNA.translate()
>>> print(protein)
MLGY*
>>> protein
Seq('MLGY*', HasStopCodon(ExtendedIUPACProtein(), '*'))

The python examples presented so far have been typed into the python terminal, but all of it could have been typed into a text file using one of many available text editors mentioned above, and run on the command line. For example, if the above code (all the stuff after the >>"}" title="\texttt{">>>"}" class="latex mathjax"> was typed into a text file called , it could be run with the command . As it stands right now, it wouldn’t print anything, highlighting one distinction between the python terminal and a python script; the variable name alone will not print its contents to the screen, but the print function will. Try putting the following code:

from Bio.Seq import Seq
DNA = Seq("ATGCTGGGATATTGA")
print(DNA)
mRNA = DNA.transcribe()
print(mRNA)
protein = mRNA.translate()
print(protein)

into a script file and run it. What is the output? Your output should look something like this:

ATGCTGGGATATTGA
AUGCUGGGAUAUUGA
MLGY*

Note that in the code above, we did not include the  that indicates the python terminal. This is because these symbols are only used in the python terminal, and are not included in ordinary python scripts. Therefore, they need to be removed when copying and pasting code from the python terminal to a script. An additional difference is we are not given the data type information that we get when using the python terminal. Printing the object just gives the actual string of that sequence.

Exercise 1

Transcribe the following DNA sequence and translate the resulting mRNA sequence. What protein sequence do you get?

>CDS
ATGCTGCGGCGAGCTCCGTCCAAAAGAAAATGGGGTTTGGTGTAAATCTGGGGGTGTAATGTTATCATAT
AAATAAAAGAAAATGTAAAAAACAAAAACAAAAACAAAAAAGCC

1.2.3 Subsequences and Reverse Complement

Subsequences for python strings work the same way as subsequences in objects. To specify a substring, we need to give a start and a stop position. The resulting substring will consist of the characters from the start to, but not including the stop. That is, the stop position is “non-inclusive”. By “non-inclusive”, it means it will contain the characters up to, but not including the stop position. This may seem strange at first, but there are some reasons for this that will become evident soon. For now, let’s just take a look at an example:

>>> from Bio.Seq import Seq
>>> DNA = Seq("ATGGCGGGAAAATGA")
>>> D1 = DNA[4:7]
>>> D1
Seq('CGG', Alphabet())

So, by specifying the range , we are including positions , , and . The resulting substring is essentially a concatenation of the characters at these positions. What if we create a substring that is the whole string itself? The result is like this:

>>> x = Seq('ACGTACGTACGT')
>>> y = x[0:len(x)]
>>> y
Seq('ACGTACGTACGT', Alphabet())

There are a couple of things going on here. First we create the string that is of length . We then create a substring that starts at position and goes up to, but not including, the length of , which is specified by . Since the indices of go from to , the substring is the same as .

If we want to print ever other character, we add to our range the step size, in this case .

>>> from Bio.Seq import Seq
>>> x = Seq('ACGTACGTACGT')
>>> y = x[0:len(x):2]
>>> y
Seq('AGAGAG', Alphabet())

It’s important to note that we don’t need to enter the start and stop position when we are traversing the whole sequence from to . We can accomplish the same thing with because the default option is to traverse the whole sequence, so with the start and stop missing it is implied that it is from to .

>>> x = Seq('ACGTACGTACGT')
>>> y = x[::2]
>>> y
Seq('AGAGAG', Alphabet())

Finally, we can reverse a sequence by using a negative step like . By using the default options on start and stop, we can get the whole sequence in reverse

>>> x = Seq('ACGTACGTACGT')
>>> x[::-1]
Seq('TGCATGCATGCA', Alphabet())

When using starts and stops in reverse, keep in mind that the same policy of going to, but not including, the stop position is used.

Because DNA is typically double stranded, we can think of the DNA sequences as being double stranded with the second strand implied. That is, even though we are given one strand of DNA, it is typically implied that there is an associated second strand that is the reverse complement of the first. We typically use the object to represent the forward strand of our sequence, but we may also want to know what the reverse strand would look like. We can produce the reverse complement very easily with the built-in Biopython method

>>> x = Seq('ACGTACGTACGT')
>>> x.reverse_complement()
Seq('ACGTACGTACGT', Alphabet())

In this case, we see that the reverse complement of our sequence is the same as . Such DNA sequences that are equal to their reverse complements are called “palindromes”. Regular palindrome phrases such as “race car” simply use the reverse, but DNA palindromes use the reverse complement.

Exercise 2

You can concatenate two sequences (attach them to each other to form one longer sequence) with the operator, such as . How can you use this to create a very long palindromic sequence?

1.3 Sequences File Formats

1.3.1 FASTA

Probably the most commonly used file format for sequences, and in fact one of the most common file formats of any kind in bioinformatics, is the FASTA file format. The FASTA file format has its origins in the program FAST, used for sequence alignment [3]. The File format is simply defined as a flat text file with one or more entries consisting of one line with a  symbol followed by a unique identifying definition line, or “defline”, and one or more lines of sequence data. Here is an example. Consider a plain text file called with the following as its only content:

>a
ACGCGTACGTGACGACGATCG
>b
ATTTCGCGACTCTGCCTACGCTAC
>c
GGGAAACCTTTTTTT

The requirement that the file is “plain text”, or without formatting and with a limited character set, is important. All too often, beginners to bioinformatics store sequence data in richer formats such as a word processing document. These types of files have fonts, font sizes, and other formatting markup that you can’t see when viewing in the word processing application. Therefore, these files are best dealt with in text editors like the ones mentioned above. To view a FASTA file from the command line without editing it, try the applications and with commands like:

$ less sequences.fa

To exit from simply type for “quit”. There aren’t really any restrictions on the sequence of the defline in the standard format, except that the defline should immediately follow the  without any intervening spaces. Any kind of sequence or set of sequences can be put into a FASTA file; all 3 billion base pairs in the complete human genome, a collection of protein sequences with zinc finger domains, or the DNA sequence of the promoter of your favorite mouse gene are all valid examples.

We can use methods found in existing libraries as part of Biopython to read the contents of a FASTA file pretty easily. For example, the following code

from Bio import SeqIO

fastaFile = "sequences.fa"
sequences = SeqIO.parse(open(fastaFile),'fasta')

will read the contents of the file and store it into the object . Now that we’ve read in the sequences, let’s try and print them back out. To do this, we will need to introduce a “for” loop.

from Bio import SeqIO

fastaFile = "sequences.fa"
sequences = SeqIO.parse(open(fastaFile),'fasta')
for record in sequences:
print(record.id, record.seq)

You’ll notice here that we have also introduced another object called here, and each sequence in our FASTA file is stored into such an object. Calling the methods and for the returns the defline text and a object. Now, let’s get back to the “for loop”. This is our first “control statement”, or a programming statement that controls when a command is to be executed. In general, a “for loop” defines a sequence of actions to be performed. In this case, the sequence of actions is to print the and for each entry in the FASTA file to standard out, which typically means print to the screen. Two other relevant pieces are the fact that the for loop statement has a colon at the end of the line, and the subsequent lines defined within the loop are indented. These are both required in python, but other programming languages may not have this restriction.

Exercise 3

Write a Biopython script that reads in a FASTA file, and prints a new FASTA file with the reverse complement of each sequence.

Write a script to read a FASTA file and print the reverse complement of each sequence. Print the GC content of each sequence.

1.3.2 FASTQ

Another very common sequence file format is “FASTQ”, which is commonly used in reporting data in high-throughput sequencing experiments. Each entry of the FASTQ file format begins with the “@” symbol, followed by a unique sequence identifier. Often this identifier encodes information about the machine and sequencing run from which this data was produced. This is followed by the corresponding read sequence, and then either just the “+” symbol, or the “+” symbol followed by the same unique identifier. Then finally there is a quality string. The quality string is the exact same length as the read sequence, and each character of it encodes a quality score. Here is an example of a single entry from a FASTQ file:

@SRR993731.910 HWI-ST880:148:D1F64ACXX:2:1101:18790:2441 length=51
TGTCCTCTGGCTCCAGGTCTCATGATGAAAAAATTTATGGAGTCCTGGACA
+SRR993731.910 HWI-ST880:148:D1F64ACXX:2:1101:18790:2441 length=51
;?<D;,2BA=3A+2AE;<<A<EEF>E@F<FFIA?DDCD<D;D9B9?##

The quality score is a character-encoded PHRED Score, defined as

(1.1)

where is the probability of a base-call error for that position [4]. This probability is computed by software as part of the base-calling pipeline during sequencing. Therefore, a high PHRED Score corresponds to a low probability of error. The value of the PHRED score corresponds to the ascii value of the character minus 33. In other words,

(1.2)

where is a python function that returns the ASCII value of an input ASCII character. Similarly, the function will return an ASCII character corresponding to the input ASCII value. The reason for the – 33 is because the ASCII characters prior to 33 are non-printable characters, such as “space”, “tab”, and “return”. Therefore, subtracting 33 sets the first printable character’s value to 0.

The PHRED score is our first example of how probabilities are used to define scoring systems in bioinformatics. In many ways, probabilistically defined scoring systems are preferred because it makes the scores more easily interpretable, and suggests a formalism with which to compute them.

We can read a FASTQ file much the same way that two read in FASTA files. The method is able to parse this format with essentially the same command as above, but with the string “FASTQ” as the second parameter to specify file format:

from Bio import SeqIO
fastqFile = "reads.fastq" # FASTQ file name
data = SeqIO.parse(fastqFile,"fastq") # parse the FASTQ file

We expect the object to be a list containing many records from the FASTQ file, which can be quite large. These records can be retrieved with the following, such as the loop:

>>> from Bio import SeqIO
>>> fastqFile = "reads.fastq"
>>> data = SeqIO.parse(fastqFile,"fastq")
>>> for record in data:
...     print(record)
...
ID: DB775P1:316:C4AGUACXX:2:1101:1748:1985
Name: DB775P1:316:C4AGUACXX:2:1101:1748:1985
Description: DB775P1:316:C4AGUACXX:2:1101:1748:1985 1:N:0:CAGATC
Number of features: 0
Per letter annotation for: phred_quality
Seq('TTTTTGTTGGGAANCTTGTGAGATTTTTGTAAATGATCGCAGTCACTTGNGCCT...GAT', SingleLetterAlphabet())
ID: DB775P1:316:C4AGUACXX:2:1101:1864:1989
Name: DB775P1:316:C4AGUACXX:2:1101:1864:1989
Description: DB775P1:316:C4AGUACXX:2:1101:1864:1989 1:N:0:CAGATC
Number of features: 0
Per letter annotation for: phred_quality
Seq('CATACACAACATANATTTGCTCATTTAGTTCCTCAAGGAACACCCGCTANTCTT...ACC', SingleLetterAlphabet())
ID: DB775P1:316:C4AGUACXX:2:1101:2323:1990
Name: DB775P1:316:C4AGUACXX:2:1101:2323:1990
Description: DB775P1:316:C4AGUACXX:2:1101:2323:1990 1:N:0:CAGATC
Number of features: 0
Per letter annotation for: phred_quality
Seq('TATGGACTACGCCGTCGAGACGGCTCACTTTGGTCTGTTCTTTAACATGNGCCA...ACG', SingleLetterAlphabet())

We can see that the function will return the PHRED quality score for each letter/character. This is returned as a list of integers, and we’ll evaluate that in a moment. To get the object from each of the these records, we can use the method:

>>> data = SeqIO.parse(fastqFile,"fastq")
>>> for record in data:
...     print(record.seq)
...
TTTTTGTTGGGAANCTTGTGAGATTTTTGTAAATGATCGCAGTCACTTGNGCCTTCAGTTGNANTCTCGATTNATNTGGAAGGTTTCAGCC
CATACACAACATANATTTGCTCATTTAGTTCCTCAAGGAACACCCGCTANTCTTATACCTTNTNTAGTATGTTTTNAAACTATTAGAAATA
TATGGACTACGCCGTCGAGACGGCTCACTTTGGTCTGTTCTTTAACATGNGCCAGTGCTGCNGNGCAGGATCTCGNACTTTCGTGGAGGAC

Consider the task of computing the average error probability, , from Equation 1.1 as a function of position. We’ll need to rearrange 1.1 and compute in terms of the .

(1.3)

One major “gotcha” when computing this kind of quantity is integer division. For python 2.7, plugging this equation into the python terminal will round the result, but in python 3, this issue is fixed. For example, let’s try two ways of computing the probability in python:

>>> 10**(-(ord("b")-33)/10)
1e-07
>>> 10**(-(ord("b"-33)/10.0)
3.162277660168379e-07

In the first case, we are dividing the exponent by , and the result is rounded. In the second case, we are using Equation , and we get the result using floating point arithmetic, which provides the desired answer.

Fortunately, the method returns a list of the values of , as if they were computed from 1.2. Therefore, we have the following

(1.4)

In practice, we’ll have to open the FASTQ file as above, and loop through the records, and we’ll need to store the probabilities to a list. Before we jump into that, let’s consider a simple example. Here is some basic python code to work with a list object. We can declare a list, in our case a list of probabilities, by the command . Next, we can add items to our list with the command with a value in the parentheses.

>>> p = []
>>> p.append(0.01)
>>> p.append(0.001)
>>> p.append(0.005)
>>> print(p)
[0.01, 0.001, 0.005]
>>> print(sum(p)/len(p))
0.00533333333333

As you can see, we append three values to our list of probabilities, and then print the average probability, by simply computing the arithmetic mean. Let’s put all these ideas together, along with a new function , that will plot a curve for us when given some values:

import sys
from Bio import SeqIO
import matplotlib as mpl
mpl.use('Agg')
import matplotlib.pyplot as pyplot

# this section takes care of reading in data from user
if len(sys.argv) != 2 or "-h" in sys.argv or "--help" in sys.argv:
print("Usage: printAverageQualityScores.py")
sys.exit()

# read in the FASTQ file name
fastq = sys.argv[1]
# parse the FASTQ file
data = SeqIO.parse(fastq,"fastq")

# initialize a list of probabilities
sum_p = [0] * 200
N = [0] * 200}$

for record in data:
for i,Q in enumerate(record.letter_annotations["phred_quality"]):
# convert the PHRED score to a probability
p_err = 10**(-float(Q)/10.0)
# append this specific probabity to array for this sequence
sum_p[i] += p_err
N[i] += 1

# now for plotting. Initialize x and y arrays to plot:
x = []
y = []

# add the average probabilties to the y values:
for i in range(len(N)):
if N[i] > 0:
pAvg = sum_p[i]/N[i]
x.append(i)
y.append(pAvg)

# plot the x and y values:
pyplot.plot(x,y)
pyplot.xlabel('position (nt)')
pyplot.ylabel('Averge error probability')
pyplot.savefig('quality.png')

The result of running this script on some sample data, in this case the reads for some ChIP-seq data for the Heat Shock Factor (HSF) in Drosophila melanogaster [5]. Note that running this script can take a long time for typical high-throughput sequencing experiments as they can be quite large. This result suggest that the probability of an error increases as the position gets larger. To some degree, errors in sequencing accumulate as once proceeds across the sequence. Repeated occurrences of the same character, or “homopolymers”, are a common source of errors in high-throughput sequencing.

Figure 1.2: The average error probability as a function of position for ChIP-seq data for Heat Shock Factor (HSF) in Drosophila melanogaster

Exercise 4

What probability of error corresponds to the following quality score characters:

- !

- A

- a

- %

1.3.3. GNU/Linux and Sequence Files

Another useful GNU/Linux command is (short for “word count”) to count the number of lines in a file. If a FASTA file was created such that it is one line per sequence, such as the example given above in section 1.3.1, then it is easy to count the number of sequences in the file. For example,

$ wc sequences.fasta
6 6 72 sequences.fasta

has returned three numbers. The first number is the number of lines in the file, which is . The second and third numbers are the number of words and number of characters, respectively. For a FASTA file with the restriction that each sequence occupies only one line in the file, then the number of sequences is the number of lines divided by 2. Similarly, the number of sequences in a FASTQ file are the number of lines divided by . For the more general FASTA file, when a sequence could in theory occupy any number of lines in the file, we could count with the deflines’ “>” symbol. We can grab out those lines with the program grep, which works by printing matches to a pattern with the syntax "grep < pattern > < file >" . So we could use the “>” symbol as our pattern, and then pipe it into “wc” with a command like:

$ grep ">" sequences.fasta | wc
3 3 9

This command produces the number of sequence entries in the file as the first number in the output (the 3). The pipe “” symbol sends the output of one program into the input of another program, as with this example we are sending the output of into for counting. Word of warning: be careful to ensure that the “>” symbol is enclosed in quotes in this command. Failure to include the quotes will result in over-writing the file, because of the meaning of the “>” symbol on the command line. Normally, a “>” symbol on the command line redirects the output of a program to a file, and over-writes the file. For example, we could even re-direct the output of the previous example like so:

$ grep ">" sequences.fasta | wc > numSequences.txt

thereby printing nothing to the screen, but storing the information to a file. Here is a table summarizing all the GNU/Linux commands discussed in this chapter.

ls

list the files in the current directory

ls -a

list all files including hidden files

ls -l

long formatted list of files

cd dir

change directory to dir

cd

change to home

pwd

show current directory

mkdir dir

create a directory dir

rm file

delete file

cp file1 file2

copy file1 to file2

mv file1 file2

rename or move file1 to file2; if file2 is a directory, move file1 into directory file2

cat file

print the contents of a file to the screen

more file

output the contents of file

less file

output the contents of file

head file

output the first 10 lines of file

tail file

output the last 10 lines of file

tail -f file

output the contents of file as it grows, starting with the last 10 lines

grep pattern file

print the lines matching pattern in the file.

wc -l file

print the number of lines in a file.

1.4 Lab 1: Introduction to GNU/Linux and Fasta files

Open up a terminal application and connect your server. Mac users can find the terminal under ““. Windows users can try MobaXterm . You will need to know your “username”, “password”, and “hostname” of the server you are connecting to. Some applications like the Mac Terminal and MobaXterm will let you navigate through files on your own computer as well using the standard Linux commands described in Table 1.3.3.

Once logged into the server, within your home directory, create a directory called “Lab1”. You can then into this directory.

$ mkdir Lab1
$ cd Lab1

You can type now, but you won’t see anything. Now create a scripts directory. It is helpful to create a directory for scripts, to keep them organized and not cluttered with data files.

$ mkdir Scripts

Next, let’s create a simple python script example. First, let’s open emacs to create a file:

$ emacs Scripts/helloworld.py

Next, we’ll type into the file the following command:

print("Hello, world!")

To save the file, let’s simply type . Then to close emacs, type . Once you’ve successfully saved and closed the file, we can run it with the following command:

$ python Scripts/helloworld.py

If you haven’t programmed anything before, this is your first program!

Now let’s download a script file, and a FASTA file. You can do this with the command , which can be run on the command line to directly download any file from the web:

$ wget http://hendrixlab.cgrb.oregonstate.edu/teaching/readFastaFile.py
$ wget http://hendrixlab.cgrb.oregonstate.edu/teaching/sequences.fasta

Try printing the contents of the FASTA file to the screen:

$ mv readFastaFile.py Scripts/.

Next, let’s move the script file “” into your scripts directory. You can do this with the command

$ mv readFastaFile.py Scripts/.

Please note the dot “.” at the end after the slash. Actually, this command will work without the dot, but it’s good to get in the habit of using it. The dot means move the file, and keep the name unchanged. If you were renaming the file, you could replace the dot with another name. You can run this script by simply typing the following:

$ python Scripts/readFastaFile.py sequences.fasta

The output of the script is just the names and sequences that are in the file. Now modify the script to print out a new FASTA file containing the reverse complement of each of the input sequences.

1.5 Biological Sequence Database

There are a tremendous number of sequence databases online today. New databases are published every month. In fact, there are journals devoted to databases! In what follows we will present some of the most commonly used databases for biological sequences. In later chapters, we will introduce databases that are relevant to that chapter.

Each of these databases have an associated webpage allowing data to be downloaded. Although these are useful, they are frequently updated and any book or notes on these interfaces would be quickly out of date. Furthermore, most people can probably easily figure out how to use their interfaces effectively. Therefore, we’ll focus on how to make use of Biopython to write scripts to download data from these resources on the command line, and learn what differentiates these databases.

1.5.1 NCBI

The National Center for Biotechnology Information isn’t a single database, but rather a very large national resource for biomedical and genomic information [6]. Their website can be accessed at and consists of a tremendous amount of biological sequence data in various forms from full genomes, to proteins, to single nucleotide polymorphisms (SNPs), to uploaded high-throughput sequence data.

Genbank

Genbank contains most of the world’s known DNA, RNA, and protein sequences, and stores bibliographic information for the sequences as well [7]. The current release of Genbank, as of October 2015, contains sequences, comprising nucleotides. You can access Genbank entries from their main page or through NCBI Nucleotide at . Genbank is archival in nature, so it can contain redundant entries.
Genbank is a database, but it is also a very detailed file format. We can use the Biopython module to retrieve a file in Genbank format. If we already know a gi accession number for the gene, then we can retrieve the data directly from NCBI. Here is an example, with portion of the output (the full Genbank file is much too large to present here).

>>> from Bio import Entrez
>>> Entrez.email = "example@oregonstate.edu"
>>> p_handle = Entrez.efetch(db="protein", id='4507341', rettype="gb", retmode="text")
>>> print(p_handle.read())
LOCUS       NP_003173    129 aa    DNA_input linear  PRI 28-NOV-2015
DEFINITION  protachykinin-1 isoform beta precursor [Homo sapiens].
ACCESSION   NP_003173
VERSION     NP_003173.1 GI:4507341
DBSOURCE    REFSEQ: accession NM_003182.2
KEYWORDS    RefSeq.
SOURCE      Homo sapiens (human)
ORGANISM  Homo sapiens
Eukaryota; Metazoa; Chordata; Craniata; Vertebrata; Euteleostomi;
Mammalia; Eutheria; Euarchontoglires; Primates; Haplorrhini;
Catarrhini; Hominidae; Homo.
REFERENCE 1 (residues 1 to 129)
AUTHOR    Agaeva,G.A., Agaeva,U.T. and Godjaev,N.M.
TITLE     [Particularities of Spatial Organization of Human Hemokinin-1
and Mouse/Rat Hemokinin-1 Molecules]
JOURNAL   Biofizika 60 (3), 457-470 (2015)
PUBMED    26349209
REMARK    GeneRIF: The spatial structures of human, mouse, and rat
hemokinin-1 protein isoforms have been presented.
...
FEATURES    Location/Qualifiers
source    1..129
/organism="Homo sapiens"
/db_xref="taxon:9606"
/chromosome="7"
/map="7q21-q22"
Protein   1..129
/product="protachykinin-1 isoform beta precursor"
/note="neuropeptide gamma; neuropeptide K; tachykinin,
precursor 1 (substance K, substance P, neurokinin 1,
neurokinin 2, neuromedin L, neurokinin alpha, neuropeptide
K, neuropeptide gamma); tachykinin 2; protachykinin;
preprotachykinin; neurokinin A; protachykinin-1; PPT"
...
CDS       1..129
/gene="TAC1"
/gene_synonym="Hs.2563; NK2; NKNA; NPK; TAC2"
/coded_by="NM_003182.2:247..636"
/note="isoform beta precursor is encoded by transcript variant beta"
/db_xref="CCDS:CCDS5649.1"
/db_xref="GeneID:6863"
/db_xref="HGNC:HGNC:11517"
/db_xref="HPRD:08876"
/db_xref="MIM:162320"
ORIGIN
1 mkilvalavf flvstqlfae eiganddlny wsdwydsdqi keelpepfeh llqriarrpk
61 pqqffglmgk rdadssiekq vallkalygh gqishkrhkt dsfvglmgkr alnsvayers
121 amqnyerrr

Genbank format is meant to be human-readable to some extent, but the files have formatting that allows them to be parsed by software tools, such as that of Biopython.

RefSeq

RefSeq (Reference Sequence) is a comprehensive and non-redundant database of biological sequences including DNA, RNA, and proteins [8]. In addition, the proteins and nucleic acid sequences corresponding to a specific gene are linked in the database. RefSeq currently maintains approximately sequence entries, each with a unique ID associated with it. You can recognize a RefSeq ID because they all have the same format:

- followed by numerical digits corresponds to RNA sequences that are not messenger RNAs.

- followed by numerical digits corresponds to messenger RNAs

- followed by numerical digits corresponds to proteins

Known and confirmed sequence IDs

- followed by numerical digits corresponds to RNA sequences that are not messenger RNAs.

- followed by numerical digits corresponds to messenger RNAs

- followed by numerical digits corresponds to proteins

We can retrieve a biological sequence for a RefSeq sequence with the Biopython module . This allows us to download sequences directly from NCBI, and print out a FASTA file.

>>> from Bio import Entrez
>>> Entrez.email = "example@oregonstate.edu"
>>> rec = Entrez.read(Entrez.esearch(db="protein", term="NP_003173"))
>>> print(rec['IdList'])
['4507341']
>>> p_handle = Entrez.efetch(db="protein", id=rec["IdList"][0], rettype="fasta";)
>>> print(p_handle.read())
>gi|4507341|ref|NP_003173.1| protachykinin-1 isoform beta precursor [Homo sapiens]
MKILVALAVFFLVSTQLFAEEIGANDDLNYWSDWYDSDQIKEELPEPFEHLLQRIARRPKPQQFFGLMGK
RDADSSIEKQVALLKALYGHGQISHKRHKTDSFVGLMGKRALNSVAYERSAMQNYERRR

Other Databases at NCBI

NCBI is a tremendous resource with numerous useful databases. Here are a few that may be useful:

- Entrez – GQuery NCBI Global Cross-database search

- Gene – Gene is a database of genes that integrates many species. The records include nomenclature, RefSeq, maps, pathways, variations, phenotypes, links to genomes

- UniGene – An attempt to computationally identify unique transcripts from the same locus, and analyze expression by tissue, age, heath status and a number of factors. It also reports related proteins and clone resources.

- GEO – Gene Expression Omnibus. A huge database of curated gene expression data, and high-throughput sequence data.

- OMIM – Comprehensive compendium of heritable traits. Includes genetic phenotypes, and allelic variants if they have been identified. This database is also associated with “morbidMap”, which identifies traits associated with diseases.

- Taxonomy – A curated database of classification and nomenclature for all organisms in the sequence database . According to NCBI, this represents 10% of the species on the planet.

1.5.2 Ensembl

The Ensembl database is a European database of biological sequences and other data. This database is produced by the European Bioinformatics Institute (EBI), which is part of the European Molecular Biology Laboratory (EMBL). This databae can be accessed at http://ensembl.org, and provides a huge database of biological data, in many ways similar to NCBI/GenBank. Ensembl also provides BioMart, which is a user friendly interface to retrieve genes, transcripts, and protein sequences.

1.5.3 UCSC Genome Bioinformatics

The University of California at Santa Cruz (UCSC) has a long history in genomics and bioinformatics research. Their website and associated databases called “UCSC Genome Bioinformatics” is a tremendous resource for data and tools for doing genomics, primarily on animal models systems. Their website at has a genome browser, download page, and software such as BLAT for rapid alignment of sequences. BLAT results are integrated into the genome browser for alignment visualization.

1.5.4 Uniprot

The Universal Protein Resource (UniProt, www.uniprot.org) collects and provides a thorough database of protein sequences. UniProt is a collaboration between EMBL, the Swiss Institute of Bioinformatics, and the Protein Information Resource (PIR). UniProt actually consists of multiple databases. For example, UniRef is a database of clustered sets of sequences at varying levels of sequence identity. UniRef100 is a collection of sequences that are identical and at least or more residues in length. This set comprises such clusters as of right now. Similarly, UniRef90 contains clusters that are at least % identical and contains clusters.

Exercise 5

What are the differences between and ? Name at least two.

1.6 Lab 2: FASTQ and Quality Scores

In this lab we will review the FASTQ file and quality scores.
Once logged into the server, within your home directory, create a directory called “Lab2”. You can then into this directory.

$ mkdir Lab2
$ cd Lab2

You can get in the habit of creating scripts directories for each lab project once you have entered that lab’s specific directory.

$ mkdir Scripts

Next, please download the python script and a FASTQ file with these commands.

$ wget http://hendrixlab.cgrb.oregonstate.edu/teaching/printAverageQualityScores.py
$ wget http://hendrixlab.cgrb.oregonstate.edu/teaching/reads.fastq

Move the script into your scripts folder with a “mv” command.

$ mv printAverageQualityScores.py Scripts/.

Next, let’s count the number of lines in the FASTQ file. Given the result, how many reads do you think are contained in this file?

$ wc -l reads.fastq

Next, let’s check what happens when you run the script without any input, or with the “-h” flag.

$ python Scripts/printAverageQualityScores.py
$ python Scripts/printAverageQualityScores.py -h

This is a usage statement printed out, making it clear how to run the script. This is very helpful for writing scripts and is good practice. Now let’s run the script with the input file, which will plot the average error probability as a function of position:

$ python Scripts/printAverageQualityScores.py reads.fastq

The script will create a file called . How does this plot look? What does this say about errors in sequencing data?
Now modify the script to instead of printing the average probability, to print the average quality score. How does this plot compare to the probability plot?
Now let’s use a quality trimmer called to trim our reads based on a quality threshold of 20. First, let’s see how this program is run.

$ fastq_quality_trimmer -h

usage: fastq_quality_trimmer [-h] [-v] [-t N] [-l N] [-z] [-i INFILE] [-o OUTFILE]
Part of FASTX Toolkit 0.0.13 by A. Gordon (gordon@cshl.edu)

[-h] = This helpful help screen.
[-t N] = Quality threshold - nucleotides with lower
quality will be trimmed (from the end of the sequence).
[-l N] = Minimum length - sequences shorter than this (after trimming)
will be discarded. Default = 0 = no minimum length.
[-z] = Compress output with GZIP.
[-i INFILE] = FASTQ input file. default is STDIN.
[-o OUTFILE] = FASTQ output file. default is STDOUT.
[-v] = Verbose - report number of sequences.
If [-o] is specified, report will be printed to STDOUT.
If [-o] is not specified (and output goes to STDOUT),
report will be printed to STDERR.

The help menu doesn’t tell you what value to subtract for some reason. This is done with the flag. The full command is this:

$ fastq_quality_trimmer -t 20 -i reads.fastq -o reads.qc.fastq -Q 33

Note: the script will no longer work for data of varying length. Now that we’ve trimmed the reads, let’s see how the resulting length distribution looks. This can be achieved by a second modification of the script to simply print the lengths of the trimmed reads.

$ wget http://hendrixlab.cgrb.oregonstate.edu/teaching/printFastqReadLengths.py

Move this script into your scripts folder by typing something like this:

$ mv printFastqReadLengths.py Scripts/.

Where the dot indicates to move the file and keep the same name. This script can be run by the following command, similar to the first script:

$ python Scripts/printFastqReadLengths.py reads.qc.fastq

Please take a look at this script and see how it works. This script will also produce a pdf file, showing the resulting read length distribution. The difference is, it is on a log-scale on the y-axis using the line:

pyplot.yscale('log')

Question: How does the read quality data compare to the length data after trimming? Is the result reasonable?

Chapter 2: Sequence Motifs

2.1 Introduction to Motifs

A biological motif, broadly speaking, is a pattern found occurring in a set of biological sequences, such as in DNA or protein sequences. A motif could be an exact sequence, such as , or it could be a degenerate consensus sequence, allowing for ambiguous characters, such as for or . Motifs can also be described by a probabilistic model, such as a position-specific scoring matrix (PSSM) or weight matrix.

2.2 String Matching

Frequently we want to search for an exact string or pattern within a larger sequence. For example, when using a restriction enzyme to cut a larger sequence at sequence-specific sites that match a particular pattern, we would like know where this pattern occurs in the larger sequence. This task can be called string matching. There are numerous algorithms that have been developed to make this task more efficient, such as the Knuth-Morris-Pratt algorithm [9] and the Burrows-Wheeler transform [10]. These approaches are beyond the scope of this course, but definitely worth mentioning. Let’s examine how to utilize the Biopython function as part of the module. We can use the function as follows to search for the short pattern .

>>> from Bio.Seq import Seq
>>> from Bio import SeqUtils
>>> pattern = Seq("ACG")
>>> sequence = Seq("ATGCGCGACGGCGTGATCAGCTTATAGCCGTACGACTGCTGCAACGTGACTGAT")
>>> results = SeqUtils.nt_search(str(sequence),pattern)
>>> print(results)
['ACG', 7, 31, 43]

You’ll note that the function takes in two arguments. The first is the sequence to search, but it is not a object, but the basic python string. The python function “str()” converts the object into a string. The second argument is the pattern or string that we are searching, but this argument clearly can be the object, which it is in this example.

Typically, we consider the DNA sequence that we are searching as double stranded, hence we want to search the forward strand and its reverse complement. In many cases for bioinformatics, such as searching an entire chromosome, it is easier to reverse complement the pattern rather that the sequence to search. In this example, this looks like

>>> results_rc = SeqUtils.nt_search(str(sequence),pattern.reverse_complement())
>>> print(results_rc)
['CGT', 11, 28, 44]

In this example, we have searched the forward strand of the DNA sequence. Alternatively, a biologist might want to know where the patterns occurs with positions defined along the reverse complement of the sequence we are searching, or the reverse strand if the DNA is double stranded. In the example of the restriction enzyme, positions defined along the reverse complement of the larger sequence are more useful for predicting the length of the resulting fragments after treatment by the restriction enzyme. This is performed as follows.

>>> results_rc = SeqUtils.nt_search(str(sequence.reverse_complement()),pattern)
>>> print(results_rc)
['ACG', 7, 23, 40]

The results are clearly different in terms of the positions, but they represent the same information. Because the larger sequence has a length of nt, we can see that we can transform between the two sets of position with the following equation:

(2.1)

Where is the position resulting from searching the reverse complement of the pattern over the forward strand of the sequence, and is the position of searching the pattern over the reverse strand of the sequence.

Exercise 6

Find the number of occurrences and locations of within the following sequence:

>sequence
AGCGATCTAGCATACTTATACGCGCGCAGCTATCGATCACTTGTGCTAGTAAAGTGCGCGCCGCA
TTAAAGTGCTAGCTAGCTACTTAGCTAGCTAGTCG

2.3 Consensus Sequences

A consensus sequence is a string of either nucleotide or protein characters along with “degenerate characters”, which specify a subset of characters. These degenerate characters can act as “wild cards”, such as , which can refer to any character, and are summarized in Table 2.1 [11]. They can also specify a more specific subset, such as , which specifies the pyrimidines and .

Symbol

Meaning

Mnemonic

R

A, G

puRine

Y

C, T

pYrimidine

W

A, T

Weak (weaker basepairs, fewer hydrogen bonds)

S

G, C

Strong (stronger basepairs, more hydrogen bonds)

K

G or T

Keto (both have a keto group)

M

A or C

aMine (both have an amine group)

B

C, G, T

not A (B comes after A)

D

A, G, T

not C (D comes after C)

A

A, C, T

not G (H comes after G)

V

A, C, G

not T or U (V comes after T and U)

N

A, C, G, T

aNy base

Table 2.1: IUPAC codes for nucleotides. In this table, everywhere that T applies, U applies as well.

We can combine the standard nucleotides and IUPAC codes to form a “consensus sequence” describing a motif. For example, “Downstream Promoter Element” (DPE) in Drosophila occurs near position +28 from the TSS of many genes, and has the consensus sequence V [12].

2.3.1 Searcing Consensus Sequences with Biopython

Biopython can also be used to search for a consensus sequence. The function is built to recognize patterns that include the wild-card characters presented in Table 2.1. If we are to search for the consensus sequence for the DPE, we would use something like this:

>>> from Bio import SeqUtils
>>> consensus = "RGWYV"
>>> sequence = "CGTAGCTAGCTCAGAGCAGGGACACGTGCTAGCAACAGCGCT"
>>> SeqUtils.nt_search(sequence,consensus)
['[AG]G[AT][CT][ACG]', 19]

As before, the results contain the pattern searched and the positions of the instances, but the wild-chard characters of the pattern are now represented in a more traditional regular expression format, with sets of character within square brackets.

2.4 Motif Finding

Motif finding can be described as the process of discovering patterns within collections of sequences. In many cases, we don’t know what the pattern looks like, and this task is a “needle in a haystack” challenge, that involves sifting through many -mers that aren’t part of the pattern, but occur frequently.

2.4.1 Sequence Complexity

Before we can start looking for motifs, we’ll need to consider things are frequently occur in biological sequence datasets, in particular DNA sequences. The most frequent -mer in the human genome is “AAAAAAAA”. Such a sequence can be called a “low complexity” sequence, and along with simple repeats, are commonly occurring sequences that can confound motif finding and sequence alignment. Sequences like “ATATATATAT” are also frequently occurring in virtually any sequence dataset, and would be discovered by a motif search if not filtered out.
The Wooton-Federhen complexity is a score that quantifies the complexity of a sequence [13]. Put simply, the WF complexity quantifies the number of possible sequences that could be generated using the same number of As , number of Cs , number of Gs  and number of Ts . This can be computed from a multinomial coefficient, which finds itself in a log in the equation:

Note that the is base , which is the size of the alphabet. A protein complexity score can be computed in an analogous fashion using base for amino acids. The factor of , where is the number of characters in the sequence, is to ensure the value is between and .

There is a program called (R. Tatusov and D.J. Lipman unpublished) that can mask sequence of low complexity. It can be run by specifying the sequence and a threshold.

$ dust sequences.fasta 30

Here we are specifying a threshold of , where the default is . Clearly, is specifying a complexity score different from WF complexity because it can be greater than .

2.4.2 Weight Matrices

Weight matrices are the most general representation of a motif. It is a probabilistic model that we will see can be used to compute the log-likelihood of a string being an instance of the motif compared to a “background” model.

Probabilistic Models of Motifs

The concept of the Position Specific Scoring Matrix (PSSM), also known as a weight matrix, was developed by Stormo et al [14].
Let’s begin by first by defining what isn’t a motif, using what is called a background model. A background model is a probabilistic representation of what a typical sequence looks like. The simplest background model is defined by single nucleotide probabilities . Under such a model, the probability of a sequence is computed as the product of the individual frequencies in the sequence. This can be expressed as

(2.2)

Consider the nucleotide frequencies for the human genome. These frequencies vary from position to position, and from chromosome to chromosome. Moreover, these frequencies will vary when comparing the promoters of genes and intergenic regions.
Our motif can be described by a matrix of probabilities . For a motif of length , we have a matrix like:

With this matrix, we can compute the probability (likelihood) of the sequence given that it is an instance of the motif by

(2.3)

Entropy and Information Content

One helpful way of describing such a model is the Shannon entropy [15]. Shannon entropy is a measure of the uncertainty of a model, in the sense of how unpredictable a sequence generated from such a model would be. For the single-nucleotide background model, the entropy is

Note that while Shannon entropy is typically denoted H, this is not to be confused with enthalpy, which is also represented with H. The entropy is maximized when each nucleotide is equally likely, that is if for all . It is intuitive that such a model would have the highest uncertainty, for example, compared to a model where and all other frequencies very low. Therefore, the maximum entropy of our background model is:

Since , we have

When the logarithms are base , the units for such a quantity is called “bits”, as is with BLAST scores (See Section 3.3). When using natural logs, the units are “nits”. We can think of this value of bits as the information content associated with knowing a particular nucleotide. A bit of information can also be understood as the number of questions necessary to unambiguously determine an unknown nucleotide. You could ask, “Is it a purine?” If the answer is “no”, you could then ask is it ? The answer to the second question always guarantees, non-canonical nucleotides aside, the nucleotide’s identity.
We can then compute the entropy at each position of our motif’s probability matrix by the expression

The Information Content of a motif at each position can be defined as the reduction in entropy. That is, the the motif provides information inasmuch as it reduces the uncertainty compared to the background model. If the change in entropy is , then the information content at position  is

Exercise 7

What is the entropy of a set of sequences with the nucleotide frequencies given by and ?

2.4.3 Relative Entropy

Another useful concept for quantifying the information in a motif is the “relative entropy”, which is also known as the Kullback-Leibler divergence and as “information gain” [16]. For biological motifs, the expression for relative entropy at a particular position is defined by the equation

(2.4)

This expression quantifies how different two discrete probability distributions are, in this case the background frequencies and one position of our motifs’s probability matrix

2.4.4 Building a Weight Matrix

A weight matrix can be built or defined when one has a collection of sequences that are determined to be instances of the motif. It is essential that the motif’s instances be precisely aligned, usually without gaps, such that each position of each sequence corresponds to the same position of the motif. Consider a set of sequences . We consider that each sequence corresponds to an instances of our motif. For example, the following are instances of the binding site for the Drosophila transcription factor Giant:

One useful contract is the Kronecker delta function. For this application, consider the expression , which returns when , and otherwise. For example, we can check if a particular nucleotide in sequence is equal to the nucleotide with the statement . Under this representation, this function is defined as:

(2.5)

This function is useful for connecting sequences to equations. For example, we can generate a matrix of counts , such that the elements contain the number of occurrences of nucleotide at position in instances of our motif:

This count matrix can then be normalized to represent the frequency of each nucleotide at each position by

Then, using equations 2.2 and 2.3 we can define a score as the log likelihood ratio of the probabilities of being an instance of the motif to being a random sequence. The expression of this score for a particular sequence is

As another example of the utility of the indicator matrix, this score can be represented as

In the second part of the equation, we have defined the weight matrix with terms given by

So, for a particular sequence, , this equation would pick out the appropriate terms necessary to compute the score

Note that for this score to make sense, you need to test a sequence that is the same length as the motif.

2.4.5 Biopython Motifs

The module

Now let’s see how we can go about building a weight matrix using Biopython [17]. First, we’ll need to import modules as usual. Next, we’ll need to specify a set of instances of our motif. In theory, this would be best done as reading in a FASTA file of instances. Next, we can create our motif from the list of instances. Here’s how this looks:

>>> from Bio import motifs
>>> from Bio.Seq import Seq
>>> instances = [Seq("CAGTT"),Seq("CATTT"),Seq("ATTA"),Seq("CAGTA"),Seq("CAGTT"),Seq("CAGTA")]
>>> motif = motifs.create(instances)
>>> print(motif.degenerate_consensus)
CAKTW

In this example, we have created a motif from a simple list of instances, each of which are objects, typed in using the method. The method is useful for creating a concise sequence representation. We lose some information by only keeping this consensus sequence, and not keeping the matrix. By creating this motif object, we can also print out the count matrix for this motif:

>>> print(motif.counts)
0 1 2 3
A: 0.00 0.00 3.00 0.00
C: 6.00 0.00 0.00 2.00
G: 0.00 6.00 3.00 2.00
T: 0.00 0.00 0.00 2.00

JASPAR sites

We often want to create a motif from external data sources. One such database of motifs is JASPAR http://jaspar.genereg.net [18]. On this webpage we can browse through various motifs. For example, in the “JASPAR CORE Insecta” section, we can see the motif MA0447.1 for the Giant binding motif http://jaspar.genereg.net/cgi-bin/jaspar_db.pl?ID=MA0447.1&rm=present&collection=CORE. After downloading the “sites” file as a FASTA file in the lower left, we can see from the file “MA0447.1.sites” that it is similar to FASTA, but contains some additional information:

>MA0447.1        gt        1
tttctgttttggcgtaTTTATGTGATgc
>MA0447.1        gt        2
ggtggcactaccctGTTACGCAATat
>MA0447.1        gt        3
tTTAATATAACgcttctatctttgttta
>MA0447.1        gt        4
gttgttacgcgtGTTACATAATgcttcg
>MA0447.1        gt        5
aaccactgtaaagctCCGGCGTATTggc
...

We can see that there is additional information encoded in the capitalized letters, indicating where the binding site is. Secondly, this format doesn’t work for most programs as a FASTA file because the string directly after the “” not unique for each line, but instead the unique information comes after with the numbering. We therefore need a special method to read in this information. The module has a method for reading this information. Let’s create a logo while we’re at it:

>>> from Bio import motifs
>>> motif = motifs.read(open("MA0447.1.sites"),"sites")
>>> print(motif.counts)
0 1 2 3 4 5 6 7 8 9
A: 9.00 0.00 1.00 29.00 0.00 5.00 0.00 30.00 33.00 0.00
C: 4.00 1.00 0.00 0.00 29.00 0.00 3.00 2.00 1.00 8.00
G: 18.00 1.00 2.00 4.00 0.00 30.00 1.00 1.00 0.00 4.00
T: 4.00 33.00 32.00 2.00 6.00 0.00 31.00 2.00 1.00 23.00

>>> motif.weblogo("giant_LOGO.pdf",format="pdf")

This last command produces a LOGO image and requires an internet connection. It actually sends the data to the website http://weblogo.berkeley.edu/ [19]. The resulting LOGO image looks like this, which is similar to the motif logo seen on the JASPAR database:

Figure 2.1: Sequence LOGO created with weblogo for the giant motif

We’ve already seen how to score a particular sequence with a weight matrix, and how to build a weight matrix form a collection of instances. In practice, we often have longer sequences where we don’t know precisely where the instances are, but know that the sequences contain subsequences that are instances of the motif.

2.5 Promoters

Broadly speaking, the promoter of a gene is the set of genomic DNA sequences that direct the gene’s transcription. This can include the core promoter that is the most proximal region around the gene’s transcription start site, and can also include other regions of genomic DNA that are involved in transcription initiation.

2.5.1 Core Promoters

The core promoter is a collection of binding sites, capable of directing the initiation of transcription, and located within from the TSS. In eukaryotes, this can include the TATA-box and the Initiator, as well as other motifs that are more organism-specific.

2.5.2 Databases of Promoters/TSSs

There are a number of databases that are devoted to cataloging the locations of the promoters and transcription start sites of genes.

EPD: Eukaryotic Promoter Database

The Eukaryotic Promoter Database (EPD) is a resource of experimentally validated promoter locations for animals, plants, and fungi. This database contains non-redundant entries of promoter locations for Humans, mouse, D. melanogaster, zebrafish, C. elegans, Arabidopsis thaliana, S. cerevisiae, and S. pombe. The database can be accessed at http://epd.vital-it.ch/.

2.6 De novo Motif Finding

2.6.1 Gibbs Sampling

Gibbs sampling has many applications in statistical physics in general [20]. Is sometimes called a Markov-chain Monte Carlo (MCMC) approach. Briefly, Gibbs sampling starts with a set of sequences in which there is a motif of unknown sequence content. First, random initial positions of the motif are selected. Next, a weight matrix is built from those instances. The weight matrix is then used to re-score the positions of the input sequences, to find a new optimal set of instances, and the weight matrix is subsequently updated. This process is repeated until convergence. Perhaps surprisingly, this method is capable of finding motifs quite well.
You might try as a challenge exercise, to implement Gibb’s Sampling with Biopython. For example, if you have a list of sequences , you can generate random positions for the initial motif instances for using the python module , and extract subsequences corresponding to these positions.

>>> from Bio.Seq import Seq
>>> import random
>>> seqs = [Seq('GTCGATCGATCGTACGTACGTACGTACGATGCTAGCTACGTACC'),
Seq('GGTTCGAGTCGAGCAAGAGCTAGCTAGCGACGTACTAC'),
Seq('GCTGATCATGCTAGCGCGTAGCTACGATCGTACGTACGATGAGCTAGCTACGTCTACGTACGTGCACA')]
>>> K = 7
>>> instances = []
>>> for seq in seqs:
...     instances.append(seq[j:j+K])
...
>>> print(instances)
[Seq('CGTACGT', Alphabet()), Seq('AAGAGCT', Alphabet()), Seq('CGCGTAG', Alphabet())]

You could create a motif from the instances as before. This will serve as the initial motif, and it would be updated repeatedly for some number of steps.

>>> from Bio import motifs
>>> motif = motifs.create(instances)
>>> print(motif.counts)
0    1    2    3    4    5    6
A: 1.00 1.00 0.00 2.00 0.00 1.00 0.00
C: 2.00 0.00 1.00 0.00 1.00 1.00 0.00
G: 0.00 2.00 1.00 1.00 1.00 1.00 1.00
T: 0.00 0.00 1.00 0.00 1.00 0.00 2.00
>>> weightMatrix = motif.pssm
>>> print(weightMatrix)
0    1    2    3    4    5    6
C: 1.42 -inf 0.42 -inf 0.42 0.42 -inf
G: -inf 1.42 0.42 0.42 0.42 0.42 0.42
T: -inf -inf 0.42 -inf 0.42 -inf 1.42

This weight matrix could in theory be used to collect the instances of the motif, and you could collect the top scoring instances for each of the input sequences, using something like this for each of the input sequences:

>>> for position,score in weightMatrix.search(seq):
...     print(position, score)

These updated instances could be used to create new motif, and this process could be repeated. I leave this as a challenge for the reader

2.6.2 MEME and the EM Algorithm

One of the most widely used software tools for motif discovery is MEME: Multiple EM for Motif Elicitation. MEME uses Expectation Maximization to find the best binding sites [21]. The algorithm computes the parameters that maximize the expected value of the (log-)likelihood pf the model. The EM algorithm attempts to learn the “missing data”, in this case the positions of the instances of the motif. Once the instances of the motif are identified, an updated weight matrix can be computed.
MEME starts with some initialization of motifs built from the enriched K-mers. Essentially this list of initial motifs is achieved through a heuristic method, which applies one iteration of EM on all K-mers. The algorithm computes an expression for the expected value of the log-likelihood. It then computes the value of the motif membership variable Z that maximizes the expected value of the log-likelihood. This process is repeated until convergence.
There are numerous options for running MEME. For example, one can select how the instances of the motifs are distributed in the input sequences. They are specified by the “-mod” option and can be:

1. oops – Exactly one instance per sequence.

2. zoops – One or zero instances per sequence.

3. anr – Any number of instances per sequence.

To find a motif in the input file “sequences.fa” that has a maximum motif length (specified with the -maxw option) and requiring exactly one occurence per sequence, we would use the following command

$ meme sequences.fa -dna -mod oops -maxw 12

Similarly, to find different motifs, each with a maximum length (width) of , and such that each sequence has zero or one instance, we would use the following command:

$ meme sequences.fa -dna -mod zoops -maxw 10 -nmotifs 2

For a full list of options, try typing on the command line.

2.7 Lab 3: Introduction to Motifs

In this lab, we are going to learn about various ways to produce motifs from sequence data. The purpose is to show what works and what doesn’t. Let’s start by creating a directory called , and into this directory, similarly as before.

2.7.1 Part 1: Building a motif and LOGO image

First, let’s download a list of sequences that contain a motif at random positions and and some flanking sequence.

$ wget http://hendrixlab.cgrb.oregonstate.edu/teaching/gt.fasta

How many lines are in this file? How many sequences does that mean are in the file?
Let’s try and create a motif from this data. What do you expect to happen?

>>> from Bio import SeqIO
>>> from Bio import motifs
>>> from Bio.Seq import Seq
>>> sequences = SeqIO.parse("gt.fasta","fasta")
>>> instances = []
>>> for data in sequences:
...     instance = data.seq
...     instances.append(instance)
...
>>> motif = motifs.create(instances)

Now the motif is created in the “motif” object. What does the count matrix look like? You can print it out by this command:

>>> print(motif.counts)

Create a motif LOGO from these counts. We can go through the web, or we can create one right now with Biopython:

>>> motif.weblogo("giant_1.pdf",format="pdf")
>>> quit()

Does this look like a good motif? How can you tell?

2.7.2 Part 2: JASPAR Database and “sites” format.

Now go to the website JASPAR at http://jaspar.genereg.net. Go to “JASPAR CORE Insecta”. Then find the motif for “giant”, labeled “gt” (motif ID MA0447.1). Click the motif LOGO image to get a pop-up of the motif details (you may need to enable pop-ups on your browser). How does the motif LOGO compare to what we found just by building a LOGO from the FASTA file of all the sites? (If you have issues downloading the file, you can get it here.)
Now let’s download the sites in JASPAR format, which is a variation on FASTA, with a little more information. Click the link that says “(Show me all binding sites)..as fasta file” in the lower left portion of the pop-up. Transfer the downloaded file to crick by using the folder icon for file transfer. By looking at the file, what differences do you see?

$ less MA0447.1.sites

Now, let’s load this data using a special function that is part of the motifs module:

>>> from Bio import motifs
>>> motif = motifs.read(open("MA0447.1.sites"),"sites")
>>> print(motif.counts)

How does this compare to what you found previously? Now create a motif LOGO and exit:

>>> motif.weblogo("giant_2.pdf",format="pdf")
>>> quit()

How does this motif compare to the image on the JASPAR webpage

2.7.3 Part 3: Running MEME on the command line

The reason why our first method didn’t work is because we just looked at sequences without performing motif finding. The JASPAR sites file encodes this information in the capital letters, letting it know where the motif instances are. Now let’s try motif finding on the first set of data to see if we can recover it.

$ meme gt.fasta -dna -mod oops -nmotifs 1 -maxw 12 -o gt_3

What do each of these parameters mean? Take a look at the output file. Note the “information content” information, as well as probability scores (p-values) for the instances.

$ less gt_3/meme.txt

Finally, by using the file transfer icon, let’s look at the meme.html output file and associated LOGO image generated by MEME. How does this compare to the actual motif in JASPAR?

Bonus question: What happens when you add the “-revcomp” flag to MEME? What is the difference in the output and why do you think this is?

$ meme gt.fasta -dna -revcomp -mod oops -nmotifs 1 -maxw 12 -o gt_4

Chapter 3: Sequence Alignments

Biological sequences evolve through a process of mutation and natural selection. By comparing two sequences, we can determine whether two sequences have a common evolutionary origin if their similarity is unlikely to be due to chance. Before we get into how this is done, we must also consider that there are many types of evolutionary relationships among sequences.

First, there is similarity, which fits the intuitive meaning of the degree of resemblance between two sequences. We might use the term identity to refere more exact situations, such the state of possessing the same subsequence. One often quantifies the percent identity between two sequences. The term homology refers to the state of sharing a common evolutionary origin. We say two sequences are homologous if they have a common ancestor. There are two types of homology. First, orthology refers to the state of being homologous sequences that arose from a common ancestral gene during speciation. Second, paralogy refers to the state of being homologous sequences that arose from a common ancestral gene from gene duplication.

Sequence alignment is the process of arranging the characters of a pair of sequences such that the number of matched characters is maximized. We can describe the alignment between two sequences with the following notation:

GCGTAACACGTGCG--
| ||| ||||||
AC--AACCCGTGCGAC

The vertical bars , or pipes, represent matching characters. Gaps, indicated by the dash are inserted in between characters in place of missing characters to optimize the number of matches. It is critical that sequence alignments are viewed in a monospace font, such as Courier, so that the width of characters don’t offset the alignment.

3.1 Alignment Algorithms and Dynamic Programming

One of the first attempts to align two sequences was carried out by Vladimir Levenstein in 1965, called “edit distance”, and now is often called Levenshtein Distance. The edit distance is defined as the number of single character edits necessary” to change one word to another. Initially, he described written texts and words, but this method was later applied to biological sequences. One of the most commonly used algorithms for computing the edit distance is the Wagner-Fischer algorithm, a Dynamic Programming algorithm.

Dynamic Programming optimally phrases the full problem as the optimal solution to the smaller pieces (sub-problems). The overall problem can then be expressed as a composition of the sub-problems. In addition to the Wagner-Fischer algorithm, numerous other dynamic programming algorithms have been developed for aligning biological sequences including the Needleman-Wunsch [22] and Smith-Waterman Algorithms [23].

3.1.1 Needleman-Wunsch Algorithm

The Needleman-Wunsch Algorithm is a global alignment algorithm, meaning the result always aligns the entire input sequences [22]. Later on in section 8.1 we will define a scoring matrix for protein alignment, but for nucleotide sequences, we often use a simpler scoring matrix such as

(3.1)

In addition to a scoring matrix, we also need to define penalties for gaps. The most common gap penalty is the linear gap penalty, defined as

which is just proportional to the length of the gap by a parameter . A more complicated approach is an “affine gap penalty”, which penalizes opening a gap by one parameter, and extending the gap by another parameter. For example, such a gap penalty can by defined by

which includes a gap open parameter and a gap extension parameter . In practice, an affine gap penalty is much more difficult to compute.

Dynamic programming for sequence alignments begins by defining a matrix or a table, to compute the scores. For example, let’s consider aligning the nucleotide sequences and . For Needleman-Wunsch, let’s define a matrix , such that the terms correspond to the score of aligning the subsequences and . We proceed from the upper left of this matrix at , and fill in the matrix as we move from left to right and from top to bottom. Here the rows of will correspond to the positions of , and the columns will correspond to the positions of .

When computing the terms of the matrix , we need to define a set of boundary conditions, namely that the score at the boundaries is associated with the penalty all the way up to that position. This is achieved by setting and for and .

We compute the terms of the matrix using a “recurrence relation”, such that the terms of a given cell of the matrix are defined in terms of the neighboring cells. Needleman-Wunsch uses the following recurrence relation:

(3.2)

Let’s consider the result of computing the matrix using the scoring matrix in 3.1, and using a linear gap penalty . The result is presented in Table 3.1.1. In this matrix, each term then corresponds to the score up to the character at that and position of the sequences and respectively. The rows will correspond to positions in the sequence , and the columns will correspond to positions of .

Figure 3.1: Each term of the matrix is computed using Equation 3.2 from the term above, to the left, and diagonally above-left.

The terms of the matrix can be filled out as is done with the following matrix, with each cell computed using the recursion relation in Equation 3.2, as depicted in Figure 3.1. The optimal path is shown in blue.

Table 3.1: A dynamic programming matrix to compute the score for Needleman-Wunsch Alignment.

–

C

C

A

T

A

C

G

A

–

0

-1

-2

-3

-4

-5

-6

-7

-8

C

-1

1

0

-1

-2

-3

-4

-5

-6

A

-2

0

0

1

0

-1

-2

-3

-4

G

-3

-1

-1

0

0

-1

-2

-1

-2

C

-4

-2

0

-1

-1

-1

0

-1

-2

T

-5

-3

-1

-1

0

-1

-1

-1

-2

A

-6

-4

-2

0

-1

1

0

-1

0

G

-7

-5

-3

-1

-1

0

0

1

0

C

-8

-6

-4

-2

-2

-1

1

0

0

G

-9

-7

-5

-3

-3

-2

0

2

1

In addition to the matrix, it is common to keep track of a traceback matrix , that keeps track of from where each term was computed from, in other words the maximum term in Eq 3.2. Table 3.1.1 demonstrates such a traceback matrix. One key complication is dealing with ties. One strategy is to favor adjacent matched characters as much as possible; therefore, we would favor diagonal terms before above or to the left.

Table 3.2: A traceback matriculates for Needleman-Wunsch.

–

C

C

A

T

A

C

G

A

–

·

·

·

·

·

·

·

·

·

C

·

×

↖

←

←

←

↖

←

←

A

·

↑

↖

↖

←

↖

←

←

↖

G

·

↑

↖

↑

↖

↖

↖

↖

←

C

·

↖

↖

←

↖

↖

↖

←

↖

T

·

↑

↑

↖

↖

←

↑

↖

↖

A

·

↑

↑

↖

←

↖

←

←

↖

G

·

↑

↑

↑

↖

↑

↖

↖

←

C

·

↖

↖

↑

↖

↑

↖

←

↖

G

↑

↑

↑

↖

↑

↑

↖

←

Finally, we have the result of the alignment. Here is the result of the Needleman-Wunsch alignment. Because it is a global alignment, the full sequence is included and the alignment ends on the first and last positions. There are, however, gaps at the first and last positions as this example illustrates.

-CAGCTAGCG-
|| || ||
CCA--TA-CGA

3.1.2 Smith-Waterman

n most applications we are only interested in aligning a small portion of the sequence to produce a local alignment. Furthermore, we don’t necessarily want to force the first and last residues to be aligned. Smith-Waterman is an alignment algorithm that has these properties [23].
We can define a set of boundary conditions for the scoring matrix , namely that the score is at the boundaries so that for and . Define the recurrence relation:

(3.3)

In addition to the different boundary conditions, a key difference between Needleman-Wunsch (global alignment) and Smith-Waterman (local alignment) is that whereas with the global alignment we start tracing back from the lower right term of the matrix, for the local alignment we start at the maximum value. This value corresponds to the last matched character of the optimal alignment.

Table 3.3: A matrix for Smith-Waterman, with an optimal path labeled in blue.

–

C

C

A

T

A

C

G

A

–

0

0

0

0

0

0

0

0

0

C

0

1

1

0

0

0

1

0

0

A

0

0

0

2

1

1

0

0

1

G

0

0

0

1

1

0

0

1

0

C

0

1

1

0

0

0

1

0

0

T

0

0

0

0

1

0

0

0

0

A

0

0

0

1

0

2

1

0

1

G

0

0

0

0

0

1

1

2

1

C

0

1

1

0

0

0

2

1

1

G

0

0

0

0

0

0

1

3

2

The optimal score corresponds to the in the last row, but second to last column. The optimal path results in an alignment with four matching positions. The traceback matrix can be built while computing the alignment matrix, and all paths are halted when a score of zero is reached.

Table 3.4: A traceback matrix for Smith-Waterman

–

C

C

A

T

A

C

G

A

–

·

·

·

·

·

·

·

·

·

C

·

↖

↖

←

·

·

↖

←

·

A

·

↑

↖

↖

←

↖

↑

↖

↖

G

·

·

·

↑

↖

↖

↖

↖

↑

C

·

↖

↖

↑

↖

↖

↖

↑

↖

T

·

↑

↖

↖

×

←

↑

↖

·

A

·

·

·

↖

↑

↖

←

←

↖

G

·

·

·

↑

↖

↑

↖

↖

←

C

·

↖

↖

←

·

↑

↖

↑

↖

G

·

↑

↖

↖

·

·

↑

↖

←

For Smith-Waterman, we typically report just the sub-alignment corresponding to the positive scores. We can report an alignment consisting of just the two sequences.

TAGCG
|| ||
TA-CG

3.1.3 Comparison

Although there are some similarities, there are a couple of key differences between Needleman-Wunsch and Smith-Waterman Algorithms. Here is a summary:

Needleman-Wunsch Algorithm

1. Computes the optimal global alignment in O(nm)

2. Backtracking begins in lower right: global adjustment

3. Allows negative scores

Smith-Waterman Algorithm

1. Computes optimal local alignment in O(nm)

2. Backtracking begins at largest value (not necessarily lower right)

3. Negative scores are zeroed out

3.1.4 Aligning DNA vs Proteins

When performing sequence alignments it is important to realize some of the key differences between aligning nucleic acid sequences and aligning protein sequences. We’ve seen that proteins can have substitution matrices, such as BLOSUM and PAM, that incorporate probabilistic models. That said, in Chapter 4 we will get into some probabilistic models of nucleotide substitution that could be incorporated into a scoring system. By building substitution matrices from curated alignments that record evolutionary changes that occur in nature, the protein substitution matrices encode the chemical similarity between amino acids. For example, scores are better for substituting between two polar amino acids compared to mutating from polar to non-polar. Furthermore, when inside the coding region of a gene, the third position of codons is more mutable because this position can typically change without changing the amino acid that it encodes.

3.2 Alignment Software

3.2.1 BLAST: Basic Local Alignment Search Tool

The BLAST algorithm (Basic Local Alignment Search Tool) developed by Altschul (1990) combines indexing of a database of sequences, and heuristics to approximate Smith-Waterman alignment, but is faster. The approach of BLAST is to index a search database using -mers, subsequences of length , for each of the sequences in the database. A query sequence is input to the program to search for similar sequences in the database. After low-complexity sequences are removed, all -mers of the query sequence are listed, and possible matches in the database are identified that would have an alignment score as good as , a predefined score threshold. The matching -mers are extended into stretches of matching -mers, that are called High-scoring Segment Pairs (HSPs), resulting in matches that are longer than . Two or more of these HSPs are combined to form a longer alignment. Ultimately Smith-Waterman alignment is performed on just these strongly matching sequences, and this is what is reported.
In summary, the approach is as follows:

1. Remove low0complexity regions or sequence repeats in the query sentence.

2. Make -mer word list of the query sequence (Proteins often  = 3)

3. List the possible matching words with a scoring matrix

4. Reduce the list of word matches with threshold T

5. Extend the exact matches to High-scoring Segment Pairs (HSPs)

6. List all HSPs and evaluate significance

7. Combine two or more HSPs into a longer alignment

8. Report the gapped Smith-Waterman local alignments of the query and each of the matched database sequences.

Figure 3.2: Make -mer word list of the query sequence (Proteins often )

3.3 Alignment Statistics

When evaluating a BLAST score, it is important to have a statistical framework for evaluating the significance of a “BLAST hit”. Here we present such a system where we consider our score as a random variable. Because BLAST identifies the maximum scoring alignment, we can describe the cumulative distribution of BLAST scores with the Generalized Extreme Value (GEV) distribution:

The parameter is the location parameter of the GEV, and is expressed here in terms of the length of the query sequence, and the length of the entire database. here is a constant that is particular to and computed from a particular database.

Figure 3.3: -mers that match with a score above are extended form High-scoring Segment Pairs (HSPs).

The p-value is the probability of a score greater than or equal to S due to chance, and is given by:

This equation comes from the Poisson distribution. If we define E-value (expected number of hits at this score or greater due to chance) as:

The p-value can then be simplified as:

After a linear transformation, the score S’ can be computed in terms of bits.

The updated equation for E-value is much simpler:

Figure 3.4: Extend the exact matches to High-scoring Segment Pairs (HSPs)

3.3.1 Running BLAST from the command line

BLAST can be run on the command line pretty easily. To do this, you need a sequence, or set of sequences to align, and a database to align to. First, let’s create the database to align to. This can be created using a FASTA file of sequences. For example, if we have the FASTA file for the human genome , we can format the database with using the following command:

$ makeblastdb -in hg38.fa -input_type fasta -title hg38 -dbtype nucl

In this command, most of the terms make sense. The input file is , the input file format is “fasta”, and the title used is . The last term specifies that the input data is nucleic acid sequences.
On some systems, an older version is installed using . With this program, the database can be created using this command:

$ formatdb -p F -t hg38 -n hg38 -i hg38.fa

In this command, the command indicates that this is a nucleotide sequence, and not a protein sequence. Specifically, the specifies protein, and the says that this is “false”, specifying that the input data is not protein. The second two commands give the database the title and name . The name specified by the command provides a basename for the output files used in the database, and also gives a label to be used when referring to the database in BLAST. Finally, the command specifies the input file, which is the FASTA file for the genome.

Next, we can run BLAST using the command . This tool allows you to run different versions of BLAST, specified by the command. To run a nucleotide query against a nucleotide database, we use . The full command is as follows:

$ blastall -p blastn -i sequences.fa -d hg38 -o sequences_hg38_blast.txt

Here we specify the input sequences, the query, with the command. Then we specify the database that we are aligning to, using the flag, referring to the database that we just created with . Finally, we specify an output file to write the results to, using the flag.

3.4 Short Read Mapping

The growth of high-throughput sequencing has led to a parallel growth of software applications for rapidly aligning short reads. Although BLAST was designed for fast alignment, these new tools are even faster for the alignment of short sequence reads. We will discuss these methods further in Chapter 9.

3.5 Lab 4: Using BLAST on the command line

In this lab, we will learn how to run BLAST on the command line. As usual, you should create and enter a directory.

3.5.1 Part 1: BLASTing to a protein database

Let’s first build a database. Please download the Swissprot database from NCBI with the following command:

$ wget ftp://ftp.ncbi.nih.gov/blast/db/FASTA/swissprot.gz

and then unzip the downloaded file with the following command:

$ gunzip swissprot.gz

Although there is no file extension, the file is a FASTA file. Let’s rename it so that we know it is a FASTA file.

$ mv swissprot swissprot.fa

Next, let’s build a database with the following command:

$ makeblastdb -in swissprot.fa -input_type fasta -title swissprot -dbtype prot

Not all of these options are required. Can you figure out which options are required by the help message printed with you run this command?

$ makeblastdb -help

Download the protein sequence infomation for human BRCA1 and create a fasta file for the sequence (https://www.ncbi.nlm.nih.gov/protein/1698399?report=fasta). Save it to a file called . Copy the sequence, and paste it into a file after opening it with nano:

$ nano brca1_pep.fasta

To save with nano, type Ctrl-X, then type Y. Next, we can BLAST the brca1 pep.fasta file we created.

$ blastp -query brca1_pep.fasta -db swissprot.fa > brca1_swissprot

Do the top hits make sense to you? You can search NCBI Protein for some of the IDs.

$ less brca1_swissprot

What are the best hits? Do the order of the sequence hits make sense in terms of what you know of the biology?

You can also BLAST the sequence to the “non-redundant” database “nr” by pasting it to the NCBI BLAST web tool: https://blast.ncbi.nlm.nih.gov/Blast.cg. Note that you could do theoretically do this by specifying “nr” for the database, but many servers don’t have this downloaded (it’s a very big file!).

3.5.2 Biopython and BLAST (optional)

You could also analyze your blast hits using Biopython. To do this, you need to set the output format to XML with the following command.

$ blastp -query brca1_pep.fasta -db swissprot -outfmt 5 > brca1_swissprot.xml

The XML can be difficult to read, but can be parsed easily. For example, you can print the alignment for each BLAST hit in the results with something like this:

>>> from Bio.Blast import NCBIXML
>>> result_handle = open("brca1_swissprot.xml")
>>> blast_record = NCBIXML.read(result_handle)
>>> for alignment in blast_record.alignments:
...     for hsp in alignment.hsps:
...         if hsp.expect < 1e-10:
...             print('sequence:', alignment.title)
...             print('length:', alignment.length)
...             print('e value:', hsp.expect)
...             print(hsp.query)
...             print(hsp.match)
...             print(hsp.sbjct)

Variations on this method could allow one to parse the BLAST output file, and extract the alignments as well.

3.5.3 Part 2: BLASTing to a genome

UCSC provides a wealth of genomic resources. You can download the Drosphila genome version dm3 at this link: http://hgdownload.soe.ucsc.edu/goldenPath/dm3/bigZips/. Download with the command at the terminal:

$ wget http://hgdownload.soe.ucsc.edu/goldenPath/dm3/bigZips/chromFa.tar.gz

Unzip the file with the command:

$ tar xvfz chromFa.tar.gz

Combine all the chromosome FASTA files into one genome file:

$ cat chr*fa > dm3.fa

In this case, the asterisk is used as a wild-card, that specifies all files with anything between a and a . You can clean up the directory by moving the files into a directory. First create the directory:

$ mkdir genome

Move the chromosome files into the directory with this command:

$ mv chr*fa genome/.

Build a blast database:

$ makeblastdb -in dm3.fa -title dm3 -dbtype nucl

Download the transcript sequence for human BRCA1 and create a FASTA file for the sequence NCBI human BRCA1 here: https://www.ncbi.nlm.nih.gov/nuccore/1147602?report=fasta.
BLAST The sequence to the genome:

$ blastn -query brca1.fa -db dm3.fa > brca1_dm3.blast

Download the RefSeq mRNA annotations with the command at the terminal:

$ wget http://hgdownload.soe.ucsc.edu/goldenPath/dm3/bigZips/refMrna.fa.gz

gunzip the file with the command:

$ gunzip refMrna.fa.gz

Create a database for the RefSeq annotations:

$ makeblastdb -in refMrna.fa -title refMrna -dbtype nucl

BLAST The sequence to the refMrna database:

$ blastn -query brca1.fa -db refMrna.fa > brca1_refMrna.blast

Discussion questions: the difference between the two results? How do you explain the difference? What is chrUextra anyway? To answer this you should look at the BLAST output with “less” in the same way you looked at other BLAST output above.

Going beyond: How does this e-value compare to BLASTing to mouse through the NCBI website? Can you find a gene in human that has a significant hit to the E. coli genome?

Chapter 4: Multiple Sequence Alignments, Molecular Evolution, and Phylogenetics

4.1 Multiple Sequence Alignment

A multiple sequence alignment is an alignment of more than 2 sequences. It turns out that this makes the problem of alignment much more complicated, and much more computationally expensive. Dynamic programming algorithm such as Smith-Waterman can be extended to higher dimensions, but at a significant computing cost. Therefore, numerous methods have been developed to make this task faster.

4.1.1 MSA Methods

There have been numerous methods developed for computing MSAs do make them more computationally feasible.

Dynamic Programming

Despite the computational cost of MSA by Dynamic Programming, there have been approaches to compute multiple sequence alignments using these approaches [24,25]. The programs MSA [25] and MULTALIN [26] use dynamic programming. This process takes computations for aligning sequences of length . The Carrillo-Lipman Algorithm uses pairwise alignments to constrain the search space. By only considering regions of the multi-sequence alignment space that are within a score threshold for each pair of sequences, the search space can be reduced.

Progressive alignments

Progressive alignments begin by performing pairwise alignments, by aligning each pair of sequences. Then it combines each pair and integrates them into a multiple sequence alignment. The different methods differ in their strategy to combine them into an overall multiple sequence alignment. Most of these methods are “greedy”, in that they combine the most similar pairs first, and proceed by fitting the less similar pairs into the MSA. Some programs that could be considered progressive alignment include T-Coffee [27], ClustalW and its variants [28], and PSAlign [29].

Iterative Alignment

Iterative Alignment is another approach that improves upon the progressive alignment because it starts with a progressive alignment and then iterates to incrementally improve the alignment with each iteration. Some programs that could be considered iterative alignment include CHAOS/DIALIGN [30], and MUSCLE [31].

Full Genome Alignments

Specialized multiple sequence alignment approaches have been developed for aligning complete genomes, to overcome the challenges associated with aligning such long sequences. Some programs that align full genomes include MLAGAN (using LAGAN) [32], MULTIZ (using BLASTZ) [33], LASTZ [34], MUSCLE [31], and MUMmer [35, 36].

4.1.2 MSA File Formats

There are several file formats that are specifically designed for multiple sequence alignment. These approaches can differ in their readability by a human, or are designed for storing large sequence alignments.

Multi-FASTA Format

Probably the simplest multiple sequence alignment format is the Multi-FASTA format (MFA), which is essentially like a FASTA file, such that each sequence provides the alignment sequence (with gaps) for a given species. The deflines can in some cases only contain information about the species, and the file name, for example, could contain information about what sequence is being described by the file. For short sequences the mfa can be human readable, but for very long sequences it can become difficult to read. Here is an example file that shows the alignment of a small (28aa) Drosophila melanogaster peptide called Sarcolamban (isoform C) with its best hits to .

D.melanogaster
--------------------------------------------------
--------------------------------MSEARNLFTTFGILAILL
FFLYLIYA---------------------VL---------------
>D.sechellia
--------------------------------------------------
--------------------------------MSEARNLFTTFGILAILL
FFLYLIYAPAAKSESIKMNEAKSLFTTFLILAFLLFLLYAFYEAAF
>D.pseudoobscura
MSEAKNLMTTFGILAFLLFCLYLIYASNNSKRWPTFCGEAEFRSENSESQ
LLRAFSYERLEQCPNKKYPPKQPTTTTTKPIKMNEARSLFTTFLILAFLL
FLLYAFYEA--------------------AF---------------
>D.busckii
--------------------------------------------------
--------------------------------MNEAKSLVTTFLILAFLL
FLLYAFYEA--------------------AF---------------

Clustal

The Clustal format was developed for the program [37], but has been widely used by many other programs [28, 38]. This file format is intended to be fairly human readable in that it expresses only a fixed length of the alignment in each section, or block. Here is what the Clustal format looks like for the same Sarcolamban example:

CLUSTAL W (1.83) multiple sequence alignment

D.melanogaster    ------------------------------------------------------------
D.sechellia       ------------------------------------------------------------
D.pseudoobscura   MSEAKNLMTTFGILAFLLFCLYLIYASNNSKRWPTFCGEAEFRSENSESQLLRAFSYERL
D.busckii         ------------------------------------------------------------

D.melanogaster    ----------------------MSEARNLFTTFGILAILLFFLYLIYA------------
D.sechellia       ----------------------MSEARNLFTTFGILAILLFFLYLIYAPAAKSESIKMNE
D.pseudoobscura   EQCPNKKYPPKQPTTTTTKPIKMNEARSLFTTFLILAFLLFLLYAFYEA-----------
D.busckii         ----------------------MNEAKSLVTTFLILAFLLFLLYAFYEA-----------
*.**:.*.*** ***:***:** :*

D.melanogaster    ---------VL---------------
D.sechellia       AKSLFTTFLILAFLLFLLYAFYEAAF
D.pseudoobscura   ---------AF---------------
D.busckii         ---------AF---------------
:

The clustal format has the following requirements, which can make it difficult to create one manually. First, the first line in the file must start with the words “” or ““. Other information in the first line is ignored, but can contain information about the version of that was used to create it. Next, there must be one or more empty lines before the actual sequence data begins. The rest of the file consists of one or more blocks of sequence data. Each block consists of one line for each sequence in the alignment. Each line consists of the sequence name, defline, or identifier, some amount white space, then up to 60 sequence symbols such as characters or gaps. Optionally, the line can be followed by white space followed by a cumulative count of residues for the sequences. The amount of white space between the identifier and the sequences is usually chosen so that the sequence data is aligned within the sequence block. After the sequence lines, there can be a line showing the degree of conservation for the columns of the alignment in this block. Finally, all this can be followed by some amount of empty lines.

MAF

The Multiple Alignment Format (MAF) can be a useful format for storing multiple sequence alignment information. It is often used to store full-genome alignments at the UCSC Genome Bioinformatics site. The file begins with a header beginning with and information about the version and scoring system. The rest of the file consists of alignment blocks. Alignment blocks start with a line that begins with the letter and a score for the alignment block. Each subsequent line begins with either an , a , or an indicating what kind of line it is. The lines beginning with contain sequence information. Lines that begin with typically follow each -line, and contain information about what is occurring before and after the sequences in this alignment block for the species considered in the line. Lines beginning with contain information about empty parts of the alignment block, for species that do not have sequences aligning to this block. For example, the following is a portion of the alignment of the Human Genome (GRCh38/hg38) with 99 vertebrates.

##maf version=1 scoring=roast.v3.3
a score=49441.000000
s hg38.chr22 10514742 28 + 50818468 acagaatggattattggaacagaataga
s panTro4.chrUn_GL393523 96163 28 + 405060 agacaatggattagtggaacagaagaga
i panTro4.chrUn_GL393523 C 0 C 0
s ponAbe2.chrUn 66608224 28 - 72422247 aaagaatggattagtggaacagaataga
i ponAbe2.chrUn C 0 C 0
s nomLeu3.chr6 67506008 28 - 121039945 acagaatagattagtggaacagaataga
i nomLeu3.chr6 C 0 C 0
s rheMac3.chr7 24251349 14 + 170124641 --------------tggaacagaataga
i rheMac3.chr7 C 0 C 0
s macFas5.chr7 24018429 14 + 171882078 --------------tggaacagaataga
i macFas5.chr7 C 0 C 0
s chlSab2.chr26 21952261 14 - 58131712 --------------tggaacagaataga
i chlSab2.chr26 C 0 C 0
s calJac3.chr10 24187336 28 + 132174527 acagaatagaccagtggatcagaataga
i calJac3.chr10 C 0 C 0
s saiBol1.JH378136 10582894 28 - 21366645 acataatagactagtggatcagaataga
i saiBol1.JH378136 C 0 C 0
s eptFus1.JH977629 13032669 12 + 23049436 ----------------gaacaaagcaga
i eptFus1.JH977629 C 0 C 0
e odoRosDiv1.KB229735 169922 2861 + 556676 I
e felCat8.chrB3 91175386 3552 - 148068395 I
e otoGar3.GL873530 132194 0 + 36342412 C
e speTri2.JH393281 9424515 97 + 41493964 I
e myoLuc2.GL429790 1333875 0 - 11218282 C
e myoDav1.KB110799 133834 0 + 1195772 C
e pteAle1.KB031042 11269154 1770 - 35143243 I
e musFur1.GL896926 13230044 2877 + 15480060 I
e canFam3.chr30 13413941 3281 + 40214260 I
e cerSim1.JH767728 28819459 183 + 61284144 I
e equCab2.chr1 43185635 316 - 185838109 I
e orcOrc1.KB316861 20719851 245 - 22150888 I
e camFer1.KB017752 865624 507 + 1978457 I

The lines contain 5 fields after the at the beginning of the line. First, the source of the column usually consists of a genome assembly version, and chromosome name separated by a dot “.”. Next is the start position of the sequence in that assembly/chromosome. This is followed by the size of the sequence from the species, which may of course vary from species to species. The next field is a strand, with “+” or “-“, indicating what strand from the species’ chromosome the sequence was taken from. The next field is the size of the source, which is typically the length of the chromosome in basepairs from which the sequence was extracted. Lastly, the sequence itself is included in the alignment block.

4.2 Phylogenetic Trees

A phylogenetic tree is a representation of the evolutionary history of a character or sequence. Branching points on the tree typically represent gene duplication events or speciation events. We try to infer the evolutionary history of a sequence by computing an optimal phylogenetic tree that is consistent with the extant sequences or species that we observe.

4.2.1 Representing a Phylogenetic Tree

A phylogenetic tree is often a “binary tree” where each branch point goes from one to two branches. The junction points where the branching takes place are called “internal nodes”. One way of representing a tree is with nested parentheses corresponding to branching. Consider the following example

((A,B),(C,D));

where the two characters and are grouped together, and the characters and are grouped. The semi-colon at the end is needed to make this tree in proper “newick” tree format. One of the fastest ways to draw a tree on the command line is the “ASCII tree”, which we can draw by using the function . To use this, we’ll need to save our tree to a text file that we can read in. We could read the tree as just text into the python terminal (creating a string), but that would require loading an additional module to use the function . Therefore, it might be just as easy to save it to a text file called that we can read in. Putting this together, we can draw the tree with the following commands:

>>> from Bio import Phylo
>>> tree = Phylo.read("tree.txt","newick")
>>> Phylo.draw_ascii(tree)
___________________________________ A
___________________________________|
|                                   |___________________________________ B
_|
|                                    ___________________________________ C
|___________________________________|
|___________________________________ D

This particular tree has all the characters at the same level, and does not include any distance or “branch length” information. Using real biological sequences, we can compute the distances along each brach to get a more informative tree. For example, we can download 18S rRNA sequences from NCBI Nucleotide. Using clustalw, we can compute a multiple sequence alignment, and produce a phylogenetic tree. In this case, the command

$ clustalw -infile=18S_rRNA.fa -type=DNA -outfile=18S_rRNA.aln

will produce the output file , which is a tree in newick tree format. The file contains the following information

(fly:0.16718,(mouse:0.00452,human:0.00351):0.05186,chicken:0.24654);

You’ll note, this is the same format as the simple example above, but rather than a simple label for each character/sequence, there is a label and a numerical value, corresponding to the brach length, separated by a colon. In addition, each of the parentheses are followed by a colon and a numerical value. In each case, the value of the branch length corresponds to the substitutions per site required to change one sequence to another, a common unit of distance used in phylogenetic trees. This value also corresponds to the length of the branch when drawing the tree. The command to draw a tree image is simply , which will allow the user to save the image.

>>> from Bio import Phylo
>>> tree = Phylo.read('18S_rRNA.dnd','newick')
>>> Phylo.draw(tree)

The resulting image can be seen in Figure 4.1, and visually demonstrates the branch lengths corresponding to the distance between individual sequences. The x-axis in the representation corresponds to this distance, but the y-axis only separates taxa, and the distance along the y-axis does not add to the evolutionary distance between sequences.

Figure 4.1: A phylogenetic tree computed with Clustal for 18S rRNA sequences for Drosophila melanogaster, Homo sapiens, Mus musculus, and Gallus gallus.

4.2.2 Pairwise Distances

Phylogenetic trees are often computed as binary trees with branch lengths that optimally match the pair-wise distances between species. In order to compute a phylogenetic tree, we need a way of defining this distance. One strategy developed by Feng and Doolittle is to compute a distance from the alignment scores computed from pair-wise alignments. The distance is defined in terms of an “effective”, or normalized, alignment score for each pair of species. This distance is defined as

so that pairs of sequences that have high scores will have a small distance between them. The effective score is defined as

Where in this expression, is the observed pairwise similarity between sequences from species and . The value is the average of the two scores when you align species and to themselves, which represents the score corresponding to aligning “identical” sequences, the maximum possible score one could get. is the average pairwise similarity between randomized, or shuffled, versions of the sequences from species and . After this normalization, the score ranges from 0 to 100.

4.3 Models of mutations

Evolution is a multi-faceted process. There are many forces involved in molecular evolution. The process of mutation is a major force in evolution. Mutation can happen when mistakes are made in DNA replication, just because the DNA replication machinery isn’t 100% perfect. Other sources of mutation are exposure to radiation, such as UV radiation, certain chemicals can induce mutations, and viruses can induce mutations of the DNA (as well as insert genetic material). Recombination is a genetic exchange between chromosomes or regions within a chromosome. During meiosis, genes and genomic DNA are shuffled between parent chromosomes. Genetic drift is a stochastic process of changing allele frequencies over time due to random sampling of organisms. Finally, natural selection is the process where differences in phenotype can affect survival and reproduction rates of different individuals.

4.3.1 Genetic Drift

Because genetic drift is a stochastic process, it can be modeled as a “rate”. The rate of nucleotide substitutions can be expressed as

where is the substitutions per site across the genome, and is the time of divergence of the two (extant) species to their common ancestor. The factor of two can be understood by the fact that it takes a total of to go from to , stopping at along the way. The mutations that occur over time separating and can be viewed as distributed over a time .

Figure 4.2: For two extant species and diverged for a time from a common ancestor , the mutation rate can be expressed as a time separating and .

In most organisms, the rate is observed to be about to mutations per generation. Some viruses have higher mutation rates mutations per generation. The generation times of different species can also affect the nucleotide substitution rate . Organisms with shorter generation times have more opportunities for meiosis per unit time.

When mutation rates are evaluated within a gene, positional dependence in nucleotide evolution is observed. Because of the degeneracy in the genetic code, the third position of most codons has a higher substitution rate. Some regions of proteins are conserved domains, hence the corresponding regions of the gene have lower mutation rates compared to other parts of the gene. Other genes such as immunoglobulins have very high mutation rates and are considered to be “hypervariable”.

Noncoding RNAs have functional constraints to preserver hairpins, and may have sequence evolution that preserves base pairing through compensatory changes on the paired nucleotide. The result is that many noncoding RNAs such as tRNAs have very conserved structure, but vary at the sequence level.

4.3.2 Substitution Models

The branches of phylogenetic trees can often represent the the expected number of substitutions per site. That is, the distance along the branches of the phylogenetic tree from the ancestor to the extant species correspond to the expected number of substitutions per site in the time it takes to evolve from the ancestor to the extant species.

A substitution model describes the process of substitution from one set of characters to another through mutation. These models are often neutral, in the sense that selection is not considered, and the characters mutate in an unconstrained way. Furthermore, these models are typically considered to be independent from position to position.

Substitution models typically are described by a rate matrix with terms that describe mutating from character to for terms where . The diagonal terms of the matrix are defined so that the sum of the rows are zero, so that

(4.1)

In general, the matrix can be defined as:

(4.2)

where the diagonal terms are defined such that it is consistent with Equation 4.1. The rate matrix is associated with a probability matrix , which describes the probability of observing the mutation from to in a time by the terms . We want these probabilities to be multiplicative, meaning that . The mutations associated with amounts of time and applied successively can be understood as the same as the mutations associated with . Furthermore, the derivative of the equation can be expressed as

(4.3)

The solution to this equation is the exponential function. The rate matrix itself can be exponentiated to compute the probability of a particular mutation in an amount of time , which can be computed using the Taylor series for the exponential function.

(4.4)

Each such model also assumes an equilibrium frequencies , which describe the probability of each nucleotide after the system has reached equilibrium.

4.3.3 Jukes-Cantor 1969 (JC69)

The simplest substitution model was proposed by Jukes and Cantor (JC69), and describes equal rates of evolution between all nucleotides [39]. The JC69 model defines a constant mutation rate , and equilibrium frequencies such that . The equilibrium frequencies describe the frequencies of each nucleotide that result after the system has evolved under this model for a “very long time”. The rate matrix for the Jukes-Cantor model is then given by:

(4.5)

It can be shown that the full expression for computing is:

(4.6)

Therefore, we can solve this system to get the probabilities , which can be expressed as

(4.7)

Finally, we note that the sum of the terms of a row of the matrix that correspond to mutation changes give us the expected value of the distance in substitutions per site. For the Jukes-Cantor model, this corresponds to . Substituting this into Equation 4.7 for the terms involving change (), gives , which can be solved for to give:

}

(4.8)

This formula is often called the Jukes-Cantor distance formula, and it gives a way to relate the proportion of sites that differ to the evolutionary distance , which stands for the expected number of substitutions in the time for a mutation rate . This formula corrects for the fact that the proportion of sites that differ, , does not take into account sites that mutate, and then mutate back to the original character.

4.3.4 Kimura 1980 model (K80)

The Jukes-Cantor model considers all mutations to be equally likely. The Kimura model (K80) considers the fact that transversions, mutations involving purine to pyrimidines or vice versa, are less likely than transitions, which are from purines to purines or pyrimidines to pyrimidines [40]. Therefore, this model has two parameters: a rate for transitions, and a rate for transversions.

(4.9)

Applying a similar derivation as was done for the JC69 model, we get the following distance formula for K80:

(4.10)

where is the proportion of sites that show a transition, and is the proportion of sites that show transversions.

4.3.5 Felsenstein 1981 model (F81)

The Felsenstein model essentially makes the assumption that the rate of mutation to a given nucleotide has a specific value equal to its equilibrium frequency , but these value vary from nucleotide to nucleotide [41]. The rate matrix is then defined as:

(4.11)

4.3.6 The Hasegawa, Kishino and Yano model (HKY85)

The Hasegawa, Kishino and Yano model takes the K80 and F81 models a step further and distinguishes between transversions and transitions [42]. In this expression, the transitions are weighted by an additional term .

(4.12)

4.3.7 Generalized Time-Reversible Model

This model has 6 rate parameters, and 4 frequencies hence, 9 free parameters (frequencies sum to 1) [43]. This method is the most detailed, but also requires the most number of parameters to describe.

4.3.8 Building Phylogenetic Trees

Now that we have a probabilistic framework with which to describe phylogenetic distances, we need some methods to build a tree from a set of pair-wise distances. Here are two basic approaches to building Phylogenetic Trees.

Unweighted Pair Group Method with Arithmetics Mean (UPGMA) Algorithm

UPGMA is a phylogenetic tree building algorithm that uses a type of hierarchical clustering [44]. This algorithm builds a rooted tree by creating internal nodes for each pair of taxa (or internal nodes), starting with the most similar and proceeding to the least similar. This approach starts with a distance matrix for each taxa and . When branches are built connecting and , an internal node is created, which corresponds to a cluster containing and . Distances are updated such that the distance between a cluster (internal node) and a leaf node is the average distance between all members of the cluster and the leaf node. Similarly, the distance between clusters is the average distance between members of the cluster.

Neighbor Joining Algorithm

One of the issues with UPGMA is the fact that it is a greedy algorithm, and joins the closest taxa first. There are tree structures where this fails. To get around this, the Neighbor joining algorithm normalizes the distances and computes a set of new distances that avoid this issue [45]. Using the original distance matrix , a new distance matrix is computed using the following formula:

Begin with a star tree, and a matrix of pair-wise distances d(i,j) between each pair of sequences/taxa, an updated distance matrix that normalizes compared to all distances is created. The updated distances are computed. The closest taxa using this distance are identified, and an internal node is created such that the distance along the branch connecting these two nearest taxa is the distance . This process s repeated using the new (internal) node as a taxa and the distances are updated.

4.3.9 Evaluating the Quality of a Phylogenetic Tree

Maximum Parsimony

Maximum Parsimony makes the assumption that the best phylogenetic tree is that with the shortest branch lengths possible, which corresponds to the fewest mutations to explain the observable characters [46 47].
This method begins by identifying the phylogenetically informative sites. These sites would have to have a character present (no gaps) for all taxa under consideration, and not be the same character for all taxa. Then trees are constructed, and characters (or sets of characters) are defined for each internal node all the way up to the root. Then each tree has a cost defined, corresponding to the total number of mutations that need to be assumed to explain that tree. The tree with the shortest total branch length is typically chosen. The length of the tree is defined as the sum of the lengths of each individual character (or column of the alignment) , and possibly using a weight for different characters (often just for all columns, but could weight certain positions more).

In this expression, the length of a character can be computed as the total number of mutations needed to explain this distribution of characters given the topology of the tree.

Maximum Likelihood

Maximum likelihood is an approach that computes a likelihood for a tree, using a probabilistic model for each tree [48 49]. The probabilistic model is applied to each branch of the tree, such as the Jukes-Cantor model, and the likelihood of a tree is the product of the probabilities of each branch. Generally speaking, we seek to find the tree that has the greatest likelihood given the .

These probabilities could be computed from an evolutionary model, such as Jukes-Cantor model. For example, if we have observed characters and , and an unknown ancestral character , and lengths of the branches from to be and respectively, we could compute the likelihood of this simple tree as

where we are summing over all possible ancestral characters , and computing the probability of mutating along the branches using a probabilistic model. This probabilistic model can be the same terms of the probability matrices discussed in Equation 4.6.

4.3.10 Tree Searching

In addition to computing the quality of a specific tree, we also want to find ways of searching the space of trees. Because the space of all possible trees is so large, we can’t exhaustively enumerate them all in a practical amount of time. Therefore, we need to sample different trees stochastically. Two such methods are Nearest Neighbor Interchange (NNI) and Subtree Pruning and Re-grafting (SPR).

Nearest Neighbor Interchange

Figure 4.3: Nearest Neighbor Interchange operates by swapping out the locations of subgraphs within a tree. For a tree with four sub-trees, there are only 2 possible interchanges.

For a tree toplogy that contains four or more taxa, the Nearest Neighbor Interchange (NNI) exchanges subtrees within the larger tree. It can be shown that for a tree with four subtrees there are only 3 distinct ways to exchange subtrees to create a new tree, including the original tree. Therefore, each such application produces two new trees that are different from the input tree [50 51].

Subtree Pruning and Regrafting

Figure 4.4: A depiction of SPR Tree searching method. A. One subtree of the larger tree structure is selected. B. An attachment point is selected. C. The subtree is then “grafted” or attached to the attachment point.

The Subtree Pruning and Regrafting (SPR) method takes a subtree from a larger tree, removes it and reattaches it to another part of the tree [52 53].

4.4 Lab 5: Phylogenetics

In this lab, we will learn some basic commands for computing phylogenetic trees, and some python commands that will draw the tree. Let’s create a new directory called to work in.

4.4.1 Download Sequences from NCBI

Download sequences in FASTA format for a gene of interest from NCBI nucleotide (https://www.ncbi.nlm.nih.gov/nuccore). Build a FASTA file containing each sequence and a defline. Shorten the defline for each species to make it easier to read later. As an example, I downloaded 18S rRNA for Human (Homo sapiens), Mouse (Mus musculus), Rat (Rattus norvegicus), Frog (Xenopus laevis), Chicken (Gallus gallus), Fly (Drosophila melanogaster) and Arabidopsis (Arabidopsis thalian). You can download these 18S rRNA sequences with the following command:

$ wget http://hendrixlab.cgrb.oregonstate.edu/teaching/18S_rRNAs.fasta

4.4.2 Create a Multiple Sequence Alignment and Phylogenetic Tree with Clustal

First, use to align the sequences, and output a multiple sequence alignment and dendrogram file.

$ clustalw2 -infile=18S_rRNAs.fasta -type=DNA -outfile=18S_rRNAs.aln

The dendogram file, indicated by the “” suffix, can be used to create an image of a phylogenetic tree using Biopython: First, enter the python terminal by typing “” and then the following:

>>> import matplotlib as mpl
>>> mpl.use('Agg')
>>> import matplotlib.pyplot as pyplot
>>> from Bio import Phylo
>>> tree = Phylo.read('18S_rRNAs.dnd','newick')
>>> Phylo.draw(tree)
>>> pyplot.savefig("myTree.png")
>>> quit()

How does the resulting tree compare to what you expect given these species?

4.4.3 Create a Multiple Sequence Alignment and Phylogenetic Tree with phyML

The program phyML provides much more flexibility in what sort of trees it can compute. To use it, we’ll need to convert our alignment file to phylip. We can do this with the Biopython module AlignIO

>>> from Bio import AlignIO
>>> alignment = AlignIO.parse("18S_rRNAs.aln","clustal")
>>> AlignIO.write(alignment,open("18S_rRNAs.phy","w"),"phylip")
1
>>> quit()

You can run phyml in the simplest way by simply typing “phyml” and then typing your alignment file:

$ phyml

Enter the sequence file name > 18S_rRNAs.phy

The program will give you a set of options, and you can optionally change them. To change the model, type “+” and then “M” to toggle through the models. Finally, hit “enter” to run the program.

>>> import matplotlib as mpl
>>> mpl.use('Agg')
>>> import matplotlib.pyplot as pyplot
>>> from Bio import Phylo
>>> tree = Phylo.read('18S_rRNAs.phy_phyml_tree.txt','newick')
>>> Phylo.draw(tree)
>>> pyplot.savefig("myTreeML.png")
>>> quit()

How does the tree created from compare to the tree created using ?

Chapter 5: Genomics

In this chapter, we will learn how to work with genomic data and genome annotations and associated file formats. It is hoped that this chapter will serve as a basic introduction to genomics, with the understanding that it is a far broader field than the types of genome annotations presented here. Many of the chapters, such as the next chapter “Transcriptomics”, and many before, could be considered as a subject matter as part of the field of genomics.

5.1 The Three Fundamental “Gotchas” of Genomics

Whereas many sciences have three or four fundamental laws that describe them, research in genomics often encounters exceptions to rules rather than universal principles. However, when learning genomics for the first time, researchers often encounter the same problems that cause the same strange things to happen in their data. These errors are so common for beginners in genomics, they can be thought of as the “fundamental gotchas” of genomics:

1. Different genome assemblies

2. Different chromosome deflines

3. 0 vs 1-based coordinates

5.1.1 Different Genome Assemblies/Annotations

The first fundamental gotcha involves the fact that for any model organism, there are typically many different genome assembly versions. For example, one might download two different datasets, such as a set of genomic positions, and try to compare them, only to realize that they were compared against different assemblies of the genome. In many cases, such as human genome assembly versions, each being incremental improvements in the accuracy of the assembly.

5.1.2 Different Chromosome Deflines

Depending on from where you downloaded the genome annotation, you could possibly encounter an assembly that uses one defline for “Chromosome 1”, and another that uses a different one. For example, some annotations such as from UCSC Genome Bioinformatics, would use “chr1”, and other databases may just use “1”. You definitely need to double check these values when comparing different formats.

5.1.3 0 vs 1-based coordinates

The third fundamental gotcha involves comparing data from different file formats, and not being aware that some file formats use -base positions, and others use -based positions. Whereas -based coordinates consider the first position of a chromosome to be position , -based coordinates consider the first position of a chromosome to be position . In many ways, is more natural for computer science, because typically programming languages use -based coordinates to describe strings. Biologists may be more accustomed to counting positions starting at . There are many different file formats for storing the positions of genomic locations, each use one of these two position systems, so it’s important to know which is which, and be aware that there is a difference.

Table 5.1: Different file formats use different position systems.

Format

Position system

GFF/GTF

1-based

BLAST results

1-based

BLAT results

1-based

maf

0-based

BAM

0-based

SAM

1-based

BED

0-based

wig

1-based

5.2 Genomic Data and File Formats

5.2.1 Formats for Genomic Locations

One of the most common form of genome annotation is that of the genome location. Often we want to annotate and label a binding site for a protein, or a gene by the genomic locations of the exons. Therefore, we need to define file formats that specify genomic locations.
Often, we can describe a genomic location by four variables that specify the “genomic coordinate”. The “genomic coordinate” can be specified by a chromosome, a start position, a stop position, and a strand. Clearly, the fundamental gotchas are important when evaluating these four values. For strand, it is most common to see “+” and “-” to refer to the forward and reverse strand. However, if you look at enough databases you will find things like “F” and “R” to denote strand.

BED Files

BED files are a simple file format to describe genomic locations. The BED file was developed by UCSC Genome Bioinformatics, and is a common format found when downloading data from that site. These descriptions are derived from the descriptions at UCSC (https://genome.ucsc.edu/FAQ/FAQformat.html). The first three required BED fields are:

1. chrom: The name of the sequence (e.g. chromosome or scaffold) that the feature is within.

2. chromStart: The beginning position of the defined genomic region, in 0-based positions.

3. chromEnd: The ending position of the defined genomic region, in 1-based positions. The chromEnd base is not included in the display of the feature, so it behaves much like ranges and substrings in python.

So using this core required format, the BED file can be a useful streamlined representation of genomic positions where strand is not important. We can expand the file format to include other information, for example when we want to store data for gene models. The next columns we could add are:

4. name: Defines the name or ID of the feature or genomic region.

5. score: A numerical value defining a score from 0 to 1000. Could be replaced with a dot.

6. strand: The strand of the genomic region, represented by “+” for forward strand or “-” for reverse complement.

7. thickStart: The beginning position of the portion of the genomic region to be drawn larger than the rest, for example a CDS to stand out from the rest of the gene.

8. thickEnd: The ending position of the portion of the genomic region to be drawn larger.

9. itemRgb: This is used to set the color of the genomic feature for some genome browsers, defined as a value for red, green, and blue as a number from 0 to 250 for each. For example, green would be (0,255,0), or blue would be (0,0,255), or many shades in between, such as a type of blue-green might be (0,100,100).

10. blockCount: For features that can be in multiple pieces, such as eons, this sets the number of blocks.

11. blockSizes: For features that can be in multiple pieces, such as exons, this comma-separated string of lengths, sets the length of each piece.

12. blockStarts: This comma-separated string of positions defines the start position of the exons or blocks.

GFF Files

GFF format is a great way to store gene annotation information. In many ways the columns of GFF are designed for genes, but by including a dot “.” for the columns that don’t apply to your annotation, you can also store more simple annotations such as the locations of a motif instance.

1. seqname – The name of the sequence. Must be a chromosome or scaffold.

2. source – The program that generated this feature.

3. feature – The name of this type of feature. Some examples of standard feature types are “CDS”, “start_codon”, “stop-codon”, and “exon”.

4. start – The starting position of the feature in the sequence. The first base is numbered 1.

5. end – The starting position of the feature (inclusive).

6. score – A score between 0 and 1000. If the track line useScore attribute is set to 1 for this annotation data set, the score value will determine the level of gray in which this feature is displayed (higher numbers= darker gray). If there is no score value, enter “.”.

7. strand – Valid entries include “+”, “-“, or “.” when strand doesn’t apply.

8. frame – If the feature is a coding exon, frame should be a number between 0-2 that represents the reading frame of the first base. If the feature is not a coding exon, the value should be “.”.

9. group – All lines with the same group are linked together into a single item.

GTF Files

A GTF file is very similar to a GFF file, but with a few different specifications. The first eight GTF fields are the same as GFF. The group field has been expanded into a list of attributes. Each attribute consists of a type/value pair. Attributes must end in a semi-colon, and be separated from any following attribute by exactly one space. The attribute list must begin with the two mandatory attributes:

1. gene_id value – A globally unique identifier for the genomic source of the sequence.

2. transcript_id value – A globally unique identifier for the predicted transcript.

BAM/SAM

Col

Field

Type

Brief description

1

QNAME

String

Query template NAME

2

FLAG

Int

bitwise FLAG

3

RNAME

String

Reference sequence NAME

4

POS

Int

1-based leftmost mapping POSition

5

MAPQ

Int

MAPping Quality

6

CIGAR

String

CIGAR string

7

RNEXT

String

Ref. name of the mate read

8

PNEXT

Int

Position of the mate/next read

9

TLEN

Int

observed Template LENgth

10

SEQ

String

segment SEQuence

11

QUAL

String

ASCII of Phred-scaled base QUALity+33

A BAM file is essentially a binary, typically indexed version of SAM. The program can allow quick retrieval of reads from a genomic region. SAM files is an optional output format for many alignment algorithms, but many people prefer to convert them to BAM files because of faster retrieval of reads for a particular genomic position due to indexing.

Quantitive Tracks

We often would like to annotate quantitative data on a genome browser; hence, it is critical to have file formats devoted to this kind of data. Consider, for example, plotting the GC content as a function of position throughout the genome.

BedGraph

track type=bedGraph
chrom1 chromStart1 chromEnd1 dataValue1
chrom2 chromStart2 chromEnd2 dataValue2

Wiggle, and BigWig

Wiggle file format is a common way to display quantitative tracks. The format is relatively simple, but takes on two different version: , and . Each of these are specified at the top of the file. For the fixed step, we have a fixed number of bases between positions presented in the file:

fixedStep chrom=chrN start=position step=stepInterval [span=windowSize]
dataValue1
dataValue2
dataValue3
...

Because the start position is specified, and the step size is specified, positions don’t need to be specified in the file. For the variable width, we’ll need to specify the position for each value:

variableStep chrom=chrN [span=windowSize]
chromStart1 dataValue1
chromStart2 dataValue2
chromStart3 dataValue3
...

The optional parameter “span” does not include the brackets when used in practice, and here only indicates that it is optional. The span essentially indicates that a value can be specified as applying to a range of positions, starting at the given chromStart value.

5.3 Genome Browsers

Genome browsers provide an interactive way to navigate the data associated with a genome in a visual way. Much like a web browser or an interactive map application. The file formats given above are exactly the kind of files that a genome browser would read and present visually.

5.3.1 IGV

The Integrated Genome Viewer (IGV) is a powerful desktop genome browser that allows you to relatively easily add and remove tracks and modify them. IGV allows you to export to various image formats, including scalable vector formats such as SVG files. In some cases the available RAM on one’s computer may be a limitation for loading too many tracks into IGV.

5.3.2 UCSC Genome Browser

The UCSC genome browser is a web-based genome browser. You can download, install, and host a version of the UCSC genome browser on your own computer, but you can also add tracks to the genome browser hosted at https://genome.ucsc.edu/cgi-bin/hgGateway.

5.3.3 Gbrowse

Gbrowse is probably one of the most common genome browsers out there. Many databases such as Flybase, or Wormbase have Gbrowse integrated in to the database for users to navigate the genomic data presented. Gbrowse allows for the display to be exported into multiple file formats, including scalable file formats.

5.3.4 JBrowse

Jbrowse is very similar to Gbrowse, but allows for asynchronous queries to the database, effectively making for a faster experience. One can scroll the position rapidly and have features immediately presented to the user without having to reload the page. Jbrowse allows the display to be exported in PNG, but not in a scalable file format.

5.4 Lab 6: Genome Annotation Data

In this lab, we will examine a script that will take a GFF annotation as input, and will extract protein-coding exons, concatenate them, and then translate them into a protein sequence. As usual, let’s work on this project in its own directory.

5.4.1 Part I: Storing a Genome to a Dictionary

A dictionary is a data structure with key-value pairs. Dictionaries are indexed by a “key”, which can be any string. This could be a useful concept for dealing with many sequences from a FASTA file, with the defline as the key.

>>> dict = {}
>>> dict["R2D2"] = "droid"
>>> dict["Vader"] = "Sith"
>>> dict["Yoda"] = "Jedi"
>>> print(dict)
{'Vader': 'Sith', 'R2D2': 'droid', 'Yoda': 'Jedi'}
>>> print(dict["R2D2"])
droid

Note how the dictionary looks when you print it out. You see a list of key-value pairs, separated by a colon “:”. Moreover, whereas the list is enclosed by square brackets when printed out, the dictionary is enclosed by curly braces when printed out. Nevertheless, both of them allow you to access a particular term with square brackets such as or .
It turns out, it is a pretty good way to store a genome, namely with the chromosome names (deflines) as the key, and the sequence as the value.

genome = {}
sequences = SeqIO.parse(genomeFasta,"fasta")
for record in sequences:
genome[record.id] = record.seq

5.4.2 Part II: Storing a GFF to a list

Next, let’s review creating a list and appending values to it. In this case, we will append (chrom,start,stop) genomic locations corresponding to CDS exons:

coords = []

GFF = open(gffFile,'r')
for line in GFF:
if "#" not in line:
chrom,source,seqtype,start,stop,score,strand,frame,attributes = line.strip().split("\t")
if "CDS" in seqtype and name in attributes:
coords.append((chrom,int(start),int(stop)))

A couple of things to note. First, we have a gene name in mind, stored in the variable “name”. We expect this to be a transcript ID because a gene ID could end up printing too many lines (one for each transcript associated with that gene). If we want to extract coding regions, it should be on a per-transcript basis. Second, we note that we are converting the positions to integers upon reading in. By default they are stored as strings, so the conversion is important. Third, note that we are excluding lines with a hash because these are typically comments that do not contain the required number of columns.

5.4.3 Part III: Find a gene of interest in Drosophila melanogaster

Let’s use NCBI Protein or NCBI Gene to find a your favorite gene. Let’s then BLAST it to Drosophila (by setting Drosophila as the taxa), and find the best Drosophila ortholog.

Next, let’s find a for one of the isoforms of the gene by navigating the flybase database and Gbrowse from flybase. Note that a looks like this: FBtr001828
Go to Flybase (http://flybase.org/) and search for a gene of interest, then paste your gene name into the “Jump to Gene” box on the upper left.

$ wget http://hendrixlab.cgrb.oregonstate.edu/teaching/flybase_r5.56.gff3

To use this script, you’ll need to download a flybase annotation. Download this one with wget:

$ wget http://hendrixlab.cgrb.oregonstate.edu/teaching/flybase_r5.56.gff3

Next, let’s look at the gff file. Let’s use less, but turn off word-wrap with the -S option.

$ less -S flybase_r5.56.gff3

Next, you can create a symbolic link to the genome file you created in Lab 4 (section 3.5) to your current directory with a command like this, but you may need to updated it based on where the file is:

$ ln -s ../Lab4/dm3.fa genome.fa

In this command, the first file “../Lab4/dm3.fa” is your original genome file. The second file name “genome.fa” is the name of the symbolic link you will be creating (but you can name it what you want). Now, let’s take a look at the script.

import sys
import re
from Bio import SeqIO
from Bio.Seq import Seq

# this section takes care of reading in data from user
usage = "Usage: " + sys.argv[0] + " "
if len(sys.argv) != 4:
print(usage)
sys.exit()

# read the input files/args.
genomeFasta = sys.argv[1]
gffFile = sys.argv[2]
name = sys.argv[3]

# read the fasta file into a dictionary.
genome = {}
sequences = SeqIO.parse(genomeFasta,"fasta")
for record in sequences:
genome[record.id] = record.seq

# initialize some variables
revcomp = False
coords = []

# read through the gffFile
GFF = open(gffFile,'r')
for line in GFF:
if "#" not in line:
chrom,source,seqtype,start,stop,score,strand,frame,attributes = line.strip().split("\t")
if "CDS" in seqtype and name in attributes:
coords.append((chrom,int(start),int(stop)))
if strand == "-":
revcomp = True

# reverse the order of the exons if on "-" strand
coords.sort(reverse = revcomp)

# collect the coding exons of the transcript
ORF = Seq('')
for chrom,start,stop in coords:
CDS = genome[chrom][start-1:stop]
if revcomp:
CDS = genome[chrom][start-1:stop].reverse_complement()
# concatenate the coding sequence
ORF += CDS

# transcribe,translate, and print
RNA = ORF.transcribe()
Protein = RNA.translate()
print(Protein)

Finally, after you copy the text of the above script into a file called “extractGeneAndProtein.py”, and put into a scripts directory you can run the script using the transcript ID for your favorite gene. Note, here we are using the genome FASTA file created in a previous lab, Lab 4 (section 3.5), but with the symbolic link created above.

$ python scripts/extractGeneAndProtein.py genome.fa flybase_r5.56.gff3 FBtr0089362

Chapter 6: Transcriptomics

Broadly speaking, Transcriptomics is the study of transcriptomes, the sum total of all transcripts in a cell. Transcriptomics seeks to build transcriptome annotations, and to measure differential expression of transcripts from different tissue types or treatments.

6.1 High-throughout Sequencing (HTS)

High-throughput sequencing (also known as deep sequencing) is a technology that has been developed in the late 20th century and continues to improve today. High-thoughtput sequencing has many applications, and most relevant for transcriptomics is deep sequencing of RNA, called RNA-seq. The word “deep” in deep sequencing refers to the depth of sequencing, characterized by:

where the depth is computed from the number of reads , the length of the reads , and the size of the transcriptome . The size of the transcriptome can be thought of the length of the union of all transcripts for a particular system. A depth of means that on average a location in the genome would have reads mapping to that location, assuming a uniform distribution of reads. This equation assumes that the reads are uniformly distributed, which is almost never true. Nevertheless it serves as a good approximation.

High-throughput sequencing can produce hundreds of millions of reads per sequencing lane, and in many cases the lane is multiplexed to include multiple samples per lane. This technology has enabled scientists to study biological phenomena at a genome-wide scale, and has enabled the discovery of a number of properties of transcription.

6.2 RNA Deep Sequencing

RNA deep sequencing is a method where a cDNA library is created for an RNA sample, and is sequenced using high-throughput sequencing, producing hundreds of millions of reads. Notably, there are different types of RNA-seq data sets. First, single-end reads involve the sequencing of one read per cDNA fragment, typically in the 5′ to 3′ direction. Paired-end reads have two reads per fragment, with the two paired-reads called “mates”. Often the first mate is sequenced in the direction of transcription, and the second mate is sequenced in the opposite 3′ to 5′ direction. This, however, can vary on the sequencing technology used. The manual for tophat 2 (https://ccb.jhu.edu/software/tophat/manual.shtml) provides the information on Figure 6.4.

Figure 6.1: A representation of single-end vs paired-end RNA sequencing. While single-end returns one read per fragment, paired-end returns two reads, sequenced from opposite ends of larger fragments.

6.2.1 Single-end Sequencing

Single-end Sequencing produces one read per fragment, so it can be good for transcript quantification, but may not resolve differences in expression across splice variants or different isoforms of the same gene. Therefore, it can be good for quantifying small RNA expression, or expression at the gene-level when splice variants are not a concern.

6.2.2 Paired-end Sequencing

Paired-end Sequencing produces two reads per fragment, and where typically a fragment size distribution is known or estimated. The result is the information of both reads in a pair, often called “mates”, can be helpful in transcriptome assembly and more precise quantification of different splice variants.

Important: to get the most out of paired-end sequencing, the fragment size should be larger than the combined read length (sum of both reads).

Typically with paired end data, one receives two files labeled and . The reads in each file correspond to pairs if they have the same read ID, excluding the possibility of the reads to be labeled and or possibly and . In practice, the library type can be determined by aligning paired reads from both the R1 and R2 files to the genome, and examining the relative orientation of the reads and overlapping transcripts.

6.3 Small RNA sequencing

For small RNA sequencing, one typically uses single-end sequencing, which results in one read per cDNA fragment, typically in the 5′ to 3′ direction.

Figure 6.2: A representation of small RNA sequencing. Reads may contain parts of the 3′ adapter.

6.3.1 Adapter Trimming

Small RNA sequencing uses size-selected small RNA samples and high-throughput sequencing. Such a protocol can be used to sequence RNA species such as microRNAs and piRNAs whose endogenous mature nucleotide sequences can be shorter than the read length used to sequence them. A challenge presented here is that the 3′ adapter sequence needs to be removed before aligning these sequences. For example, the program can be used to trim adapter sequences from an input FASTQ file.

cutadapt -a ATCTCGTATGCCGTCTTCTGCTTG reads.fastq > reads_trimmed.fastq

where we are trimming off the Illumina TruSeq Small RNA 3′ adapter sequence.

6.3.2 Alignment of small RNA reads

In general, for both small RNAs and large RNAs, RNA-seq read alignment typically takes a FASTQ file as input, aligns to the genome, and produces a BAM or SAM file as the output.

One method for aligning reads such as small RNA reads is . As one of the first methods for aligning deep sequence data, it does not allow gaps (not until bowtie2) but has other functionality such as colorspace read mapping. As a simple example, you would need to build an index of the genome that you’ll be aligning to (here we are using human genome version hg38) using . The index of the genome, much like a BLAST index created with provides a quick look-up of genomic locations to assist with the alignment.

$ bowtie-build -f hg38.fasta hg38

Then you would align the reads (a FASTQ file) to the genome index to create an alignment (a SAM file).

$ bowtie -m 50 -l 20 -n 2 -S -q hg38 smallRNA_reads.fastq smallRNA_hg38.sam

This command takes in several options. In general, the command is the following:

$ bowtie [options]

Where the square bracket terms are optional, and the angle brackets are required (or highly recommended for most practical purposes). In this example, we have a few options specified:

$ bowtie -m 50 -l 20 -n 2 -S -q hg38 smallRNA_reads.fastq smallRNA_hg38.sam

The specifies reads to have no more than 50 hits to the genome. The specifies a seed length that will be used for matching to the genome. The number of mismatches to the seed is specified as in the command. The specifies SAM output, and specifies FASTQ input.

6.3.3 Colorspace

As mentioned, has the ability to align colorspace reads. Colorspace reads typically come in a file like this:

>311_5120_1770
T322302111212131102211023200
>311_5120_1780
T303223021112121311011230122

The numbers correspond to the “colors”, which are not literal colors, but rather groups of dimers designated by the groups in Table 6.3.3 and in Figure 6.3.

AA

CC

GG

TT

0=blue

AC

CA

GT

TG

1=green

AG

CT

GA

TC

2=yellow

AT

CG

GC

TA

3=red

The quality scores (PHRED score values) will be in an accompanying file:

>read1
31 31 29 31 28 29 33 22 32 28 33 32 28 30 32 30 31 ...
>read2
31 28 29 32 28 22 32 28 33 32 28 30 32 30 31 31 32 ...

Consider the sequence . The sequence is determined through a process called “double interogation”. Each nucleotide is determined by two color characters. Therefore, the first character, called the “primer base”, gets the process going. The primer base is removed along with the adjacent color character before alignment. The rest is determined as follows. First number is after the , leading to the only dimer in the row starting with , which is . Next, the number is after the from before. The only dimer available starting with is . This process is continued across the sequence giving the complete sequence .

A. Colorspace designations can be determined on the following matrix, where the row gives the first nucleotide, and the column gives the second. B. Each dimer in a color space group is closed under the three nucleic acid transformations: reverse, complement, and reverse complement.

Color space also has the property that the reverse colorspace sequence is the reverse complement of the original. In the previous example, we saw that the colorspace sequence encoded the DNA sequence
. The reverse colorspace sequence would be , which encodes the reverse complement . This rule holds accept for the the first character, which is the priming base. Note that both colorspace sequences also had an additional number at the beginning ( for the first one, and for its reverse complement) which does not work with this reverse/reverse complement rule.

6.3.4 Quantifying small RNA Expression

For small RNA sequencing, the length is less of a concern. For example, with mature microRNAs, whether it is 18 or 25nt, you typically sequence one read per fragment. This fact suggest that small RNA reads don’t need to be normazlied by the length. Therefore, for a microRNA with reads mapping to it out of total mapped reads, we would compute the expression as

(6.1)

With the module , you can load a BAM file with the following command:

import pysam
bam = pysam.AlignmentFile("smallRNA_hg38.bam", "rb")

In this command, the “rb” term specifies “read” and “bam”. If it was a SAM file, you would just use “r”. You can retrieve the reads overlapping a genomic location with a command like this:

miR_count = bam.count(chrom,start,end)

Where are a predefined set of variables that could have been read in from a GFF file similar to Lab 6. For this function, the chromosome is a string, and the start ad stop are both ints. These ints behave like normal python boundaries, with start 0-based and the end being 0-based non-inclusive.

A better method of counting the microRNAs might also check the strand of the read.

import pysam
bam = pysam.AlignmentFile("smallRNA_hg38.bam", "rb")

chrom,start,end,strand = location

miR_count = 0
for read in bam.fetch(chrom,start,end):
if read.is_reverse:
if strand == '-':
miR_count += 1
else:
if strand == '+':
miR_count += 1

6.4 Long RNA sequencing

The sequencing of long RNAs, such as mRNAs and long non-coding RNAs (lncRNAs) requires long read sequencing, or simply called “RNA-seq”. Single-end Sequencing produces one read per fragment, so it can be good for transcript quantification, but may not resolve differences in expression across splice variants or different isoforms of the same gene. Therefore, it can be good for quantifying small RNA expression, or expression at the gene-level when splice variants are not a concern.

6.4.1 Quality Trimming

It is good practice to trim off low-quality nucleotides from the reads based on the PHRED score. Similar to adapter trimming, this will improve the alignment.

$ fastq_quality_trimmer -t 30 -l 18 -i input_reads.fastq -o output_trimmed_reads.fastq

Is the PHRED score threshold specified by a good one? The threshold specifies that trimmed reads less than 18nt should be discarded.

6.4.2 Paired-end Library Types

The different arrangements of paired-end reads depticted in Figure 6.4. In practice, the library type can be determined by aligning paired reads from both the R1 and R2 fastq file to the genome, and examining the relative orientation of the reads and overlapping transcripts.

A representation of an the different library types found in RNA-seq data used in , , and . The most common is fr-secondstrand, with the read in the direction of transcription.

6.4.3 Methods of Read Alignment

The first step of an RNA-seq analysis is aligning the reads to the genome. Actually, there are approaches to RNA-seq without a genome discussed in section 6.4.6. For organisms with an assembled genome, the typical approach is to align the reads to the genome. There are many different programs to do this, with different features. For eukaryotic organisms, it is important that the alignment allows for gaps corresponding to introns. Another important consideration is the number of mismatches allowed in the alignment. For a highly polymorphic genome, for example, one would want to increase the number of allowed mismatches under the assumption that the transcripts sequenced could be transcribed from a slightly different genome.

RNA-seq read alignment typically takes a FASTQ file as input, aligns to the genome, and produces a BAM or SAM file as the output.
Methods for aligning RNA-seq reads include and , although is considered the most updated software. These programs allow gaps and both are derived from .

Alignment with

The first step of aligning is to have an index of the genome. As with bowtie, you would need to build an index of the genome that you’ll be aligning to (here we are using human genome version hg38) using . Note that the index files created here won’t be compatible with and vice versa.

$ hisat2-build hg38.fasta hg38

Then you would align the reads (a FASTQ file) to the genome index to create an alignment (a SAM file) using .

$ hisat2 -q --rna-strandness F -x hg38 -u RNAseq_reads.fastq -S RNAseq_hg38.sam

This command takes in several options. In general, the command is the following:

$ hisat2 [options] -x -u -S

Where the square bracket terms are optional, and the angle brackets are required (or highly recommended for most practical purposes). In this example, we have a few options specified:

$ hisat2 -q --rna-strandness F -x hg38 -u RNAseq_reads.fastq -S RNAseq_hg38.sam

The specifies that the reads are sequenced in the direction of transcription. It could also be or if not stated would assume unstranded reads by default. The specifies FASTQ input, but this is the default so not strictly needed. Other options could include if you wanted to keep intronic gaps less than 100kb.

Alignment with

As an example of aligning reads to the genome, let’s consider , which is part of the Tuxedo suite. For the other steps of the RNA-seq analysis pipeline we will also examine parts of the tuxedo suite pipeline.

In the case of , which uses bowtie2 for the alignment, the creation of the index is based on the Burrows-wheeler transform. The index of the genome consists of multiple files, but each has the same prefix, or “file base” as an input. For example, with the mouse genome that has a file base , we will have the six files , , , , , .

So, using default parameters, aligning with would be performed by the following command:

$ tophat -o reads_mm10_tophat mm10 reads.fastq

where the option species the output directory where the output files are going to go. Probably the most important of the output files in this directory is the BAM file called . This file contains a sorted record of the alignment information. There are other files created such as log files that could provide useful information for tracking down issues if something goes wrong. Also among these output files is , which provides a BED file of all splice junctions identified, and the number of reads that span that particular junction.

Of course, we don’t always want to run with default parameters. For example, we might want to align with a reference annotation, such that only gaps are considered that correspond to gaps in an input GTF file. To achieve this, we would use the following command:

$ tophat -o reads_mm10_tophat --no-novel-juncs -G ensembl_mm10.gtf mm10 reads.FASTQ

Here we have specified to not include any novel junctions with the flag. Required with this flag is to specify a GTF file with the flag. In this case we have specified an annotation from Ensembl for mm10.
We may consider adding to this a command to specify the number of mismatches (the default is 2). This can be changed with the flag.

Quantifying Expression

Gene expression studies measure the expression levels of genes by attempting to quantify the number of mRNA transcripts per gene. There are a number of ways to compute gene expression, with varying accuracy and different pros and cons.

6.4.4 Quantifying Expression with Microarrays

Microarrays allow for high-throughput measurement of gene expression through the use of hybridization probes. The probes are DNA sequences that are complementary to the transcripts that you would like to measure (or more precisely the cDNA created from the transcripts). Because of the requirement of probes to detect the gene expression, one or more probes for each gene needs to be designed. This puts a constraint on the genes that one can detect. Also, different hybridization energies can complicate the analysis in some cases, as this would need to be controlled for.

6.4.5 Quantifying Expression with RNA-seq

When quantifying expression with RNA-seq, it is important to consider a normalization of the data that best controls for the fact that sequencing depth may vary from experiment to experiment, and gene lengths are highly variable. Longer genes generate more fragments, and hence result in more reads per physical mRNA molecule. One method that has these properties is RPKM, which stands for Reads Per Kilobase of gene length per Million reads mapped. For a gene of lenxth , with reads mapping to it out of total mapped reads, we would compute the expression with

(6.2)

Taking this idea further, another measure of expression is the FPKM, or Fragments Per Kilobase of gene length per Million reads mapped. Here, the fragment refers to the cDNA fragment from which the read was sequenced. For example with paired-end reads, we have two reads per fragment. If one of the mates was of poor quality or did not map, the reads would underestimate the number of fragments mapping to this location. For example, when aligning with , there is an option , which would only allow reads where both mates properly map the the same chromosome or scaffold. When this option is not used, there are cases where one fragment corresponds to one read (either or ) or both reads. When the option is used, each fragment corresponds to two reads.

Measuring Differential Expression with RNA-seq

The program is part of the software package, and it computes expression values of genes and identifies significantly differentially expressed genes. with default parameters, we can run cuffdiff to compare the expression of two RNA-seq experiments with the following command (input on one line):

$ cuffdiff -o treatment_vs_control -L treatment,control \\
ensembl_mm10.gtf treatment_mm10/accepted_hits.bam control_mm10/accepted_hits.bam

So we note that in this case it is required that we specify a GTF annotation file to use to use such that the expression of the genes annotated in this file are quantified. Furthermore, we need two BAM files, presumably both aligned with Tophat using the same parameters. The parameter specifies “labels” for the BAM file compared, such that the output files will use these in the file header. For each BAM file, we could also give SAM files, or a comma-separated list of BAM files corresponding to individual replicates for that experiment. The command specifies the output directory in which the results will be print. Lastly, the is used to split the command on two lines here, but is not necessary to type.

The results of a cuffdiff computation is a set of many files. Among them are differential expression files called files. There are files for the genes and for the isoforms for the genes given by and files respectively. These files contain the information on the genomic location of the transcript, as well as FPKM values for the samples compared. Furthermore, it has the log fold change of these expression values, as well as output from a statistical test to identify significantly differentially expressed genes.

Intersection Count, Union Count, and the True Expression

When counting reads that overlap a gene model, there are different approaches to quantifying the expression of a gene with multiple transcript isoforms. For example, the intersection count would take the expression of all exons that are common to each isoform in a particular gene. The union count take the perspective of counting the reads that overlap the union of all exons for a particular gene. As pointed out in Trapnell et al. 2013, there are issues with both of these models. Figure 6.5 demonstrates a scenario with both the intersection and union count models fail to accurately depict the expression of a gene.

Figure 6.5: A depiction of the intersection and union models, and how they can differ from the actual expression. This diagram is inspired from Figure 1 from Trapnell et al. 2013. In each case, the expression per length is depicted for the two isoforms where each exon of the longer isoform is of length , and both for the shorter (red) isoform is .

6.4.6 Transcriptome Assembly

Transcriptome assembly is the process by which the aligned reads are compiled into transcript models that best represent the read data. At minimum, a transcript model must consist of reads that map to that location. Introns must be supported by a gapped alignment that spans the intronic region, typically further requiring that the donor and acceptor sites are present at the start and end of the intron (GU on the 5′ end and AG on the 3′ end of the intron). We may optionally want to assemble a transcriptome for an RNA-seq experiment when we don’t have a transcriptome annotation available. We may also want to perform this step when none is available. A schematic of how this process is performed is depicted in Figure 6.6

Figure 6.6: A schematic representation of an RNA-seq pipeline for computation of transcript models.

Stringtie

In general, running uses commands of the following form:

$ stringtie [options] RNAseq_reads.bam -G human_RefSeq.gtf -o transcriptome_v2.gtf

Where the square bracket terms are options. In fact, the is also optional, but builds the transcriptome assembly using the provided GTF file as a guide. Other options include for example would require 2 reads to build a splice junction, rather than the default value of .

You can also use stringtie to combine and merge several GTF files into one. There are two ways to do this.

$ stringtie --merge [options] gtf_file_list.txt

Where is a list of GTF files. Alternatively, you can specify the GTF files in a series after the command:

$ stringtie --merge [options] human_transcripts1.gtf human_transcripts2.gtf

Cufflinks

The next step of the Tuxedo suite that corresponds to transcriptome assembly is . This can be run with default parameters with the following command:

$ cufflinks reads_mm10_tophat/accepted_hits.bam

where we have specified the BAM file produced from the previous commands. We may want to add a reference annotation to our command. This may seem counter-intuitive since the purpose is to create our own new transcriptome assembly, but doing so allows cufflinks to assign gene names corresponding to known gene names in the annotation file. We can do this with the updated command:

$ cufflinks -g ensembl_mm10.gtf reads_mm10_tophat/accepted_hits.bam

Note that we are specifying the same GTF file as before, but cufflinks uses a lower-case “-g”, whereas uses an upper-case “-G”. We may want to consider is to make the alignment run “quietly”, by reducing the number of printed messages with the option. Another option worth considering is the option (capital i), which specifies the maximum intron length, or gap-length to consider. Restricting this can occasionally remove artifacts since the default is bp.

Transcriptome Assembly without a Reference Genome

There are approaches to assemble a transcriptome and quantify gene/transcript expression without a reference genome. Prominent among them is the suite called “Trinity”, which consists of three pieces of software, as the name suggests [54].

Figure 6.7: A schematic representation of a non-reference-guided transcriptome assembly, Trinity (Grabherret al.2013).

The software components of Trinity work in series to create a transcriptome assembly in three steps (Figure 6.7). First, clusters reads by common -mers (subsequences of length ). Next, clusters contigs into components consisting of matched -mers, and builds a de Bruijn graph from these components. In the last step, assembles transcript models from the de Bruijn graph [54].

6.5 Single-Cell Transcriptomics

While traditional RNA-seq measures gene expression on many, typically heterogeneous cells, single-cell RNA-seq (scRNA-seq) is a method for measuring gene expression in hundreds of individual cells [55]. While these methods do reduce heterogeneity compared to bulk measurements, they do have a large amount of technical noise [56].

6.6 Transcription Initiation

Transcription Initiation is a remarkable process, because it is very much a needle-in-a-haystack phenomena. Somehow, the cell is able to pinpoint specific single-nucleotide positions out of the genome, and use these specific positions as the start of a transcript. This first nucleotide that is transcribed is called the Transcription Start Site (TSS).

In many cases, transcription initiation begins with the recognition of the core promoter, which is a collection of binding sites that direct the initiation of transcription. In eukaryotes, RNA polymerase II does not directly recognize the core promoter sequences. Instead, a collection of proteins called transcription (activation) factors (TAFs) mediate the binding of RNA polymerase and the initiation of transcription. Only after certain transcription factors are attached to the promoter does the RNA polymerase bind to it. The completed assembly of transcription factors and RNA polymerase bind to the promoter, forming a transcription initiation complex.

In prokaryotes, transcription begins with the binding of RNA polymerase to the promoter in DNA. At the start of initiation, the core enzyme is associated with a sigma factor that aids in finding the appropriate -35 and -10 base pairs upstream of promoter sequences. When the sigma factor and RNA polymerase combine, they form a holoenzyme.

6.6.1 Methods of Mapping Transcription Start Sites (TSSs)

There are many methods that have been developed for the mapping of transcription start sites. Most of these methods rely on the biochemistry of the 5′ end of transcribed RNA.

Cap Analysis of Gene Expression (CAGE)

CAGE sequencing begins by capping the 5′ ends of small transcript fragments, which are then extracted, reverse transcribed to DNA, PCR amplified and sequenced. CAGE reads tend to map at the transcription start site of genes , but in practice do map at and around these locations. The CAGE technique precedes that of high-throughput sequencing and was used before then, but has benefited from the depth associated with new sequencing techniques.

Rapid Amplification of edna Ends (RACE)

In RACE, a cDNA copy of the mRNA is produced through reverse transcription, PCR amplification of the cDNA copies. There are different methods for 5′ RACE and 3′ RACE, resulting in the proper mapping of both the 5′ end and the 3′ end of genes. 3′ RACE uses the poly(A) tail of genes for priming during the reverse transcription. 5′ RACE uses a gene-specific primer that recognizes a known region of the gene of interest.

Database of Transcription Start Sites (DBTSS)

The Database of Transcription Start Sites (DBTSS) is a resource of experimentally mapped transcription start sites, through a method called TSS-seq, and is primarily focused on human and mouse, but contains other species. This database now contains 491 million TSS tag sequences collected from 20 different human tissues and 7 different human cell cultures. The DBTSS can be accessed at http://dbtss.hgc.jp/.
The TSS-seq method ligates the Illumina sequencing adapter to the 5′ cap site of the mRNA, performs full-length cDNA creation, and high-throughput sequencing. The uniquely mapping sequencing reads are then mapped to the genome, and clustered into bp bins. TSCs overlapping internal exons are then removed, and each cluster is either associated with the most likely RefSeq transcript, or labeled as intergenic

6.7 Transcription

Transcription is carried out by RNA Polymerase, which makes a RNA copy of the template strand of a gene’s DNA sequence. Measurements of the binding of RNA Polymerase II (Pol II) to genes has revealed that not all genes have the same binding profiles over their promoter regions and gene bodies.

6.7.1 Measuring RNA Polymerase binding: Pol II ChIp-seq

Pol II ChIP-seq provides a signal the binding of Pol II to the genomic DNA. It measures the recruitment of Polymerase to the genome, but the Ser5 of the Pol II complex needs to be phosphorylated for transcription to begin. Therefore, many of these protocols use antibodies that specifically recognize the Ser5-Phosphorylated form of Pol II as well as other forms.

Briefly, chromatin samples are crosslinked and sonicated into fragments. Chromatin is immunoprecipitated with a RNA Pol II antibody, or cocktail of antibodies that recognize different forms of Pol II. DNA fragments are then isolated and sequenced with HTS. These reads are then mapped to the genome, and the number of reads that map to a particular genomic region is indicative of the time that Pol II spends binding to that region.

6.7.2 RNA Polymerase II Stalling

One of the phenomena observed using Pol II ChIP-seq is stalling (Pol II Stalling) [57]. Stalled Pol II is capable of quickly responding to input signals because the transcript is already initiated, and hence is called “poised” [58] for rapid activation. Pol II Stalling is associated with a greater proportion of Pol II ChIP-seq reads at and around the TSS of genes compared to the gene body. For a Pol II signal given by for position of the genome, The state of being stalled is defined by a large stalling index:

where indicates the maximum signal value within some distance of the TSS of the gene, and indicates the median signal value within the gene body [57].

Figure 6.8: The Pol II signal over a gene can indicate its transcriptional status. A. Genes with a large proportion of Pol II binding at the TSS are “Stalled”. B. Genes with roughly uniformly high Pol II binding for the promoter compared to the gene body are called active. C. Genes with little or no Pol II binding are categorized as “No Pol II”.

6.8 Elongation

Elongation is characterized by the release of Pol II from the promoter, and the production of RNA transcripts. In mammalian genes, the Pol II elongation rates are about 0.5kb per minute in the first few kilobases, and increases to 2-5kb per minute after 15kb [59].

6.8.1 Measuring Nascent Transcription: GRO-seq

While Pol II ChIP-seq measures DNA-associated Pol II, Global Run-On Sequencing (GRO-seq) measures transcriptionally engaged Pol II, by detecting nascent RNA transcripts. GRO-seq involves a nuclear run-on experiment on a genome-wide scale, where engaged Pol II incorporates bromo-tagged nucleotides, and is followed by RNA isolation and deep sequencing [60].

6.8.2 Divergent Transcription

One of the phenomena observed using Global Run-On sequencing (GRO-seq) is that of divergent transcription [61]. Divergent transcription is characterized by nascent transcripts associated with both the sense and anti-sense strand of the gene.

Figure 6.9: GRO-seq can identify cases of divergent transcription, where nascent transcription is observed both sense and antisense to a genes TSS.

6.9 Lab 7: RNA-seq

The data for this project was taken from a paper by Kurasz et al., which describes a study of gene expression in Salmonella enterica serovar Typhimurium and the effect of nucleic acid damage induced by mitomycin (MMC), which is a chemotherapy drug, and other compounds [62]. We will look at the MMC data in this project. I’m assuming you’re working in a new directory called .
Due to the large amount of data generated, I have some streamlined versions of these files to make it possible to complete this lab in less than an hour.

6.9.1 Step 1: Download the genome and build a hisat2 index

First thing we should do is download the genome for Salmonella. For reference, I found this by searching NCBI in the “Assembly” section and found this link here: https://www.ncbi.nlm.nih.gov/assembly/?term=Salmonella+Typhimurium. I renamed the file and made it available at this link:

$ wget http://hendrixlab.cgrb.oregonstate.edu/teaching/salmonella_genome.fasta

Next, we should build a index with the following command:

$ hisat2-build salmonella_genome.fasta sal

In this command, I have used the filebase “sal”, which is specified as the second argument to above, but you could call it what you want and change the subsequent commands accordingly. Note that several files are created as a result:

$ ls sal.*
sal.1.ht2 sal.2.ht2 sal.3.ht2 sal.4.ht2 sal.5.ht2 sal.6.ht2 sal.7.ht2 sal.8.ht2

6.9.2 Step 2: Download the FASTQ files and align with hisat2

First, we will download the FASTQ file. I found this experiment on GEO here: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE82167}. As a side note, what I did was go to the “Samples” section, expanded that section, clicked the links to find the “SRR…” IDs for samples. I ran a command like this (but please don’t run this because the original files are too big!):

$ ~/sratoolkit.2.9.4-centos_linux64/bin/fastq-dump -A SRR3621123 >& SRR3621123.err
$ ~/sratoolkit.2.9.4-centos_linux64/bin/fastq-dump -A SRR3621125 >& SRR3621125.err

These commands assume you have installed the SRA toolkit in your home directory, and are using version 2.9.4, and are using CentOS; please alter these as needed. To make this more feasible in class, I have created the following shortened version of the resulting files:

$ wget http://hendrixlab.cgrb.oregonstate.edu/teaching/salmonella_ctrl.fastq
$ wget http://hendrixlab.cgrb.oregonstate.edu/teaching/salmonella_MMC.fastq

Question: How many lines are in these files and how many reads does it have?

We might want to quality trim these files, as well as remove any potential adapters. You can trim the files with using the commands:

$ skewer -x AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -q 30 -Q 30 salmonella_ctrl.fastq -o salmonella_ctrl
$ skewer -x AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -q 30 -Q 30 salmonella_MMC.fastq -o salmonella_MMC

For this command, it is requiring a PHRED score of 30 for both the 3′ end as well as the average score. Meaning, nucleotides will be removed from the 3′ end if they are lower quality than 30, and will discard reads with an average score less than 30.

The alignment itself is done with using mostly default options to create a SAM file as the output.

$ hisat2 -x sal -U salmonella_ctrl-trimmed.fastq -S salmonella_ctrl.sam
$ hisat2 -x sal -U salmonella_MMC-trimmed.fastq -S salmonella_MMC.sam

6.9.3 Step 3: Convert SAM to BAM, sort, and index the files

In order for cuffdiff to run properly, we need to sort and index the files. SAM will work actually, but it needs to be sorted. Perhaps the easiest and fastest thing is do do this with samtools, which is like a “Swiss Army knife” for SAM and BAM files. First, convert to BAM using samtools

$ samtools view -b -S -o salmonella_ctrl.bam salmonella_ctrl.sam
$ samtools view -b -S -o salmonella_MMC.bam salmonella_MMC.sam

Next, we sort and index these files:

$ samtools sort salmonella_ctrl.bam salmonella_ctrl.sort
$ samtools index salmonella_ctrl.sort.bam
$ samtools sort salmonella_MMC.bam salmonella_MMC.sort
$ samtools index salmonella_MMC.sort.bam

6.9.4 Step 4: Identify differentially expressed genes with Cuffdiff.

To do this, you’ll need to download the GFF file for the organism.

$ wget http://hendrixlab.cgrb.oregonstate.edu/teaching/salmonella.gff

Please read more information on the options of cuffdiff, check here: http://cole-trapnell-lab.github.io/cufflinks/cuffdiff/index.html. We are going to run a basic version here, and the only option we are setting is to create a label for the output file columns with meaningful labels.

$ cuffdiff -L WT,MMC salmonella.gff salmonella_ctrl.sort.bam salmonella_MMC.sort.bam

Do you find any differentially expressed genes? These will correspond to those significant at an FDR of 0.05, by default. They make this easy to putting a “yes” on the last column for genes that are significant.

$ grep yes gene_exp.diff

Do these genes make sense given the paper, and can you find information about their function at NCBI?

You can look for enriched functional annotations using the “GO Enrichment” tool at http://geneontology.org/. GO Terms are part of a controlled vocabulary for gene functions. Salmonella is an available organism on this page, which enables us to study gene enrichment associated with the differentially expressed genes from this experiment. By extracting the gene IDs, you could paste them into this GO Enrichment tool, and compute enriched functional categories.

Chapter 7: Noncoding RNAs

It was once believed that most genes encode proteins and most transcripts are messenger RNAs (mRNAs). According to some studies on the matter, the number of genes that do not encode proteins outnumbers those that do [63]. Unlike mRNAs, which are translated into proteins by the ribosome, these noncoding RNAs (ncRNAs) act functionally as an RNA transcript. NcRNAs that are greater than 200nt are called long noncoding RNAs (lncRNAs), and those that are less than 50nt can be called small noncoding RNAs (sncRNAs). The remaining intermediate size range from 50-199nt can be called short RNAs, although these exact length thresholds are a bit arbitrary.

7.1 Small Noncoding RNAs (srcRNAs)

As we learned in Chapter 6, the transcriptome includes small transcripts. These sncRNAs are involved in gene regulation and gene silencing, regulation of splicing, and possibly other epigenetic roles throughout the cell. Many sncRNAs act as guides to direct the sequence-specific activity of enzymes to other RNAs such as mRNAs or other transcripts. Other more recent studies have implicated sncRNAs in directing the binding to DNA through triplex interactions.

7.1.1 microRNAs

Figure 7.1: The microRNA biogenesis pathway.

MicroRNAs (miRNAs, miRs) are small endogenous regulatory RNA molecules (18-25nt) that are involved in post-transcriptional gene silencing [64]. miRNAs are processed from longer primary transcripts. Hairpins are cut at the base by Drosha, producing a microRNA precursor (pre-miR). These hairpin precursors are shuttled out of the nucleus and cut at the loop to produce double-stranded RNA. Typically one of the products is incorporated into the RNA-induced silencing complex (RISC), and the other is degraded.
Target sites for microRNAs are in the 3′ UTR of protein coding genes, although they have been observed in protein-coding regions as well. The target sites are around 7 nt in length, complementary to positions 2-8 of the mature microRNA.

7.1.2 piRNAs

Figure 7.2: The piRNA “ping pong” biogenesis pathway. Ago3 cuts the piRNA primary transcript by using a transposon-derived small RNA as a guide, and Aub cuts the transposon transcript using a mature piRNA as a guide. Both small RNAs are truncated by exonuclease digestion followed by methylation.

Piwi-interacting RNAs (piRNAs, piRs) are small endogenous RNAs molecules (21-31nt) that are part of the genome’s natural defense against transposons, and may be involved in post-transcriptional gene silencing and other epigenetic functions [65]. piRNAs are the most abundant class of small RNAs [66]. piRNAs are processed from longer primary transcripts based on complementary interactions with transposon transcripts. These small RNAs do not arise from hairpins, but rather the duplex interactions with the transposon transcripts. The enzymes Aubergene (Aub) and Argonaute 3 (Ago3) for processing [67]. Ago3 cuts the piRNA primary transcript by using a transposon-derived small RNA as a guide, and Aub cuts the transposon transcript using a mature piRNA as a guide.

7.1.3 Duplex RNA interactions

Many small RNAs function through the RNA duplex formation due to complementary sequence interactions. Let’s discuss how we can describe these types of interactions.

Consider the complementary RNA sequences and that could form a base-paired (heterodimeric) RNA duplex .

5' ACGCAUU 3'
|||||||
3' UGCGUAA 5'

The reaction for this duplex formation could be expressed as follows:

This is under the assumption that there are no structures formed by and on their own that need to unfold for the reaction to proceed forward, and that the individual molecules and don’t form homodimers on their own. In many applications, this equation is depicted the other way, and the process is described as a “melting” process whereby the duplex becomes unpaired.

Table 7.1: Base pair energies for the Nearest Neighbor model (source https://rna.urmc.rochester.edu/NNDB/turner04/wc.html). In each example, the dimers are given in the 5’ to 3’ direction, but it should benoted that they pair anti-parallel.

Sequence

(kcal/mol)

AA:UU

-0.93

AU:AU

-1.10

UA:UA

-1.33

CA:UG

-2.11

GU:AC

-2.24

CU:AG

-2.08

GA:UC

-2.35

CG:CG

-2.36

GC:GC

-3.42

GG:CC

-3.26

GC-init.

2.045

AU-init

2.495

symm

0.43

The energy required to drive such a reaction is set by the equilibrium constant for this reaction

And the free energy required to melt the duplex is given by

The temperature at which the paired and unpaired nucleotides become equally abundant is defined as the melting temperature

In our example above, the RNA duplex with sequences and forms a duplex. Since this system is completely paired, we can describe it as a sequence of base pairs , corresponding to the base pair at position $i$. In this example, using 1-based coordinates, we have . The nearest neighbor model describes the basepair energy between a duplex as the sum of dimer interactions, plus initialization terms on the ends. The full expression for the duplex energy is then

Calculating Duplex Energy with the Nearest Neighbor Model

Let’s consider an example. and . In this case, the

The summation terms can be expanded as follows

Adding up all the terms gives the final solution.

7.1.4 RNA/DNA Hybrids

RNA/DNA hybrids are structures formed by a “hybrid duplex” consisting of one strand of RNA and one strand of DNA [68,69]. If was assume that the RNA has no structure (i.e. is a sncRNA) and that the DNA doesn’t fold into a structure when the double helix is opened, the energy parameters of such an interaction can be computed using the difference between the DNA/DNA duplex energy and the RNA/DNA hybrid energy. The can be computed from the difference of these two values. Therefore, the formula is given by:

Figure 7.3: The steps in forming an RNA/DNA hybrid can be viewed as opening the DNA, followed by the RNA hybridizing with the DNA. Each step has a thermodynamic energy.

DNA/DNA duplex

first →
second ↓

G

C

A

T

G

-1.84

-2.17

-1.28

-1.45

C

-.2.24

-1.84

-1.44

-1.30

A

-1.3

-1.45

-1.00

-0.58

T/U

-1.44

-1.28

-0.88

-1.00

RNA/DNA hybrid

G

C

A

U

-2.2

-1.2

-1.4

-1.0

-2.4

-1.5

-1.6

-0.8

-1.4

-1.0

-0.4

-0.3

-1.5

-0.9

-1.0

-0.2

RNA/DNA – DNA/DNA

first →
second ↓

G

C

A

T

G

-0.36

0.97

-0.12

0.445

C

-0.16

0.34

-0.16

0.5

A

-0.1

0.45

0.6

0.28

T/U

-0.06

0.38

-0.12

0.8

7.1.5 Triplexes

Triple-helix (triplex) structures between a double-stranded nucleic acid sequence and a single-stranded nucleic acid sequence can form for specific triplets of nucleotides. For example, triplexes can form between a single-stranded RNA sequence and double-stranded DNA. The energy parameters of these interactions are not known, but interactions can be predicted with a program like : http://bioinformatics.org.au/tools/triplexator/

$ triplexator -l 15 -e 20 -o output.txt -ss transcripts.fasta -ds promoters.fasta

The output, stored to the file specified by the flag, is a tab-delimited file containing information on the locations of triplex forming oligonucleotides (TFOs) [70].

7.2. Long Noncoding RNAs

7.2.1 Quantifying coding potential

In order to identify noncoding RNAs, perhaps it is easier to characterize what is coding, and label the remaining as noncoding. Many approaches have been developed to quantify a transcripts potential to encode a protein ranging from comparative approaches, to machine learning, to experimental methods.

Coding Potential and Homology

The first approach is to use homology, such as by BLAST or using protein domains on predicted ORFs. If there is significant homology between the translation of an ORF and a known protein, then there is a good chance that the gene in question is also protein coding. Similarly, if an ORF translates to a sequence that has a known protein domain, then it is likely protein coding.
However, the caveat here is to make sure that the gene has a valid, functional ORF. If a mutation has caused the ORF to be lost (e.g. introduce an early stop codon or several frame-shifts) then there is a good chance it is a “pseudogene”, which could be considered noncoding.

Evolutionary Models

Evolutionary methods for identifying coding potential rely on the fact that protein coding genes have well defined sequences of codons, which constitute an open reading frame. Open reading frames (ORFs) are defined by a start codon positioned at a multiple of three nucleotides away from a stop codon. When considered comparatively across different species in a multiple alignment, one rarely sees nucleotide insertions or deletions that cause a “frame shift”. Frame shifts are when the phase of the codons changes, giving rise to a transcript that encodes a completely different amino acid sequence, if at all. Similarly, multiple sequence alignments of protein coding transcripts are depleted for mutations that cause a non-synonymous change, or mutations that alter what amino acid a codon encodes. The software PhyloCSF does exactly this, and quantifies the likelihood of a given transcripts being coding or noncoding based on its multiple alignment [71].

Machine Learning Approaches to Coding Potential

Machine learning approaches define a set of “features” that describe the phenomena that the program seeks to learn. For identifying protein coding genes, the features describe aspects of the transcript being evaluated. For example, the program (Coding Potential Assessment Tool) quantifies the frequencies of -mers in the transcript, as well as the length of the transcript, and is able to distinguish protein coding genes from noncoding fairly well [72].

Ribosome Binding Profiles

Figure 7.4: While protein coding genes exhibit a sharp decrease in ribosome binding after the stop codon, noncoding RNAs are more stochastic over the length of untranslated ORFs.

Ribosomal profiling quantifies the binding of the Ribosome to transcripts at each position. It has been determined that protein coding genes have a distinct ribosomal profile over the open reading frame, characterized by a sharp drop-off in binding immediately after the stop codon. Noncoding RNAs have a more stochastic binding profile around trinucleotides matching the stop codon at the end of an (untranslated) ORF [73]. Figure 7.4 depicts the kind of data produced in these experiments. To quantify the potential for coding, the authors developed a score called a Ribosomal Release Score (RRS). The equation for the RRS is

(7.1)

For protein coding genes, we expect this ribosome ratio to be high, since the Ribosome read count should be much greater than that of for such genes.

Detecting Peptides with Mass Spectrometry

Some approaches have used Mass Spectrometry to detect small peptides, and validate that a particular open reading frame was real. In these cases, small peptides can be detected that were predicted to be translated from a small ORF [74].

7.3 RNA Structure Prediction

Many ncRNAs operate by forming a structure, consisting of basepairs between nucleotides (formed by hydrogen bonds), as well as regions of unpaired nucleotides. An RNA structure can be viewed at many different levels, with the “primary structure” consisting of the sequence itself. The secondary structure consists of a mapping of which nucleotides are paired, and can be represented on a 2D surface. The tertiary structure consists of the three dimensional configuration of the molecule. The quaternary structure comprises the inter-molecular interactions that these molecules can participate in.

While the primary structure of the RNA can be described by a sequence of the characters , such that corresponds the the nucleotide at position of the molecule the sequence represents. The secondary structure can be represented by a set of ordered pairs such that each pair of positions that form a base pair constitutes an ordered pair , and the corresponding nucleotides and are complementary.

One could remove restriction that i<j, then you could consider the basepairs unordered pairs .

7.3.1 The Nussinov Algorithm

The Nussinov Algorithm predicts the number of basepairs, and ultimately the structure for an RNA sequence using dynamic programming [75]. This is done by understanding that whether or not a pair of nucleotides will pair depends on the intervening nucleotides from to .

Let’s consider an RNA sequence with nucleotides . Consider the possible base pair formed between positions and . This will of course depend on whether and are complementary or wobble pairs capable of forming hydrogen bonds.

Figure 7.5: A visual representation of the recursion relation in the Nussinov algorithm, which can be broken up into four possibilities.

Let’s build a scoring matrix such that the terms of this matrix correspond to the number of basepairs between these two positions.

(7.2)

The first three options are probably pretty clear. Either it adds a basepair at , or it doesn’t. If it adds a basepair, then increment to the previous value , that corresponds to the best structure prediction for the intervening positions.
If it doesn’t add a pair, then simply extend a previously computed number of basepairs to .

The fourth term is the most complex. In this condition, there is an intermediate nucleotide , such that the structures produced from to and from to are optimal. After going through all positions, the optimal number of basepairs for a sequence of length corresponds to the value of

7.4 Destabilizing energies

This nearest neighbor model so far only considers the stabilizing energy (negative terms) of a base pair stem, and does not consider the destabilizing energy of structures such as hairpin loops that contributed positive, destabilizing energy to a structure. The Nussinov Algorithm only considers the base pairs, and does not consider that not all basepairs contribute the same energy, as seen in the tables above. Furthermore, it does not take into account the destabilizing energy of the loops, bulges, and other unpaired structures. For example, this image shows some of the common forms of destabilizing energies in RNA structures. Fortunately, there are other algorithms that do compute the minimum free energy including the contribution of these types of structures, and have already been implemented in relatively easy to use software.

Table 7.3: The destabilizing energy of hairpin loops with closing C·G basepair, with A·U basepair, and for a bulge loop in kcal/mol.

length (nt)

C·G-closed loop

A·U-closed loop

bulge loop

1

999

999

2.8

2

999

999

3.9

3

8.4

8.0

4.5

4

5.9

7.5

5.0

5

4.1

6.9

5.2

6

4.3

6.4

5.3

7

4.5

6.6

5.5

8

4.6

6.8

5.6

9

4.8

6.9

5.7

10

4.9

7.0

5.8

12

5.0

7.1

5.9

14

5.2

7.3

6.1

16

5.3

7.4

6.2

18

5.4

7.4

6.3

20

5.5

7.6

6.4

25

5.7

7.7

6.5

30

5.9

7.9

6.7

The destabilizing energy of a hairpin can be summarized in Table 7.3[76]. For a loop of length , the destabilizing energy caused by the loop can be approximated by:

(7.3)

but this equation is only valid for . Other forms of destabilizing energies include hairpin loops, internal loops, bulges, and multiloops (Figure 7.6). In some cases, multiloops can gain stability from coaxial stacking, when two stems on opposite sides of a multiloop share a common axis (Figure 7.6D) [77].

Figure 7.6: Some types of destabilizing energies include A hairpin loops B internal loops C bulges and D multiloops. Figure adapted from Danaee et. al

7.4.1 RNAfold

is software that comes as part of the Vienna RNA package [79]. To compute a structure for a sequence, simply input the sequence data as a FASTA file using . The GNU/Linux program simply prints the file to the screen, but upon piping it into a file will run the program. For example, consider the microRNA hairpin sequence for miR-1 in humans:

$ cat mir-1.fasta | RNAfold --noLP -p
>hsa-mir-1-1 MI0000651
UGGGAAACAUACUUCUUUAUAUGCCCAUAUGGACCUGCUAAGCUAUGGAAUGUAAAGAAGUAUGUAUCUCA
(((((.((((((((((((((((..((((((((.....)))...))))).)))))))))))))))).))))) (-29.10)
{((((.((((((((((((((((..((((((((.....))}...))))).)))))))))))))))).))))} [-29.84]
(((((.((((((((((((((((..((((((((.....)))...))))).)))))))))))))))).))))) {-29.10 d=2.64}
frequency of mfe structure in ensemble 0.299677; ensemble diversity 3.99

The result is a prediction of the minimum free energy structure for the molecule (in dot-bracket notation), and the minimum free energy of the structure (in kcal/mol). The dot-bracket notation can be understood as containing matching parentheses that correspond to basepairs. Each left parenthesis has a corresponding right parenthesis, and their positions match which nucleotides are predicted to be paired. Dots correspond to unpaired nucleotides. In addition, by default, two postscript image files are created. The first is a secondary structure image ending in and the second is a dotplot image ending in . The resulting image for this sequence is shown in Figure 7.7.

Figure 7.7:The structure and dot plot for a microRNA produced from RNAfold.

The secondary structure image is a representation of how the structure of the RNA would look in a 2D plane. In a real 3D tertiary structure, there would be helical turns in the stem region of the structure. The dotplot figure shows a matrix where below the main diagonal demonstrates the base pairs present in the minimum free energy structure. Each term of the matrix corresponds to a pair of positions along the sequence. If the two nucleotides at positions and are paired, then that cell of the matrix is filled. Above the main diagonal shows the probability of each basepair. You will notice some faint dots around the main larger ones, corresponding to basepairs that could form in suboptimal configurations.

Next, let’s consider computing the structure of a tRNA. Suppose that the FASTA file contains a tRNA for the yeast species Saccharomyces cerevisiae. We can compute the structure as before with the command:

$ cat trna.fasta | RNAfold --noLP -p
>Saccharomyces_cerevisiae_chrXVI.trna12-PheGAA (622631-622540) Phe (GAA) 92 bp Sc: 69.79
GCGGAUUUAGCUCAGUUGGGAGAGCGCCAGACUGAAGAAAAAACUUCGGUCAAGUUAUCUGGAGGUCCUGUGUUCGAUCCACAGAAUUCGCA
((((((...((((........))))....((((((((......))))))))......((((..((((........))))..)))))))))). (-30.40)
((((((,..((((........)))),{..((((((((......))))))))..},..{(((..((((........))))..)))))))))). [-31.39]
((((((...((((........))))....((((((((......))))))))......((((..((((........))))..)))))))))). {-30.40 d=3.52}
frequency of mfe structure in ensemble 0.199156; ensemble diversity 5.29

The result is something that looks kind of like a cloverleaf structure. An important fact about tRNAs is that they are chemically modified in the cell, which results in a different structure than one might predict. In the end, the specific structure formed by a tRNA probably looks different than the predictions here due to these modifications. The images produced by are depicted in Figure 7.8.

Figure 7.8: The structure and dotplot for a tRNA produced from RNAfold.

7.5 Lab 8: RNA Structure

7.5.1 RNA Secondary Structure Prediction with RNAfold

Copy and paste the following hairpin sequence to a FASTA file called “hairpin.fasta”:

>hairpin
ACAUACUUCUUUAUAUCCAUAUGAACCUGCUAAGCUAUGGAUGUAAAGAAGUAUGU

We can compute the secondary structure of this hairpin with RNAfold:

$ cat hairpin.fasta | RNAfold

Comment on the hairpin structure dot-bracket sequence and the energy (given in parentheses). How big is the hairpin loop in terms of number of nucleotides? How many nucleotides are paired?

7.5.2 RNA Duplex Prediction with RNAduplex

Now let’s use RNAduplex to compute the energy of the RNA duplex formed by two RNA strands. Create a FASTA file called with the following 5′ and 3′ “arms” of the hairpin:

>five
ACAUACUUCUUUAUAUCCAUA
>three
UAUGGAUGUAAAGAAGUAUGU

Now compute the duplex energy:

$ cat arms.fasta | RNAduplex

Given the differences in energies between the duplex and the hairpin, How much destabilizing energy would you say the hairpin loop give?

7.5.3 The Destabilizing Energy of a Hairpin Loop

Let’s make the calculation of the hairpin loop’s destabilizing energy more rigorous. Let’s try increasing the length of the hairpin loop by inserting nucleotides to the hairpin loop sequence above in increments of two. Note: You should choose nucleotides that aren’t self-complementary to avoid creating more base pairs. Also, check the dot-bracket sequence to make sure that no new base pairs are added.
Make a table with the difference in the hairpin energies and the duplex energies for the different lengths. How does the difference in energy change with the length of the loop? How does this compare to equation 7.3? Try plotting the difference . This difference should be due in part to the positive energy “cost” associated with the loop.

Chapter 8: Proteins

Many transcripts encode Proteins, which are sequences of amino acids that fold into specific structures responsible for their functions. To a great extent, the sequence of the protein determines its structure, which in turn determines its function. Therefore, identifying similar protein sequences can suggest structures and functions. The field of protein bioinformatics is certainly extensive, and we will only be able to cover a small portion to introduce some concepts that may be useful.

To go from mRNA to protein, we need an open reading frame (ORF) defined by a start codon in frame with a stop codon. Typical mRNAs have many possible ORFs, with other signals, such as the Kozak sequence, determining which ORF is used. In most cases, the longest ORF is the coding DNA sequence (CDS), or the region that is actually translated. In general it is important to distinguish ORFs, which are quite common, from the CDS.

Once we have the CDS sequence, it is easy to compute the corresponding peptide sequence with Biopython, as we have seen in section 1.2.2.

8.1 Protein Alignment

Our previous discussion of alignment applies to proteins as well as nucleic acid sequences. However, with 20 amino acids, our scoring system can be updated to account for chemical similarity, and the fact that some amino acids are more likely to mutated to others.
Typically, when aligning two sequences, the assumption is that the two sequences have a common evolutionary origin, and the association between single characters in the alignment correspond to the evolutionary history of those characters. In order to optimize such an alignment, we need to define a scoring system used to compare two characters.

Let’s define some probabilistic models of sequence alignments so that we can describe the scoring system with such a formalism. Much of what we will describe here was developed for protein sequences, but the formalism could be used more generally. First, let’s formulate the probability of the two sequences given a random model. That is, the two sequences are just two random strings of characters that have been independently generated from a generative probabilistic model.

Consider an ungapped alignment of the sequences and :

x[1]x[2]x[3]x[4] ... x[n]
| | | | |
y[1]y[2]y[3]y[4] ... y[n]

Let’s work out a scoring system for this, and then later introduce gaps.
For a random model, each of the characters and at each position is just randomly selected according to their frequency in the database

Alternatively, if the sequences have a common ancestor, we could also define a probability as describing the frequency of occurrence of substitutions in the database between two characters and .

If we combine our two probabilistic models, we can compute a odds ratio

and by taking a logarithm, we can define a score as a log-odds ratio, which is more convenient for our purpose because it turns the product into a sum

The last terms give us our similarity matrix terms defined as

8.1.1 PAM Matrices

A Point Accepted Mutation Matrix (PAM Matrix), is was the first systematically defined matrix for scoring the similarity of peptide sequences. This scoring system was created by Margaret Dayhoff in 1978 [80]. First we define a “Mutation Matrix” that describes the probability of mutating from to .

The Mutation matrix is time reversible, meaning

Under such a time-reversible model, the observed frequency of mutating to is equal to the observed frequency of mutating to . This is a common assumption, particularly because we can’t observe, and can only infer, which character is ancestral. We will see below that we can compute the most likely or most parsimonious ancestral character, but the databases that were originally used to build these scoring matrices just used alignments; hence, we can only observe substitutions.
If our database has occurrences of amino acid , then we define our normalized frequency of that amino acid as

where is the total number of amino acids in the database. Because , the probabilities are normalized.

We can’t measure how much time has transpired between mutations in the sequences of our database. Therefore, we have to talk about how many substitutions have taken place as a percentage.
These mutation matrices are scaled such that there is amino acid change:

Therefore, this equation scales the mutation rates, and gives us restrictions on the protein alignments that we can consider in our database.
The actual PAM matrix that is used in scoring is defined as

To get more distant relationships, we can raise the matrix to a power.

Consider :
The result of this matrix product gives the probability of mutating from to in twice the amount of time defined by the matrix , and considering all intermediate amino acids .

PAM matrices of a higher order, are then defined by taking the log of this higher-order power of the original PAM matrix.

Therefore, the larger is, it describes more distant homology. For example, is an option in NCBI BLAST to specify more distance relationships.

8.1.2 BLOSUM Matrices

Another very common set of protein substitution matrices is BLOSUM [81]. BLOSUM is actually a series of matrices BLOSUM-x, where they are built from alignments that are at least identical.
These alignments come from the Blocks Database , which is still available at .
Consider a single column of one of these Block alignments

In this example, there are 3, 2 s, and 2 s. We wish to compute the number of pairs for each pair of amino acids and . When the amino acid is the same, for example is the number of pairs that can be selected from items, or , where is the number of occurrences of the amino acid . This expression can be expanded in terms of factorials as follows.

(8.1)

When the amino acids are different, such as , the value is computed as the product of the number of occurrences of each amino acid, like . To avoid double counting, we can only consider cases when latex i \le j$.
These frequencies are normalized to make a matrix, with terms defined as

Where the second sum in the denominator just goes up to to count half of the symmetric matrix to avoid double-counting pairs.
We can get the normalized frequencies of each animo acid from these quantities by

where the term is because there is a probability of that a random amino acid selected from the pairs corresponding to the pairs is . Next we want to compute the probability of selecting the amino acids and by random.

To understand this equation, imagine selecting a pair of amino acids at random from the database. If they are the same, then the probability of that happening is just . If they are different, then there are two ways that could be selected: first then , and first then . Finally, the BLOSUM matrix would be computed as a log-likelihood ratio

In practice, the values for the BLOSUM matrices are rounded to integers because when they were created computer memory and computation was expensive, and this required less space than storing decimals and performing floating point arithmetic.

8.1.3 Karlin and Altschul Generalization

Karlin and Altschul formulated a generalized scoring matrix that is similar to both PAM matrices and BLOSUM matrices, defined as

where is a positive parameter that scales the matrix [82].

8.1.4 Biopython and Substitution Matrices

We can access values of a given substitution matrix using Biopython and the module . For example, the most commonly used substitution matrix and the default for NCBI protein BLAST is BLOSUM62. We can print the values of the matrix or access a particular term by the following

>>> from Bio.SubsMat import MatrixInfo
>>> S = MatrixInfo.blosum62
>>> S['Q','Q']
5
>>> S['W','Q']
-2

If you play around with this matrix (actually it is a dictionary that takes a pair of characters as keys) you will realize that some pairs are stored, but others are not and return errors when one tries to access them. For example, is stored, but not . To get around this, an additional function can be created:

>>> def getScore(a,b,S):
...     if (a,b) not in S:
...         return S[b,a]
...     return S[a,b]
...
>>> getScore('Q','W',S)
-2

8.2 Functional Annotation of Proteins

Computational methods can be used to infer protein function when a new protein is identified, but even if the inference is very strong, experiments still need to be done to validate the prediction. That said, there are many ways in which a hypothesis about protein function can be inferred, including sequence homology, sequence motifs such as “domains”, which could be considered a type of homology, structure-based function prediction, genomic co-location (gene clusters), gene expression data, in particular co-expression.

8.2.1 Protein Evolution and Homology

The first method that one might use to infer protein function is by sequence homology. Frequently, the more fundamental molecular functions of proteins are very well conserved. There are a number of mechanisms of protein evolution that one must consider when evaluating protein homology.

In DNA replication, there are often errors that cause major changes in proteins, which can be highly disruptive to the function of the protein, but perhaps in some cases may confer new function. First, is the event of gene fusion, where two genes may be merged into one gene, combining different functional regions from the two ancestral genes. For example, a gene encoding a protein with a DNA-binding domain could fuse with a gene encoding a protein with a protein-binding domain, producing a new gene encoding a protein with both domains. In many cases, a gene fusion event is caused by a chromosomal translocation event, but such translocations can also cause other mutations in genes, such as loss of a portion of a gene. In addition, regions of chromosomes can be inverted during replication, resulting in gene fusions or deletions of portions of genes. Along the same lines, a chromosomal deletion, where a portion of a chromosome is deleted during DNA replication, can also lead to gene fusion events or loss of a portion of a gene.

Gene duplication events can happen as part of a chromosomal duplication, or as part of a local duplication. When this happens, multiple copies of a gene can be created. The two copies of the gene could be beneficial in some cases, or harmful in others. In these cases, one of the copies can be possibly silenced. When both copies remain, one of the copies can undergo less selective pressure, because the bulk of the work is carried out by one of them. In this case, the other copy can accumulate mutations over time, which could result in new function or a specialized function. This process where one copy accumulates mutations and ultimately carries out a specialized role is often called sub-functionalization of the protein [83, 84].

8.2.2 Protein Domains

A protein domain is a conserved part of the protein that has a distinct structure and function. Often domains occur as highly conserved blocks of conservation that can be clearly observed in a sequence alignment. Often distantly related genes have domains conserved beyond the rest of the gene. Domains also will fold and evolve independently of the rest of the protein, resulting in increased conservation over the domain region. The idea of a protein domain was introduced by Wetlaufer in 1973 through observing common stable units in X-ray crystallography studies [85].
A protein domain can be understood as similar to a motif, although some databases have more complex representations when necessary to describe a more complex pattern. To identify protein domains, one needs a sufficiently diverged set of proteins in order to accumulate enough mutations outside the domain so that the domain itself can be identified. However, with a set of proteins that is too diverged, once could even lose the conservation of the domain.

Many databases of protein domains exist, including Pfam, PANTHER, PROSITE, and Interpro.

Pfam Database

Pfam is a database of proteins, protein families [86], and domains, and is available at http://pfam.sanger.ac.uk. Pfam domains are a commonly used annotation of protein domains that are relatively easy to use.
For example, one can use the program HMMer (pronounced “hammer”) to scan Pfam domains downloaded as an hmm file. The command would be something like this:

$ hmmscan --domtblout domainTable.txt Pfam-A.hmm proteins.fasta

To produce a table similar to a BLAST output with locations of matches, e-values, and p-values.

8.3 Secondary Structure prediction

To get a full view of a protein’s structure, we’ll need a way to compute its tertiary structure, which is extremely difficult. In practice, we have to rely on experimental evidence such as X-ray crystallographic studies or NMR structures. Some studies will use a known structure of a protein and apply it to a homologous protein. Secondary structure for proteins can be computed, however, with reasonable accuracy. The secondary structure of a protein is a list of the positions corresponding to alpha helices and beta strands.
The software can predict secondary structural features and provide a confidence score for each position [87]. For example, the program can be run on the command line with a command like:

$ jnet -p human_catalase.fasta

to produce an output file indicating the locations of these regions:

Length = 527 Homologues = 1
RES : MADSRDPASDQMQHWKEQRAAQKADVLTTGAGNPVGDKLNVITVGPRGPLLVQDVVFTDEMAHFD
ALIGN : ---------HHHHHHHHHHHHHH--EEE---------EEEEEEE-----EEEEEEEE-----H--
CONF : 88888887525778889988874122440687775752678883477752575455201000013
FINAL : ---------HHHHHHHHHHHHHH--EEE---------EEEEEEE-----EEEEEEEE--------

where the value is the confidence, a number from to , that indicates the quality of the prediction. In this representation, the represent alpha helices, and the indicate beta-strands. Here we have run the program with one protein sequence, but it could be run with multiple sequences for greater accuracy.

8.4 Gene Ontology

Gene Ontology (GO) is a database of controlled vocabulary terms that describe gene/protein function [88]. A valuable web resource for GO terms and related data is available at their website (http://geneontology.org). These terms form a hierarchy, where certain higher-level terms “contain” the lower-level terms in conceptual scope.

Figure 8.1: GO Terms are organized into a hierarchy of concepts, with the lower terms representing more refined functions.

There are many levels of evidence for the assignment of a GO term, indicated by “evidence codes”, given by a symbol of a few letters. For example, indicates that the term was “Inferred by experiment” Others such as for “traceable author statement” indicate that there is a statement in a publication making the claim of the association.

Table 8.1: Table of some GO term evidence codes and their meaning.

TAS

Traceable Author Statement

EXP

Inferred from Experiment

IDA

Inferred from Direct assay

IPI

Inferred from Physical Interaction

IMP

Inferred from Mutant Phenotype

IGI

Inferred from Genetic Interaction

IEP

Inferred from Expression Pattern

ISS

Inferred from Sequence or Structural Similarity

ISO

Inferred from Sequence Othology

Figure 8.2: GO Terms are organized as a directed acyclic graph (DAG), which can connect nodes lower in the graph, but does not contain cycles linking back up further up the tree.

8.4.1 Multiple hypothesis testing

Frequently when performing a study on the effect of a mutation or other treatment, one looks at the GO terms of the genes that are significantly differentially expressed due to the experimental condition compared to a control. One evaluates whether the enrichment of a particular GO term in the data exceeds that of the expected rate of this term occurring in randomly sampled gene sets of the same size. This significance is typically computed from a p-value describing the enrichment.

However, in doing so, one needs to perform a multiple test correction, to account for the fact that many hypotheses were tested. For example, if you found a one-in-a-million event after examining one million events, you might not be surprised. The simplest multiple-test correction is the Bonferroni correction [89], where when testing hypotheses, we need to test for significance at level by examples satisfying the relationship:

(8.2)

Another method is the Benjamini-Hochberg procedure, which is a little bit more complex [90]. To do this procedure, one needs to sort the list of p-values in ascending order such that such that the subscript given by indicates the p-value with rank . This procedure gives us a list of significant hypotheses with a rate of false-discovery defined by the imposed , a number between and . After the sorting, all p-values with a rank less than the maximum rank such that:

(8.3)

where defines the false discovery rate, are deemed significant.

Figure 8.3: Significant p-values are determined by the largest rank such that . Note that while the p-value is sorted so that it is monotonically increasing as a function of , the term is not necessarily always increasing.

8.5 Lab 9: Proteins

In this lab, we learn some methods for predicting properties of protein sequences, such as domains and secondary structure.

8.5.1 Scanning for domains

First, let’s look at protein domains from Pfam. First, we’ll need to download a domain file from Pfam with the following command:

$ wget ftp://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/Pfam-A.hmm.gz}

This file contains a set of HMM models for protein domains defined for the Pfam database. To use this file, we must first unzip the file with gunzip.

$ gunzip Pfam-A.hmm.gz

This will produce the unzipped file , which can be used to scan for these domains using the HMMer software program hmmscan as part of the HMMer package. However, we still need to index this domain file. This can be done with the HMMer software.

$ hmmpress Pfam-A.hmm

As we have seen in this course, there are many ways to download a protein sequence. Since we will work on protein secondary structure prediction, let’s consider downloading something from PDB, the Protein Databank.
The PDZ domain is a structural domain involved in signaling in many organisms ranging from bacteria to animals. The structure consists of roughly 5 beta-strands and 2 alpha-helices, as demonstrated by this image of the tertiary structure [91].

Figure 8.4: Tertiary structure image created using WebGL at http://www.rcsb.org/3d-view/2NB4/1

You can see more detail on this structure at PDB (http://www.rcsb.org/3d-view/2NB4/1). One useful piece of data in this database entry is the secondary structure information found at the PDB database here: http://www.rcsb.org/pdb/explore/remediatedSequence.do?structureId=2NB4.
The sequence can be obtained in the upper-right hand corner of the PDB database page under “Download Files” and then under “FASTA sequence”. Here is the PDZ domain sequence:

>2NB4
PLTRPYLGFRVAVGRDSSGCTTLSIQEVTQTYTGSNG
GADLMGPAFAAGLRVGDQLVRFAGYTVTELAAFNTVV
ARHVRPSASIPVVFSRDGVVMSATIVVGELE

Download the FASTA file here PDZ sequences and compute a multiple sequence alignment for the sequences to see if it improves the domain identification.
Now we can save this sequence to a file called . We can then run to see if there are any domains in this sequence:

$ hmmscan --domtblout domains_2NB4.txt Pfam-A.hmm pdz.fasta

Does this output make sense? you can view it with the following command, using “less -S” to turn off word-wrapping.

$ less -S domains_2NB4.txt

The first few columns define which domain was found. As expected the only domains are PDZ. What is unexpected, perhaps, is the positions of the two Pfam domains pertaining to PDZ (defined by the “ali coord” entry and “from” and “to” columns.) don’t always start at the same position, suggesting that to some degree the boundaries are a bit “fuzzy” at times.

Also, note that another column provides the e-value, which is interpreted similarly to BLAST, as the expected number of occurrences of this domain by chance
The secondary structure of the protein can be computed with the following command:

$ jnet -p pdz.fasta

This produces a file format that gives the sequence and a letter designation for whether the region is an alpha-helix, designated by an “H”, or a beta-strand, designated by an “E”. How closely does the predicted structure correspond to the true structure seen here: https://www.rcsb.org/pdb/explore/remediatedSequence.do?structureId=2NB4?

You can run this program on an alignment, but you need to remove gaps from the multiple sequence alignment. To do this, first compute an alignment. I used the PDZ file above

$ clustalw2 -infile=PDZ_sequences.fasta -type=Protein -outfile=PDZ_sequences.aln

Next, you need to convert to FASTA alignment format and remove gaps. You can covert to FASTA using AlignIO (see Lab 5, section 4.4), but this will not remove gaps. To remove gaps, very industrious students might try to write a python script to do this. Note that you need to remove gaps a column of the alignment. Meaning, if one sequence has a gap, that specific position of the alignment would be removed for all sequences. Alternatively, you an use the program to do this. You can remove the gaps and convert to FASTA in one step:

$ trimal -in PDZ_sequences.aln -fasta -nogaps > PDZ_sequences.fa

Next, you can run on the resulting alignment:

$ jnet -p PDZ_sequences.fa

The problem here is whether or not the alignment with gaps removed is informative. If your sequences are too far away, a lot of gaps could be removed, rendering the positions of the resulting gap-free alignment different from the original sequences.

Chapter 9: Gene Regulation

Genes can be regulated in several different ways, and be regulated transcriptionally (at the gene-level), post-transcriptionally (at the mRNA level), translationally, and post-translationally. Splicing of mRNA can also be regulated and there are many other inputs to regulation still.

Clearly, transcriptional regulation and post-transcriptional regulation are best suited for study by nucleic acid sequence bioinformatics because techniques such as ChIP-seq and RNA-seq can be used to measure the inputs and outputs of these modes of regulation. Here we will focus on regulation by transcription factors using ChIP-seq and regulation by microRNAs using small RNA-seq and RNA-seq.

9.1 Transcription Factors and ChIP-seq

Transcription factors (TFs) are proteins that bind to DNA and are capable of regulating the transcription of genes. TFs can be activators, which activate transcription of the genes they regulate, or repressors, which block the transcription of the genes they regulate.

9.1.1 ChIP-seq Peak identification

ChIP-seq measures the genome-wide binding of a TF by sequencing DNA fragments that are bound to the TF in a controlled experiment [92–94]. Antibodies that recognize epitopes in the protein of interest are used to “pull-down” the protein. When cross-linking is used, such as by adding formaldehyde, everything sticks together, holding the protein to the DNA that it is bound to. After the antibodies and associated chromatin are precipitated out from the solution and isolated, they can be reverse-crosslinked, and then a DNA-extraction is performed. The fragments of DNA are then sequenced using deep sequencing. Despite best efforts to just isolate fragments directly bound the the TF of interest, a lot of spurious reads make it through the process as well. Therefore, we need to sequence a control sample without the antibody, either using the immunogloblulin (IgG) or just naked DNA (input), so that we can look for statistical enrichment of the reads that align to a particular genomic location from the ChIP-Sample compared to this control sample.

Figure 9.1: ChIP-seq peaks are characterized by a local enrichment of ChIP-seq reads compared to control, and are observed over the binding sites for the TF under investigation.

Short read alignment with

Once we have the the reads from the high-throughput sequencer, typically in a FASTQ file, we need to align to the genome. This can be achieved with and a bowtie index of the genome [95].

$ bowtie -m 1 -S -q hg38 TF_ChIPSeq.fastq TF_ChIPSeq_hg38_bowtie.sam
$ bowtie -m 1 -S -q hg38 DNA_input.fastq DNA_input_hg38_bowtie.sam

In this particular example, we are aligning the fastq files such that we are only retaining alignments that have one hit to the genome. This is done with the “” option. If one wanted to retain only reads that map to at most hits, then we would use “” where is some integer. Recall that is the “base” of the bowtie index of the genome. It is actually defined by a set of files, as described in section 6.4.3. While on the subject of bowtie indexes, we can either download them from the web, such as from Illumina’s iGenome datasets (https://support.illumina.com/sequencing/sequencing_software/igenome.html), or we can create one from a genome FASTA file with the following command:

$ bowtie-build -f hg38_genome.fasta hg38

Identifying Peaks with Peak Seq

The software PeakSeq was one of the first pieces of software used to identify ChIP-seq peaks using a statistical model [96]. The program divides the genome into segments of length such the number of reads aligning to each genomic segment from the ChIP-seq sample is , along with the number of reads from the control sample . The statistical significance is assessed using a binomial model. That is, the reads for a particular segment is assumed to be Bernoulli distributed random variables that take the value of with a probability and a value of with a probability . Then, it is determined whether the number of ChIP-seq reads is larger than expected due to chance.

First, a best-fit line of slope through the data for all segments is computed. The data that are fit is the x-axis has the value of and the y-axis has the value , and each data point is a segment. The best-fit line through the data gives a scale factor that can be used to scale the number of reads for the control sample, to be compared to the number of reads for the ChIP sample. When the scale factor is used, the binomial model uses a probability of when comparing to .

The p-value for a binomial variable is computed as one minus the cumulative distribution function for the binomial distribution , which gives the probability of observing less than or equal to successes in trials of a Bernoulli random variable with probability . The p-value is computed by

Each p-value is used in a Benjamini-Hochberg correction described in section 8.4.1.

Identifying Peaks with MACS

Figure 9.2: A. ChIP-seq reads are sequenced from fragments bound to the TF. B. Reads are ligated to 5′ adapters. C. Reads are sequenced from the 5′ end due to how the adapters are attached. D. This results in a bimodal distribution of reads around the location of the TF binding site.

is a software that can find ChIP-seqs peaks similar to that of [97]. The difference is that models the length distribution of the ChIP-DNA fragments, which helps it locate binding sites more accurately. There are different versions of , but probably the most commonly used is , which can be run to identify peaks with commands like

$ macs14 -t TF_ChIP_hg38.bam -c DNA_input_hg38.bam -f BAM -n TF_ChIP_vs_DNA_input -g hs -p 0.01

This program supports other input formats other than , but this is probably the most common. The command gives it the effective genome size. This can either be a number, or a two letter abbreviation for the more common model systems. Finally, the specifies a p-value threshold defining what peaks to include.

Peak Modeling with GPS

Other approaches, such as is implemented in Genome Positioning System (GPS) model the distribution of reads around the motif.

9.1.2 Chromatin Modifications

While transcription factors tend to have sharp, well-defined peaks, ChIP-seq performed on antibodies that recognize chromatin modifications tend to be enriched for broad domains that can extend for many megabases. Certain chromatin marks, such as H3K27me3 can be spread out over very large polycomb-repressed domains. Some programs, such as SICER, are specially designed to identify these kinds of broad domains [98]. Other programs such as chromHMM will train an HMM on a set of samples, and segment the genome into multiple genomic states [99].

9.2 MicroRNA regulation and Small RNA-seq

MicroRNAs are processed from endogenous primary transcripts called “pri-miRs” that can either be transcribed as a noncoding RNA or as the the intron of a protein-coding gene. The base of one or more specific hairpin sequences in the pri-miR are processed out by Drosha, resulting in a microRNA hairpin precursor called a “pre-miR”. The pre-miR is exported out of the nucleus, and then the loop of the hairpin is removed by Dicer, which is an RNase III enzyme. The result is an RNA duplex of two mature microRNA sequences. The resulting mature microRNAs (miRs) are small (18-25nt) RNA molecules that regulate gene expression throughout the eukaryotes. The word for microRNAs is typically written in lower-case as in “microRNAs” except when at the beginning of the sentence. The regulation by miRs is post-transcriptional, and can either lead to degradation of the transcript it regulates, or translational inhibition. MicroRNAs operate through complementary sequences in the 3′ UTRs of transcripts that match the “seed” sequence, defined as positions 2-8 of the mature microRNA. That is, the seed sequence of the miR is complementary to the target site in the 3′ UTR of the gene it regulates.

9.2.1 Read abundance

The first thing we might be interested in knowing about a microRNA is the extent to which its mature miR is expressed. This can be done with small RNA-seq, which is RNA-seq performed on size-selected fragments. When aligning microRNA reads, we first have to trim the 3′ adapters from the FASTQ files as explained in section 6.3.1.

$ bowtie -m 50 -l 20 -n 2 -S -q hg38 smallRNASeq.fastq smallRNA_hg38_bowtie.sam

To compare across samples, it is best to use a normalized expression count, such as Reads Per Million (RPM), which allows expression levels to be compared between samples of varying depth. For microRNA with the number of reads mapped to it being , the RPM is computed by

where is the number of reads in the sample. This equation is similar to RPKM from equation 6.2, but doesn’t normalize per length.
This is because microRNAs always have reads spanning the full locus defining the mature microRNA when mapped to the genome, so being a few nucleotides longer doesn’t give one mature miR annotation an advantage over a shorter one. Also, normalizing by length can artificially inflate differences between a short miR of length and a longer miR of length (both lengths exist in microRNA databases) even when they both have the same number of reads.

9.2.2 Argonaute CLIP-seq

To directly measure functional mature microRNAs that are actually incorporated into the RISC complex, one can perform AGO-IP, or an immunoprecipitation with an antibody for Argonaute (AGO) [100].

9.2.3 MicroRNA target prediction

There are many databases that contain information on microRNA target sites, such as http://www.microrna.org/microrna/getGeneForm.do, which allows you to retrieve alignments of 3′ UTRs and associated target sites.

The software allows one to predict target sites from a set of miR sequences and a set of UTR sequences. The format is pretty simple, but because of the information used, a simple FASTA file does not suffice. An additional column is required for the taxa ID of the species used for each example. For the microRNAs, the file format contains a tab-delimited list of miR-Family names, seed sequences, and taxa ID:

let-7/98    GAGGUAG    10090
let-7/98    GAGGUAG    10116
let-7/98    GAGGUAG    9031
let-7/98    GAGGUAG    9606
let-7/98    GAGGUAG    9615
miR-18    AAGGUGC    10090
miR-18    AAGGUGC    10116
miR-18    AAGGUGC    9031
miR-18    AAGGUGC    9606
miR-18    AAGGUGC    9615
miR-1/206    GGAAUGU    10090
miR-1/206    GGAAUGU    10116
miR-1/206    GGAAUGU    9031
miR-1/206    GGAAUGU    9606
miR-1/206    GGAAUGU    9615

For example, in this example is the taxa ID for human. These taxa IDs are designated by NCBI Taxonomy, at http://www.ncbi.nlm.nih.gov/taxonomy. The 3′ UTRs are kept in a similar file that contains sequences from a multiple sequence alignment of the UTR sequences. The UTR file also has taxa IDs for each sequence:

CDC2L6 10090 CCCACUCCCU---CU-------------------------GCUUGGCCUUGGA-----------CUCCAGCAGGGUGGUAUUUGUGUUAC---AAAGACCCCCAG...
CDC2L6 10116 CCCACUGCCC------------------------------GCUUGGCCUCGGGGAGCACAGAGCCGGCGGCAGGGUGGUUUCUGCAAUAC---AAAGAACCCCAC...
CDC2L6 9031 GAGGCUUAACAGACUGCACUGAAAAAGGAAUGGAUUAAAAGCCAGA----AGA-----------CUCCAGCAAUAUGAAGUUCGUGUUGAUGAGAAGAACCAAA-...
CDC2L6 9606 CCAGCUCCCG---UUGGGCCAGGCCAG--------CCCAGCCCAGAGCACAGG-----------CUCCAGCAAUAUG---UCUGCAUUGA---AAAGAACCAAAA...
CDC2L6 9615 CCG-CCACGG------------------------------CCCAGGGCACAGA-----------CUCCAGCAAUAUGAUGUCCGCAUUGA---CAAGAACCAAA-...
FNDC3A 10090 AAAUAUAACUUUAUUUUUAA-CACU--GUAUUACAUUUUAUUUUGUCAUGUACU-AAAAUUAUUUCUGUA---UUGCUUUUAC--AAAAUGGUGGCAUUUAGCAC...
FNDC3A 10116 AAAUAUAACUUUAUUUUUAA-CAC----UAUUACAUUUUAUUUUGUCAUGUACU-AAAAUUAUUUCUGUA---UUGCUUUUAC--AAAAUGGUGGCAUUUAGCAC...
FNDC3A 9031 AGAAACAGAUUUAUUUUAGAAUGCUGCCCAUUACAUUUUACUUUCUCAUAAUCUAAAAAAAAAUUCUGUUCUCUUGCUUUUACAAA-AACA--GGCAUUUAGCAC...
FNDC3A 9606 AAAUAUAACUUUAUUUUUUA--ACU--CUAUUACAUUUUAUUUUGUCAUGUACU-AAAAUUAUUUCUGUA---UUGCUUUUAUAAAAAACAGUGGCAUUUAGCAC...
FNDC3A 9615 AAAUAUAACUUUAUUUUUAA-UACU--GUAUUACAUUUUAUUUUGUCAUGUACU-AAAAUUAUUUCUGUA---UUGCUUUUAC-AAAAACAGUGGCAUUUAGCAC...

These two examples are in fact the first few lines if the sample files given with the TargetScan software. TargetScan is a database of conserved microRNA target sites, but also is a set of software that can be used to predict conserved target sites from two files like these.

$ perl targetscan_70.pl miR_Family_info_sample.txt UTR_Sequences_sample.txt targetscan_output.txt

9.3 Regulatory Networks

Regulatory networks are described by a graph representing the regulatory connections between genes. In some studies, just a small number of connections are described. Other genome-wide studies consider the larger network of connections as a whole. For these genome-wide approaches, the statistical trends in the data become the focus of the investigation.

9.4 Lab 10: ChIP-seq

The data for this project was taken from the paper Zhou et al. 2011, entitled “Integrated approaches reveal determinants of genome-wide binding and function of the transcription factor Pho4” [101]. This paper is studying the binding of the transcription factor Pho4, which binds to a sequence-specific motif described in the paper. We will be looking at the data corresponding to Pho4 under the No Pi conditions.
The raw data for this project is available at this GEO entry: GSE29506, but I have prepared the fastq files for download below.
Step 1: download the script from this directory with wget.

$ wget http://als2077a-unix1.science.oregonstate.edu/Pho4/bedToFasta.pl
$ wget http://als2077a-unix1.science.oregonstate.edu/Pho4/sce_INPUT_NoPi_ChIPSeq.fastq
$ wget http://als2077a-unix1.science.oregonstate.edu/Pho4/sce_Pho4_NoPi_ChIPSeq.fastq
$ wget http://als2077a-unix1.science.oregonstate.edu/sce/sce_R64_1_1.fa

You will download a perl script, 2 FASTQ files, and a FASTA file for the yeast genome. The other data you will use was downloaded in Lab 8 (section 7.5). It would be a good idea to check the file sizes with “” to confirm that the download was complete.

$ ls -l *fastq
-rw-r--r-- 1 dhendrix dhendrix 2817509732 Mar 10 23:23 sce_INPUT_NoPi_ChIPSeq.fastq
-rw-r--r-- 1 dhendrix dhendrix 746715226 Mar 10 23:31 sce_Pho4_NoPi_ChIPSeq.fastq

Step 2: Align the fastq files to the genome: You need to align both the ChIP-seq data and the INPUT (control) data.
To do this, we first need to create an index to the genome with this command:

$ bowtie-build sce_R64_1_1.fa sce_R64_1_1

Next, align the FASTQ files using bowtie and the following command.

$ bowtie -m 1 -S -q sce_R64_1_1 sce_Pho4_NoPi_ChIPSeq.fastq sce_Pho4_NoPi_ChIPSeq.sam
$ bowtie -m 1 -S -q sce_R64_1_1 sce_INPUT_NoPi_ChIPSeq.fastq sce_INPUT_NoPi_ChIPSeq.sam

The option specified by ensures that only reads with exactly one hit to the genome will be kept.
Step 3: Find the ChIP-seq peaks with macs. Note that we are using a stringent p-value threshold here to pull out only the strongest peaks.

$ macs14 -t sce_Pho4_NoPi_ChIPSeq.sam -c sce_INPUT_NoPi_ChIPSeq.sam -n Pho4_vs_INPUT -g 1.2e7 -f SAM -p 1e-10 >& macs14.err

Macs will produce both peak bed files, and summit bed files. As you may have guessed, the peaks are the full enriched regions, and the summits are the tip of the peaks, and only constitute a small window of 2bp.
Step 4: Use the bed file of the summits to create a FASTA file of of +/- 30bp around the summits:

$ perl Scripts/bedToFasta.pl Pho4_vs_INPUT_summits.bed sce_R64_1_1.fa 30

This perl script was in the directory of Pho4 stuff you downloaded above. Basically, it takes the genomic locations in the bed file, and prints them to a FASTA file. Note that in general, this program is run with the inputs as follows

$ perl bedToFasta.pl [bed file] [genome fasta] [sequence buffer]

Check out the perl script to see what it looks like. It’s similar to python in many ways, but different in many as well.
Step 5: Run MEME to find motifs in the peak region:

$ meme Pho4_vs_INPUT_summits.fasta -mod zoops -dna -nmotifs 3 -maxw 8 -revcomp

If you did all the steps correctly, the binding site for Pho4 should be one of the motifs in the MEME output.
Questions to answer for this project: How many peaks do you find at this threshold? You can count the number of lines in a bed file produced for the peaks with the following command:

$ wc -l bedfile

Which motif corresponds to the motif in the paper? Include the LOGO for the motif you found in your results. How many instances of the motif did you find? What is the information content of the motif? Even though this was run with “zoops”, you might still expect every peak to have an instance of the motif. How many of the peaks have an instance of the motif?

1

Appendix A: Mathematical Preliminaries

Some students with limited mathematical background may benefit from a brief introduction to some common mathematical notation frequently used in this text.

An efficient way to write a sum of many variables is with the sigma notation . Consider adding up the values . This can be compactly represented as:

Similarly, if we want to represent a product of many variables, we can use the following notation:

An example of such a product function is “factorial”, . This quantity is the produce of the integers from to , and represents the number of possible arrangements of distinct objects.

Logarithms are an important function to know in bioinformatics as they are commonly used in scoring systems. Logarithms are powerful because they have the following useful algebraic property (among others):

We can combine the properties of our sums and products and logs, with the following equation:

Another very useful mathematical construct is the Kronecker Delta function, and can be used for counting. It is defined as follows:

2

Appendix B: Probability

B.0.1 Probability Distributions

Probability Mass Function

A probability mass function for a discrete random variable is defined such that is the probability of the variable taking the value . These distributions are normalized, meaning summing over all possible values will equal . In other words:

Probability Density Functions

For continuous random variables , a Probability Density Function gives the value corresponding to the likelihood of the random variable taking on that value, and probabilities are computed from

Integrating over all values of will equal .

B.0.2 Conditional Probability

In probability, we often want to talk about conditional probability, which gives us the probability of an event given that we know that another event has occurred.

Bayes’ Theorem tells us that:

We can also expand a probability of an event into the conditional probabilities of each other condition, times the probability of these other events.

This expression is often called the evidence. We can express the probability of given , which is called the posterior probability by

Also in this equation, is called the likelihood, and is called the prior.

B.0.3 Bernouli Distribution

One of the simplest probability distributions we should consider is the Bernoulli Distribution. It is essentially a single event, with a probability of success denoted , and a probability of failure of .

Example 1: A coin toss. Heads or tails? Typically for a fair coin, this will be such that .

Example 2: Selecting a single nucleotide from the genome at random. Is it purine (A or G) or is it pyrimidine (C or T)? For the human genome, the GC content is about , but it varies from chromosome to chromosome.

B.0.4 Binomial Distribution

Central to the Binomial distribution is the binomial coefficient The number of ways to have successes out of trials is . The probability of observing successes out of trials is . Since , and raised to the th power is still , we can see that these probabilities are normalized (sum to one) by the equation:

Example 1: Consider tossing a fair coin 100 times, and repeating this for 1000 trials. A histogram for , number of heads in tosses is shown in Figure B.1. The expected value of is , and the variance is .

Figure B.1: A histogram of the number of heads observed in 100 tosses of a fair coin, repeated 1000 times.

B.0.5 Multinomial Distribution

The multinomial distribution deals which situations when there are more than two possible outcomes (for example, DNA nucleotides). The multinomial coefficient, , in this example describes the number of ways to have a sequence of length , with the number of occurrences of A,C,G,and T to be and respectively.

The probability of a particular set of observed counts depends on the frequencies by the expression:

B.0.6 Poisson Distribution

The Poisson Distribution describes the number of observed events given an expected number of occurrences . Consider, for example, the probability of a red car driving past a particular street corner in your neighborhood. Its probability mass function is given by:

It has one parameter , which is also the expected value of and the variance ().
It can be shown that for fixed , and for , the Binomial distribution is equivalent to the Poisson distribution.

In many applications, the parameter for the Poisson distribution is treated as a rate. In this case, the expected value of is , so that describes the rate of occurrence of the event per unit time (or per unit distance, or whatever the application is). The probability distribution is given by:

Example: Consider modeling the number of mutations observed in a length of DNA .
Consider the important case when for the Poisson Distribution.

The probability that is not zero, is therefore:

(B.1)

Now consider being replaced by a random variable .

B.0.7 Exponential Distribution

The Exponential distribution is related to the Poisson by equation . It is the probability of an event occurring in a length . The cumulative distribution function is the same as this last equation:

And the probability density function is the derivative of the cumulative distribution function:

Note that whereas the Poisson distribution is a discrete distribution, the Exponential distribution is continuous and can take on any real value such that .

B.0.8 Normal Distribution

A very important distribution for dealing with data, the normal distribution models many natural processes. We already saw that the binomial distribution looks like the normal distribution for large . Indeed, the Central Limit Theorem, the sum of random variables approaches the normal distribution for large . Here is the p.d.f.:

B.0.9 Extreme Value Distribution

The maximum value in a set of random variables is also a random variable. This variable is described by the Extreme Value Distribution, also known as the Gumbel Distribution.

For , and with a location parameter , a location parameter , and a shape parameter . In the limit that , the distribution reduces to:

In many cases, for practical purposes, one makes the assumption that .

3

Bibliography

[1] F.H.C. Crick. On protein synthesis. Symp. Soc. Exp. Biol, XII:139–163, 1956. http://profiles.nlm.nih.gov/SC/B/B/F/T/_/scbbft.pdf.

[2] F. Crick. Central dogma of molecular biology. Nature, 227(5258):561–563, Aug 1970.

[3] David J Lipman and William R Pearson. Rapid and sensitive protein similarity searches. Science, 227(4693):1435–1441, 1985.

[4] Brent Ewing and Phil Green. Base-calling of automated sequencer traces using phred. ii. error probabilities. Genome research, 8(3):186–194, 1998.

[5] Michael J Guertin and John T Lis. Chromatin landscape dictates hsf binding to target dna elements. PLoS Genet, 6(9):e1001114, 2010.

[6] T Tatusova, M DiCuccio, A Badretdin, V Chetvernin, S Ciufo, and W Li. The ncbi handbook. 2013.

[7] James Ostell and Eric W Sayers. Dennis a. benson, ilene karsch-mizrachi, karen clark, david j. lipman[j]. Nucleic acids research, 1:6, 2011.

[8] Kim D Pruitt, Tatiana Tatusova, and Donna R Maglott. Ncbi reference sequences (refseq): a cu-rated non-redundant sequence database of genomes, transcripts and proteins. Nucleic acids research, 35(suppl 1):D61–D65, 2007.

[9] Donald E. Knuth, Jr. James H. Morris, and Vaughan R. Pratt. Fast pattern matching in strings. SIAMJ. Comput., 6(2):323–350, 1977.

[10] Michael Burrows and David J. Wheeler. A block sorting lossless data compression algorithm. Technical Report 124, Digital Equipment Corporation, May 1994.

[11] Andrew D Johnson. An extended iupac nomenclature code for polymorphic nucleic acids. Bioinformatics, 26(10):1386–1389, 2010.

[12] Thomas W Burke and James T Kadonaga. Drosophila tfiid binds to a conserved downstream basal promoter element that is present in many tata-box-deficient promoters. Genes & development, 10(6):711–724, 1996.

[13] John C Wootton and Scott Federhen. [33] analysis of compositionally biased regions in sequence databases. Methods in enzymology, 266:554–571, 1996.

[14] G. D. Stormo, T. D. Schneider, L. Gold, and A. Ehrenfeucht. Use of the ’Perceptron’ algorithm to distinguish translational initiation sites in E. coli. Nucleic Acids Res., 10(9):2997–3011, May 1982.

[15] Claude Elwood Shannon. A mathematical theory of communication. ACM SIGMOBILE Mobile Computing and Communications Review, 5(1):3–55, 2001.

[16] Solomon Kullback and Richard A Leibler. On information and sufficiency. The annals of mathematical statistics, 22(1):79–86, 1951.

[17] Peter JA Cock, Tiago Antao, Jeffrey T Chang, Brad A Chapman, Cymon J Cox, Andrew Dalke, Iddo Friedberg, Thomas Hamelryck, Frank Kauff, Bartek Wilczynski, et al. Biopython: freely available python tools for computational molecular biology and bioinformatics. Bioinformatics, 25(11):1422–1423, 2009.

[18] Anthony Mathelier, Xiaobei Zhao, Allen W Zhang, Fran ̧cois Parcy, Rebecca Worsley-Hunt, David J Arenillas, Sorana Buchman, Chih-yu Chen, Alice Chou, Hans Ienasescu, et al. Jaspar 2014: an extensively expanded and updated open-access database of transcription factor binding profiles. Nucleic acids research, page gkt997, 2013.

[19] Gavin E Crooks, Gary Hon, John-Marc Chandonia, and Steven E Brenner. Weblogo: a sequence logo generator. Genome research, 14(6):1188–1190, 2004.

[20] Stuart Geman and Donald Geman. Stochastic relaxation, gibbs distributions, and the bayesian restoration of images. IEEE Transactions on pattern analysis and machine intelligence, (6):721–741, 1984.

[21] Timothy L Bailey, Charles Elkan, et al. Fitting a mixture model by expectation maximization to discover motifs in bipolymers. 1994.

[22] Saul B Needleman and Christian D Wunsch. A general method applicable to the search for similarities in the amino acid sequence of two proteins. Journal of molecular biology, 48(3):443–453, 1970.

[23] Temple F Smith, Michael S Waterman, et al. Identification of common molecular subsequences. Journal of molecular biology, 147(1):195–197, 1981.

[24] Humberto Carrillo and David Lipman. The multiple sequence alignment problem in biology.SIAMJournal on Applied Mathematics, 48(5):1073–1082, 1988.

[25] David J Lipman, Stephen F Altschul, and John D Kececioglu. A tool for multiple sequence alignment. Proceedings of the National Academy of Sciences, 86(12):4412–4415, 1989.

[26] Florence Corpet. Multiple sequence alignment with hierarchical clustering. Nucleic acids research, 16(22):10881–10890, 1988.

[27] C ́edric Notredame, Desmond G Higgins, and Jaap Heringa. T-coffee: A novel method for fast and accurate multiple sequence alignment. Journal of molecular biology, 302(1):205–217, 2000.

[28] Julie D Thompson, Toby J Gibson, and Des G Higgins. Multiple sequence alignment using clustalw and clustalx. Current protocols in bioinformatics, (1):2–3, 2003.

[29] Sing-Hoi Sze, Yue Lu, and Qingwu Yang. A polynomial time solvable formulation of multiple sequence alignment. Journal of computational biology, 13(2):309–319, 2006.

[30] Michael Brudno, Rasmus Steinkamp, and Burkhard Morgenstern. The chaos/dialign www server for multiple alignment of genomic sequences. Nucleic acids research, 32(suppl2):W41–W44, 2004.

[31] Robert C Edgar. Muscle: multiple sequence alignment with high accuracy and high throughput. Nucleic acids research, 32(5):1792–1797, 2004.

[32] Michael Brudno, Chuong B Do, Gregory M Cooper, Michael F Kim, Eugene Davydov, Eric D Green, Arend Sidow, Serafim Batzoglou, NISC Comparative Sequencing Program, et al. Lagan and multilagan: efficient tools for large-scale multiple alignment of genomic dna. Genome research, 13(4):721–731, 2003.

[33] Mathieu Blanchette, W James Kent, Cathy Riemer, Laura Elnitski, Arian FA Smit, Krishna M Roskin, Robert Baertsch, Kate Rosenbloom, Hiram Clawson, Eric D Green, et al. Aligning multiple genomic sequences with the threaded blockset aligner. Genome research, 14(4):708–715, 2004.

[34] Francesca Chiaromonte, VB Yap, and W Miller. Scoring pairwise genomic sequence alignments. In Biocomputing 2002, pages 115–126. World Scientific, 2001.

[35] Arthur L Delcher, Simon Kasif, Robert D Fleischmann, Jeremy Peterson, Owen White, and Steven LSalzberg. Alignment of whole genomes. Nucleic acids research, 27(11):2369–2376, 1999.

[36] Arthur L Delcher, Adam Phillippy, Jane Carlton, and Steven L Salzberg. Fast algorithms for large-scale genome alignment and comparison. Nucleic acids research, 30(11):2478–2483, 2002.

[37] Ramu Chenna, Hideaki Sugawara, Tadashi Koike, Rodrigo Lopez, Toby J Gibson, Desmond G Higgins,and Julie D Thompson. Multiple sequence alignment with the clustal series of programs. Nucleic acids research, 31(13):3497–3500, 2003.

[38] Fabian Sievers, Andreas Wilm, David Dineen, Toby J Gibson, Kevin Karplus, Weizhong Li, RodrigoLopez, Hamish McWilliam, Michael Remmert, Johannes S ̈oding, et al. Fast, scalable generation of high-quality protein multiple sequence alignments using clustal omega. Molecular systems biology, 7(1):539, 2011.

[39] Thomas H Jukes, Charles R Cantor, et al. Evolution of protein molecules. Mammalian protein metabolism, 3(21):132, 1969.

[40] Motoo Kimura. A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences. Journal of molecular evolution, 16(2):111–120, 1980.

[41] Joseph Felsenstein. Evolutionary trees from dna sequences: a maximum likelihood approach. Journal of molecular evolution, 17(6):368–376, 1981.

[42] Masami Hasegawa, Hirohisa Kishino, and Taka-aki Yano. Dating of the human-ape splitting by a molecular clock of mitochondrial dna. Journal of molecular evolution, 22(2):160–174, 1985.

[43] Simon Tavaré. Some probabilistic and statistical problems in the analysis of dna sequences. Lectures on mathematics in the life sciences, 17(2):57–86, 1986.

[44] Robert R Sokal. A statistical method for evaluating systematic relationship. University of Kansas science bulletin, 28:1409–1438, 1958.

[45] Naruya Saitou and Masatoshi Nei. The neighbor-joining method: a new method for reconstructing phylogenetic trees. Molecular biology and evolution, 4(4):406–425, 1987.

[46] James S Farris. Methods for computing wagner trees. Systematic Biology, 19(1):83–92, 1970.

[47] Walter M Fitch. Toward defining the course of evolution: minimum change for a specific tree topology. Systematic Biology, 20(4):406–416, 1971.

[48] Ziheng Yang, Sudhir Kumar, and Masatoshi Nei. A new method of inference of ancestral nucleotide and amino acid sequences. Genetics, 141(4):1641–1650, 1995.

[49] Jeffrey M Koshi and Richard A Goldstein. Probabilistic reconstruction of ancestral protein sequences. Journal of molecular evolution, 42(2):313–320, 1996.

[50] Joseph Felsenstein and Joseph Felenstein. Inferring phylogenies, volume 2. Sinauer associates Sunderland, MA, 2004.

[51] Benjamin L Allen and Mike Steel. Subtree transfer operations and their induced metrics on evolutionary trees. Annals of combinatorics, 5(1):1–15, 2001.

[52] David L Swofford. Phylogeny reconstruction. Molecular systematics, pages 411–501, 1990.

[53] Magnus Bordewich and Charles Semple. On the computational complexity of the rooted subtree prune and regraft distance. Annals of combinatorics, 8(4):409–423, 2005.

[54] Manfred G Grabherr, Brian J Haas, Moran Yassour, Joshua Z Levin, Dawn A Thompson, Ido Amit, Xian Adiconis, Lin Fan, Raktima Raychowdhury, Qiandong Zeng, et al. Full-length transcriptome assembly from rna-seq data without a reference genome. Nature biotechnology, 29(7):644, 2011.

[55] Serena Liu and Cole Trapnell. Single-cell transcriptome sequencing: recent advances and remaining challenges. F1000Research, 5, 2016.

[56] Philip Brennecke, Simon Anders, Jong Kyoung Kim, Aleksandra A Kolodziejczyk, Xiuwei Zhang, Valentina Proserpio, Bianka Baying, Vladimir Benes, Sarah A Teichmann, John C Marioni, et al. Accounting for technical noise in single-cell rna-seq experiments. Nature methods, 10(11):1093, 2013.

[57] Julia Zeitlinger, Alexander Stark, Manolis Kellis, Joung-Woo Hong, Sergei Nechaev, Karen Adelman, Michael Levine, and Richard A Young. Rna polymerase stalling at developmental control genes in the drosophila melanogaster embryo. Nature genetics, 39(12):1512, 2007.

[58] Ginger W Muse, Daniel A Gilchrist, Sergei Nechaev, Ruchir Shah, Joel S Parker, Sherry F Grissom, Julia Zeitlinger, and Karen Adelman. Rna polymerase is poised for activation across the genome. Nature genetics, 39(12):1507, 2007.

[59] I. Jonkers and J. T. Lis. Getting up to speed with transcription elongation by RNA polymerase II. Nat. Rev. Mol. Cell Biol., 16(3):167–177, Mar 2015.

[60] Charles G Danko, Stephanie L Hyland, Leighton J Core, Andre L Martins, Colin T Waters, Hyung WonLee, Vivian G Cheung, W Lee Kraus, John T Lis, and Adam Siepel. Identification of active transcriptional regulatory elements from gro-seq data. Nature methods, 12(5):433, 2015.

[61] Amy C Seila, J Mauro Calabrese, Stuart S Levine, Gene W Yeo, Peter B Rahl, Ryan A Flynn, Richard A Young, and Phillip A Sharp. Divergent transcription from active promoters. science, 322(5909):1849–1851, 2008.

[62] Jennifer E Kurasz, Christine E Hartman, David J Samuels, Bijoy K Mohanty, Anquilla Deleveaux, Jan Mr ́azek, and Anna C Karls. Genotoxic, metabolic, and oxidative stresses regulate the rna repair operon of salmonella enterica serovar typhimurium. Journal of bacteriology, 200(23):e00476–18, 2018.

[63] Matthew K Iyer, Yashar S Niknafs, Rohit Malik, Udit Singhal, Anirban Sahu, Yasuyuki Hosono, Terrence R Barrette, John R Prensner, Joseph R Evans, Shuang Zhao, et al. The landscape of long noncoding rnas in the human transcriptome. Nature genetics, 47(3):199, 2015.

[64] Victor Ambros. The functions of animal micrornas. Nature, 431(7006):350, 2004.

[65] Alexei A Aravin, Gregory J Hannon, and Julius Brennecke. The piwi-pirna pathway provides an adaptive defense in the transposon arms race. science, 318(5851):761–764, 2007.

[66] Anita G Seto, Robert E Kingston, and Nelson C Lau. The coming of age for piwi proteins. Molecular cell, 26(5):603–609, 2007.

[67] Julius Brennecke, Alexei A Aravin, Alexander Stark, Monica Dus, Manolis Kellis, Ravi Sachidanandam,and Gregory J Hannon. Discrete small rna-generating loci as master regulators of transposon activity in drosophila. Cell, 128(6):1089–1103, 2007.

[68] Naoki Sugimoto, Shu-ichi Nakano, Misa Katoh, Akiko Matsumura, Hiroyuki Nakamuta, Tatsuo Ohmichi, Mari Yoneyama, and Muneo Sasaki. Thermodynamic parameters to predict stability of rna/dna hybrid duplexes. Biochemistry, 34(35):11211–11216, 1995.

[69] Julian L Huppert. Thermodynamic prediction of rna–dna duplex-forming regions in the human genome. Molecular BioSystems, 4(6):686–691, 2008.

[70] Fabian A Buske, Denis C Bauer, John S Mattick, and Timothy L Bailey. Triplexator: detecting nucleic acid triple helices in genomic and transcriptomic data. Genome research, 22(7):1372–1381, 2012.

[71] M. F. Lin, I. Jungreis, and M. Kellis. PhyloCSF: a comparative genomics method to distinguish protein coding and non-coding regions. Bioinformatics, 27(13):i275–282, Jul 2011.

[72] L. Wang, H. J. Park, S. Dasari, S. Wang, J. P. Kocher, and W. Li. CPAT: Coding-Potential Assessment Tool using an alignment-free logistic regression model. Nucleic Acids Res., 41(6):e74, Apr 2013.

[73] M. Guttman, P. Russell, N. T. Ingolia, J. S. Weissman, and E. S. Lander. Ribosome profiling provides evidence that large noncoding RNAs do not encode proteins. Cell, 154(1):240–251, Jul 2013.

[74] A. A. Bazzini, T. G. Johnstone, R. Christiano, S. D. Mackowiak, B. Obermayer, E. S. Fleming, C. E. Vejnar, M. T. Lee, N. Rajewsky, T. C. Walther, and A. J. Giraldez. Identification of small ORFs in vertebrates using ribosome footprinting and evolutionary conservation. EMBO J., 33(9):981–993, May 2014.

[75] Ruth Nussinov, George Pieczenik, Jerrold R Griggs, and Daniel J Kleitman. Algorithms for loop matchings. SIAM Journal on Applied mathematics, 35(1):68–82, 1978.

[76] David H Mathews, Jeffrey Sabina, Michael Zuker, and Douglas H Turner. Expanded sequence dependence of thermodynamic parameters improves prediction of rna secondary structure. Journal of molecular biology, 288(5):911–940, 1999.

[77] Rahul Tyagi and David H Mathews. Predicting helical coaxial stacking in rna multibranch loops. Rna, 13(7):939–951, 2007.

[78] Padideh Danaee, Mason Rouches, Michelle Wiley, Dezhong Deng, Liang Huang, and David Hendrix. bprna: large-scale automated annotation and analysis of rna secondary structure. Nucleic acids research, 46(11):5381–5394, 2018.

[79] Ronny Lorenz, Stephan H Bernhart, Christian H ̈oner Zu Siederdissen, Hakim Tafer, Christoph Flamm, Peter F Stadler, and Ivo L Hofacker. Viennarna package 2.0. Algorithms for molecular biology, 6(1):26,2011.

[80] MO Dayhoff, R Schwartz, and BC Orcutt. 22 a model of evolutionary change in proteins. In Atlas of protein sequence and structure, volume 5, pages 345–352. National Biomedical Research Foundation Silver Spring MD, 1978.

[81] Steven Henikoff and Jorja G Henikoff. Amino acid substitution matrices from protein blocks. Proceedings of the National Academy of Sciences, 89(22):10915–10919, 1992.

[82] Samuel Karlin and Stephen F Altschul. Methods for assessing the statistical significance of molecular sequence features by using general scoring schemes. Proceedings of the National Academy of Sciences, 87(6):2264–2268, 1990.

[83] Arlin Stoltzfus. On the possibility of constructive neutral evolution. Journal of Molecular Evolution, 49(2):169–181, 1999.

[84] Allan Force, Michael Lynch, F Bryan Pickett, Angel Amores, Yi-lin Yan, and John Postlethwait. Preservation of duplicate genes by complementary, degenerative mutations. Genetics, 151(4):1531–1545, 1999.

[85] Donald B Wetlaufer. Nucleation, rapid folding, and globular intrachain regions in proteins. Proceedingsof the National Academy of Sciences, 70(3):697–701, 1973.

[86] S. El-Gebali, J. Mistry, A. Bateman, S. R. Eddy, A. Luciani, S. C. Potter, M. Qureshi, L. J. Richardson, G. A. Salazar, A. Smart, E. L. L. Sonnhammer, L. Hirsh, L. Paladin, D. Piovesan, S. C. E. Tosatto, and R. D. Finn. The Pfam protein families database in 2019. Nucleic Acids Res., 47(D1):D427–D432,Jan 2019.

[87] J. A. Cuff and G. J. Barton. Application of multiple sequence alignment profiles to improve protein secondary structure prediction. Proteins, 40(3):502–511, Aug 2000.

[88] M. Ashburner, C. A. Ball, J. A. Blake, D. Botstein, H. Butler, J. M. Cherry, A. P. Davis, K. Dolinski, S. S. Dwight, J. T. Eppig, M. A. Harris, D. P. Hill, L. Issel-Tarver, A. Kasarskis, S. Lewis, J. C.Matese, J. E. Richardson, M. Ringwald, G. M. Rubin, and G. Sherlock. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat. Genet., 25(1):25–29, May 2000.

[89] Charles W Dunnett. A multiple comparison procedure for comparing several treatments with a control. Journal of the American Statistical Association, 50(272):1096–1121, 1955.

[90] Yoav Benjamini and Yosef Hochberg. Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal statistical society: series B (Methodological), 57(1):289–300, 1995.

[91] Steven T Runyon, Yingnan Zhang, Brent A Appleton, Stephen L Sazinsky, Ping Wu, Borlan Pan, Christian Wiesmann, Nicholas J Skelton, and Sachdev S Sidhu. Structural and functional analysis of the pdz domains of human htra1 and htra3. Protein Science, 16(11):2454–2471, 2007.

[92] David S Johnson, Ali Mortazavi, Richard M Myers, and Barbara Wold. Genome-wide mapping of in vivo protein-dna interactions. Science, 316(5830):1497–1502, 2007.

[93] Raja Jothi, Suresh Cuddapah, Artem Barski, Kairong Cui, and Keji Zhao. Genome-wide identification of in vivo protein–dna binding sites from chip-seq data. Nucleic acids research, 36(16):5221, 2008.

[94] Dominic Schmidt, Michael D Wilson, Christiana Spyrou, Gordon D Brown, James Hadfield, and Duncan T Odom. Chip-seq: using high-throughput sequencing to discover protein–dna interactions. Methods, 48(3):240–248, 2009.

[95] Ben Langmead. Aligning short sequencing reads with bowtie. Current protocols in bioinformatics, 32(1):11–7, 2010.

[96] Joel Rozowsky, Ghia Euskirchen, Raymond K Auerbach, Zhengdong D Zhang, Theodore Gibson, Robert Bjornson, Nicholas Carriero, Michael Snyder, and Mark B Gerstein. Peakseq enables systematic scoring of chip-seq experiments relative to controls. Nature biotechnology, 27(1):66, 2009.

[97] Yong Zhang, Tao Liu, Clifford A Meyer, J ́erˆome Eeckhoute, David S Johnson, Bradley E Bernstein, Chad Nusbaum, Richard M Myers, Myles Brown, Wei Li, et al. Model-based analysis of chip-seq (macs). Genome biology, 9(9):R137, 2008.

[98] Shiliyang Xu, Sean Grullon, Kai Ge, and Weiqun Peng. Spatial clustering for identification of chip-enriched regions (sicer) to map regions of histone methylation patterns in embryonic stem cells. In Stem Cell Transcriptional Networks, pages 97–111. Springer, 2014.

[99] Jason Ernst and Manolis Kellis. Chromhmm: automating chromatin-state discovery and characterization. Nature methods, 9(3):215, 2012.

[100] Peter M Clark, Phillipe Loher, Kevin Quann, Jonathan Brody, Eric R Londin, and Isidore Rigoutsos. Argonaute clip-seq reveals mirna targetome diversity across tissue types. Scientific reports, 4:5947,2014.

[101] Xu Zhou and Erin K O’Shea. Integrated approaches reveal determinants of genome-wide binding and function of the transcription factor pho4. Molecular cell, 42(6):826-836, 2011.

4

Creative Commons License

This work is licensed by David Hendrix under a
Creative Commons Attribution 4.0 International License (CC BY)

You are free to:
Share — copy and redistribute the material in any medium or format
Adapt — remix, transform, and build upon the material

The licensor cannot revoke these freedoms as long as you follow the license terms.

Under the following terms:
Attribution — You must give appropriate credit, provide a link to the license, and indicate if
changes were made. You may do so in any reasonable manner, but not in any way that suggests
the licensor endorses you or your use.
No additional restrictions — You may not apply legal terms or technological measures that legally restrict others from doing anything the license permits.

5

Recommended Citations

APA outline:

Source from website:

- (Full last name, first initial of first name).  (Date of publication). Title of source. Retrieved from https://www.someaddress.com/full/url/

Source from print:

- (Full last name, first initial of first name).  (Date of publication). Title of source. Title of container (larger whole that the source is in, i.e. a chapter in a book), volume number, page numbers.

Examples

If retrieving from a webpage:

- Berndt, T. J. (2002). Friendship quality and social development. Retrieved from insert link.

If retrieving from a book:

- Berndt, T. J. (2002). Friendship quality and social development. Current Directions in Psychological Science, 11, 7-10.

MLA outline:

Author (last, first name). Title of source. Title of container (larger whole that the source is in, i.e. a chapter in a book), Other contributors, Version, Number, Publisher, Publication Date, Location (page numbers).

Examples

- Bagchi, Alaknanda. “Conflicting Nationalisms: The Voice of the Subaltern in Mahasweta Devi’s Bashai Tudu.” Tulsa Studies in Women’s Literature, vol. 15, no. 1, 1996, pp. 41-50.

- Said, Edward W. Culture and Imperialism. Knopf, 1994.

Chicago outline:

Source from website:

- Lastname, Firstname. “Title of Web Page.” Name of Website. Publishing organization, publication or revision date if available. Access date if no other date is available. URL .

Source from print:

- Last name, First name. Title of Book. Place of publication: Publisher, Year of publication.

Examples

- Davidson, Donald, Essays on Actions and Events. Oxford: Clarendon, 2001. https://bibliotecamathom.files.wordpress.com/2012/10/essays-on-actions-and-events.pdf.

- Kerouac, Jack. The Dharma Bums. New York: Viking Press, 1958.

6

Versioning

This page provides a record of changes made to this guide. Each set of edits is acknowledged with a 0.01 increase in the version number. The exported files for this toolkit reflect the most recent version.

If you find an error in this text, please fill out the form at bit.ly/33cz3Q1

Version

Date

Change Made

Location in text

0.1

MM/DD/YYYY