Comp 7 / Bio 40 Project: Reliability of Metagenomics Reads

Due December 12, 2016


In this class, we are studying 16S rRNA metagenomic sequencing, in which the sequences read come from one of the hypervariable regions of the 16S rRNA genes, which are generally well conserved across prokaryotes. Because of the variability in the region we are sequencing, it is possible to use computational means to try to identify the taxonomic classification of each sequence read. Of course, the sequence data can be noisy; noise may be represented by nucleotides reported as 'N' in the sequence data files. In addition, not every bacterial species has previously been sequenced, and some species may have sufficiently similar sequences even in these hypervariable regions, making it difficult to classify them exactly.

The MiSeq metagenomics pipeline uses the method of Wang et al. ( Assignment of rRNA Sequences into the New Bacterial Taxonomy. Q. Wang, G. M. Garrity, J. M. Tiedje, J. R. Cole. Appl. Environ. Microbiol. 73(16):5261, 2007) to taxonomically classify sequences. This is a probablistic method, so the classification at each level of the taxonomic hierarchy is associated with a confidence score, which corresponds roughly to an estimate of the probability that the classification is correct.

The software uses a cutoff of 80% confidence to report a result. If the classification software is less than 80% confident in its classification at any given taxonomic level, it reports that the sequence is "Unclassified" at that level. Some samples have many unclassified reads, while others have relatively few.


We are going to focus on the most specific level of taxonomic classification reported by our software, the genus level. Our hypothesis is the reads that could not be classified with greater than 80% certainty at the genus level had a higher percentage of undetermined nucleotides (Ns). In this assignment, you will analyze next-generation sequencing data to confirm or refute this hypothesis.


You will write a python program to address this question by combining data from two data files produced by the MiSeq software. We will help you design the program and will even give you the outline of the code to start with, but you will fill in each of the pieces!

In order to test your hypothesis, you will need to extract information from a data file to determine which reads are "unclassified" at the genus level, and which are not. You will also need to get the actual nucleotide sequences reported for each of those reads. These are stored in separate data files, so you will have to match the data up between two files using the read identifiers.

Matching data between two files is a very common problem in bioinformatics, and it is one that is easily solved using the dictionary data structure that you have recently learned. First we will describe the two file formats you will encounter, and then we will describe the program you will need to complete.

FASTQ File Format

Most next-generation raw sequencing data comes in the form of FASTQ files. FASTQ files list each read that the sequencer captured. Each read in a FASTQ file consists of four lines. The first line starts with the "@" character and then a unique identifier naming the read; the second line contains the nucleotide sequence of the read (with N representing reads that are too noisy for a confident nucleotide call); the third line starts with the "+" character and optionally includes the read name again, and the fourth line contains a set of symbols the same length as the read, representing the quality score for each nucleotide. A sample read from a FASTQ file is shown below.

@M01675:6:000000000-A6HK1:1:1101:14905:1750 1:N:0:39

For this analysis, you will not need to use the third or fourth lines of each set. Keep this in mind when you are writing a function to read in the FASTQ file.

Classification File

The classification file is a text file that contains the read ID and the phylogenetic classification for each read. Phylogeny classification can happen at each level of scientific classification. For example, consider the Horse and E. coli bacterium below.

These organisms can be classified at any of the following levels:

E. ferus
E. coli

At the Domain level, we could say that the organism on the left is a eukaryote, and the organism on the right is a prokaryote. However, at the Species level, we can say that the organism on the left is a horse and the organism on the right is an E. coli bacterium.

When running a 16S experiment, the software tries to identify each read to the most fine-grained level of phylogeny that it can from the sequence information. Sometimes, the software will be able to classify a read to the genus level (what you are interested in), and sometimes the software won’t have enough information to perform the classification to that level.

Each read in the classification file consists of two lines. The first line includes the read name (similar to that from the fastq file, with the differences described below), and the second line lists the classification at each taxonomic level and the corresponding probability. The classifications are semicolon-delimited, and the columns correspond to the classification at the level of kingdom, phylum, class, order, family, and genus.

Note that there is no species-level classification reported; we don't have enough sequence variability nor depth of sequencing to make accurate claims at that level.

A sample two-row entry from the classfication file looks like the following:

>39 M01675:6:000000000-A6HK1:1:1101:14905:1750
Bacteria; 0.96; Firmicutes; 0.96; Bacilli; 0.96; Lactobacillales; 0.71; Enterococcaceae; 0.61; Melissococcus; 0.61; 

This suggests 96% confidence that read M01675:6:000000000-A6HK1:1:1101:14905:1750 comes from the domain of Bacteria, but only 61% confidence of the mapping to the genus Melissococcus. In the description of the code outline below, we refer to the second line in each two-row classification file entry (the one starting with "Bacteria" in our example above) as a "classification string").

Note also that some classification strings are shorter than this and don't include data for some of the finer levels of classification. For example, you might see the following lines, which only classify that sequence down to the class level:

>39 M01675:6:000000000-A6HK1:1:1101:14975:1719
Bacteria; 0.96; Proteobacteria; 0.91; Gammaproteobacteria; 0.91;

Keep in mind that the read ID lines used by the classification file and the corresponding ID lines in the FASTQ file are slightly different and have to be modified to be compatible with one another. Specifically, a read-name line in the classification file looks like:

>39 M01675:6:000000000-A6HK1:1:1101:14064:1631
while the comparable line in the fastq file is:
@M01675:6:000000000-A6HK1:1:1101:14064:1631 1:N:0:39
We suggest you use just M01675:6:000000000-A6HK1:1:1101:14064:1631 as the read name. So in the first case, you will need to remove the '>39 ' part of the line, while in the second case you will need to remove the ' 1:N:0:39' part and the leading '@' character. Suggestions for doing this are more thoroughly documented in the skeleton code and below.

Code Introduction

We have provided skeleton code that you can use to construct your progam. We have broken the code up into functions that you need to write for this assignment. We have also written a main function that can be viewed as the high-level summary of the program.

Functions that need to be implemented in the skeleton code simply say, “pass” under the function definition. For example:

def readInClasses(input):

We have provided this skeleton code to make this project more manageable. However, if you would like to create your own program from scratch, feel free to disregard the provided code.

Skeleton Code Overview

At the top of the skeleton code, there are two constants: CLASS_FILE and READS_FILE. Add the names of the classification file and the FASTQ file in quotation marks to each of these constants, respectively. Make sure that these files are located in the same directory as the code.

The classification file is called metagen.txt and the fastq file, metagen.fastq. Both are available for download from the projects page. There should also be smaller sample files for use when you are testing your code. The code is initially set to use just the smallest of these (rand_10.txt and rand_10.fasta). Only replace these file names with the names of the larger test files and, eventually, the full sized files, after you get your code to work on the smaller data sets.

Now take a look at the main function in the skeleton code. It performs the following steps:

  1. It uses the readInClasses function to create a dictionary classes that maps read IDs to their corresponding classification strings (a string containing the line from the classification file that contains the taxonomic classifications and their confidence scores for that read).
  2. Using classes, it creates a dictionary classified, that maps the read IDs to a boolean value, either True or False, depending on whether the confidence in the classification of that read ID at the genus level is above or equal to the constant CUTOFF (set to 0.8), in which case the value for the read ID in classified is True, and False otherwise.
  3. Create one more dictionary, reads, mapping read IDs to a string containing the nucleotide sequence for that read from the FASTQ file.
  4. This step is the heart of your calculation. Using the two dictionaries, classified and reads, the findAvgNCount() function should go through each read id that is a key of classified, and calculate the percent N's for the corresponding read. Define two lists, one for the classified ids and the other for the unclassified ids. Add the float representing the percent N's to one of the two lists depending on whether the value for the read id in classified is True or False (i.e., whether the read for the id is classified at the genus level or not). Now use the helper function avgList() to average the values in each list, and return a list that contains just the two averages.
  5. Finally, the report function nicely prints the values returned by the previous step.

