Polishing with Racon
Racon is intended as a standalone consensus module to correct raw contigs generated by rapid assembly methods which do not include a consensus step. The goal of Racon is to generate genomic consensus which is of similar or better quality compared to the output generated by assembly methods which employ both error correction and consensus steps, while providing a speedup of several times compared to those methods. It supports data produced by both Pacific Biosciences and Oxford Nanopore Technologies.
Racon can be used as a polishing tool after the assembly with either Illumina data or data produced by third generation of sequencing. The type of data inputed is automatically detected.
Racon takes as input only three files: contigs in FASTA/FASTQ format, reads in FASTA/FASTQ format and overlaps/alignments between the reads and the contigs in MHAP/PAF/SAM format. Output is a set of polished contigs in FASTA format printed to stdout. All input files can be compressed with gzip (which will have impact on parsing time).
We are going to use racon to do an initial correction. The medaka documentation advises to do four rounds with racon before polishing with medaka since medaka has been trained with racon polished assemblies. We are only doing one round here.
Mapping of Nanopore reads to the assembly
In order to use racon, we need a mapping of the reads to assembly. We use bwa for this task.
First we need to create an index on our assembly, this has already been done for the pilon polishing:
#already done for pilon
bwa index ~/workdir/assembly/assembly.contigs.fasta
Then we will run the mapping. Check the usage of bwa mem:
Usage: bwa mem [options] <idxbase> <in1.fq> [in2.fq]
Algorithm options:
-t INT number of threads [1]
-k INT minimum seed length [19]
-w INT band width for banded alignment [100]
-d INT off-diagonal X-dropoff [100]
-r FLOAT look for internal seeds inside a seed longer than {-k} * FLOAT [1.5]
-y INT seed occurrence for the 3rd round seeding [20]
-c INT skip seeds with more than INT occurrences [500]
-D FLOAT drop chains shorter than FLOAT fraction of the longest overlapping chain [0.50]
-W INT discard a chain if seeded bases shorter than INT [0]
-m INT perform at most INT rounds of mate rescues for each read [50]
-S skip mate rescue
-P skip pairing; mate rescue performed unless -S also in use
-e discard full-length exact matches
Scoring options:
-A INT score for a sequence match, which scales options -TdBOELU unless overridden [1]
-B INT penalty for a mismatch [4]
-O INT[,INT] gap open penalties for deletions and insertions [6,6]
-E INT[,INT] gap extension penalty; a gap of size k cost '{-O} + {-E}*k' [1,1]
-L INT[,INT] penalty for 5'- and 3'-end clipping [5,5]
-U INT penalty for an unpaired read pair [17]
-x STR read type. Setting -x changes multiple parameters unless overridden [null]
pacbio: -k17 -W40 -r10 -A1 -B1 -O1 -E1 -L0 (PacBio reads to ref)
ont2d: -k14 -W20 -r10 -A1 -B1 -O1 -E1 -L0 (Oxford Nanopore 2D-reads to ref)
intractg: -B9 -O16 -L5 (intra-species contigs to ref)
Input/output options:
-p smart pairing (ignoring in2.fq)
-R STR read group header line such as '@RG\tID:foo\tSM:bar' [null]
-H STR/FILE insert STR to header if it starts with @; or insert lines in FILE [null]
-j treat ALT contigs as part of the primary assembly (i.e. ignore <idxbase>.alt file)
-v INT verbose level: 1=error, 2=warning, 3=message, 4+=debugging [3]
-T INT minimum score to output [30]
-h INT[,INT] if there are <INT hits with score >80% of the max score, output all in XA [5,200]
-a output all alignments for SE or unpaired PE
-C append FASTA/FASTQ comment to SAM output
-V output the reference FASTA header in the XR tag
-Y use soft clipping for supplementary alignments
-M mark shorter split hits as secondary
-I FLOAT[,FLOAT[,INT[,INT]]]
specify the mean, standard deviation (10% of the mean if absent), max
(4 sigma from the mean if absent) and min of the insert size distribution.
FR orientation only. [inferred]
Note, that there is an option for Oxford Nanopore 2D-reads:
-x STR read type. Setting -x changes multiple parameters unless overridden [null]
pacbio: -k17 -W40 -r10 -A1 -B1 -O1 -E1 -L0 (PacBio reads to ref)
ont2d: -k14 -W20 -r10 -A1 -B1 -O1 -E1 -L0 (Oxford Nanopore 2D-reads to ref)
intractg: -B9 -O16 -L5 (intra-species contigs to ref)
We use this default option for our mapping (Note that we need a sam file for racon):
cd ~/workdir/
mkdir nanopore_mapping
bwa mem -t 14 -x ont2d ~/workdir/assembly/assembly.contigs.fasta ~/workdir/basecall/basecall_trimmed.fastq.gz > ~/workdir/nanopore_mapping/mapping.sam
Run racon
Check the usage of racon:
racon --help
usage: racon [options ...] <sequences> <overlaps> <target sequences>
<sequences>
input file in FASTA/FASTQ format (can be compressed with gzip)
containing sequences used for correction
<overlaps>
input file in MHAP/PAF/SAM format (can be compressed with gzip)
containing overlaps between sequences and target sequences
<target sequences>
input file in FASTA/FASTQ format (can be compressed with gzip)
containing sequences which will be corrected
options:
-u, --include-unpolished
output unpolished target sequences
-f, --fragment-correction
perform fragment correction instead of contig polishing
(overlaps file should contain dual/self overlaps!)
-w, --window-length <int>
default: 500
size of window on which POA is performed
-q, --quality-threshold <float>
default: 10.0
threshold for average base quality of windows used in POA
-e, --error-threshold <float>
default: 0.3
maximum allowed error rate used for filtering overlaps
-m, --match <int>
default: 5
score for matching bases
-x, --mismatch <int>
default: -4
score for mismatching bases
-g, --gap <int>
default: -8
gap penalty (must be negative)
-t, --threads <int>
default: 1
number of threads
--version
prints the version number
-h, --help
prints the usage
Then we can call racon with our mapping, the read file and the assembly file. We use 14 threads to do this:
cd ~/workdir/
mkdir racon
racon -m 8 -x -6 -g -8 -w 500 -t 14 ~/workdir/basecall/basecall_trimmed.fastq.gz ~/workdir/nanopore_mapping/mapping.sam ~/workdir/assembly/assembly.contigs.fasta > racon/racon.fasta
The options:
-m 8 -x -6 -g -8 -w 500
are used, because they were also used for the training of medaka and we want to have similar error profiles of the draft.