2020 Alignment (archive)

Open

After raw sequence files are generated (in FASTQ format), quality-checked, and pre-processed in some way, the next step in most 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, and with sufficient (but not absolute - and this tradeoff is an important consideration when constructing alignment pipelines) accuracy 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. These programs vary widely in their design, inputs, outputs, and applications. 

In this section, we will primarily focus on two of the most versatile mappers: BWA and Bowtie2, the latter being part of the Tuxedo suite (e.g. transcriptome-aware Tophat2) which also includes tools for manipulating NGS data after alignment.

This should be second nature by now

You have already worked with a paired-end yeast ChIP-seq dataset, which we will continue to use here.  For the sake of uniformity, however, we will set up a new directory in your scratch area called 'alignment' and fill it with our sequencing data.  To set up your scratch area properly and move into it, execute something like:

Now you have created the alignment directory, moved into it, and created a subdirectory for our raw fastq files.  We will be using four data sets that consist of five files (since the paired-end data set has two separate files for each of the R1 and R2 reads).  To copy them over, execute something like:

We first moved ourselves into the fastq directory, then copied over all files that end in "fastq.gz" in the corral directory specified in the second line.  These are descriptions of the five files we copied:

Before we get to alignment, we need a genome to align to.  We will use four different references here: 

• the yeast genome (sacCer3)

• the human genome (hg19)

• a Vibrio cholerae genome (0395; our name: vibCho)

• the microRNA database mirbase (v20), human subset

NOTE: For the sake of simplicity, these are not necessarily the most recent versions of these references - for example, hg19 is the second most recent human genome, with the most recent called hg38.  Similarly, the most recent mirbase annotation is v21.

These are the four reference genomes we will be using today, with some information about them (and here is information about many more genomes):

Searching genomes is computationally hard work and takes a long time if done on un-indexed, linear genomic sequence.  So aligners require that references first be indexed to accelerate later retrieval.  The aligners we are using each require a different index, but use the same method (the Burrows-Wheeler Transform) to get the job done. This involves taking a FASTA file as input, with each chromosome (or contig) as a separate FASTA entry, and producing an aligner-specific set of files as output. Those output index files are then used by the aligner when performing the sequence alignment, and subsequent alignments are reported using coordinates referencing the original FASTA reference files.

hg19 is way too big for us to index here, so we're not going to do it (especially not all at the same time!). Instead, we will "point" to an existing set of hg19 index files, which are all located at:

We can index the references for the yeast genome, the human miRNAs, and the V. cholerae genome, because they are all tiny compared to the human genome. We will grab the FASTA files for yeast and human miRNAs two references and build each index right before we use them. We will also grab the special file that contains the V. cholerae genome sequence and annotations (a .gbk file), and generate the reference FASTA and some other interesting information when we get to that exercise. These references are currently at the following locations: