Introduction to Genome Assembly

Overview

question Questions
  • How do we perform a very basic genome assembly from short read data?
objectives Objectives
  • assemble some paired end reads using Velvet
  • examine the output of the assembly.
requirements Requirements

time Time estimation: 30 minutes

Genome assembly with Velvet: Background

Velvet is one of a number of de novo assemblers that use short read sets as input (e.g. Illumina Reads). The assembly method is based on the manipulation of de Bruijn graphs, via the removal of errors and the simplication of repeated regions.

comment Comment

For information about Velvet, you can check its (nice) Wikipedia page.

For this tutorial, we have a set of reads from an imaginary Staphylococcus aureus bacterium with a miniature genome (197,394 bp). Our mutant strain read set was sequenced with the whole genome shotgun method, using an Illumina DNA sequencing instrument. From these reads, we would like to rebuild our imaginary Staphylococcus aureus bacterium via a de novo assembly of a short read set using the Velvet assembler.

Agenda

In this tutorial, we will deal with:

  1. Get the data
  2. Evaluate the input reads
  3. Assemble reads with Velvet
  4. Collect some statistics on the contigs
  5. Discussion

Get the data

We will now import the data that we will use for the tutorial.

hands_on Hands-on: Getting the data

  1. Create and name a new history for this tutorial.
  2. Import from Zenodo or from the data library the files:
    • Copy the link location (Right-click on the filename then “Copy Link Address”)
    • Open the Galaxy Upload Manager
    • Select Paste/Fetch Data
    • Paste the link into the text field
    • Change the data-type to fastqsanger
    • Press Start
  3. Change the name of the files to mutant_R1 and mutant_R2.

    As a default, Galaxy uses the link as the name of the new dataset. It also does not link the dataset to a database or a reference genome.

    tip Tip: Rename a dataset

    • Click on the galaxy-pencil pencil icon for the dataset to edit its attributes
    • In the central panel, change the Name field
    • Click the Save button
  4. Inspect the content of a file.

    tip Tip: Inspecting the content of a dataset

    • Click on the galaxy-eye (eye) icon next to the relevant history entry
    • View the content of the file in the central panel

    question Questions

    1. What are four key features of a FASTQ file?
    2. What is the main difference between a FASTQ and a FASTA file?

    solution Solution

    1. Each sequence in a FASTQ file is represented by 4 lines: 1st line is the id, 2nd line is the sequence, 3rd line is not used, and 4th line is the quality of sequencing per nucleotide
    2. In a FASTQ file, not only are the sequences present, but information about the quality of sequencing is also included.

The reads have been sequenced from an imaginary Staphylococcus aureus bacterium using an Illumina DNA sequencing instrument. We obtained the 2 files we imported (mutant_R1 and mutant_R2)

question Question

Why do we have 2 files here if we only sequenced the bacteria once?

solution Solution

  1. The bacteria has been sequenced using paired-end sequencing. The first file corresponds to forward reads and the second file to reverse reads.

Evaluate the input reads

Before doing any assembly, the first questions you should ask about your input reads include:

We will evaluate the input reads using the FastQC tool. This tool runs a standard series of tests on your read set and returns a relatively easy-to-interpret report. We will use it to evaluate the quality of our FASTQ files and combine the results with MultiQC.

hands_on Hands-on: FastQC on a fastq file

  1. FastQC tool with the following parameters
    • “Short read data from your current history” to (Multiple datasets) mutant_R1.fastq and mutant_R2.fastq
  2. MultiQC tool with the following parameters
    • “Software name” to FastQC
    • “Result file” to the raw data files generated by FastQC

MultiQC generates a webpage combining reports for FastQC on both datasets. It includes these graphs and tables:

comment Comment

For a fuller discussion of FastQC outputs and warnings, see the FastQC website link, including the section on each of the output reports, and examples of “good” and “bad” Illumina data.

We won’t be doing anything to these data to clean it up as there isn’t much need. Therefore we will get on with the assembly!

Assemble reads with Velvet

Now, we want to assemble our reads to find the sequence of our imaginary Staphylococcus aureus bacterium. We will perform a de novo assembly of the reads into long contiguous sequences using the Velvet short read assembler.

The first step of the assembler is to build a de Bruijn graph. For that, it will break our reads into k-mers, i.e. fragments of length k. Velvet requires the user to input a value of k (k-mer size) for the assembly process. Small k-mers will give greater connectivity, but large k-mers will give better specificity.

