Genome Assembly (velvet) -- GVA2016

Overview

Velvet is a De Bruijn graph assembler works fairly rapidly on short (microbial) genomes. In this tutorial we will use velvet to assemble an E. coli genome from simulated Illumina reads. Genome assembly is quite difficult (though as Oxford Nanopore comes online it will likely get much easier and involve new tools). Genome assembly should only be used when you can not find a reference genome that is close to your own, if you are engaged in metagenomic projects where you don't know what organisms may be present, and in situations where you believe you may have novel sequence insertions into a genome of interest (Note that in this case however you would actually want to grab reads that do not map to your reference genome (and their pair in the case of paired end and mate-pair sequencing) rather than performing these functions on the fastq files you get from the gsaf.

Unfortunatly, at the time of this class, velvet is not available on lonestar5. While we hope and think that this may just be because lonestar5 is very new still and it just hasn't happened yet, there is obviously no guarantee that it will come online as a module or a timeframe for when that may happen. Luckily (but somewhat annoyingly) it is available as a module on stampede. As we have discussed in class, stampede compute nodes are not able to access the BioITeam repositories so you will have to make sure you have copied your files to the locations you want them before you enter an idev node.

Learning Objectives

• Run velvet to perform de novo assembly on fragment, paired-end, and mate-paired data.
• Use contig_stats.pl to display assembly statistics.
• Find proteins of interest in an assembly using Blast.

Table of Contents

Data

First log into stampede using the same log in credentials you have been using for lonestar5. Next, let's copy the fastq read files.

