1.2.1 Assembling paired-end reads.¶

join_paired_ends.py -f C01D01F_sub.fastq -r C01D01R_sub.fastq -o C01D01

The previous command works by calling the python script join_paired_ends.py and having it act on the forward and reverse reads for sample C01D01.

1.2.2 Sanity check and file inspection.¶

There are some questions you may be having: What files does this return? How many reads were successfully merged?

First let’s take a look at what the outputs of join_paired_ends.py.

cd C01D01
ls -lah

You’ll see three files here: fastqjoin.join.fastq, fastqjoin.un1.fastq, and fastqjoin.un2.fastq. fastqjoin.join.fastq contains all of the reads that were successfully merged, while fastqjoin.un1.fastq and fastqjoin.un2.fastq contain those reads which failed to merge form the forward and reverse files. Let’s have a look at our merged sequences.

head fastqjoin.join.fastq

We see that the file is still in fastq format (ie there are 4 lines for each sequence beginning with a header line). We could go through an individually count every sequence that is there, or we could make use of grep and wc to automatically count.

grep  '@HWI' fastqjoin.join.fastq | wc -l

The above command actually works in two steps. Everything before the | is searching for the string @HWI in our merged reads file fastq.join.fastq. We search for @HWI because every sequence header will start with this. The command after the | counts the number of lines from the output in the first half of the command. Therefore the command counts every single header of the sequences in the merged file which equals the total number of merged sequences. We see that the program managed to merge 7670 sequences.

We can double check our sanity by using a positive control. Let’s investigate how many sequences there were in the original unmerged files.

grep '@HWI' ../C01D01F_sub.fastq | wc -l

This returns 10,000. So our the program managed to merged about 75% of the reads. Just in case you aren’t convinced as to how this works lets go ahead and redirect the grep command to a new file instead of to the word count command.

grep '@HWI' fastqjoin.join.fastq > C01D01_MergedHeaders.txt

The > redirects the output of grep to a file instead of having every single result fly by on the screen. If you head the C01D01_MergedHeaders.txt file you’ll see the following.

Then feel free to use wc -l on this file to check the number of lines. The -l command in wc signifies that we want to count every line instead of every word.

wc -l C01D01_MergedHeaders.txt

The output is once again 7670. Hope you’re convinced now!

Next, we are going to automate the previous process so we don’t have to individually type in 54 lines of code. This means we need to clean up the directory for the sample we just did though. Move back a directory in the EDAMAME_16S/Fastq directory and execute the following command.

rm -r C01D01

1.2.3 Automate paired-end merging with a shell script.¶

We would have to execute an iteration of the join_paired_ends.py command for every pair of reads that need to be assembled. This could take a long time. So, we’ll use a shell script to automate the task.

To download the script and list onto the AMI, first navigate to the “Fastq” directory, use curl to get the files, and make a new Merged_Reads directory to put the merged reads into.

curl -O https://raw.githubusercontent.com/edamame-course/Amplicon_Analysis/master/resources/Merged_Reads_Script.sh

curl -O https://raw.githubusercontent.com/edamame-course/Amplicon_Analysis/master/resources/list.txt
mkdir Merged_Reads

Change permissions on the script to make it executable:

chmod 755 Merged_Reads_Script.sh

Execute the script from the Fastq Directory. Bonus: A great opportunity to try tmux!

./Merged_Reads_Script.sh

Bonus: A little about this shell script For those of you interested in how this script works I recommend you take a look at it either through your terminal with less or online here. The script uses a “for loop” to repeat a set of commands for every single pair of files we have. For each sample the script first performs the merging with join_paired_ends.py and saves the output to a directory with the name of the sample. Secondly, because of some weird naming conventions the join_paired_ends.py function follows, the script renames and moves the resulting merged read files to a single directory each with their own name. Finally, the script removes the original directory that join_paired_ends.py as a cleanup step. This gets rid of 54 unnecessary directories making saving our instance space and our eyes the strain of having to look at them all everytime we want to ls in the directory.

1.2.4 Sanity check #2.¶

How many files were we expecting from the assembly? There were 54 pairs to be assembled, and we are generating one assembled fastq for each. Thus, the Merged_reads directory should contain 54 files. Navigate up one directory, and then use the wc (word count) command to check the number of files in the directory.

ls Merged_Reads | wc -l

The terminal should return the number “54”. Let’s move this whole directory up one level so that we can access more easily with QIIME:

mv Merged_Reads ..

Congratulations! You’ve assembled paired-end reads!

1.5.1 Preamble¶

Picking OTUs is sometimes called “clustering,” as sequences with some threshold of identity are “clustered” together to into an OTU.

Important decision: Should I use a de-novo method of picking OTUs or a reference-based method, or some combination? (Or not at all?). The answer to this will depend, in part, on what is known about the community a priori. For instance, a human or mouse gut bacterial community will have lots of representatives in well-curated 16S databases, simply because these communities are relatively well-studied. Therefore, a reference-based method may be preferred. The limitation is that any taxa that are unknown or previously unidentified will be omitted from the community. As another example, a community from a lesser-known environment, like Mars or a cave, or a community from a relatively less-explored environment would have fewer representative taxa in even the best databases. Therefore, one would miss a lot of taxa if using a reference-based method. The third option is to use a reference database but to set aside any sequences that do not have good matches to that database, and then to cluster these de novo.

For this tutorial, we are going to use an OTU-picking approach that uses a reference to identify as many OTUs as possible, but also includes any “new” sequences that do not hit the database. It is called “open reference” OTU picking, and you can read more about it in this paper by Rideout et al. Discussion of the workflow by the QIIME developers is here.

We use the QIIME workflow command: pick_open_reference_otus.py for this step. Documentation is here. The default QIIME 1.9.1 method for OTU picking is uclust, which we use here for simplicity’s sake. However, we encourage you to explore different OTU clustering algorithms to understand how they perform. They are not created equal.

1.5.2 OTU picking¶

This next step will take about 45 minutes. Note: Use tmux!

The link explains tmux in more detail, but for now, execute the following:

tmux new -s OTUP

Before continuing, make sure you are in the “EDAMAME_16S” directory.

pick_open_reference_otus.py -i combined_seqs.fna -o uclust_openref/ -f

In the above script:

• We tell QIIME to look for the input file -i, “combined_seqs.fna”.
• We specify that output files should go in a new folder, uclust_openref/
• We tell the program to overwrite already-existing files in the folder if we are running this program more than once (-f)

Other default parameters of interest:

• Singletons are removed from the OTU table (default flag –min_otu_size)
• Alignment is performed with PyNAST
• Taxonomy is assigned with uclust
• Default workflow values can be changed using a parameter file
• We do not perform prefiltering, as per the recommendations of Rideout et al.

Finally, we need to exit tmux:

ctrl+b d

Alignment output¶

Navigate into the pynast_aligned_seq directory directory. There are four files waiting there: one file of sequences that failed to align, one of sequences that did align, one of “pre-alignment” files, and a log file. Inspect each.

If you want to build a tree with some other out-of-QIIME software, this is where you would find the rep_set alignment. The log file provides a rep-sequence by rep-sequence report. If you needed align sequences as an independent step, you would use align_seqs.py; documentation here.

How many failed alignments were there?

count_seqs.py -i rep_set_failures.fasta

We see that there were 690 rep. sequences that failed to align, and approximately 22,206 that did. (Also, notice what short-read alignments generally look like...not amazing).

Taxonomy files¶

Move up a directory and then cd into the uclust_assigned_taxonomy directory.

more rep_set_tax_assignments.txt

In the “taxonomy” directory, you will find a log file and the specific taxonomic assignments given to each representative sequence, linked to the OTU ID of that representative sequence

#####Other Very Important Files in uclust_openref/ ######OTU map

nano final_otu_map.txt

Explore this file. It links the exact sequences in each sample to its OTU ID. You should see an OTU ID (starting with the number of the first OTU that was picked) in the the left most column. After that number, there is a list of Sequence IDs that have been clustered into that OTU ID. The first part of the sequence ID is the SampleID from which it came, and the second part is the sequence number within that sample.

You will notice that some files have “mc2” appended to them. “mc2” designates that the minimum count of sequences for any particular OTU was 2. In other words, that file has singleton OTUs already filtered out.

Sanity check: How can you compare the OTUs in the full dataset versus the singletons-omitted dataset?

more rep_set.fna