Sequence Alignment

We are now ready to align our sequences to the reference. We will use SILVA, the high-quality ribosomal RNA databases, as our reference. This is an important step to improve the clustering of your operational taxonomic unit (OTUs). An OTU is an operational definition used to classify groups of closely related individuals. Here, we need to align our sequences to a reference alignment using align.seqs. Our sequences aren’t just random strings of nucleotides; they are part of a known gene (such as 16S). After alignment is complete, we will get the alignment statistics by running Summary.seqs, and perform additional data cleaning using Filter.seqs.

Step 1: Sequence alignment

Have a look at the alignment output:

After the alignment is done, have a look at the summary output (log file):

The Start and End columns tell us that the majority of reads aligned between positions 1968 and 11550, which is what we expect to find given the reference file we used. However, some reads align to very different positions, which could indicate insertions or deletions at the terminal ends of the alignments or other complicating factors.

Also notice the Polymer column in the output table. This indicates the average homopolymer length. Since we know that our reference database does not contain any homopolymer stretches longer than 8 bases, any reads containing such long stretches are likely the result of PCR errors and we would be wise to remove them.

Next we will clean our data further by removing poorly aligned sequences and any sequences with long homopolymer stretches.

Step 2: Post-alignment Data Cleaning - to ensure all our reads overlap our region of interest:

Remove poorly aligned sequences

Click the View data (eye icon) next to the Filter.seqs ... : filtered fasta result in the History panel. Your resulting alignment (filtered fasta output) should look something like this (don't click Show all or Save):

These are all our representative reads again, now with additional alignment information.

Click on the View data icon next to Filter.seqs ... : logfile, and when the results are loaded in the center panel, scroll to the very bottom of the page. You should see the following additional information:

From this log file we see that while our initial alignment was 13425 positions wide, after filtering the overhangs (trump parameter) and removing positions that had a gap in every aligned read (vertical parameter), we have trimmed our alignment down to a length of 362.

Because any filtering step we perform might lead to sequences no longer being unique, we deduplicate our data by re-running the Unique.seqs tool:

Pre-clustering

The next step in cleaning our data, is to merge near-identical sequences together. Sequences that only differ by around 1 in every 100 bases are likely to represent sequencing errors, not true biological variation. Because our contigs are ~250 bp long, we will set the threshold to 2 mismatches.

Step 3: Chimera Removal

We have now thoroughly cleaned our data and removed as much sequencing error as we can. Next, we will look at a class of sequencing artefacts known as chimeras.

Chimeras occur when the polymerase falls off one sequence and then reattaches to another during a single round of PCR. This results in a sequence that looks unique but is, in fact, just half of one and half of another e.g. not real. Needless to say, we do not want such sequencing artefacts confounding our results. We’ll do this chimera removal using the VSEARCH algorithm Rognes et al. 2016 that is called within mothur, using the Chimera.vsearch tool.

(slide credit: http://slideplayer.com/slide/4559004/ )

Chimera.vsearch will split the data by sample and check for chimeras. The recommended way of doing this is to use the abundant sequences as our reference.