2020 Alignment workflow

1 Overview
- 1.1 Stage the alignment data
  - 1.1.1 Get the alignment exercises files
2 Reference Genomes
- - 2.1.1 BWA hg19 index location
- 2.2 Exploring FASTA with grep
  - 2.2.1 grep to match contig names in a FASTA file
  - 2.2.2 Get the alignment exercises files
3 Aligner overview
4 Exercise #1: BWA global alignment – Yeast ChIP-seq
- 4.1 Overview ChIP-seq alignment workflow with BWA
- 4.2 Introducing BWA
  - 4.2.1 BWA suite usage
- 4.3 Building the BWA sacCer3 index
  - 4.3.1 bwa index usage
  - 4.3.2 Get the alignment exercises files
  - 4.3.3 Prepare BWA reference directory for sacCer3
  - 4.3.4 Build BWA index for sacCer3
  - 4.3.5 BWA index files for sacCer3
- 4.4 Performing the bwa alignment
  - 4.4.1 Prepare to align yeast data
  - 4.4.2 Start an idev session
  - 4.4.3 bwa aln commands for yeast R1 and R2
  - 4.4.4 Pairing of BWA R1 and R2 aligned reads
- 4.5 Using cut to isolate fields
  - 4.5.1 Cut syntax for a single field
  - 4.5.2 Cut syntax for multiple fields
  - 4.5.3 Grep pattern that doesn't match header
  - 4.5.4 Get contig name info with cut
  - 4.5.5 Filter contig name of * (unaligned)
  - 4.5.6 Count aligned SAM records
  - 4.5.7 Calculate alignment rate
5 Exercise #2: Basic SAMtools Utilities
- - 5.1.1 SAMtools suite usage
- 5.2 samtools view
  - 5.2.1 samtools view usage
  - 5.2.2 Get the alignment exercises files
  - 5.2.3 Convert SAM to binary BAM
  - 5.2.4 View BAM records
- 5.3 samtools sort
  - 5.3.1 samtools sort usage
  - 5.3.2 Sort a BAM file
- 5.4 samtools index
  - 5.4.1 samtools index usage
  - 5.4.2 Index a sorted bam
- 5.5 samtools flagstat
  - 5.5.1 Run samtools flagstat using tee
  - 5.5.2 samtools flagstat output
- 5.6 samtools idxstats
  - 5.6.1 Use samtools idxstats to summarize mapped reads by contig
  - 5.6.2 samtools idxstats output
6 Exercise #3: BWA PE alignment with BioITeam script
- - 6.1.1 Run multiple alignments using the TACC batch system
  - 6.1.2 Commands to run multiple alignment scripts
- 6.2 Output files
- 6.3 Verifying alignment success
  - 6.3.1 Checking the alignment log file
  - 6.3.2 Count the number of successful alignments
  - 6.3.3 Check for failed alignment tasks
- 6.4 Checking alignment statistics
  - 6.4.1 <prefix>.samstats.txt output
  - 6.4.2 <prefix>.iszinfo.txt output
  - 6.4.3 Review multiple alignment rates
- 6.5 TACC batch system considerations

Overview

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 login5.ls5.tacc.utexas.edu. This should be second nature by now

Then stage the sample datasets and references we will use.

Get the alignment exercises files

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

These are descriptions of the FASTQ files we copied:

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 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.)

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

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:

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.:

/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.:

>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.