Friday, July 11, 2014

Cool Unix Commands

I will add to this list as I discover new ones.  If you have a favorite or useful command feel free to include it in a comment on this post.


Convert a FASTQ file to FASTA (originally posted here):
sed -n '1~4s/^@/>/p;2~4p' 

NOTE: this assumes that each FASTQ entry spans only four lines as is customary.



Convert a SAM file to FASTA

awk '{OFS=""}{print $1, "\n", $10; }' file.sam > file.fasta

NOTE: You will loose a lot of information in the sam file.  You can save more of that info by adding column variables to the print statement.  Also, you may have to change the column variable numbers depending on your sam file format.  This is just a general example.



Replace spaces in file names with underscore (originally posted here)

rename ' ' '_' *

NOTE:  do NOT put spaces in file names!!  This is so annoying!



Get a histogram of sequence lengths from FASTA/Q files (from Surge Biswas)

FASTQ:  cat <fastq file> | awk '{if(NR%4==2) print length($1)}' | sort -n | uniq -c
FASTA:  cat <fasta file> | awk '{if(NR%4==0) print length($1)}' | sort -n | uniq -c



Do arithmetic operations on the bash command line

echo $((1 + 1))
echo $((1 - 1))
echo $((1 * 1))
echo $((1 / 1))
echo $(((1+3) / (1+1)))

For floating point operations you can use the bc tool.  For example

echo "scale=1; 1/2" | bc



Add a comment to a bash command on the command line

<command>; # this is a comment line

A practical example:  mv file1 old_file1; # there is now a new file1 is a more recent version

NOTE:  Do you ever have a long and complex command for which you would like to save a simple note?  You can use this little trick and the note will be saved along side your command in your history.  The next time you look through your history to rerun the command you will also see the associated note.



Count the number of bases in a FASTA file

grep -v ">" file.fasta | wc | awk '{print $3 - $1}'

(from martinghunt on SEQanswers)


Thursday, June 26, 2014

How to Learn Bioinformatics

Introduction

At least once a month someone asks me for help learning bioinformatics.  I love it when this happens because it usually means they want to take control of their own analysis thereby freeing up my time for problems that interest me.  This post is a collection of tips and resources for people wanting to learn how to do bioinformatics.

Keep These Things in Mind:
  • Learning the basics of bioinformatics is easy.  The basics as described in this post are often taught in high school.  However, don't get frustrated if you don't understand everything all at once.  Learning anything new takes time and practice no matter its difficulty.  
  •  A little bit goes a long way.  I estimate that nearly 90% of my work is occupied by simple routine procedures.  Learning how to do these tasks will substantially expand your ability to analyze and interpret your data.
  • Google it.  Google is the best resource for learning new techniques and trouble shooting problems.  If you have a question type exactly what you would say to a person into the google search bar.  When you take questions to your bioinformatics friends it's likely they won't know the answer offhand and will google it anyway.
  • Try it.  If you're not sure about something try it and see what happens.  Generally, there is very little danger is just trying a command to see if and how it works.  That being said it's a good idea to backup important files and data just incase something goes very wrong.  Every Unix programmer that I know has deleted a really important file using the rm command (which is one of the few irreversible Unix commands).  It's going to happen to you too so make a backup.

Learn the Unix Basics
  • Get on a Unix machine.  Doing is the most important aspect of learning Unix.  You will never fully understand the basic concepts if you only read about them.  Mac users have it easy because OSX is build on a unix shell.  Simply open the terminal application and you are ready to start with an online tutorial.  For non-mac users I recommend finding an old computer and installing a Linux/Unix operating system like Ubuntu.  A slightly more difficult approach would be to partition the drive of an existing computer to dual boot a Linux/Unix OS along with the existing OS.
  • Complete an online tutorial
  • Buy a book if you are a book learner.  However, the basic can pretty much all be learned using online materials.  My favorite Unix book is O'Reilly's Unix Power Tools.

Learn a Scripting Language
  • Pick a scripting language.  Scripting languages are computer languages that are not compiled (i.e. they are interpreted by the computer on the fly).  The two most popular bioinformatics scripting languages are Perl and Python.  Both languages have their strengths and weakness, but I personally prefer Perl.
  • Complete an online tutorial for your language.
  • Buy a book.  My favorite Perl book is Perl Best Practices by Damian Conway.  This book is a must have for all Perl programmers!  I don't have much experience with Python books, so I would recommend looking at book reviews before making a purchase.

Learn a Statistical/Graphing Language
  • Pick a language for doing statistical operations and building figures.  Languages like R and Matlab are prime choices for both statistics and graphics.  Both languages have their strengths and weaknesses, but I personally prefer R.  If you choose R I highly recommend using the ggplot2 library for building figures.  
  • Complete an online tutorial for your language.
  • Give up Excel.  Excel is a powerful program but lacks the flexibility of computer languages like R and Matlab.  While there is a steeper learning curve for R and Matlab, you will substantially enhance your ability to do statistical analyses and build graphics by getting away from Excel.

Learn Basic Bioinformatics Procedures and Corresponding Software Tools

For example:
This is only a small list of procedures and tools primarily focusing on DNA sequence analysis.  For a more comprehensive list see OMICtools.


Find a problem

I strongly encourage new bioinformaticians to find some real data to do meaningful science using the above principles and skills.  If you don't personally have data I recommend downloading data from a public repository (i.e. Genbank).  A similar alternative would be to choose a paper that uses a procedure you are interested in learning and recapitulate the results.  

Wednesday, March 26, 2014

Predicting Full-Length Ribosomal Gene Sequences

Introduction

The 16S ribosomal gene has been used extensively in biology for distinguishing relatedness between species.  This gene has regions of DNA that are highly conserved among almost all living organisms and other regions that have high DNA sequence variability.  The conserved regions are ideal for building PCR primers that can amplify DNA from many different organism.  The variable regions that are amplified using these conserved primers can be used to determine the relatedness between two or more organisms.  Closely related species typically have much more similar DNA sequences than distantly related species.

Typically PCR is used to amplify a portion of the 16S ribosomal gene for sequencing.  However, whole genome sequences or whole metagenome sequences also contain short DNA reads originating from the 16S gene.  These reads can be separated from the pool of other genomic reads and assembled into the entire 16S gene.  EMIRGE (Miller, et al. 2011) is an algorithm for reconstructing full-length ribosomal genes from short read DNA sequences.


EMIRGE

EMIRGE reconstructs full-length ribosomal genes from short read DNA sequences.  It first maps reads to a database of known 16S genes such as the SILVA or greengenes database.  After the initial mapping, EMIRGE estimates the probability that a given read was generated from the reference to which it mapped.  Based on these probability estimates, reference sequences are changed to reflect the 16S sequences that are likely to be represented by the set of reads.  Reads are then remapped to the adjusted 16S sequence database and the processes is repeated until an equilibrium is achieved.  The resulting database of 16S sequences reflect the likely 16S genes represented by the input set of short reads.

This software was primarily built to infer the set of 16S genes from whole metagenome reads.  However, it can also be used to infer the single 16S gene from genomic sequences from a single isolate.  Full-length 16S genes are difficult to assemble even when only reads from a single genome are considered.

In the Dangl lab, we use EMIRGE to predict full-length 16S genes from reads generated from a single genome of bacteria.  An example of the EMIRGE command we use is:

emirge.py my_output_dir -1 fwd_reads.fastq -2 rev_reads.fastq -b SSURef_NR99_115_tax_silva_formated -f SSURef_NR99_115_tax_silva_formated.fasta -i 600 -s 1000 -l 250

The descriptions of each parameter are below:

my_output_dir: the output and working directory for EMIRGE.

-1: the forward or single-end genomic sequencing reads

-2: the reverse end genomic sequencing reads

-b: the bowtie index of the 16S sequence database

-f: the fasta file of the 16S sequence database

-i: insert size of paired-end reads

-s: standard deviation of insert size for paired-end reads

-l: max length of reads


Other Details

EMIRGE uses bowtie to map reads to the reference database.  To build the bowtie index of the reference database the following command was used:

bowtie-build SSURef_NR99_115_tax_silva_formated.fasta SSURef_NR99_115_tax_silva_formated

Also, the database downloaded from SILVA had to be reformatted using this Perl script.  This script requires BioUtils.

Monday, March 24, 2014

2014 JGI Users Meeting Notes

Here are some notes from a few of the speakers at the JGI Users meeting in California.  In general the speakers were fantastic.  Some general themes of the conference include:  single-cell genomics, synthetic biology, fungal metagenomics, and metabolics.  A person take-home message for me was the need for creative biological solutions to common issues that the human race currently faces or will face in the near future.

Mark Ackermann (opening keynote) – A Single Cell Perspective on Bacterial Interactions
- Focused on phenotypic heterogeneity, when identical cells have different functional profiles.
- Most genes don’t have clonal variation but in the ones that do how is that heterogeneity important for the community.
- Salmonella is an example of phenotypic heterogeneity.  One cell type causes inflammation and one uses the inflammation response to reproduce and cause full infection.
- Different cell types survive better in different environmental conditions.
- Another example of phenotypic heterogeneity is in alpine lakes where there are generally large amounts of ammonium that bacteria can use as a nitrogen source.  However, there are some cells that fix their own nitrogen in the event that ammonium runs out.
- preliminary data show that neighboring cells are more likely to be of the same cell type.

Mary Berbee – Pectinases link Early Fungal Evolution to the Land Plant Lineage
- Sequenced early divergent fungal groups.
- The relationship between the early branching groups is still poorly resolved.
- Showed some cool trees where she had overlaid two trees to highlight difference between the two.  I would like to know what software she used to do this.
- Her trees were based on whole genomes but I’m not sure how she built them.

Rytas Vilgalys – Understanding the Forest Microbiome:  A Fungal Perspective
- Oak and pine share many fungi while populus has more different fungi.
- Soils from the same region are likely to share the same fungi.
- Populus of different genotypes do not assembly different fungi.  At least not nearly as different as fungi from different regions.
- They have isolated ~1,800 fungal isolates.  These isolate represent only ~15% of the isolates that are likely populus endophytes.
- Many fungal isolates stimulate plant growth.
- They are re-inoculating these isolates to confirm they are endophytic.
- Mortierella elongata is an isolate that stimulates plant growth in populus and Arabidopsis thaliana.
- M. elongata also harbors bacterial symbionts (Glomeribacter which are known to affect lipid fermentation and is a sister to Burkholderia.  These bacteria cannot be cultured possibly because they rely so heavily on the host for nutrients). 
- M. elongata migrate to the roots.
- Different genes are expressed in M. elongata grown in culture than those sampled from the rhizosphere.
- Different genes are expressed in M. elongata inoculated on different hosts.

Eddy Rubin
- Bacterial genes are typically ~900bp.
- In a couple of sequenced genomes they saw average bacterial gene lengths as low as 200bp.  However, when they adjust the codon table by replacing one of the stop codons to code for a glycine predicted genes have an average length of 900bp!  Some bacteria use different codon translations! 
- Natalia Ivanova is a gene annotation specialist they consulted for help in this analysis.
- They found evidence of recoding in lots of other bacteria by looking at sequenced isolates.
- Didn’t find evidence of recoding in archea. 
- They show that phages which use different codon profiles can circumvent host cell machinery to match their codon profile!
- CRISPR regions in bacterial cells often contain phage elements that correspond to different codon profiles.  This is further evidence that phages with different codon profiles can infect cells with canonical codon profiles. 

Nicole Dublier –Metagenomics and Metaproteomic Analyses of Symbioses between Bacteria and Gutless Marine Worms
- Bacteria can use hydrogen to produce more energy than methane.  Nature 2011
- They discovered key genes able to metabolize hydrogen.
- The second half of the talk was about gutless worms living in shallow water.  They completely dependent on bacterial symbionts for feeding and waste excretion. 
- There are species specific symbionts.
- Her proteomics data yield more obvious features than comparative genomics.  As an example she shows how one isolate contains a protein that does the function of 3 different proteins in the canonical Calvin Cycle.  DNA sequencing confirmed this observation but would have been a “needle-in-a-haystack” for a comparative genomics project.  This work published in PNAS.

