Before you start the alignment and analysis processes, it us useful to perform some initial quality checks on your raw data. You may also need to pre-process the sequences to trim them or remove adapters. Here we will assume you have paired-end data from one of GSAF's Illumina sequencers.
Use our summer school reservation (CoreNGS-Wed) when submitting batch jobs to get higher priority on the ls6 normal queue today:
|
The first order of business after receiving sequencing data should be to check your data quality. This often-overlooked step helps guide the manner in which you process the data, and can prevent many headaches.
FastQC is a tool that produces a quality analysis report on FASTQ files.
Useful links:
First and foremost, the FastQC "Summary" should generally be ignored. Its "grading scale" (green - good, yellow - warning, red - failed) incorporates assumptions for a particular kind of experiment, and is not applicable to most real-world data. Instead, look through the individual reports and evaluate them according to your experiment type.
The FastQC reports I find most useful, and why:
For many of its reports, FastQC analyzes only the first ~100,000 sequences in order to keep processing and memory requirements down. Consult the Online documentation for each FastQC report for full details. |
Make sure you're in an idev session. If you're in an idev session, the hostname command will display a name like c455-021.ls6.tacc.utexas.edu. But if you're on a login node the hostname will be something like login2.ls6.tacc.utexas.edu. If you're on a login node, start an idev session like this:
|
FastQC is available as part of BioContainers on ls6. To make it available:
# Load the main BioContainers module then load the fastqc module module load biocontainers # make take a while module load fastqc |
It has a number of options (see fastqc --help | more) but can be run very simply with just a FASTQ file as its argument.
|
# make sure you're in your $SCRATCH/core_ngs/fastq_prep directory cds cd core_ngs/fastq_prep fastqc small.fq |
Exercise: What did FastQC create?
ls -l shows two new items.
|
Let's unzip the .zip file and see what's in it.
unzip small_fastqc.zip |
What was created?
ls -l shows one new item, the small_fastqc directory (note the "d" in "drwxrwxr-x")
ls -l small_fastqc shows the directory contents:
|
You can't run a web browser directly from your "dumb terminal" command line environment. The FastQC results have to be placed where a web browser can access them. One way to do this is to copy the results back to your laptop, for example by using scp from your computer (read more at Copying files from TACC to your laptop).
For convenience, we put an example FastQC report at this URL:
https://web.corral.tacc.utexas.edu/BioinformaticsResource/CoreNGS/yeast_stuff/Sample_Yeast_L005_R1.cat_fastqc/fastqc_report.html
Exercise: Based on this FastQC output, should we trim this data?
The Per base sequence quality report does not look good. The data should probably be trimmed (to 40 or 50 bp) before alignment. |
Newer versions of FastQC have slightly different report formats. See this example:
https://web.corral.tacc.utexas.edu/BioinformaticsResource/CoreNGS/reports/wcaar_mqc_report.html
FastQC reports are all well and good, but what if you have dozens of samples? It quickly becomes tedious to have to look through all the separate FastQC reports, including separate R1 and R2 reports for paired end datasets.
The MultiQC tool helps address this issue. Once FastQC reports have been generated, it can scan them and create a consolidated report from all the individual reports.
Whats even cooler, is that MultiQC can also consolidate reports from other bioinformatics tools (e.g. bowtie2 aligner statistics, samtools statistics, cutadapt, Picard, and may more). And if your favorite tool is not known by MultiQC, you can configure custom reports fairly easily. For more information, see this recent Byte Club tutorial on Using MultiQC.
Here we're just going to create a MultiQC report for two paired-end ATAC-seq datasets – 4 FASTQ files total. First stage the data:
mkdir -p $SCRATCH/core_ngs/multiqc/fqc.atacseq cd $SCRATCH/core_ngs/multiqc/fqc.atacseq cp $CORENGS/multiqc/fqc.atacseq/*.zip . |
You should see these 4 files in your $SCRATCH/core_ngs/multiqc/fqc.atacseq directory:
50knuclei_S56_L007_R1_001_fastqc.zip 5knuclei_S77_L008_R1_001_fastqc.zip 50knuclei_S56_L007_R2_001_fastqc.zip 5knuclei_S77_L008_R2_001_fastqc.zip |
Now make the BioContainers MultiQC accessible in your environment.
Make sure you're in an idev session. If you're in an idev session, the hostname command will display a name like c455-020.ls6.tacc.utexas.edu. But if you're on a login node the hostname will be something like login1.ls6.tacc.utexas.edu. If you're on a login node, start an idev session like this:
|
# Load the main BioContainers module if you have not already module load biocontainers # may take a while # Load the multiqc module and ask for its usage information module load multiqc multiqc --help | more |
|
Even though multiqc has many options, it is quite easy to create a basic report by just pointing it to the directory where individual reports are located:
cd $SCRATCH/core_ngs/multiqc multiqc fqc.atacseq |
Exercise: How many reports did multiqc find?
Based on its execution output, it found 4 reports
|
Exercise: What was created by running multiqc?
One file was created (multiqc_report.html) and one directory (multiqc_data). |
You can see the resulting MultiQC report here: https://web.corral.tacc.utexas.edu/BioinformaticsResource/CoreNGS/reports/atacseq/multiqc_report.html.
And an example of a MultiQC report that includes both standard and custom plots is this is the Tag-Seq post-processing MultiQC report produced by the Bioinformatics Consulting Group: https://web.corral.tacc.utexas.edu/BioinformaticsResource/CoreNGS/reports/mqc_tagseq_trim_JA21030_SA21045_mouse.html
There are two main reasons you may want to trim your sequences:
There are a number of open source tools that can trim off 3' bases and produce a FASTQ file of the trimmed reads to use as input to the alignment program.
The FASTX Toolkit provides a set of command line tools for manipulating both FASTA and FASTQ files. The available modules are described on their website. They include a fast fastx_trimmer utility for trimming FASTQ sequences (and quality score strings) before alignment.
Make sure you're in an idev session. If you're in an idev session, the hostname command will display a name like c455-021.ls6.tacc.utexas.edu. But if you're on a login node the hostname will be something like login3.ls6.tacc.utexas.edu. If you're on a login node, start an idev session like this:
|
FASTX Toolkit is available as a BioContainers module.
module load biocontainers # takes a while module spider fastx module load fastxtools |
Here's an example of how to run fastx_trimmer to trim all input sequences down to 50 bases.
Where does fastx_trimmer read its input from? And where does it write its output? Ask the program for its usage.
# will fastx_trimmer give us usage information? fastx_trimmer --help # no, it wants you to use the -h option to ask for help: fastx_trimmer -h |
The usage: its help information
fastx_trimmer [-h] [-f N] [-l N] [-t N] [-m MINLEN] [-z] [-v] [-i INFILE] [-o OUTFILE] |
Because the [-i INFILE] [-o OUTFILE] options are shown in brackets [ ], reading from a file and writing to a file are optional. That means that by default the program reads its input data from standard input and writes trimmed sequences to standard output:
|
# make sure you're in your $SCRATCH/core_ngs/fastq_prep directory cd $SCRATCH/core_ngs/fastq_prep zcat Sample_Yeast_L005_R1.cat.fastq.gz | fastx_trimmer -l 50 -Q 33 > trim50_R1.fq |
Exercise: compressing fastx_trimmer output
How would you tell fastx_trimmer to compress (gzip) its output file?
Type fastx_trimmer -h (help) to see program documentation |
You could supply the -z option like this:
Or you could gzip the output yourself.
See the 3x+ difference in file sizes when the output is compressed with ls -lh trim* |
Exercise: other fastx toolkit programs
What other FASTQ manipulation programs are part of the FASTX Toolkit?
Type fastx_ then tab twice (completion) to see their names. |
The FASTX Toolkit also has programs that work on FASTA files. To see them, type fasta_ then tab twice (completion) to see their names.
Data from RNA-seq or other library prep methods that result in short fragments can cause problems with moderately long (50-100bp) reads, since the 3' end of sequences can be read into (or even through) to the 3' adapter at different read offsets . This 3' adapter contamination can cause the "real" insert sequence not to align because the adapter sequence does not correspond to the bases at the 3' end of the reference genome sequence.
Unlike general fixed-length trimming (e.g. trimming 100 bp sequences to 50 bp), specific adapter trimming removes differing numbers of 3' bases depending on where the adapter sequence is found.
You must tell any adapter trimming program what your R1 and R2 adapters look like. The GSAF website describes the flavors of Illumina adapter and barcode sequences in more detail: Illumina - all flavors (USE with Caution, this is outdated but can be useful for a basic understanding of the adapters, the GSAF primarily only uses UDI's for all projects). |
The cutadapt program, available in BioContainers, is an excellent tool for removing adapter contamination.
Make sure you're in an idev session. If you're in an idev session, the hostname command will display a name like c455-021.ls6.tacc.utexas.edu. But if you're on a login node the hostname will be something like login3.ls6.tacc.utexas.edu. If you're on a login node, start an idev session like this:
|
module load biocontainers module spider cutadapt module load cutadapt cutadapt --help |
A common application of cutadapt is to remove adapter contamination from RNA library sequence data. Here we'll show that for some small RNA libraries sequenced by GSAF, using their documented small RNA library adapters.
When you run cutadapt you give it the adapter sequence to trim, and the adapter sequence is different for R1 and R2 reads. Here's what the options look like (without running it on our files yet).
cutadapt -m 22 -O 4 -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC <fastq_file> |
cutadapt -m 22 -O 4 -a TGATCGTCGGACTGTAGAACTCTGAACGTGTAGA <fastq_file> |
Notes:
Figuring out which adapter sequence to use when can be tricky. Your sequencing provider can tell you what adapters they used to prep your libraries. For GSAF's adapter layout, please refer to Illumina - all flavors (USE with Caution, this is outdated but can be useful for a basic understanding of the adapters, the GSAF primarily only uses UDI's for all projects) (you may want to read all the "gory details" below later).
The top strand, 5' to 3', of a read sequence looks like this.
The -a argument to cutadapt is documented as the "sequence of adapter that was ligated to the 3' end". So we care about the <Read 2 primer> for R1 reads, and the <Read 1 primer> for R2 reads. The "contaminent" for adapter trimming will be the <Read 2 primer> for R1 reads. There is only one Read 2 primer:
The "contaminent" for adapter trimming will be the <Read 1 primer> for R2 reads. However, there are three different Read 1 primers, depending on library construction:
Since R2 reads are the reverse complement of R1 reads, the R2 adapter contaminent will be the RC of the Read 1 primer used. For ChIP-seq libraries where reads come from both DNA strands, the TruSeq Read 1 primer is always used.
For RNAseq libraries in this class, we use the small RNA sequencing primer as the Read 1 primer.
|
Exercise: other cutadapt options
The cutadapt program has many options. Let's explore a few.
How would you tell cutadapt to trim trailing N's?
Then, in the less pager, type /trim <enter> to look for the first occurrence of the string "trim", then n to look for subsequent occurrences. |
The relevant option is --trim-n |
How would you control the accuracy (error rate) of cutadapt's matching between the adapter sequences and the FASTQ sequences?
Use the less pager to search for terms like "error" or "accuracy". |
Then, in the less pager, type /error <enter> to look for the first occurrence of the string "error", then n to look for subsequent occurrences. |
The relevant option is -e <floating point error rate> or --error-rate=<floating point error rate>:
|
Suppose you are processing 100 bp reads with 30 bp adapters. By default, how many mismatches between the adapter and a sequence will be tolerated?
cutadapt's default error rate is 0.1 (10%) |
Up to three mismatches will be tolerated when the whole 30 bp adapter is found (10% of 30). If only 20 of the 30 adapter bases are found, up to two mismatches will be tolerated (10% of 20). |
How would you require a more stringent matching (i.e., allowing fewer mismatches)?
| Providing --error-rate=0.05 (or -e 0.05) as an option, for example, would specify a 5% error rate, or no more than 1 mismatching base in 20. |
Let's run cutadapt on some real human miRNA (micro-RNA) data.
First, stage the data we want to use. This data is from a small RNA library where the expected insert size is around 15-25 bp.
mkdir -p $SCRATCH/core_ngs/fastq_prep cd $SCRATCH/core_ngs/fastq_prep cp $CORENGS/human_stuff/Sample_H54_miRNA_L004_R1.cat.fastq.gz . cp $CORENGS/human_stuff/Sample_H54_miRNA_L005_R1.cat.fastq.gz . |
Exercise: How many reads are in these files? Is it single end or paired end data?
|
Looking at the FASTQ file names, we see this is two lanes of single-end reads (L004 and L005). The data from lane 4 has 2,001,337 reads, the data from lane 5 has 2,022,237 reads. |
Exercise: How long are the reads?
You could just Look at the size of the actual sequence on the 2nd line of any FASTQ entry and count the characters.... But you're experts now! So challenge yourself. Use a combination of tail and head to extract the 2nd line of the .gz file. Then use the wc program, but not with the -l option (check wc --help). |
|
These are 101-base reads. wc -c counts the "invisible" newline character, so subtract 1 from the character count it returns for a line. Here's a way to strip the trailing newline characters from the quality scores string before calling wc -c to count the characters. We use the echo -n option that tells echo not to include the trailing newline in its output. We gemerate that text using sub-shell evaluation (an alternative to backtick evaluation) of that zcat ... command:
|
Adapter trimming is a rather slow process, and these are large files. So to start with we're going to create a smaller FASTQ file to work with.
# Remember, FASTQ files have 4 lines per read zcat Sample_H54_miRNA_L004_R1.cat.fastq.gz | head -2000 > miRNA_test.fq |
Now execute cutadapt like this:
|
cutadapt -m 20 -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC miRNA_test.fq \ 2> miRNA_test.cuta.log \ | gzip > miRNA_test.cutadapt.fq.gz |
Notes:
You should see a miRNA_test.cuta.log log file when the command completes. How many lines does it have?
|
Take a look at the first 15 lines.
|
It will look something like this:
This is cutadapt 1.18 with Python 3.7.1 Command line parameters: -m 20 -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC miRNA_test.fq Processing reads on 1 core in single-end mode ... Finished in 0.06 s (113 us/read; 0.53 M reads/minute). === Summary === Total reads processed: 500 Reads with adapters: 492 (98.4%) Reads that were too short: 64 (12.8%) Reads written (passing filters): 436 (87.2%) Total basepairs processed: 50,500 bp Total written (filtered): 10,909 bp (21.6%) |
Notes:
The cutadapt --help output describes its usage as follows:
From this we see that the
And this says that input reads can also be provided on standard input, if that argument is a hyphen ( - ). So input data can come:
What about cutadapt output (the trimmed reads)? The brackets around the usage -o option indicate that the resulting trimmed FASTQ can be written to a file, but is not by default. This implies that cutadapt by default writes its results to standard output. So output can go
Finally, as we've seen, cutadapt also writes diagnostic output. Where does it go? The usage line doesn't say anything about diagnostics explicitly. But in the Output section of cutadapt --help:
Careful reading of this suggests that:
|
Special care must be taken when removing adapters for paired-end FASTQ files.
Now we're going to run cutadapt on the larger FASTQ files, and also perform paired-end adapter trimming on some yeast paired-end RNA-seq data.
First stage the 4 FASTQ files we will work on:
mkdir -p $SCRATCH/core_ngs/cutadapt cd $SCRATCH/core_ngs/cutadapt cp $CORENGS/human_stuff/Sample_H54_miRNA_L004_R1.cat.fastq.gz . cp $CORENGS/human_stuff/Sample_H54_miRNA_L005_R1.cat.fastq.gz . cp $CORENGS/custom_tracks/Yeast_RNAseq_L002_R1.fastq.gz . cp $CORENGS/custom_tracks/Yeast_RNAseq_L002_R2.fastq.gz . |
Instead of running cutadapt on the command line, we're going to submit a job to the TACC batch system to perform single-end adapter trimming on the two lanes of miRNA data, and paired-end adapter trimming on the two yeast RNAseq FASTQ files.
Paired end adapter trimming is rather complicated, so instead of trying to do it all in one command line we will use one of the handy BioITeam scripts that handles all the details of paired-end read trimming, including all the environment setup.
The BioITeam has an a number of useful NGS scripts that can be executed by anyone on ls6. or stampede2. They are located in the /work/projects/BioITeam/common/script/ directory. For groups that participate in BRCF pods, the scripts are available in /mnt/bioi/script on any compute server. |
The name of the script we want is trim_adapters.sh. Just type the full path of the script with no arguments to see its help information:
/work/projects/BioITeam/common/script/trim_adapters.sh |
You should see something like this:
trim_adapters.sh 2020_04_20
Trim adapters from single- or paired-end sequences using cutadapt. Usage:
trim_adapters.sh <in_fq> <out_pfx> [ paired min_len adapter1 adapter2 ]
Required arguments:
in_fq For single-end alignments, path to input fastq file.
For paired-end alignemtts, path to the the R1 fastq file
which must contain the string 'R1' in its name. The
corresponding 'R2' must have the same path except for 'R1'
out_pfx Desired prefix of output files.
Optional arguments:
paired 0 = single end alignment (default); 1 = paired end.
min_len Minimum sequence length after adapter removal. Default 32.
adapter1 3' adapter. Default GATCGGAAGAGCACACGTCTGAACTCCAGTCAC (NEB).
Specifiy 'illumina' for AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC
(standard Illumina TruSeq3 indexed adapter).
adapter2 5' adapter. Default TGATCGTCGGACTGTAGAACTCTGAACGTGTAGA (NEB).
Specifiy 'illumina' for AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTA
(standard Illumina TruSeq universal adapter).
Environment variables:
show_only 1 = only show what would be done (default not set)
keep 1 = keep intermediate file(s) (default 0, don't keep)
cuta_args other cutadapt options (e.g. '--trim-n --max-n=0.25')
Examples:
export cuta_args='-O 5'; trim_adapters.sh my.fastq.gz h54_b1 1 40
trim_adapters.sh my_fastq.gz yeast_b3 1 28 Illumina Illumina |
Based on this information, here are the 3 cutadapt commands we want to execute:
/work/projects/BioITeam/common/script/trim_adapters.sh Sample_H54_miRNA_L004_R1.cat.fastq.gz H54_miRNA_L004 0 20 /work/projects/BioITeam/common/script/trim_adapters.sh Sample_H54_miRNA_L005_R1.cat.fastq.gz H54_miRNA_L005 0 20 /work/projects/BioITeam/common/script/trim_adapters.sh Yeast_RNAseq_L002_R1.fastq.gz yeast_rnaseq 1 |
Let's put these command into a cuta.cmds commands file. But first we need to learn a bit about Editing files in Linux.
Exercise: Create cuta.cmds file
Use nano or emacs to create a cuta.cmds file with the 3 cutadapt processing commands above. If you have trouble with this, you can copy a pre-made commands file:
cd $SCRATCH/core_ngs/cutadapt cp $CORENGS/tacc/cuta.cmds . |
Or use this "cat to MARKER" trick, also known as an heredoc. The MARKER tag can be anything; below it is EOL.
cd $SCRATCH/core_ngs/cutadapt cat > cuta.cmds << EOL /work/projects/BioITeam/common/script/trim_adapters.sh Sample_H54_miRNA_L004_R1.cat.fastq.gz H54_miRNA_L004 0 20 /work/projects/BioITeam/common/script/trim_adapters.sh Sample_H54_miRNA_L005_R1.cat.fastq.gz H54_miRNA_L005 0 20 /work/projects/BioITeam/common/script/trim_adapters.sh Yeast_RNAseq_L002_R1.fastq.gz yeast_rnaseq 1 EOL |
When you're finished you should have a cuta.cmds file that is 3 lines long (check this with wc -l).
Next create a batch submission script for your job and submit it to the normal queue with a maximum run time of 2 hours.
Since batch jobs can't be submitted from an idev session, make sure you are back on a login node (just exit the idev session).
cd $SCRATCH/core_ngs/cutadapt launcher_creator.py -j cuta.cmds -n cuta -t 01:00:00 -a OTH21164 -q normal sbatch --reservation=CoreNGS-Wed cuta.slurm showq -u # or, if you're not on the reservation: launcher_creator.py -j cuta.cmds -n cuta -t 01:00:00 -a OTH21164 -q development sbatch cuta.slurm showq -u |
How will you know your job is done?
Your cuta job will no longer be displayed in the showq -u output.
|
All our BioITeam scripts, if they complete without errors, will write a line to their logfile that includes the words "completed successfully!". So another way of checking that each command completed is to search for that text in the logfiles. Here we use the powerful grep (general regular expression processor) tool:
|
You should see several log files when the job is finished:
Take a look at the first part of the yeast_rnaseq.acut.pass1.log log file:
|
It will look something like this:
This is cutadapt 1.18 with Python 3.7.1 Command line parameters: -m 32 -a GATCGGAAGAGCACACGTCTGAACTCCAGTCAC --trim-n --paired-output yeast_rnaseq_R2.tmp.cuta.fastq -o yeast_rnaseq_R1.tmp.cuta.fastq Yeast_RNAseq_L002_R1.fastq.gz Yeast_RNAseq_L002_R2.fastq.gz Processing reads on 1 core in paired-end legacy mode ... WARNING: Legacy mode is enabled. Read modification and filtering options *ignore* the second read. To switch to regular paired-end mode, provide the --pair-filter=any option or use any of the -A/-B/-G/-U/--interleaved options. Finished in 151.54 s (24 us/read; 2.55 M reads/minute). === Summary === Total read pairs processed: 6,440,847 Read 1 with adapter: 3,875,741 (60.2%) Read 2 with adapter: 0 (0.0%) Pairs that were too short: 112,847 (1.8%) Pairs written (passing filters): 6,328,000 (98.2%) |
The corresponding yeast_rnaseq.acut.pass2.log file looks like this:
This is cutadapt 1.18 with Python 3.7.1 Command line parameters: -m 32 -a TGATCGTCGGACTGTAGAACTCTGAACGTGTAGA --paired-output yeast_rnaseq_R1.cuta.fastq -o yeast_rnaseq_R2.cuta.fastq yeast_rnaseq_R2.tmp.cuta.fastq yeast_rnaseq_R1.tmp.cuta.fastq Processing reads on 1 core in paired-end legacy mode ... Finished in 83.64 s (13 us/read; 4.54 M reads/minute). === Summary === Total read pairs processed: 6,328,000 Read 1 with adapter: 90,848 (1.4%) Read 2 with adapter: 0 (0.0%) Pairs that were too short: 0 (0.0%) Pairs written (passing filters): 6,328,000 (100.0%) Total basepairs processed: 1,198,172,994 bp Read 1: 639,128,000 bp Read 2: 559,044,994 bp Total written (filtered): 1,197,894,462 bp (100.0%) |
Exercise: Verify that both adapter-trimmed yeast_rnaseq fastq files have 6,328,000 reads
For more on printf, which is available in most programming languages, see https://alvinalexander.com/programming/printf-format-cheat-sheet/ |