Quality Control¶
Shell Variables¶
Assign the variables in this notebook.
[1]:
set -u
# We used these in the last notebook
DATA_BASE="/data/hts_2019_data"
RAW_FASTQS="$DATA_BASE/hts2019_pilot_rawdata"
# The following are new for this notebook
CUROUT=$HOME/work/scratch/bioinf_intro
QC=$CUROUT/qc_output
Making New Directories¶
The directories do carry over between notebooks, they are a function of the server, so we only need to make the directories that are new in this notebook
[2]:
mkdir -p $QC
Now let’s check to be sure that worked. We will run ls and check that these directories now exist in the $OUTPUT directory.
[3]:
ls $CUROUT
count_out igv qc_output stuff_for_igv_shorter_intron.tgz trimmed_fastqs
genome myinfo star_out stuff_for_igv.tgz
Fastqc¶
Now that we have had a quick look at our fastq file, let’s run some basic analysis using a program called fastqc. This will take less than a minute to run. Remember that while “long” jobs are running it will say “In [*]:” in the left margin, once it is done, a number will replace the asterisk.
[4]:
fastqc -h
FastQC - A high throughput sequence QC analysis tool
SYNOPSIS
fastqc seqfile1 seqfile2 .. seqfileN
fastqc [-o output dir] [--(no)extract] [-f fastq|bam|sam]
[-c contaminant file] seqfile1 .. seqfileN
DESCRIPTION
FastQC reads a set of sequence files and produces from each one a quality
control report consisting of a number of different modules, each one of
which will help to identify a different potential type of problem in your
data.
If no files to process are specified on the command line then the program
will start as an interactive graphical application. If files are provided
on the command line then the program will run with no user interaction
required. In this mode it is suitable for inclusion into a standardised
analysis pipeline.
The options for the program as as follows:
-h --help Print this help file and exit
-v --version Print the version of the program and exit
-o --outdir Create all output files in the specified output directory.
Please note that this directory must exist as the program
will not create it. If this option is not set then the
output file for each sequence file is created in the same
directory as the sequence file which was processed.
--casava Files come from raw casava output. Files in the same sample
group (differing only by the group number) will be analysed
as a set rather than individually. Sequences with the filter
flag set in the header will be excluded from the analysis.
Files must have the same names given to them by casava
(including being gzipped and ending with .gz) otherwise they
won't be grouped together correctly.
--nofilter If running with --casava then don't remove read flagged by
casava as poor quality when performing the QC analysis.
--extract If set then the zipped output file will be uncompressed in
the same directory after it has been created. By default
this option will be set if fastqc is run in non-interactive
mode.
-j --java Provides the full path to the java binary you want to use to
launch fastqc. If not supplied then java is assumed to be in
your path.
--noextract Do not uncompress the output file after creating it. You
should set this option if you do not wish to uncompress
the output when running in non-interactive mode.
--nogroup Disable grouping of bases for reads >50bp. All reports will
show data for every base in the read. WARNING: Using this
option will cause fastqc to crash and burn if you use it on
really long reads, and your plots may end up a ridiculous size.
You have been warned!
-f --format Bypasses the normal sequence file format detection and
forces the program to use the specified format. Valid
formats are bam,sam,bam_mapped,sam_mapped and fastq
-t --threads Specifies the number of files which can be processed
simultaneously. Each thread will be allocated 250MB of
memory so you shouldn't run more threads than your
available memory will cope with, and not more than
6 threads on a 32 bit machine
-c Specifies a non-default file which contains the list of
--contaminants contaminants to screen overrepresented sequences against.
The file must contain sets of named contaminants in the
form name[tab]sequence. Lines prefixed with a hash will
be ignored.
-a Specifies a non-default file which contains the list of
--adapters adapter sequences which will be explicity searched against
the library. The file must contain sets of named adapters
in the form name[tab]sequence. Lines prefixed with a hash
will be ignored.
-l Specifies a non-default file which contains a set of criteria
--limits which will be used to determine the warn/error limits for the
various modules. This file can also be used to selectively
remove some modules from the output all together. The format
needs to mirror the default limits.txt file found in the
Configuration folder.
-k --kmers Specifies the length of Kmer to look for in the Kmer content
module. Specified Kmer length must be between 2 and 10. Default
length is 7 if not specified.
-q --quiet Supress all progress messages on stdout and only report errors.
-d --dir Selects a directory to be used for temporary files written when
generating report images. Defaults to system temp directory if
not specified.
BUGS
Any bugs in fastqc should be reported either to simon.andrews@babraham.ac.uk
or in www.bioinformatics.babraham.ac.uk/bugzilla/
[5]:
# fastqc --threads 2 --quiet $DEMUX/r1.8A_pilot.fq.gz -o $QC
fastqc --extract $RAW_FASTQS/21_2019_P_M1_S21_L002_R1_001.fastq.gz -o $QC
Started analysis of 21_2019_P_M1_S21_L002_R1_001.fastq.gz
Approx 5% complete for 21_2019_P_M1_S21_L002_R1_001.fastq.gz
Approx 10% complete for 21_2019_P_M1_S21_L002_R1_001.fastq.gz
Approx 15% complete for 21_2019_P_M1_S21_L002_R1_001.fastq.gz
Approx 20% complete for 21_2019_P_M1_S21_L002_R1_001.fastq.gz
Approx 25% complete for 21_2019_P_M1_S21_L002_R1_001.fastq.gz
Approx 30% complete for 21_2019_P_M1_S21_L002_R1_001.fastq.gz
Approx 35% complete for 21_2019_P_M1_S21_L002_R1_001.fastq.gz
Approx 40% complete for 21_2019_P_M1_S21_L002_R1_001.fastq.gz
Approx 45% complete for 21_2019_P_M1_S21_L002_R1_001.fastq.gz
Approx 50% complete for 21_2019_P_M1_S21_L002_R1_001.fastq.gz
Approx 55% complete for 21_2019_P_M1_S21_L002_R1_001.fastq.gz
Approx 60% complete for 21_2019_P_M1_S21_L002_R1_001.fastq.gz
Approx 65% complete for 21_2019_P_M1_S21_L002_R1_001.fastq.gz
Approx 70% complete for 21_2019_P_M1_S21_L002_R1_001.fastq.gz
Approx 75% complete for 21_2019_P_M1_S21_L002_R1_001.fastq.gz
Approx 80% complete for 21_2019_P_M1_S21_L002_R1_001.fastq.gz
Approx 85% complete for 21_2019_P_M1_S21_L002_R1_001.fastq.gz
Approx 90% complete for 21_2019_P_M1_S21_L002_R1_001.fastq.gz
Approx 95% complete for 21_2019_P_M1_S21_L002_R1_001.fastq.gz
Analysis complete for 21_2019_P_M1_S21_L002_R1_001.fastq.gz
Once fastqc is done running we can view the results by finding the output in the Jupyter browser, it should be in:
[6]:
echo $QC
/home/jovyan/work/scratch/bioinf_intro/qc_output
MultiQC¶
FastQC is a useful tool, but it has one problem: it generates one report for each FASTQ file. When you have more than a handful of FASTQs (as most projects will), it is tedious to look at each one, and there is no simple way to compare them.
MultiQC is a solution to this problem. It mines the results from FastQC (and other HTS analysis tools) and generates reports that combine and summarize results for all the FASTQs analyzed.
Run FastQC on Multiple FASTQs¶
Let’s pick a selection of FASTQs to look at - we will run FastQC on the Lane 1 results for samples 10 through 19 (because it is easy to specify these files with a simple command).
[7]:
ls $RAW_FASTQS/1?_2019_*_L001_R1_001.fastq.gz
/data/hts_2019_data/hts2019_pilot_rawdata/10_2019_P_M1_S10_L001_R1_001.fastq.gz
/data/hts_2019_data/hts2019_pilot_rawdata/11_2019_P_M1_S11_L001_R1_001.fastq.gz
/data/hts_2019_data/hts2019_pilot_rawdata/12_2019_P_M1_S12_L001_R1_001.fastq.gz
/data/hts_2019_data/hts2019_pilot_rawdata/13_2019_P_M1_S13_L001_R1_001.fastq.gz
/data/hts_2019_data/hts2019_pilot_rawdata/14_2019_P_M1_S14_L001_R1_001.fastq.gz
/data/hts_2019_data/hts2019_pilot_rawdata/15_2019_P_M1_S15_L001_R1_001.fastq.gz
/data/hts_2019_data/hts2019_pilot_rawdata/16_2019_P_M1_S16_L001_R1_001.fastq.gz
/data/hts_2019_data/hts2019_pilot_rawdata/17_2019_P_M1_S17_L001_R1_001.fastq.gz
/data/hts_2019_data/hts2019_pilot_rawdata/18_2019_P_M1_S18_L001_R1_001.fastq.gz
/data/hts_2019_data/hts2019_pilot_rawdata/19_2019_P_M1_S19_L001_R1_001.fastq.gz
--threads 2 tells fastqc to run using 2 cores (of the 16 that we all share). In general you should not do this because it is potentially inconsiderate of other users. I am using it so we don’t have to wait too long for this to run during class.
[8]:
fastqc --quiet --threads 2 --extract $RAW_FASTQS/1?_2019_*_L001_R1_001.fastq.gz -o $QC
Run MultiQC¶
[9]:
multiqc -h
/home/jovyan/.local/lib/python3.5/site-packages/multiqc/utils/config.py:45: YAMLLoadWarning: calling yaml.load() without Loader=... is deprecated, as the default Loader is unsafe. Please read https://msg.pyyaml.org/load for full details.
configs = yaml.load(f)
/home/jovyan/.local/lib/python3.5/site-packages/multiqc/utils/config.py:51: YAMLLoadWarning: calling yaml.load() without Loader=... is deprecated, as the default Loader is unsafe. Please read https://msg.pyyaml.org/load for full details.
sp = yaml.load(f)
Usage: multiqc [OPTIONS] <analysis directory>
MultiQC aggregates results from bioinformatics analyses across many
samples into a single report.
It searches a given directory for analysis logs and compiles a HTML
report. It's a general use tool, perfect for summarising the output from
numerous bioinformatics tools.
To run, supply with one or more directory to scan for analysis results. To
run here, use 'multiqc .'
See http://multiqc.info for more details.
Author: Phil Ewels (http://phil.ewels.co.uk)
Options:
-f, --force Overwrite any existing reports
-d, --dirs Prepend directory to sample names
-dd, --dirs-depth INTEGER Prepend [INT] directories to sample names.
Negative number to take from start of path.
-s, --fullnames Do not clean the sample names (leave as full
file name)
-i, --title TEXT Report title. Printed as page header, used
for filename if not otherwise specified.
-b, --comment TEXT Custom comment, will be printed at the top
of the report.
-n, --filename TEXT Report filename. Use 'stdout' to print to
standard out.
-o, --outdir TEXT Create report in the specified output
directory.
-t, --template [default_dev|geo|simple|sections|default]
Report template to use.
--tag TEXT Use only modules which tagged with this
keyword, eg. RNA
--view-tags, --view_tags View the available tags and which modules
they load
-x, --ignore TEXT Ignore analysis files (glob expression)
--ignore-samples TEXT Ignore sample names (glob expression)
--ignore-symlinks Ignore symlinked directories and files
--sample-names PATH File containing alternative sample names
-l, --file-list Supply a file containing a list of file
paths to be searched, one per row
-e, --exclude [module name] Do not use this module. Can specify multiple
times.
-m, --module [module name] Use only this module. Can specify multiple
times.
--data-dir Force the parsed data directory to be
created.
--no-data-dir Prevent the parsed data directory from being
created.
-k, --data-format [json|tsv|yaml]
Output parsed data in a different format.
Default: tsv
-z, --zip-data-dir Compress the data directory.
-p, --export Export plots as static images in addition to
the report
-fp, --flat Use only flat plots (static images)
-ip, --interactive Use only interactive plots (HighCharts
Javascript)
--lint Use strict linting (validation) to help code
development
--pdf Creates PDF report with 'simple' template.
Requires Pandoc to be installed.
--no-megaqc-upload Don't upload generated report to MegaQC,
even if MegaQC options are found
-c, --config PATH Specific config file to load, after those in
MultiQC dir / home dir / working dir.
--cl-config, --cl_config TEXT Specify MultiQC config YAML on the command
line
-v, --verbose Increase output verbosity.
-q, --quiet Only show log warnings
--version Show the version and exit.
-h, --help Show this message and exit.
[10]:
multiqc $QC --outdir $QC
/home/jovyan/.local/lib/python3.5/site-packages/multiqc/utils/config.py:45: YAMLLoadWarning: calling yaml.load() without Loader=... is deprecated, as the default Loader is unsafe. Please read https://msg.pyyaml.org/load for full details.
configs = yaml.load(f)
/home/jovyan/.local/lib/python3.5/site-packages/multiqc/utils/config.py:51: YAMLLoadWarning: calling yaml.load() without Loader=... is deprecated, as the default Loader is unsafe. Please read https://msg.pyyaml.org/load for full details.
sp = yaml.load(f)
[INFO ] multiqc : This is MultiQC v1.7
[INFO ] multiqc : Template : default
[INFO ] multiqc : Searching '/home/jovyan/work/scratch/bioinf_intro/qc_output'
Searching 221 files.. [####################################] 100%
[INFO ] fastqc : Found 11 reports
[INFO ] multiqc : Compressing plot data
[WARNING] multiqc : Previous MultiQC output found! Adjusting filenames..
[WARNING] multiqc : Use -f or --force to overwrite existing reports instead
[INFO ] multiqc : Report : ../../scratch/bioinf_intro/qc_output/multiqc_report_1.html
[INFO ] multiqc : Data : ../../scratch/bioinf_intro/qc_output/multiqc_data_1
[INFO ] multiqc : MultiQC complete
Once multiqc is done running we can view the results by finding the output in the Jupyter browser, it should be in a file named multiqc_report.html in :
[11]:
echo $QC
/home/jovyan/work/scratch/bioinf_intro/qc_output