Erin Nuccio – Mapping Soil Carbon from Cradle to Grave:  Using Omics and Isotope Analyses to Identify the Microbial Blueprint for Root-enhanced Decomposition of Organic Matter.
- The general question is how do microbes transform and stabilize root carbon in soil.
- Carbon can affect nitrogen rates.
- Plants fix carbon for microbes in the soil.
- Looking at the rhizosphere over time it gradually deviates from bulk soil in carbon levels at time points of 3, 6, 9, and 12 weeks.
- Some preliminary data show that bacteria prefer carbon excreted by plant over as an energy source over nitrogen liter material (ie material artificially added to the system).

Michael Fischbach – A Gene-to-Molecule Approach to the Discovery and Characterization of Natural Products
- Discovers natural gene products.  By gene products I think he means functional protein units.
- Undiscovered gene products are often coded by clusters of genes.
- Has some type of algorithm to computationally discover these clusters that may produce unknown gene products.
- Lots of his most interesting clusters were found on human associated microbes.
- Discovered several oligosaccharide clusters.  These bacteria were very difficult to work with but these clusters and the functions they provide to the human host are of high interest.
- The general observation of this study was that microbes in our gut are making products for which we have no idea what they are or how they function.  It’s like taking several prescription drugs for your entire life!  We need to figure out what is going on in there. 

Kelly Matzen – Genetic Control of Mosquitoes
- In the 50’s DDT was used to control mosquito populations and subsequently mosquito born disease such as dengue.  However, DDT is know to be detrimental to the environment in several ways and therefore is being used much less.  We are starting to see diseases like dengue make a comeback in places like Florida and of course in places like Central and South America.
- Right now the most effective control is pesticides. 
- They are releasing massive numbers of sterile male mosquitoes to control (ie reduce) mosquito populations.  This technique has been successfully used before in the United States to control populations of other insects many years ago.
- This technique seems to be working in the small field studies they have been conducting. 
- There is some push back from legislators but in general it seems like good solution.

Cameron Coates – Characterization of Cyanobacterial Hydrocarbon composition and Distribution of Biosynthetic Pathways
- Cyanobacteria produce over 30% of the earth’s oxygen.
- They are very diverse and live in all sorts of habitats on earth.
- They can produce hydrocarbons where are relevant of use of biofuels.  However, they don’t produce large amounts of hydorcarbons.
- They looked at the evolution of cyanobacteria hydrocarbon pathways.  There are two main pathways.  Several clades have both pathways suggesting a large amount of horizontal gene transfer. 
- This work was published in PLOS ONE.

June Medford – Making Better Plants:  Synthetic Approaches in Plant Engineering
- They created a biological input/output system.  This allows for some external factor to cause a reaction that can be observed in the plant.
- They use a pariplasmic binding protein as the input signal because it can quickly defuse through the cell wall and are then translocated to the nucleus to transcriptionally regulate some response. 
- They can theoretically use this system as a flag for pollutants or other dangers that we currently use very expensive technology to detect.
- They are currently developing a system to detect TNT where the response signal of the plant is to turn white.  This system can detect traces of TNT 10x smaller than a dog!  There are still some kinks to work through like response time.  But looks like a very promising system.  This idea has countless unexplored applications!

Kankshita Swaminathan – Genome Biology of Miscanthus
- Miscanthus is in the same clade as sugar cane, corn, and sorghum.  These plants have been amenable to breading.
- The genomic sequence of sorghum is very close to Miscanthus except that Miscanthus has had a whole genome duplication event.
- In the winter all the nutrients migrate to the rhizome leaving only the stalk above ground.  The stalk is the most important element for biofuels and can be harvested without significantly depleting soil nutrients.

