-
Notifications
You must be signed in to change notification settings - Fork 7
1. BioLinux and quality control
- BioLinux and crash course to command line
- First look to the OSD102 metagenome
- Clipping the Illumina adapters
- Merging reads
- Quality trimming
First we will open the BioLinux virtual machine. You should have already an installed Virtual Box following these instructions here.
After booting BioLinux, we will unzip the files needed for the course located in the Desktop. Either just double click on the file microbeco_course.tar.gz or use the tar command after changing the directory to Desktop to extract it. For example like this:
cd ~/Desktop
tar xvfz microbeco_course.tar.gzWe can explore the contents of the folder using the file browser or using the tree command on the command line:
tree microbeco_coursemicrobeco_course
├── analysis
│ ├── 0-pre-processing
│ │ ├── 0-adapter-clipping
│ │ ├── 1-read-merging
│ │ │ ├── bbmerge
│ │ │ └── pear
│ │ └── 2-quality-trimming
│ │ ├── merged
│ │ └── paired
│ ├── 1-read-level-analysis
│ │ ├── 0-gene-prediction
│ │ ├── 1-rRNA-identification
│ │ └── 2-sequence-annotation
│ └── 2-assembly-level-analysis
│ ├── 0-metagenomic-assembly
│ │ ├── OSD102_assembly-idba
│ │ └── OSD102_assembly-megahit
│ │ └── intermediate_contigs
│ ├── 1-assembly-evaluation
│ │ ├── 0-metaquast-evaluation
│ │ └── 1-contig-read-coverage
│ │ ├── 0-megahit-assembly
│ │ └── 1-idba_ud-assembly
│ ├── 2-assembly-annotation
│ │ ├── OSD102_contigs_idba_annotation
│ │ └── OSD102_contigs_megahit_annotation
│ └── 3-assembly-gene-abundances
├── data
│ ├── DB
│ └── fastq
└── scripts
Some of the computational steps already have been calculated and the results are included in this folder structure.
The raw data is a compressed FASTQ file. Details about the FASTQ format can be found in the original publication or check wikipedia for a quick description. Below you can find an example of a single read in FASTQ format:
@HWI-M02024:109:000000000-ABR79:1:1101:18881:2628 1:N:0:CCTTGAAA
GTGTTGTCGAAGGCGACCAGGGGCGCACCTTTCTCGACCGCCGCGCCGTCTTCGGCCAGCCATTGCACCGACAGGTTCCACTCGTCGGTCTTGGGCACAGAGAACTCGGTCGAGACGAGCGTGTCGACTTCCCCGGTCAAGAACAGGTGCTGCTCCAAAGCGCCGTGCGTGACGGCGATGACGTCCCGCTTGGCGTCGGGGGATTGCAGCGGCGAGCACCCTGCCAGCGCCAGCAGCACCGACGCCCAGGCAACGCGCCGGCTAGGACTCCGCATGATGCTCCTGCGCGTCAGGGGCCCT
+
CCCCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG5EGGGGGGGGGGDEGGF>EGGGGGGGGGGGGGGGGGGDGGGGGGGGGGGGGGGEGGGGGG6EDEGGGGGGGCFFGFGDDDGGGGG3C<FGG<FGDGDFFFFFF=F?@?FBFB@EEF?>BF>DF
For Illumina derived reads we can find differences depending of the CASAVA version used (the CASAVA software has been recently discontinued by Illumina). Those differences are based in how the quality values are reported, the so called PHRED quality scores. Phred quality scores represent the estimated probability of an error. Phred scores are integer values (0 to 41) encoded using ASCII characters. Depending on the technology and the base caller used, it could be encoded using ASCII 33 or ASCII 64 (that means which ASCII character is used as start point for the encoding, in case of phred33, the ASCII character 33 is ! and for phred64, the ASCII character 64 is @).
In the case of the OSD reads, CASAVA version is 1.8+, and are encoded using phred33. Here you can find a nice explanation of the rationale behind the quality scores. High-throughput sequencing can present different pathologies, a nice read is this paper where several of this pathologies are described.
To have an idea of how our raw reads look, we will use FASTQC to inspect the quality of the reads. We can start FASTQC using the BioLinux dash or in the command line:
fastqc &Once FASTQC is started we can open both OSD fastq files found in ~/Desktop/microbeco_course/data/fastq
Have a look at the report. Can you spot problems on the raw reads? In case the sequences present problems, can you think on how we can solve them?
Illumina libraries (we will only describe TruSeq used for the OSD) are constructed by ligating adapters that will bind to the flowcell to short fragments of DNA. Here and here you can find a more detailed description of the Illumina chemistry.
As you have observed in the FASTQC report, we have a high presence of Illumina universal adapters. That will be a problem during the downstream analyses and submit for the submission of the reads to ENA or NCBI SRA where technical sequences must be removed before submission. To remove those adapters there are plenty of programs that perform this task like Trimmomatic, Cutadapt, Sickle, BBTools and many others. For the course we will use bbduk from BBTools. It is implemented in Java and therefore architecture independent (can run on Windows, Linux or Mac OS X). It fast and very accurate in clipping adapters. In the command line we need to type:
cd ~/Desktop/microbeco_course/analysis/0-pre-processing/0-adapter-clipping
bbduk.sh -Xmx1g in1=~/Desktop/microbeco_course/data/fastq/OSD102_R1.fastq.gz in2=~/Desktop/microbeco_course/data/fastq/OSD102_R2.fastq.gz out1=OSD102_R1.ac.fastq.gz out2=OSD102_R2.ac.fastq.gz ref=/usr/local/bioinf/bbmap/resources/truseq.fa.gz ktrim=r k=23 mink=11 hdist=1 tpe tboBBtools already contains a file with a collection of adapters in fasta format, in our case we use the TruSeq found in the file /usr/local/bioinf/bbmap/resources/truseq.fa.gz. In case you are using different adapters you should add them to a file. The different parameters used with bbduk are ref for the adpater file, ktrim=r specifies right-trimming (3'-adapters), k=23 the k-mer size to use, kmin=11 minimum k-mer size, hdist=1 allow one mismatch, tbo to trim adapters based on pair overlap detection and tpe to trim both reads to the same length. For more detailed information you can execute bbduk.sh without any parameter.
We can check the effect of the adapter clipping with FASTQC on the new files OSD102_R1.ac.fastq.gz and OSD102_R2.ac.fastq.gz
Does the adapter clipping fixed all the problems we observed on the raw reads?
In our example we have paired-end reads that can be merged if overlapping region exist. Merging reads will result in a longer read that can be useful for different downstream analyses like annotation, reducing the number of errors or simply reducing the size of the dataset. There are many different programs for read merging, some of them are FLASH, PANDASeq, COPE, SeqPrep, PEAR, BBmerge... We will use PEAR and BBmerge.
For BBmerge we do in the command line:
cd ~/Desktop/microbeco_course/analysis/0-pre-processing/1-read-merging/bbmerge
bbmerge.sh in1=~/Desktop/microbeco_course/analysis/0-pre-processing/0-adapter-clipping/bbduk/OSD102_R1.ac.fastq.gz in2=~/Desktop/microbeco_course/analysis/0-pre-processing/0-adapter-clipping/bbduk/OSD102_R2.ac.fastq.gz out=OSD102.me.bbmerge.fastq outu=OSD102.R1.bbmerge.fastq outu2=OSD102.R2.bbmerge.fastqHere, this forum thread explains what BBmerge does with the quality scores in the overlapping region.
if b1 equals b2:
q=max(q1, q2)+min(q1, q2)/4
(capped at 41 by default, though this can be overridden)
if b1 does not equal b2:
q=max(q1, q2)-min(q1, q2)
For PEAR:
cd ~/Desktop/microbeco_course/analysis/0-pre-processing/1-read-merging/pear
pear -f ~/Desktop/microbeco_course/analysis/0-pre-processing/0-adapter-clipping/bbduk/OSD102_R1.ac.fastq.gz -r ~/Desktop/microbeco_course/analysis/0-pre-processing/0-adapter-clipping/bbduk/OSD102_R2.ac.fastq.gz -o OSD102_meAfter PEAR is done, we rename to something easier and removed discarded
mv OSD102_me.assembled.fastq OSD102.me.pear.fastq
mv OSD102_me.unassembled.forward.fastq OSD102.R1.pear.fastq
mv OSD102_me.unassembled.reverse.fastq OSD102.R2.pear.fastq
rm OSD102_me.discarded.fastqHow many reads have been assembled by each program? Check the resulting merged files in FASTQC, do we see anything special? Why do we merge before quality filtering?
After the reads have been merged, we can trim them based on the PHRED quality scores. We will use a Q-value of 20 to retain high quality reads. We will use BBduk for quality trimming. BBduk uses the phred algorithm (Mott) to trim the sequences, here you can find a detailed explanation of the algorithm. First we will quality trim the merged reads, and afterwards the ones that remain single. In a command line:
cd ~/Desktop/microbeco_course/analysis/0-pre-processing/2-quality-trimming/merged
bbduk.sh -Xmx1g in=~/Desktop/microbeco_course/analysis/0-pre-processing/1-read-merging/bbmerge/OSD102.me.bbmerge.fastq out=OSD102.me.qc.fastq trimq=20 qtrim=rl minlen=100 maxns=0
bbduk.sh -Xmx1g in=~/Desktop/microbeco_course/analysis/0-pre-processing/1-read-merging/bbmerge/OSD102.R1.bbmerge.fastq in2=~/Desktop/microbeco_course/analysis/0-pre-processing/1-read-merging/bbmerge/OSD102.R2.bbmerge.fastq outs=OSD102.SR.qc.fastq out=OSD102.PE.qc.fastq trimq=20 qtrim=rl minlen=100 maxns=0The different parameters used by bbmerge.sh are trimq=20 to trim to a quality of 20 with the phred algorithm, qtrim=rl to trim left and right sides, minlen=100 to discard reads shorter than 100nt and maxns=0 to remove reads with ambiguous bases. After all reads have been quality trimmed and filtered we can concatenate them in a file:
cat OSD102.me.qc.fastq OSD102.PE.qc.fastq OSD102.SR.qc.fastq > OSD102.qc.fastqAdditionally, we will quality trim the original paired reads before merging. We will use them for assembly later. In the command line:
bbduk.sh -Xmx1g in=~/Desktop/microbeco_course/analysis/0-pre-processing/0-adapter-clipping/bbduk/OSD102_R1.ac.fastq.gz in2=~/Desktop/microbeco_course/analysis/0-pre-processing/0-adapter-clipping/bbduk/OSD102_R2.ac.fastq.gz outs=OSD102.SR.qc.fastq out=OSD102.R1.qc.fastq out2=OSD102.R2.qc.fastq trimq=20 qtrim=rl minlen=100 maxns=0 How many reads are discarded? Check the resulting reads with FASTQC, do we see any improvement?