Sequencing Informatics

Next Generation DNA Sequencing Resources

//

Alignment for NGS

Alignment for NGS: Tutorial

NGS Sequencing Informatics

Stuart Brown, Stratos Efstathiadis

NYU Center for Health Informatics & Bioinformatics

Alignment

Most NGS experiments involve sequence reads collected from known genomes: gene expression, variants, ChIP-seq, etc.
First step in data analysis is to align reads to the reference genome.
This is computationally demanding due to the huge volume of NGS data and the large size of the reference genome (human = 3.2 Gb)

Tutorial Data

We can use a small data set for this alignment tutorial that was originally created by the Galaxy team for their ChIP-seq tutorial.
It contains a small number of reads that map just to one part of Chromosome 19 on the Mouse genome.

https://main.g2.bx.psu.edu/u/james/p/exercise-chip-seq
The small size of this data set will make all of the compute operations go more quickly. A real data set will take MUCH LONGER at every step, and use more
storage space.

Unix

For this tutorial, we will assume that you have access to a computer with a Unix or Linux operating system installed plus the standard software compilers
and libraries needed to install and run software.
This can be done on any Macintosh computer, your own Unix server, or a “cloud” server such as the Amazon Elastic Computer Cloud (EC2)
See our Unix and HPC tutorials for more details on setting up your system and basic Unix commands

NGS Alignment

Pre-NGS alignment algorithms such as Smith-Waterman and BLAST are both too slow and too sensitive (able
to align sequences with many mismatches)
NGS reads should be near-perfect matches to the reference genome.
Need to process hundreds of millions of reads in a few hours.
- time is critical because alignment is used for QC of NGS experiments, lab may be waiting
  before running next sample

There are many NGS aligners

Alignment Software

Illumina provides the ELAND alignment tool
- fast, max of 2 mismatches in 32 base ‘seed’
- uses a hash index of the NGS reads (query data), which uses more RAM (and runs slower) as machines produce more data
- constantly improving its handling of indels
BWA
and Bowtie are popular open source tools based on the Burrows Wheeler Transformation
- fast
- only works well on sequences with <3% difference from the reference genome
- better gapped alignments (indels) than ELAND
- provides a “mapping quality” score for the alignment of each read

Burrows Wheeler Transformation

Invented by David Wheeler in 1983 (Bell Labs).

A Block Sorting Lossless Data Compression Algorithm
. Systems Research Center Technical Report No 124. Palo Alto, CA: Digital Equipment Corporation, Burrows M, Wheeler DJ. 1994
Originally developed for compressing large files (bzip2)
Lossless, fully reversible
Align reads to a BW transformed suffix array index of the reference genome, which enables fast exact matching
Short-read algorithm: alters the read sequence to match the reference exactly.
Results in great gains in speed, efficiency, and reduces RAM requirement vs. other types of indexes

BWA Features

Compresses the genome as a BWT index
Identical sequences (genome repeats) only occur once in the index, reads that match the genome in multiple locations get a mapping quality of 0
Fast and moderate memory footprint (<4GB)
Gapped alignment for both SE and PE reads
Effective pairing to achieve high alignment accuracy; suboptimal hits considered in pairing.
Non-unique reads are placed randomly with a mapping quality of 0; all hits can be output in a concise format.
The default configuration works for most typical input.
Automatically adjust parameters based on read lengths and error rates.

Install BWA software

Prepare Reference Genome Index

Obtain a file that contains the reference genome of your organism.
Download from UCSC Genome Browser
http://hgdownload.soe.ucsc.edu/downloads.html
[Note: to align the ch19 tutorial data, you will need to download the mouse mm9 genome, or just mouse chromosome 19]
Unzip files for individual Ref Geneome chromosomes and join them together:

tar zvfx chromFa.tar.gz

cat *.fa > ref.fa
Index the reference genome :

bwa index -a bwtsw ref.fa

Align your FASTQ file

Your data should be in FASTQ format
Align to the index:

bwa aln ref.fa data.fastq.gz > data.sai
Export the alignment in SAM format

bwa samse ref.fa data.sai data.fastq.gz > data.aln.sam
For paired end data, align two FASTQ files to index, then use sampe to export one SAM file

bwa sampe ref.fa read1.sai read2.fq.gz > pe_aln.sam

Alignment Output Format

SAM format is a human readable text file
Shows each read, quality string for bases, alignment position on the genome, alignment quality score, and any mismatches or indels of the read vs reference
(using CIGAR code).
All current NGS aligners now produce output in the BAM format, a compressed binary file, indexed by genome position.

SAM format: aligned to reference genome

SAMTools commands

Download and install SAMTools

http://sourceforge.net/projects/samtools/files/samtools/
Make the BAM file:

samtools view –bt ref.fa –o data.bam data.sam
Sort the BAM file:

samtools sort data.bam data.sorted.bam
Index the BAM file:

samtools index data.sorted.bam data.st_index.bam

Visualization

BAM file format also simplifies visualization of NGS data.
Index allows viewer to quickly find and load only read data for a specific genomic interval.
Integrative Genomics Viewer (IGV) from Broad Institute is our current favorite.
Use the Java Webstart to download and run IGV (use the 1.2 GB version): http://www.broadinstitute.org/igv/

Note: if you are using a remote server (your own or Cloud), then you will need to download the BAM and .bai files from the server to your desktop/laptop in
order to visualize.
IGV is easy to learn to use.
It has a nice Users Guide:

http://www.broadinstitute.org/software/igv/UserGuide
•Very briefly, load the BAM file and choose the genome (mouse mm9 for the tutorial data)