Annalee Newitz (closing keynote) - How Humans Will Survive a Mass Extinction
- Humans have a very good chance of surviving a mass extinction because we are very adaptable.  However, our focus should be how we can preserve the diversity of the earth as it is now.
- A mass extinction is when greater than 70% of the earth's species are killed.
- Five mass extinctions have occurred in the history of the earth.  Perhaps the largest was caused by cyanobacteria because they released large amounts of oxygen into the atmosphere.  Close to 90% of species became extinct as a result.
- Climate change is inevitable regardless of wither or not humans are the cause.
- The questions we should be asking are:  how can we respond to these changing climates and what can we do to preserve the world as we know it.
- Space travel seems like an important step in human survival.  



Tuesday, February 25, 2014

The Pan/Core/Accessory Genome

Introduction
The term "pan genome" was coined in 2005 by Tettelin in a paper describing the genomes of eight pathogenic Streptococcus strains.  The pan genome is the set of all unique genes from a set of genomes (ie gene union).  The core genome is the set of genes found in each genome (ie gene intersection).  The accessory genome is the genes unique to a particular genome (ie strain specific genes).  

Previous Studies
Read et al. (2012) discusses the pan and core genome of phytoplancton.  They estimate the size of the E. huxleyi pan genome to be large because there are several thousand genes in the reference that are missing from all of the three well-sequenced isolates.  This is definitely more of a core genome paper (the analysis of which is easy when you have a reference).  Perhaps a better way of showing the diversity of the pan genome is rarefaction curves on the number of homologous genes or perhaps k-mer content.  The lead investigator, Igor Grigoriev, is the fungal genomics lead investigator at the JGI.

Pan and Core Genome Dynamics
In general, as more genomes are added to the analysis set, the core-genome shrinks and the pan-genome grows.  Collins (2012) describe these dynamics in their Molecular Biology and Evolution publication by using an infinitely many genes (IMG) model.  In summary, the Collins IMG model is based on the idea of three types of gene classes:  core, shell, and cloud genes.  Core genes are those found in all genomes, shell genes are gained and lost from genomes at a relatively slow rate, and cloud genes are rapidly gained and lost from genomes.  Empirical data from a set of Bacillaceae genomes support the Collins IMG model.  The Collins IMG model can be used to predict the size of core- and pan-genomes.  In the Dangl lab this has become a question of substantial interest for determining which relevant clades require more isolate genomes for more robust functional genomics analyses.

Core vs Accessory
Given a set of genes from various genomes what gene set is most interesting?  Of course this question will most heavily depend on previous knowledge of the genomes in question.  In the Dangl lab we are interested in microbes that inhabit the endosphere (inner plant root).  Because the set of microbes living inside these roots exhibit a different profile than surrounding soil, one could hypothesize that a single or small set of genes are responsible for a microbes ability to inhabit the inner root.  Under this hypothesis the core-genome would be of particular interest because it should contain these genes.  However, in practice the core genomes is primarily composed of common cellular functions ubiquitous to all bacteria.

Because the core-genome is likely to reveal nothing of specific interest the accessory-genome seems most interesting.  The genes in this set alludes to what makes a particular genome functionally distinct and interesting.  They get at questions like:  "What functions does a particular bacteria provide to the community."

Software for pan-genome analysis
The primary aim in any pan-genome analysis is grouping orthologous genes from different genomes.  To do this many pan-genome pipelines utilize algorithms and databases such as GO, COG, KEGG, eggNOG, Pfam, etc.

