Overview

The fastQC tool was presented in the second tutorial on the first day of the class as the go to tool for quality control analysis of fastq files, but there is an underlying issue that checking each fastq file is quite daunting and evaluating each file individually can introduce its own set of artifacts or biases. The MultiQC tool represents a tool which works directly on fastQC reports to quickly generate summary reports to both identify samples that are different among a group and to make global decisions about how to treat a set of files.

Learning Objectives

In this tutorial, we will:

1. work with some simple bash scripting from the command line (for loops) to generate multiple fastqc reports simultaneously and look at 272 plasmid samples.
2. work with MultiQC to make decisions about read preprocessing.
3. identify outlier files that are clearly different from the group as a whole and determine how to deal with these files.

Installing multiqc

Hopefully by now you have had enough experience installing packages via conda that it is second nature to think of going to https://anaconda.org/ searching for multiqc, finding  https://anaconda.org/bioconda/multiqc, and using the information there to install the package.

conda activate GVA2021
conda install -c bioconda multiqc

Found conflicts! Looking for incompatible packages.
This can take several minutes.  Press CTRL-C to abort.
failed                                                                                                                                                                                                                

UnsatisfiableError: The following specifications were found to be incompatible with each other:

Output in format: Requested package -> Available versions

As we have seen several times now, when we have difficulties getting programs installed together, and need them to interact with each other, the easiest solution is often to create a new environment, specifying all desired packages at the same time.

Get some data and verify access to fastqc

Copy the plasmid sequencing files found in the BioITeam directory gva_course/plasmid_qc/ to a new directory named GVA_multiqc. There are 2 main ways to do this particularly since there are so many files (544 total). 

cp -r $BI/gva_course/plasmid_qc/ $SCRATCH/GVA_multiqc

mkdir $SCRATCH/GVA_multiqc 
cp $BI/gva_course/plasmid_qc/* $SCRATCH/GVA_multiqc

fastqc --version

Generating FastQC commands

Throughout the first part of the course we focused on working with a single sample and thus were able to type commands 1 at a time. We further only had a few input files that we were dealing with in an individual tutorial thus tab completion and ls are very useful. Here we are dealing with 544 files which is more than the total number of files we dealt with in all the required tutorials combined, and nobody wants to type out 544 commands 1 at a time. Therefore, we are going to construct a single commands file with 544 lines that we can use to launch all commands without having to know the name of any single file.  To do so we will use the bash 'for' command.

cd $SCRATCH/GVA_multiqc
ls

For loops on the command line have 3 parts:

1. A list of something to deal with 1 at a time. Followed by a ';'
   1. for f in *.gz; in the following example
2. Something to do with each item in the list. this must start with the word 'do'
   1. do echo "fastqc -o fastqc_output $f "; in the following example
3. The word "done" so bash knows to stop looking for more commands.
   1. done in the following example, but we add a final redirect (>) so rather than printing to the screen the output goes to a file (fastqc_commands in this case)

for f in *.gz; do echo "fastqc -o fastqc_output $f"; done  > fastqc_commands

Use the head and wc -l to see what the output is and how much there is of it respectively.

Next we need to make the output directory for all the fastqc reports to go into and send the fastqc_commands file to the queue to execute. Like our breseq tutorial, this involves the use of a .slurm file.

mkdir fastqc_output
cp /corral-repl/utexas/BioITeam/gva_course/GVA2021.launcher.slurm multiqc.slurm
nano multiqc.slurm

Again while in nano you will edit the same lines you edited in the in the breseq tutorial. Note that most of these lines have additional text to the right of the line. This commented text is present to help remind you what goes on each line, leaving it alone will not hurt anything, removing it may make it more difficult for you to remember what the purpose of the line is

The changes to lines 22 and 23 are optional but will give you an idea of what types of email you could expect from TACC if you choose to use these options. Just be sure to pay attention to these 2 lines starting with a single # symbol after editing them.

Line 27 does not change this time as we are still working with the GVA2021 environment. It is included in the table to help remind you of its importance, and needing to list the correct environment on this line. 

Again use ctl-o and ctl-x to save the file and exit.

mkdir fastqc_output

sbatch fastqc.slurm

Run MultiQC

Hopefully by now the job we submitted with all the fastqc commands has finished. Use the showq -u command to check. If the showq -u command still shows your job in the 'waiting' section it has not started, once it does start i would expect it to finish in ~4-5 minutes. It is possible that an error occurs and does not finish so we need to check for output from the files.

If you ls fastqc_output directory you are hit in the face with more files and directories than you have seen in any directory during this class. You immediately can notice that there is a directory and a compressed version of each of those directories, but in order to know things worked correctly, we need to make sure that we have 2 files for each of our 544 samples. The easiest way to do that in my opinion is to pipe that output to the wc -l command to count the total number of lines.

ls fastqc_output; ls fastqc_output | wc -l

Assuming you see both directories and their associated compressed files,  and a total count of 1088 run the 'multiqc -h' command to look through the options and see if you can figure out how to guild the command. 

cd fastqc_output
multiqc .

In this case (much as is the case with FastQC) while there are are a reasonable number of options that can be used, none are truly needed for evaluating fastq files. The only requirement is that you specify where the FastQC output that you want to generate a single report for is located. In the example above you changed into the directory containing those results and then specified multiqc should look in you current directory for the files. It would have been comparable to stay in the existing directory, and instead use a command of "multiqc fastqc_output". 

Evaluate MultiQC report

As the multiqc_report.html file is a html file, you will need to transfer it back to your laptop to view it. Hopefully, by now you have learned how to do this without needing the scp tutorial open to help you. If not, consider getting my attention on zoom so i can try to help clear up any confusion you may be having.

Once you have the report back on your local computer, open it and begin looking at it. Unlike the FastQC report, the multiqc report comes with detailed help information for each section you can access with the different ?help icons on the right as well as a video at the top of the page. If anything is not clear, and you'd like help clearing it up, let me know.

Optional Exercise

1. Using information in the MultiQC report, modify the bash loop used to create the fastqc_commands file above to create a cutadapt_commands file that could modify all 544 files at once.
2. Move over to the trimmomatic tutorial and come back to trim all adapter sequences from all files and rerun fastqc/multiqc to see what a difference trimming makes on overall quality.

Return to the Genome Variant Analysis Course 2020 Home Page