Principles and Applications of Modern
DNA Sequencing
EEEB GU4055
Session 4: Scientific Python
Today's topics
1. review of topics thus far.
2. file I/O
3. fastq to fasta
4. genome annotation
5. introducing numpy and pandas
Python Dictionaries
A look-up table. Store values associated with look-up keys. Very efficient data structure for storing and retrieving data.
# make a dict using curly bracket format
adict = {
'key1': 'val1',
'key2': 'val2',
}
# request a value using a key
print(adict['key1'])
val1
Comments are a key part of writing good code
# import the random library
import random
# create a list with 1000 random numbers between 0-10
integer_list = [random.randint(0, 10) for i in range(1000)]
# create an empty dictionary
counter = {}
# iterate over elements of the integer list
for item in integer_list:
# conditional True if item is not already in the dict keys
if item not in counter:
# set the value to 1 for this key
counter[item] = 1
# item is already in dict keys
else:
# increment value by 1 for this key
counter[item] += 1
Python Advanced
You have now learned all of the core object types in Python. From these simple objects more complex Python object can be built. Thousands of complex software tools have been developed from creatively combining these objects with Python coding routines.
# The core Python object types
a_integer = 8
b_float = 0.2345
c_string = "a string"
d_list = ["a", "list", "of", "strings"]
e_tuple = ("a", "tuple", "of", "strings")
f_dict = {"a key": ["and value in a dictionary"]}
File I/O: bash to Python
%%bash
wget http://eaton-lab.org/data/40578.fastq.gz -q -O 40578.fastq.gz
import os
import gzip
import requests
# the URL to file 1
url1 = "https://eaton-lab.org/data/40578.fastq.gz"
# open a file for writing and write the content to it
with gzip.open("40578.fastq.gz", 'wb') as ffile:
ffile.write(requests.get(url1).content)
File I/O: file objects
Reading and writing files in Python is done through File objects. You first create an object, then use it is some way (functions), and finally close it.
# open a file object in write-mode
ofile = open("./datafiles/helloworld.txt", 'w')
# write a string to the file
ofile.write("hello world")
# close the file object
ofile.close()
File I/O: file objects
The term 'with' creates context-dependence within the indented block. The object can have functions that are automatically called when the block starts or ends. For an open file object the block ending calls .close(). This often is simpler better code.
# a simpler alternative: use 'with', 'as', and indentation
with open("./helloworld.txt", 'w') as ofile:
ofile.write("hello world")
File I/O: file objects
To reiterate, these two code blocks are equivalent.
# open a file object in write-mode
ofile = open("./helloworld.txt", 'w')
# write a string to the file
ofile.write("hello world")
# close the file object
ofile.close()
# a simpler alternative: use 'with', 'as', and indentation
with open("./helloworld.txt", 'w') as ofile:
ofile.write("hello world")
File I/O: file objects
Compression or decompression is as simple as writing or reading using a File object from a compression library (e.g., gzip or bz2).
import gzip
# open a gzipped file object in write-mode to expect byte strings
ofile = gzip.open("./helloworld.txt", 'wb')
# write a byte-string to the file
ofile.write(b"hello world")
# close the file object
ofile.close()
File I/O: reading
Open a file and call .read() to load all of the contents at once to a string or bytes object.
# open a file
fobj = open("./data.txt", 'r')
# read data from this file to create a string object
fdata = fobj.read()
# close the file
fobj.close()
# print part of the string 'fdata'
print(fdata[:500])
File I/O: reading
Open a file and call .read() to load all of the contents at once to a string or bytes object.
# open a gzip file
fobj = gzip.open("./data.fastq", 'r')
# read compressed data from this file and decode it
fdata = fobj.read().decode()
# close the file
fobj.close()
# print part of the string 'fdata'
print(fdata[:500])
File I/O: reading
Once again, the 'with' context style of code is a bit more concise.
# open a gzip file and read data using 'with'
with gzip.open("./data.fastq", 'r') as fobj:
fdata = fobj.read().decode()
print(fdata[:500])
File I/O: file paths
File paths are an important concept to master in bioinformatics. Be aware of the absolute path to where you are, and where the files you want to operate on are located, and understand how relative paths can also point to these locations.
# os.path.abspath() prints the absolute path from a relative path
os.path.abspath(".")
# os.path.join() combines two or more parts of paths together with /
os.path.join("/home/deren", "file1.txt")
# os.path.basename() and os.path.dirname() return the dir and filename
os.path.basename("/home/deren/file.txt")
# check if a file exists before trying to do something with it
os.path.exists("badfilename.txt")
The FASTA file format
The fasta format is commonly used to store sequence data. We learned about it in our first notebook assignment and also saw some empirical examples representing full genomes. The delimiter ">" separates sequences. Files typically end in .fasta, .fna, (DNA specific) or .faa (amino acids specific).
>mus musculus gene A
AGTCAGTCAGCGCTAGTCATAACACGCAAGTCAATATATACGACAGCAGCTAGCTACTTCGACA
CAGTCGATCAGCTAGCTGACTACTATATATTTTTATATGTAAAAAAAACATATGCGCGCTTTTG
GGGGAGTATTTTATGCATATCATGCAGCATATAGGTAGCTGTGCATGCTGCTAGCACGATCGTA
CATGCTAGCTAGCTAGCTAGCTAGCTAGCTGACTAGCTAGTGCTAGCTAGCTATATATATATAT
>mus musculus gene B
ACGTACGTACGTACGTAGCTAGCTACATGCTAGCATGCATGCTAGCTAGCTATATATAGCCCCC
CAGCGGGGGGCGTCATGCATAAAAAAAAAAAGCATCATGCCGCGCCCCTAGCGCGTATTTTCTT
...
The FASTA file format
Challenge: Write code to combine a fasta header (e.g., "> sequence name") and a random sequence of DNA to a create valid fasta data string. Then write the data to a file and save it as "datafiles/sequence.fasta".
# a function to return a random DNA string
def random_dna(length):
dna = "".join([random.choice("ACGT") for i in range(length)])
return dna
# format dna string to fasta format
dna = random_dna(20)
fasta = "> sequence name\n" + dna
# write it to a file
with open("./datafiles/sequence.fasta", 'w') as out:
out.write(fasta)
The FASTQ file format
The fastq format is commonly used to store sequenced read data. It differs from fasta in that it contains quality (confidence) scores. Each sequenced read represented by four lines, and a single file can contain many millions of reads.
@SEQ_ID
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65
FASTQ quality scores
Quality scores are encoded using ASCII characters, where each
character can be translated into an integer score that is
log10 probability the base call is incorrect:
$$ Q = -10 * log_{10}(P) $$
# load a phred Q score as a string:
phred = "IIIIIGHIIIIIHIIIIFIIIDIHGIIIBGIIFIDIDI"
# ord() translates ascii to numeric
q_scores = [ord(i) for i in phred]
# values are offset by 33 on modern Illumina machines
q_scores = [ord(i) - 33 for i in phred]
print(q_scores)
[40, 40, 40, 40, 40, 38, 39, 40, 40, 40, 40, 40, 39, 40, 40, ...
FASTQ quality scores
Quality score is an integer log10 probability the base call
is incorrect:
$$ Q = -10 * log_{10}(P) $$
# Q=30 means 3 decimal places in the probability of being wrong (0.001)
import math
print(-10 * math.log10(0.001))
# print the probability associated with the first few q_scores
probs = [10 ** (q / -10) for q in q_scores]
print(probs)
30.0
[0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.00015848931924611142, ...
FASTQ conversion
Now that you understand reading and writing files, working with string and list objects, and the format of fastq and fasta file formats, you are prepared to write a function to convert from one to the other.
def fastq2fasta(in_fastq, out_fasta):
"""
(1) Write a function;
(2) read 'datafiles/40578.fastq.gz' from disk;
(3) convert to fasta format; and
(4) write result to a file
"""
# 2. read in the fastq file
with gzip.open(in_fastq, 'rb') as indata:
fastq = indata.read().decode()
# 3. convert to fasta: start with an empty list
fastadata = []
# split file into separate reads on \n@ delimiter
reads = fastq.split("\n@")
for read in reads:
# split read into 4 lines
lines = read.split("\n")
# join ">", line[0], and line[1] to make fasta formatted sequence
fastadata.append(">" + lines[0] + "\n" + lines[1])
# join strings in the list back together into a string
fasta = "\n".join(fastadata)
# write to a file
with open(out_fasta, 'w') as out:
out.write(fasta)
Dictionaries in action: translation
GENCODE = {
'ATA': 'I', 'ATC': 'I', 'ATT': 'I', 'ATG': 'M',
'ACA': 'T', 'ACC': 'T', 'ACG': 'T', 'ACT': 'T',
'AAC': 'N', 'AAT': 'N', 'AAA': 'K', 'AAG': 'K',
'AGC': 'S', 'AGT': 'S', 'AGA': 'R', 'AGG': 'R',
'CTA': 'L', 'CTC': 'L', 'CTG': 'L', 'CTT': 'L',
'CCA': 'P', 'CCC': 'P', 'CCG': 'P', 'CCT': 'P',
'CAC': 'H', 'CAT': 'H', 'CAA': 'Q', 'CAG': 'Q',
'CGA': 'R', 'CGC': 'R', 'CGG': 'R', 'CGT': 'R',
'GTA': 'V', 'GTC': 'V', 'GTG': 'V', 'GTT': 'V',
'GCA': 'A', 'GCC': 'A', 'GCG': 'A', 'GCT': 'A',
'GAC': 'D', 'GAT': 'D', 'GAA': 'E', 'GAG': 'E',
'GGA': 'G', 'GGC': 'G', 'GGG': 'G', 'GGT': 'G',
'TCA': 'S', 'TCC': 'S', 'TCG': 'S', 'TCT': 'S',
'TTC': 'F', 'TTT': 'F', 'TTA': 'L', 'TTG': 'L',
'TAC': 'Y', 'TAT': 'Y', 'TAA': '_', 'TAG': '_',
'TGC': 'C', 'TGT': 'C', 'TGA': '_', 'TGG': 'W',
}
GENCODE["CTA"]
L
How are genes identified?
Genome annotation is the process of labeling genomic elements,
including genes and their parts. Your reading by Yandell and Ence
introduces the concepts of genome annotation, a process that
has evolved rapidly over the last decade.
We used a translation dictionary in Python to search a string
of DNA for start and stop codons to an open reading frame (ORF)
: a region that could be translated. This was a relatively crude approach.
Modern approaches use RNA-seq, in which RNA is extracted and reverse-transcribed into cDNA
-- the subset of the genome matching to coding genes -- that
is then sequenced and mapped to a reference genome. Thus we only
examine a subset of the genome for annotation.
Genome annotations in practice
Annotated genomes of model organisms, like humans and
Drosophila, are works in progress. However, they are considered
to be highly accurate in comparison to most other genomes that
have been sequenced recently. This is because their genomes are
assembled better (more contiguous), and because they are able to
build upon decades of experimental work to infer the function
of genes.
Your next assignment will cover assembly statistics, like N50,
what it means, how to calculate it, and what type of values
represent good versus poor genome assemblies.
Genome annotations in practice
1. Repeat mask the genome to prevent data mapping to
repetitive regions (e.g., transposons), which also have
open reading frames (ORFs) and can interfere with the
identification of other ORFs associated with genes.
2. Map RNA-seq reads to a reference (e.g., using TopHat) and assemble
transcript contigs from overlapping mapped reads (e.g., using
Cufflinks).
3. ab initio gene prediction: little or no external evidence. Find the most likely coding sequence (CDS) in ORFs. Model parameters (e.g., GC content, intron lengths) affect accuracy, and vary among organisms. Additional evidence includes RNA-seq and Protein sequence data to identify introns/exons and UTRs.
4. Annotation: combining evidence from multiple data types or analyses.
Two more Python object: Arrays and DataFrames
numpy and pandas are the reason Python is popular for data science.
# numpy arrays are similar to lists, but also very different.
import numpy as np
arr = np.array([0, 5, 2, 1, 14])
# pandas dataframes are tables {column name: data}
import pandas as pd
df = pd.DataFrame({"column1": arr})
[ 0 5 2 1 14]
column1
0 0
1 5
2 2
3 1
4 14
A short numpy intro
numpy arrays are super efficient data structures. All data in an array is of the same type (int8, int64, float64) which makes computations fast. In addition, is performs broadcasting on arrays and includes many numerical functions.
# create an array from a list
arr = np.array([0, 5, 2, 1, 14])
# create a 1-d array of all zeros
arr = np.zeros(100)
# create a 2-d array of all zeros
arr = np.zeros((10, 10))
# create a 2-d array of all zeros of integer type data
arr = np.zeros((10, 10), dtype=np.int8)
A short numpy intro
numpy arrays are super fast and efficient. All data is of the same 'type' (int8, int64, float64). Functions can be 'broadcast', and many numerical functions are supported.
# array of 1-10 reshaped to be 5 rows 2 columns
narr = np.arange(10).reshape((5, 2))
# print the shape of the array
print(narr.shape)
# arrays can be indexed just like lists
print(narr)
(5, 2)
array([[0, 1],
[2, 3],
[4, 5],
[6, 7],
[8, 9]])
A short numpy intro
numpy arrays are super fast and efficient. All data is of the same 'type' (int8, int64, float64). Functions can be 'broadcast', and many numerical functions are supported.
# two arrays can be added together value-by-value
xarr = np.arange(10) + np.arange(10)
# to lists are concatenated when added, not processed per-value
larr = list(range(10)) + list(range(10))
print(xarr)
print(larr)
[ 0 2 4 6 8 10 12 14 16 18]
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
A short numpy intro
numpy arrays are super fast and efficient. All data is of the same 'type' (int8, int64, float64). Functions can be 'broadcast', and many numerical functions are supported.
# arrays can be indexed and sliced just like lists but in more dimensions
xarr = np.zeros(1000).reshape((10, 10, 10))
print(xarr[0])
[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]
A short numpy intro
Many useful functions/modules in numpy, like .random.
# sample random values with numpy
np.random.normal(0, 1, 10)
array([-0.35549967, -2.1416518 , -0.49230544, 1.47456753, 1.31386496,
-1.38097489, -0.2578635 , -1.60208958, -0.45677291, 0.91109757])
A short pandas intro
DataFrames are tables that can be selected by row or column names or indices.
data = pd.DataFrame({
"randval": np.random.normal(0, 1, 10),
"randbase": np.random.choice(list("ACGT"), 10),
"randint": np.random.randint(0, 5, 10),
})
randval randbase randint
0 0.602565 A 1
1 -0.657427 G 1
2 -0.907259 C 0
3 0.775811 T 0
4 0.601185 T 4
5 2.155603 G 2
...
A short pandas intro
Read/write data to and from tabular formats (e.g., CSV).
# load comma or tab-separated data from a file
df = pd.read_csv("datafile.csv", sep="\t")
# load a datafile from a URL
df = pd.read_csv("https://eaton-lab.org/data/iris-data-dirty.csv", header=None)
0 1 2 3 4
0 5.1 3.5 1.4 0.2 Iris-setosa
1 4.9 3.0 1.4 0.2 Iris-setosa
2 4.7 3.2 1.3 0.2 Iris-setosa
3 4.6 3.1 1.5 0.2 Iris-setosa
4 5.0 3.6 1.4 0.2 Iris-setosa
.. ... ... ... ... ...
A short pandas intro
Indexing rows, columns, or cells. Similar to numpy or Python lists you can use indexing with .iloc
# select row 0, all columns
data.iloc[0, :]
# select column 0, all rows
data.iloc[:, 0]
# select a range of columns and rows
data.iloc[:4, :4]
# select a specific cell
data.iloc[3, 2]
A short pandas intro
To index by row or column names use .loc.
# select row 0, all columns
data.loc[0, :]
# select column "randint", all rows
data.loc[:, "randint"]
# create a boolean mask (True, False) where
mask = data.loc[:, "randint"] > 1
# apply mask to select all rows where mask is True
data.iloc[mask, :]
Assignment
Your assignment will introduce numpy and pandas for operating
on data related to genome annotations. You will calculate genome
N50 using numpy, and you will read and operate on a GFF file using
pandas. Revisit the Yandell paper as needed.
Two assigned short chapters: Chapters 2,3 of the Python Data Science Handbook.
One assigned primary reading: “OrthoDB: A Hierarchical Catalog of Animal, Fungal and Bacterial Orthologs.” Nucleic Acids Research 41 (Database issue): D358–65. https://doi.org/10.1093/nar/gks1116