Here are some notes on various pipelines developed for pan-genome analyses.
  • GET_HOMOLOGUES (my recommendation)
    • This is my program of choice for pan/core/accessory genome analyses.  
    • Fantastic documentation
    • Options for bidirectional blast hit (BDBH), COGtriangle, and/or orthoMCL algorithms for building clusters of orthologous groups.
    • Builds clear figures
    • Several options for powerful downstream analyses. 
    • Parallelization options
  • Panseq (Laing, 2010)
    • Nice web-based interface.  
    • Seems to work (as opposed to some of the following programs)
    • Output formats are not as user friendly or as concise get_homologues
    • Job completion email never comes in.  Be sure to save the link somewhere.
    • More suited for small, quick analyses
  • PGAT (Brittnacher, 2011)
  • PGAP (Zhao, 2011) 
    • Did not install because it requires the old version of blast (blastall)
  • PanFunPro (Lukjancenko, 2013)
    • Still in the development stage.   I had a quick look at the source code and there were some things didn't make sense.  I'll give the developers a little more time to work out the kinks.  
    • The installation can take some time because of some large dependencies (eg. InterProScan).  Furthermore, the installation for the PanFunPro Perl scripts could be streamlined using a tool like Module::Build.  However, the installation instructions for it's dependencies are well written making installation remarkable easy. 
  • PanOCT (Fouts, 2012)
    • Primarily an algorithm for determining homology between a set of genes from 2 or more eukaryote genomes.  
    • Considers conservation of neighboring genomic regions for determining homology.  The basic idea is that two genes are truly homologous (as opposed to paralogous) will be situated in the same genomic location.

Wednesday, February 12, 2014

A Brief Introduction to Sequence Assembly

Assembly Background
Sequence assembly is one of the overarching challenges in bioinformatics.  To understand the assembly problem it helps to understand some basics of DNA sequencing.  Consider a bacterium having a genome comprised of a single 5 megabase (5 million base pairs) chromosome.  Ideally, sequencing machines would start at the beginning of the chromosome and read each of the 5 million base pairs until arriving at the end.  Unfortunately, the current technology is limited to reading sequences between 30 and ~10,000+ bases.  The assembly problem is to take these short segments of DNA called reads and overlap them in such a way to recreate the original 5Mb chromosome.

To illustrate this consider the set of character strings below that come from a quote by Theodore Roosevelt (spaces have been replaced with "_" for clarity).  Can you put the pieces together to find out what it says?


You should end up with something that looks like this:


This is more or less what assembly programs attempt to do with DNA.  Some things to notice in the above example:
  • Repeats can be problematic during assembly.  Notice that the word "you" is used twice in this sentence.  Looking at the two character strings "Believe_yo" and "you're_ha" you may have incorrectly merged them to form the character string "Believe_you're_ha" which is an incorrect assembly.  DNA repeats are common in genomes and can fragment assemblies or cause assembly mistakes.
  • Longer reads can help with the repeat problem.  For example, given only two long character strings "Believe_you_can_and_you're" and "can_and_you're_halfway_there" it is much easier to unambiguously assembly the quotation despite the fact that the word "you" is used twice.
  • Sequencing errors complicate assembly.  For example, if the character string "halfway" was sequenced as "calfway" there would be no way to finish the assembly correctly because "you're_ha" does not overlap with "calfway."  
  • Coverage (i.e. the number of times a character is represented in the assembly) helps distinguish sequencing errors.  For example, if we sample from the quote (or in the genome in the case of DNA) many times and see the character string "halfway" 100 times and the character string "calfway" only once we can assume that "calfway" was incorrectly sequenced.  

De Novo vs Mapping Assembly
De novo is a latin expression meaning "from the beginning" (Wikipedia).  De novo sequence assemblies are build with no external information beyond the raw sequencing reads.  First pass de novo assembles are called "draft" assemblies because the genome remains fragmented (i.e. discontinuous) and may contain assembly errors.  Extensive resequencing and curation is generally required to "complete" or "finish" a genome assembly.

Alternatively, mapping assembly uses a reference sequence as an anchor to orient sequenced reads.  After reads are ordered based on their location in the reference sequence a consensus sequence is generated from all the mapped reads.  The consensus sequence can differ from the reference sequence but differences are generally single base differences scattered throughout the genome.  Mapping assembly is most useful when the reference sequence and sequenced organism are closely related.

In some cases (e.g. metagenomics), a combination of de novo and mapping assemblies may be advantageous.  However, hybrid assembly algorithms and protocols are not well explored.

De Novo Assembly Algorithm Classes
There are two main classes of assembly algorithms used in de novo assembly:  overlap-layout-consensus (OLC) and de bruijn graph (DBG).  Similar to the above example, OLC first finds reads with overlapping ends, builds a layout graph based on these overlaps, and lastly generates a consensus sequence as the graph is traversed.  OLC was the first assembly method developed and works well with long-read, low-coverage sequencing technologies like Sanger (and possibly PacBio).

