Page Comparison

Tip

title	Reservations

Use our summer school reservation (

CoreNGSday4

CoreNGS) when submitting batch jobs to get higher priority on the ls6 normal queue

today:

sbatch --reservation=CoreNGSday4 <batch_file>.slurm

Table of Contents

Overview

Image Removed

After raw sequence files are generated (in FASTQ format), quality-checked, and pre-processed in some way, the next step in many NGS pipelines is mapping to a reference genome.

For individual sequences it is common to use a tool like BLAST to identify genes or species of origin. However a normal NGS dataset will have tens to hundreds of millions of sequences, which BLAST and similar tools are not designed to handle. Thus a large set of computational tools have been developed to quickly align each read to its best location (if any) in a reference.

Even though many mapping tools exist, a few individual programs have a dominant "market share" of the NGS world. In this section, we will primarily focus on two of the most versatile general-purpose ones: BWA and Bowtie2 (the latter being part of the Tuxedo suite which includes the transcriptome-aware RNA-seq aligner Tophat2 as well as other downstream quantifiaction tools).

Stage the alignment data

First connect to ls6.tacc.utexas.edu and start an idev session. This should be second nature by now

Code Block

language	bash
title	Start an idev session

idev -p normal -m 180 -A OTH21164 -N 1

Then stage the sample datasets and references we will use.

Code Block

language	bash
title	Get the alignment exercises files