Implementation Strategy

We suggest that you start by getting some of the smaller and simpler functions to work first. Look at the function descriptions below. Helper functions like avgList and getPctNs should be easy to write and debug on their own. You must write some code to test these, although you won't need to run all the tests in your final version. Then try something like readInFastq. Test each of these out before you go on to another function.

Important note: For a project of this size, it is essential that you only write and test one function at a time! Once one function is working and thoroughly tested, then go on to the next one. If you come to us having written the whole program but never tested any of it and claiming it doesn't work, we will not be able to help you without removing everything you wrote!

Function Details

We have listed and described each of the functions you need to implement below:

def readInClasses(inputfile):
This function reads in the classification file (inputfile) into a dictionary so that for each read, the read ID is the key and the classification string is the value. (That is, this dict has the form {id: classification, ....}.) Before adding the key-value pair to the dictionary, use the helper function prepareKey() to remove unnecessary parts of the line to create the read ID. This function returns the dictionary created here.

def partitionClasses(classes):
This function uses the dictionary returned from readInClasses (classes) and creates a dictionary, classified, where each read id is the key and the value is either True or False. True indicates a read whose confidence score for genus classification was greater than or equal to the CUTOFF constant at the top of the file, and False indicates a read whose confidence score for genus classification was strictly less than the CUTOFF. The confidence score for genus classification can be obtained by calling the getScore() function. This function should return the created dictionary, classified.

def getScore(classString):
Takes in a classification string (classString) and finds the percent certainty of the classification at the genus level. If the read was classified to the genus level, the list generated from splitting the classString on semicolons will have 12 elements, and the element with index 11 will be the percent certainty at the genus level. If the list has fewer than 12 elements, then the read was not classified to the genus level, and the function should return 0.0. Make sure that this function always returns a float.

def prepareKey(idName):
Gets rid of the >, number, and white space at the beginning of each ID (idName). This function returns the remaining string as a read ID.

def readInFastq(input):
Reads in the FASTQ file (input) into a dictionary so that for each read, the read ID is the key and the nucleotide string is the value. Each identifier contains a second block of text that should be disregarded. Consider using the split() function to discard this second block of text. The 3rd and 4th lines of each fastq entry should be discarded. Consider using the modulo operator ("%") to determine which lines you are going to use and which lines you will skip. This function should return the dictionary mapping read IDs to nucleotide sequence strings.

def findAvgNCount(classified, reads):
Determines the average percent N content of the reads that are classified and the reads that are not classified. The classified dictionary contains the read ids as keys and True or False as values. The reads dictionary contains the read ids as keys and the nucleotide strings as values. This function should construct two lists of the percent N content, one for the classified reads and another for the unclassified reads. A helper function getPctNs() can be used to calculate the percent N content of a nucleotide string passed in as an argument. This function should also use the helper function avgList() to compute the average of each of the two lists. Store the two averages in a list called averagedPercent, where averagedPercent[0] is the average percent N content of the classified reads and averagedPercent[1] is the average percent N content of the unclassified reads. Return this list containing the two averaged values.

def getPctNs(read):
Determines the percent N content of nucleotide string passed in (read). It is good defensive programming to make sure that all of the characters in read are capitalized. This can be accomplished in Python by using the upper function. This function should also convert the decimal value to the percentage. For example, if 0.002 of the nucleotides are N, then 0.2% of the nucleotides are N. This function returns a float representing the percentage of the nucleotides in the input parameter read that are N.

def avglist(listToAverage):
Calculates the average of listToAverage. Assumes that all elements in the list are either ints or floats. This function returns a float.

def report(avgNCount):
avgNCount is the list returned by the function findAvgNCount. This function prints a description of the average percent N content in the classified reads and the unclassified reads. Before printing these values, round these percentages to 4 decimal places using the round function (look it up in Learning Python).