DBG based assemblers convert the set of reads into a set of k-mers (i.e. short DNA sequences of length k).  These k-mers are then used to build a de bruijn graph from which the genomic sequence is inferred.  DBG assemblers work well with high-coverage sequencing methods like Illumina and Ion Torrent.  


Useful Links

Friday, November 15, 2013

New SPAdes Assembler Used on MiSeq Reads for 6 Burkholderia Genomes

Introduction

This post introduces a practical assembly strategy using the new SPAdes assembler with Illumina MiSeq reads.  Utilizing the sequence analysis tools sickle, FLASH, SPAdes, and in-house Perl scripts we assemble 6 Burkholderia genomes to draft status using both paired-end (PE) and mate-pair (MP) reads.  This exercise demonstrates the practicality of using the Illumina MiSeq for small scale assembly projects.

The SPAdes Assembler

Many popular de novo assemblers, including SPAdes, rely on a computational data structure called a de Bruijn graph.  SPAdes uses a multisized de Bruijn graph to balance trade-offs between small and large k-mer sizes.  A smaller k causes more repeat regions to collapse into the same node tangling the graph.  However, larger k values can fragment the graph in low coverage regions.  Multisized de Bruijn graphs allow the size of k to vary based on regional coverage depth.  SPAdes also has a unique method for handling PE reads.  Typically PE information is incorporated after the initial assembly by mapping PE reads back to assembled contigs.  When corresponding pairs map to the ends of different contigs those contigs can be merged into a single scaffold.  SPAdes improves on this by incorporating PE distance information directly into the de Bruijn graph.  Because of these features SPAdes can be run with any number of single-end (SE), PE, or MP read files.

Assembling Burkhoderia Genomes with SPAdes

Burkholderia is a genus of bacteria with strains having a variety of environmental effects ranging from plant growth promoting activity to human and plant pathogenesis.  Because our lab is particularly interested in phenotypes of association between bacteria and plants, we sequenced 6 strains of burkholderia isolated from the endophyte compartment (i.e. inner root) of A. thaliana using two runs on our Illumina MiSeq machine.  The first run was done using a standard PE library, and the second was done using a MP library.


The pipeline used for assembly is outlined and described below.  



FLASH
FLASH is a software package the merges PE reads into a single read.  In PE sequencing there is a distribution of DNA fragment lengths from which ends are sequenced (Supp Fig 1).  As reads continue to lengthen more PE reads can potentially overlap.  The main advantage to merging the overlapping reads is the error correction step.  Furthermore, longer reads tend to yield better assemblies because of their ability to span repetitive regions that fragment assemblies.  While longer merged reads may span repetitive regions it is still possible that spanned repeats are incorrectly assembled.  In terms of draft assembly this is a minor point.

After experimenting with various combinations of parameters I decided on the following FLASH parameters:

-m 25 -r 250 -f 500 -s 100

Adapter Trimming
A preprocessing step is required to trim MP reads at the Illumina junction adapter.  For details on how to do this and why it is necessary see this post.  

Sickle
Graph assemblers can be sensitive to erroneous k-mers produced from sequencing errors which convolute the graph with false nodes and connections.  Rigorous quality filtering prior to assembly reduces the number of false k-mers thereby improving assembly.  Sickle is a software package for quality filtering which relies on a sliding window algorithm over the read quality scores.  If the average quality of a window drops below a given minimum the remaining bases are trimmed off.  Any trimmed or untrimmed reads shorter than a given length are removed.  Sickle also produces a file containing high quality reads where the corresponding pair failed the quality check.  These reads can be used in downstream analysis as single-end reads.  

The parameter setting that I used when running Sickle are:

-t sanger -l 127 -q 20