mkdir -p $SCRATCH/core_ngs/references/fasta
cp $CORENGS/references/*.fa     $SCRATCH/core_ngs/references/fasta/

mkdir -p $SCRATCH/core_ngs/alignment/fastq
cp $CORENGS/alignment/*fastq.gz $SCRATCH/core_ngs/alignment/fastq/
cd $SCRATCH/core_ngs/alignment/fastq

These are descriptions of the FASTQ files we copied:

...

Reference Genomes

Before we get to alignment, we need a reference to align to. This is usually an organism's genome, but can also be any set of names sequences, such as a transcriptome or other set of genes.

Here are the four reference genomes we will be using today, with some information about them. These are not necessarily the most recent versions of these references (e.g. the newest human reference genome is hg38 and the most recent miRBase annotation is v21. (See here for information about many more genomes.)

...

Searching genomes is computationally hard work and takes a long time if done on linear genomic sequence. So aligners require that references first be indexed to accelerate lookup. The aligners we are using each require a different index, but use the same method (the Burrows-Wheeler Transform) to get the job done.

Building a reference index involves taking a FASTA file as input, with each contig (contiguous string of bases, e.g. a chromosome) as a separate FASTA entry, and producing an aligner-specific set of files as output. Those output index files are then used to perform the sequence alignment, and alignments are reported using coordinates referencing names and offset positions based on the original FASTA file contig entries.

We can quickly index the references for the yeast genome, the human miRNAs, and the V. cholerae genome, because they are all small, so we'll build each index from the appropriate FASTA files right before we use them.

hg19 is way too big for us to index here so we will use an existing set of BWA hg19 index files located at:

Code Block

language	bash
title	BWA hg19 index location

/work/projects/BioITeam/ref_genome/bwa/bwtsw/hg19

...

The BioITeam maintains a set of reference indexes for many common organisms and aligners. They can be found in aligner-specific sub-directories of the /work/projects/BioITeam/ref_genome area. E.g.:

Code Block

language	bash

/work/projects/BioITeam/ref_genome/
   bowtie2/
   bwa/
   hisat2/
   kallisto/
   star/
   tophat/

Exploring FASTA with grep

It is often useful to know what chromosomes/contigs are in a FASTA file before you start an alignment so that you're familiar with the contig naming convention – and to verify that it's the one you expect. For example, chromosome 1 is specified differently in different references and organisms: chr1 (USCS human), chrI (UCSC yeast), or just 1 (Ensembl human GRCh37).

A FASTA file consists of a number of contig name entries, each one starting with a right carat ( > ) character, followed by many lines of base characters. E.g.:

Code Block

>chrI
CCACACCACACCCACACACCCACACACCACACCACACACCACACCACACC
CACACACACACATCCTAACACTACCCTAACACAGCCCTAATCTAACCCTG
GCCAACCTGTCTCTCAACTTACCCTCCATTACCCTGCCTCCACTCGTTAC
CCTGTCCCATTCAACCATACCACTCCGAACCACCATCCATCCCTCTACTT

How do we dig out just the lines that have the contig names and ignore all the sequences? Well, the contig name lines all follow the pattern above, and since the > character is not a valid base, it will never appear on a sequence line.

We've discovered a pattern (also known as a regular expression) to use in searching, and the command line tool that does regular expression matching is grep (general regular expression parser). Read more about grep here: Advanced commands: grep.

Regular expressions are so powerful that nearly every modern computer language includes a "regex" module of some sort. There are many online tutorials for regular expressions, and several slightly different "flavors" of them. But the most common is the Perl style (http://perldoc.perl.org/perlretut.html), which was one of the fist and still the most powerful (there's a reason Perl was used extensively when assembling the human genome). We're only going to use simple regular expressions here, but learning more about them will pay handsome dividends for you in the future.

Here's how to execute grep to list contig names in a FASTA file.

Code Block

language	bash
title	grep to match contig names in a FASTA file

# If you haven't staged the fasta files
cds
mkdir -p core_ngs/references/fasta
cd core_ngs/references/fasta
cp $CORENGS/references/*.fa .

cd $SCRATCH/core_ngs/references/fasta
grep -P '^>' sacCer3.fa | more

Notes:

The -P option tells grep to Perl-style regular expression patterns.
- This makes including special characters like Tab ( \t ), Carriage Return ( \r ) or Linefeed ( \n ) much easier that the default POSIX paterns.
- While it is not required here, it generally doesn't hurt to include this option.
'^>' is the regular expression describing the pattern we're looking for (described below)
sacCer3.fa is the file to search.
- lines with text that match our pattern will be written to standard output
- non matching lines will be omitted
We pipe to more just in case there are a lot of contig names.

Now down to the nuts and bolts of the pattern: '^>'

First, the single quotes around the pattern – this tells the bash shell to pass the exact string contents to grep.

As part of its friendly command line parsing and evaluation, the shell will often look for special characters on the command line that mean something to it (for example, the $ in front of an environment variable name, like in $SCRATCH). Well, regular expressions treat the $ specially too – but in a completely different way! Those single quotes tell the shell "don't look inside here for special characters – treat this as a literal string and pass it to the program". The shell will obey, will strip the single quotes off the string, and will pass the actual pattern, ^>, to the grep program. (Note that the shell does look inside double quotes ( " ) for certain special signals, such as looking for environment variable names to evaluate. Read more about Quoting in the shell.)

So what does ^> mean to grep? We know that contig name lines always start with a > character, so > is a literal for grep to use in its pattern match.

We might be able to get away with just using this literal alone as our regex, specifying '>' as the command line argument. But for grep, the more specific the pattern, the better. So we constrain where the > can appear on the line. The special carat ( ^ ) metacharacter represents "beginning of line". So ^> means "beginning of a line followed by a > character".

Exercise: How many contigs are there in the sacCer3 reference?

mkdir -p $SCRATCH/

.

Expand

title	Setup (if needed)

Code Block

language	bash
title	Get the alignment exercises files

Code Block

language	bash
title	Request an interactive (idev) node

# Request a 180 minute interactive node on the normal queue using our reservation
idev -m 120 -N 1 -A OTH21164 -r CoreNGS
idev -m 120 -N 1 -A TRA23004 -r CoreNGS

# Request a 120 minute idev node on the development queue 
idev -m 120 -N 1 -A OTH21164 -p development
idev -m 120 -N 1 -A TRA23004 -p development

Code Block

language	bash
title	Submit a batch job

# Using our reservation
sbatch --reseservation=CoreNGS <batch_file>.slurm

Note that the reservation name (CoreNGS) is different from the TACC allocation/project for this class, which is OTH21164.

File Name	Description	Sample
Sample_Yeast_L005_R1.cat.fastq.gz	Paired-end Illumina, First of pair, FASTQ	Yeast ChIP-seq
Sample_Yeast_L005_R2.cat.fastq.gz	Paired-end Illumina, Second of pair, FASTQ	Yeast ChIP-seq
human_rnaseq.fastq.gz	Paired-end Illumina, First of pair only, FASTQ	Human RNA-seq
human_mirnaseq.fastq.gz	Single-end Illumina, FASTQ	Human microRNA-seq
cholera_rnaseq.fastq.gz	Single-end Illumina, FASTQ	V. cholerae RNA-seq

Reference	Species	Base Length	Contig Number	Source	Download
hg19	Human	3.1 Gbp	25 (really 93)	UCSC	UCSC GoldenPath
sacCer3	Yeast	12.2 Mbp	17	UCSC	UCSC GoldenPath
mirbase v20	Human subset	160 Kbp	1908	miRBase	miRBase Downloads
vibCho (O395)	Vibrio cholerae	~4 Mbp	2	GenBank	GenBank Downloads

alignment type	aligner options	pro's	con's
global with bwa	single end reads: bwa aln <R1> bwa samse paired end reads: bwa aln <R1> bwa aln <R2> bwa sampe	simple to use (take default options) good for basic global alignment	multiple steps needed
global with bowtie2	bowtie2	extremely configurable can be used for RNAseq alignment (after adapter trimming) because of its many options	complex (many options)
local with bwa	bwa mem	simple to use (take default options) very fast no adapter trimming needed good for simple RNAseq analysis the secondary alignments it reports can provide splice junction information	always produces alignments with secondary reads must be filtered if not desired
local with bowtie2	bowtie2 --local	extremely configurable no adapter trimming needed good for small RNA alignment because of its many options	complex – many options

Option	Effect
-l	Specifies the length of the seed (default = 32)
-k	Specifies the number of mismatches allowable in the seed of each alignment (default = 2)
-n	Specifies the number of mismatches (or fraction of bases in a given alignment that can be mismatches) in the entire alignment (including the seed) (default = 0.04)
-t	Specifies the number of threads

Versions Compared

Old Version 188

New Version Current

Key

Overview

Stage the alignment data

Reference Genomes

Exploring FASTA with grep

Overview

Stage the alignment data

Reference Genomes

Exploring FASTA with grep

Exercise #1: BWA global alignment – Yeast ChIP-seq

Overview ChIP-seq alignment workflow with BWA

Introducing BWA

Aligner overview

Aligner overview

Exercise #1: BWA global alignment – Yeast ChIP-seq

Overview ChIP-seq alignment workflow with BWA

Introducing BWA

Building the BWA sacCer3 index

Building the BWA sacCer3 index

Performing the bwa alignment

Performing the bwa alignment

samtools flagstat

Tip
	Make sure your output files are not empty
Double check that output was written by doing ls -lh and making sure the file sizes listed are not 0.