hands_on Hands-on: Assemble the reads

  1. FASTQ interlacer tool with the following parameters
    • “Type of paired-end datasets” to 2 separate datasets
    • “Left-hand mates” to mutant_R1.fastq
    • “Right-hand mates” to mutant_R2.fastq

    Currently our paired-end reads are in 2 files (one with the forward reads and one with the reverse reads), but Velvet requires only one file, where each read is next to its mate read. In other words, if the reads are indexed from 0, then reads 0 and 1 are paired, 2 and 3, 4 and 5, etc. Before doing the assembly per se, we need to prepare the files by combining them.

  2. velveth tool with the following parameters
    • “Hash Length” to 29
    • “Input Files”: click on Insert Input Files
    • “file format” to fastq
    • “read type” to shortPaired reads
    • “Dataset” to the pairs output of FASTQ interlacer

    The tool takes our reads and break them into k-mers.

  3. velvetg tool with the following parameters
    • “Velvet Dataset” to the output of velveth
    • “Using Paired Reads” to Yes

    This last tool actually does the assembly.

Two files are generated:

Collect some statistics on the contigs

question Question

  1. How many contigs have been built?
  2. What is the mean, min and max length of the contigs?

solution Solution

  1. 190
  2. To compute this information, we can use the Datamash tool on the 2nd columns (length). Be careful with the first line, the header. As a result, we obtain: 597.82 as mean, 1 as min and 12904 as max. It would mean that the smallest contig has a length of 1 bp, even smaller than k. The length on the 2nd column corresponds to length of the contig in k-mers. This means that the smallest contig has a length of 1k = 29. So to obtain the real length, we need to add k-1 to the length. We then obtain a mean contig length of 625.82 bp, a min contig of 29 bp and a max contig of 12,932 bp.

This table is limitted, but we will now collect more basic statistics on our assembly.

hands_on Hands-on: Collect fasta statistics on our contigs

  1. Quast tool with
    • “Contigs/scaffolds output file” to the output of velvetg
    • “Type of data” to contig
    • “Reference File” to wildtype.fna
    • “Type of organism” to Prokaryotes
    • “Lower Threshold” to 500
    • “Thresholds” to 0,1000

This tool generates 5 output files, but we will focus on the HTML report and the Icarus viewer.

question Question

  1. What is represented in the Icarus viewer?

solution Solution

  1. Icarus is a novel genome visualizer for accurate assessment and analysis of genomic draft assemblies. It draws contigs ordered from longest to shortest, highlights N50, N75 (NG50, NG75) and long contigs larger than a user-specified threshold

The HTML report reports many statistics computed by QUAST to assess the quality of the assembly:

question Question

  1. How many contigs have been constructed?
  2. Which proportion of the reference genome do they represent?
  3. How many misassemblies have been found?
  4. Has the assembly introduced mismatches and indels?
  5. What are N50 and L50?
  6. Is there a bias in GC percentage induced by the assembly?

solution Solution

  1. 190 contigs have been constructed, but only 47 have a length > 500 bp.
  2. The contigs represents 87.965% of the reference genome.
  3. 1 misassembly has been found: it corresponds to a relocation, i.e. a misassembly event (breakpoint) where the left flanking sequence aligns over 1 kbp away from the right flanking sequence on the reference genome.
  4. 8.06 mismatches per 100 kbp and 4.03 indels per 100 kbp are found.
  5. N50 is the length for which the collection of all contigs of that length or longer covers at least half an assembly. In other words, if contigs were ordered from small to large, half of all the nucleotides will be in contigs this size or larger. And L50 is the number of contigs equal to or longer than N50: L50 is the minimal number of contigs that cover half the assembly.
  6. The GC % in the assembly is 33.64%, really similar to the one of the reference genome (33.43%).

Discussion

hands_on (Optional) Hands-on: Rerun for values k ranging from 31 to 101

  1. velveth tool with the same parameters as before except
    • “Hash Length” to a value between 31 and 101
  2. velvetg tool with the same parameters as before
  3. Quast tool with the same parameters as before

We have completed an assembly on this data set for a number of k values ranging from 29 to 101. A few of the assembly metrics appear below.

contigs
Figure 2: Number of contigs in the assembly for various k-mer sizes
largest_contig
Figure 3: Largest contig in each of the assemblies by k-mer size
total_bp
Figure 4: Total number of base pairs in all the contigs for each assembly by k-mer size
n50
Figure 5: N50 metric for each of the assemblies by k-mer size

question Questions

  1. Are there any distinct features in the charts?
  2. Does it look like one assembly might be better than some of the others?

The reasons for these patterns will be discussed in detail in the De Bruijn graph assembly slides and tutorial.

keypoints Key points

  • We assembled some Illumina fastq reads into contigs using a short read assembler called Velvet
  • We showed what effect one of the key assembly parameters, the k-mer size, has on the assembly
  • It looks as though there are some exploitable patterns in the metric data vs the k-mer size.

Useful literature

Further information, including links to documentation and original publications, regarding the tools, analysis techniques and the interpretation of results described in this tutorial can be found here.

congratulations Congratulations on successfully completing this tutorial!