Pre-processing of sequencing reads
Once sequencing reads are obtained from the sequencing machine, they need to be pre-processed for further analysis. This step includes the quality control of initial reads and read trimming that includes removing adapter sequences, filtering out low quality reads and trimming reads off low quality base pairs.
QC of sequencing reads
FastQC calculates statistics about the composition and quality of raw sequences, while Fastq Screen looks for possible contaminations.
mkdir QC $RUN fastqc resources/A549_25_3chr10_*.fastq.gz -o ./QC/ Started analysis of A549_25_3chr10_1.fastq.gz Approx 5% complete for A549_25_3chr10_1.fastq.gz Approx 10% complete for A549_25_3chr10_1.fastq.gz Approx 15% complete for A549_25_3chr10_1.fastq.gz Approx 20% complete for A549_25_3chr10_1.fastq.gz ... Approx 85% complete for A549_25_3chr10_2.fastq.gz Approx 90% complete for A549_25_3chr10_2.fastq.gz Approx 95% complete for A549_25_3chr10_2.fastq.gz Analysis complete for A549_25_3chr10_2.fastq.gz
We can display the results with a browser; e.g., Firefox, for each file individually or all files with one command:
firefox QC/A549_25_3chr10_1_fastqc.html firefox QC/*.html
Below is an example of a poor quality dataset. As you can see, the average quality drops dramatically towards the 3’-end.
FastQ Screen requires a number of databases to be installed. Note, the program is included in the Docker and Singularity images but the databases are not. A number of databases prepared by the authors of the program can be downloaded by using the following command (DO NOT LAUNCH IT IN THE CLASS because it will take a while!!! But you will need to do it when completing the Final project):
$RUN fastq_screen --get_genomes
This will download Bowtie indexes for 11 genomes (arabidopsis, drosophila, E. coli, human, lambda, mouse, mitochondria, phiX, rat, worm and yeast) and 3 collection of sequences (adapters, vectors, rRNA). The files will be downloaded in the FastQ_Screen_Genomes folder. The file fastq_screen.conf will be also downloaded in this folder. To use the tool, you will have to modify the fastq_screen.conf by providing the full path to the Bowtie2 executable (/usr/local/bin/bowtie2 if you use our singularity image) and full paths to the downloaded folders with genome index files. Here we show where to change the executables:
# This is a configuration file for fastq_screen ########### ## Bowtie # ########### ## If the bowtie binary is not in your PATH then you can ## set this value to tell the program where to find it. ## Uncomment the line below and set the appropriate location ## #BOWTIE /usr/local/bin/bowtie/bowtie BOWTIE2 /usr/local/bin/bowtie2 ...
FastQ Screen runs check on a random subset of 100,000 reads (that can be changed using option –subset).
To execute FastQ Screen:
$RUN fastq_screen --conf FastQ_Screen_Genomes/fastq_screen.conf \ resources/A549_0_1chr10_1.fastq.gz \ --outdir ./QC/A549_0_1 Using fastq_screen v0.13.0 Reading configuration from 'fastq_screen.conf' Aligner (--aligner) not specified, but Bowtie2 path and index files found: mapping with Bowtie2 Adding database Human Adding database Mouse Adding database Rat Adding database Drosophila Adding database Worm Adding database Yeast Adding database Arabidopsis Adding database Ecoli Adding database rRNA Adding database MT Adding database PhiX Adding database Lambda Adding database Vectors Adding database Adapters Using 7 threads for searches Option --subset set to 100000 reads Processing A549_0_1chr10_1.fastq.gz Counting sequences in A549_0_1chr10_1.fastq.gz Making reduced sequence file with ratio 711:1 ...
Below is an example of the FastQ Screen results for A549_0_1_1.fastq.gz which we prepared.
wget https://biocorecrg.github.io/RNAseq_course_2019/precomp_res/A549_0_1_fastq_screen.tar.gz tar -zvxf A549_0_1_fastq_screen.tar.gz A549_0_1_fastq_screen/ A549_0_1_fastq_screen/A549_0_1chr10_1_screen.html A549_0_1_fastq_screen/A549_0_1chr10_1_screen.txt A549_0_1_fastq_screen/A549_0_1chr10_1_screen.png mv A549_0_1_fastq_screen QC firefox QC/A549_0_1_fastq_screen/A549_0_1chr10_1_screen.html
Initial processing of sequencing reads
Before mapping reads to the genome/transcriptome or performing a de novo assembly, the reads has to be pre-processed, if needed, as follows:
- Demultiplex by index or barcode (it is usually done in the sequencing facility)
- Remove adapter sequences
- Trim reads by quality
- Discard reads by quality/ambiguity
- Filter reads by k-mer coverage (recommended for the de novo assembly)
- Normalize k-mer coverage (recommended for the de novo assembly)
As shown before, both the presence of low quality reads and adapters are reported in the fastqc output.
Adapters are usually expected in small RNA-Seq because the molecules are typically shorter than the reads, and that makes an adapter to be present at 3’-end. Let’s run FastQC on a fastq file for small RNA-Seq.
$RUN fastqc resources/subsample_to_trim.fq.gz -o ./QC Started analysis of subsample_to_trim.fq.gz Approx 5% complete for subsample_to_trim.fq.gz Approx 10% complete for subsample_to_trim.fq.gz Approx 15% complete for subsample_to_trim.fq.gz Approx 20% complete for subsample_to_trim.fq.gz Approx 25% complete for subsample_to_trim.fq.gz ...
Once you got the FastQC report (above), how to figure out the sequence(s) of the adapter(s) that needs to be trimmed?
Let’s use skewer to trim the Illumina small RNA 3’ adapter.
$RUN skewer resources/subsample_to_trim.fq.gz -x TGGAATTCTCGGGTGCCAAGG -o QC/subsample_to_trim .--. .-. : .--': :.-. `. `. : `'.' .--. .-..-..-. .--. .--. _`, :: . `.' '_.': `; `; :' '_.': ..' `.__.':_;:_;`.__.'`.__.__.'`.__.':_; skewer v0.2.2 [April 4, 2016] Parameters used: -- 3' end adapter sequence (-x): TGGAATTCTCGGGTGCCAAGG -- maximum error ratio allowed (-r): 0.100 -- maximum indel error ratio allowed (-d): 0.030 -- minimum read length allowed after trimming (-l): 18 -- file format (-f): Sanger/Illumina 1.8+ FASTQ (auto detected) -- minimum overlap length for adapter detection (-k): 3 Thu Apr 18 17:51:18 2019 >> started |=================================================>| (100.00%) Thu Apr 18 17:51:25 2019 >> done (6.789s) 1000000 reads processed; of these: 30171 ( 3.02%) short reads filtered out after trimming by size control 2220 ( 0.22%) empty reads filtered out after trimming by size control 967609 (96.76%) reads available; of these: 958360 (99.04%) trimmed reads available after processing 9249 ( 0.96%) untrimmed reads available after processing log has been saved to "QC/subsample_to_trim-trimmed.log".
We can look at the read distribution after the trimming of the adapter by inspecting the log-file or re-launching FastQC.
$RUN fastqc QC/subsample_to_trim-trimmed.fastq -o QC Started analysis of subsample_to_trim-trimmed.fastq Approx 5% complete for subsample_to_trim-trimmed.fastq Approx 10% complete for subsample_to_trim-trimmed.fastq Approx 15% complete for subsample_to_trim-trimmed.fastq Approx 20% complete for subsample_to_trim-trimmed.fastq Approx 25% complete for subsample_to_trim-trimmed.fastq Approx 30% complete for subsample_to_trim-trimmed.fastq ...
Let’s explore the tool skewer in more detail, using “skewer –help” command.
- Which parameter indicates the minimum read length allowed after trimming? And what is its default value?
- Which parameter indicates the threshold on the average read quality to be filtered out?
- Using skewer filter out reads in “subsample_to_trim.fq-trimmed.fastq” that have average quality below 30 and trim them on 3’-end until the base quality is reached 30. How many reads were filtered out and how many remained?