cds
mkdir velvet_tutorial
cp $BI/ngs_course/velvet/data/*/* velvet_tutorial
cp $BI/bin/contig_stats.pl .  # because we are on stampede not lonestar5 we need to copy this file before we start an idev session
cd velvet_tutorial

Now we have a bunch of Illumina reads. These are simulated reads. If you'd ever like to simulate some on your own, you might try using Mason.

login1$ ls
paired_end_2x100_ins_1500_c_20.fastq  paired_end_2x100_ins_400_c_20.fastq  single_end_100_c_50.fastq
paired_end_2x100_ins_3000_c_20.fastq  paired_end_2x100_ins_400_c_25.fastq
paired_end_2x100_ins_3000_c_25.fastq  paired_end_2x100_ins_400_c_50.fastq

There are 4 sets of simulated reads:

Note that these fastq files are "interleaved", with each read pair together one-after-the-other in the file. The #/1 and #/2 in the read names indicate the pairs.

login1$ head paired_end_2x100_ins_1500_c_20.fastq
@READ-1/1
TTTCACCGTTGACCAGCACCCAGTTCAGCGCCGCGCGACCACGATATTTTGGTAACAGCGAACCATGCAGATTAAATGCACCTGCGGGAGCGAGCTGCAA
+
*@A+<at:var at:name="55G" />T@@I&+@A+@@<at:var at:name="II" />G@+++A++GG++@++I@+@+G&/+I+GD+II@++G@@I?<at:var at:name="I" />@<at:var at:name="IIGGI" /><at:var at:name="A4" />6@A,+AT=<at:var at:name="G" />+@AA+GAG++@
@READ-1/2
TTAACACCGGGCTATAAGTACAATCTGACCGATATTAACGCCGCGATTGCCCTGACACAGTTAGTCAAATTAGAGCACCTCAACACCCGTCGGCGCGAAA
+
I@@H+A+@G+&+@AG+I>G+I@+CAIA++$+T<at:var at:name="GG" />@+++1+<at:var at:name="GI" />+ICI+A+@<at:var at:name="I" />++A+@@A.@<G@@+)GCGC%I@IIAA++++G+A;@+++@@@@6

Often your read pairs will be "separate" with the corresponding paired reads at the same index in two different files (each with exactly the same number of reads).

Velvet Assembly

Now let's use Velvet to assemble the reads.

First, you will need to get an idev node and load the velvet module.

 module load velvet

Using velvet consists of a sequence of two commands:

1. velveth - analyzes kmers in the reads in preparation for assembly
2. velvetg - constructs the assembly and filters contigs from the graph

Look at the help for each program.

The <hash_length> parameter of velveth is the kmer value that is key to the assembly process. Choosing it controls the tradeoff between sensitivity (lower hash_length, more reads included, longer contigs) and specificity (higher hash length, less chance of misassembly, more reads ignored, shorter contigs)). There is more discussion about choosing an appropriate kmer value in the Velvet manual and in this blog post.

Velvet has an option of specifying the insertion size of a paired read file (-ins_length). If no size is given, Velvet will guess the insertion size. We'll just have Velvet guess the size.

Velvet also has an option to specify the expected coverage of the genome, which helps it choose how to resolve repeated sequences (-exp_cov). We set this parameter to estimate this from the data. A common problem with using Velvet is that you have many very short contigs and the last line of output tells you that it used 0 of your reads. This is caused by not setting this parameter. The default is NOT auto.

We'll need to create a commands file and submit it to TACC. Let's make the commands file say:

velveth single_out 61 -fastq single_end_100_c_50.fastq && velvetg single_out -exp_cov auto -amos_file yes
velveth pairedc20_out 61 -fastq -shortPaired paired_end_2x100_ins_3000_c_20.fastq paired_end_2x100_ins_1500_c_20.fastq paired_end_2x100_ins_400_c_20.fastq && velvetg pairedc20_out -exp_cov auto -amos_file yes
velveth pairedc25_out 61 -fastq -shortPaired paired_end_2x100_ins_3000_c_25.fastq paired_end_2x100_ins_400_c_25.fastq && velvetg pairedc25_out -exp_cov auto -amos_file yes
velveth pairedc50_out 61 -fastq -shortPaired paired_end_2x100_ins_400_c_50.fastq && velvetg pairedc50_out -exp_cov auto -amos_file yes

If you are assembling large genomes or have high coverage depth data in the future, you will probably need to submit your jobs to the "largemem" queue.

Velvet Output

Explore each output directory that was created by Velvet.

Checking the tail of the Log files in each of the output folders, we see lines like the following:

Median coverage depth = 2.657895
Final graph has 9748 nodes and n50 of 191, max 1427, total 1865207, using 281499/2314900 reads

Median coverage depth = 11.131337
Final graph has 265 nodes and n50 of 127102, max 397974, total 4558511, using 1464201/2314900 reads

Median coverage depth = 11.109244
Final graph has 203 nodes and n50 of 698134, max 1032531, total 4585717, using 1465818/2314900 reads

Median coverage depth = 13.353287
Final graph has 202 nodes and n50 of 698626, max 1139610, total 4602729, using 1758595/2777880 reads

With better read pairs that link more distant locations in the genome, there are fewer contigs, and contigs are are longer, giving us a more complete picture of linkage across the genome.

The complete E. coli genome is about 4.6 Mb. Why weren't we able to assemble it, even with this "perfect" data?

More assembly statistics: contig_stats.pl

The output file stats.txt contains information about every contig in the assembly, but it isn't sorted and can be a bit overwhelming.

You can generate some summary stats and graphs about each assembly using the contig_stats.pl script that we copied from  $BI/bin before luanching the idev session. You probably want to change into the directory of results for a specific assembly before running this command, since it generates several output files in the current working directory. Since you are on $SCRATCH and in a directory that is not in your path simply typing the command will not work. As such you need to explicitly tell the command line to launch the file like you have done with "commands" files in other tutorials.  You will need to copy the PNG files back to your computer to view it using scp.

./contig_stats.pl -plot -tool Velvet contigs.fa

Transfer back several of the different outputs into their own direcotry and being comparing them to determine which library and set of parameters seemed to work best.

What do I do now?

Many choices:

1. Get a better assembly: maybe add a different library size, or go into a detailed genome completion project (commonly called "finishing") using a sequence assembly editor like consed or gap4 or AMOS. (Be careful though, the amount of data in NGS data sets can be very difficult for these programs to deal with, since many were designed for Sanger sequencing reads.) You may have a lot of PCR products to make to close gaps and/or to order and orient scaffolds. consed in particular makes this pretty easy, but it may still consume a lot more time and money than the initial shotgun assembly. You can identify some misassemblies by mapping the original reads to the assembly and then viewing them in IGV to look for discordant mate pairs, for example.
2. Look for things: If you're just after a few homologs, an operon, etc. you're probably done. Most assemblers will be able to take 2x100 data and give you full gene sequences since these are non-repetitive and so assemble well. You can turn the contigs.fa into a blast database (formatdb or makeblastdb depending on which version of blast you have) and start blasting away.

Further Reading

• Examples of using single-end, paired-end, and mate-paired data in Velvet at the GenomeFactory.
• Another velvet tutorial