-
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 most computational steps have been already calculated and the results included on the folder structure.
The raw data is a compressed FASTQ file. Details for the FASTQ format can be found on the original publication or for a quick description we can look at wikipedia. Below you can find an example of a 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 of 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-througput 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 opened we can open both OSD fastq files found in ~/Desktop/microbeco_course/data/fastq
Have a look to 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 as are the ones 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 or if we want to submit 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 perfom this task like Trimmomatic, Cutadapt, Sickle, BBTools and many otheres. For the course we will use bbduk from BBTools, is a java programmed fast tool, architecture indpendent (can run in Windows, Linus or Mac OS X) and very accurate on clipping adapters. In the command line we:
cd ~/Desktop/microbeco_course/analysis/0-pre-processing/0-adapter-clipping
bbduk.sh -Xmx1g in1=~/Desktop/microbeco_course/data/OSD102_R1_shotgun_raw.fastq.gz in2=~/Desktop/microbeco_course/data/OSD102_R2_shotgun_raw.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 exists an overlapping region. Merging reads will result in a longer read that can be useful for different downstream analyses like annotation, reducing the number errors or simply reducing the size of the dataset. There are many differents 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 a terminal:
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 is explained 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 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 trim the reads. 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, abd after the ones that remain singel. In a terminal:
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 ambiguos 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.fastqAlso we will quality trim the original paired reads before merging, later we will use them for assembly
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?