The '-t' parameter specifies the quality value type.  The '-l' parameter specifies the minimum read length and '-q' is the minimum average window quality score.  At this stage it is important to consider the '-l' parameter if running SPAdes is your next step.  The SPAdes parameter '-k' takes a list of k-mers from which to build a multi-sized graph.  For long Illumina reads the SPAdes manual advises '-k' be set to '21,33,55,77,99,127.'  Because the longest k-mer is 127bp reads in you dataset should be at least that long.  To allow for shorter reads when running Sickle simply lower the '-l' parameter and ensure the largest '-k' value is shorter than or equal to that value.  Because many of the MP reads are already short I ran sickle with the -l parameter set to 75 and the max SPAdes -k parameter set to 75.

I should note that SPAdes has a built-in correction algorithm--hammer.  Hammer looks at k-mer frequency to identify potential false k-mers.  However, I recommend running Sickle first for 2 reasons:
  1. Sickle is a much faster correction algorithm than hammer.  Hence, limiting the number of false k-mers that hammer must identify speeds up SPAdes.
  2. Sickle uses quality values.  While quality values are not perfect, they do correlate with the correctness of bases.  If reads have any low-quality k-mers it is possible by chance that some of these k-mers will be the same causing hammer to classify them as a true k-mer.  However, it is likely that one error inside a high quality k-mer will not be trimmed by Sickle, but will be identified by hammer.  This intuition says that a combination of Sickle and hammer is the best regime for quality filtering.
SPAdes
I ran the SPAdes assembler twice for each genome--first with both PE and MP reads and then with only the PE reads.  The second run shows that improvements were made when MP reads are included.  The following command was used to run SPAdes with both PE and MP reads:

spades.py -k 21,33,55,75 -t 5 --careful -s merged_PE_reads.fastq -1 fwd_PE_sickled_reads.fastq -2 rev_PE_sickled_reads.fastq -s PE_sickled_single_reads.fastq --mp1-1 fwd_MP_sickled_reads.fastq --mp1-2 rev_MP_sickled_reads.fastq --mp1-s singles_MP_sickled_reads.fastq --mp1-rf -o spades_out

Using only PE reads (i.e. the merged reads, fwd high quality reads, reverse high quality reads, and single high quality reads) SPAdes was run using the following command:

spades.py -k 21,33,55,77,99,127 -t 5 --careful -s out.extendedFrags.fastq -1 out.notCombined_1_qc.fastq -2 out.notCombined_2_qc.fastq -s out.notCombined_singles_qc.fastq   -o spades_out

QUAST
When the assembly is complete QUAST can be used to compare and view assembly results.  For examples see the Results section of this post.


Results

Assembly
All 6 genomes assembled very well when both PE and MP reads were used.


Does merging PE reads help?
To illustrate how merging reads prior to assembly can improve assembly, I assembled one Burkholderia genome using three different input sets:
  1. Merged and unmerged reads
  2. Raw PE reads with no merging
  3. Only reads that merge

Under the assumption that improved assembly statistics (e.g. number of contigs, N50, etc.) indicate a better assembly, a combination of merged and unmerged reads (option 1) yields the best assembly.

Does adding MP reads help?
When assembling with only PE reads the resulting assemblies are worse than when MP reads are included.


To test if the assembly improvements were primarily a factor of increased coverage or information incorporated in MP reads, I reassembled the PE+MP reads for B1/BMP1 treating the MP reads as SE reads.  Despite MP reads substantial increase in coverage (52 to 103), when treated as SE reads only minor assembly improvements are observed.  This suggests that distance information provided by MP reads can have a greater effect on improving assemblies than additional coverage.


Conclusions

Using a combination of FLASH, Sickle, SPAdes, and in-house Perl scripts we assemble 6 Burkholderia genomes to draft status.  This exercise demonstrates the practicality of the Illumina MiSeq for small scale sequencing projects consisting of a handful of bacterial genomes.  A single PE MiSeq run may be sufficient to assembly a small number of bacterial genomes, however the additional coverage and pairing information of MP reads can improve PE assemblies.  Merging overlapping PE reads and quality trimming before assembly can increase N50 and other assembly statistics.


Supplemental Figures
Figure 1 -- PE insert size distribution for one sample (B1)