cgeb2001's DokuWiki!

You've loaded an old revision of the document! If you save it, you will create a new version with this data.

====== Polishing your MinION assembly ======
Documentation by Jon Jerlstrom Hultqvist and Shelby Williams

The following are tools for polishing Oxford Nanopore assemblies, to help eliminate some of the error associated with the long-read data. Each polisher uses a mapping step prior to correcting the assemblies. BWA MEM is used here as the mapping tool, but can be replaced with another mapper.

----

**Racon** [[https://github.com/isovic/racon]]

Racon is a super quick option for polishing assemblies. It is made to work with assemblies made with miniasm (another quick option) but can easily be implemented for assemblies created by Abruijn, Canu, ect. Racon works best when it is used in iteration; an assembly should be polished more than once for best results. In this way, Racon can improve upon polishing that was attempted by assemblers, such as Abruijn.

Racon uses three input files: the assembly file you wish to polish (in fasta format), the original raw reads of that assembly (in fastq format), and an overlap file (in PAF or MHAP). This overlap file can be easily created using Minimap, and this step of the polishing process is already integrated into the shell script:

Shell script:
{{ :bman_racon_1sh.doc |}}

Formatting for this script:
<code>
minimap -t8 assembly.fasta /scratch2/user/path/to/raw_reads.fastq | racon /scratch2/user/path/to/raw_reads.fastq - assembly.fasta output_assembly.polish1.fasta
</code>
Note: for each polish performed, the script must be changed to reflect the new file names.

----

**Pilon** [[https://github.com/broadinstitute/pilon/wiki]]

Pilon uses mapped illumina reads to improve the assembly quality in several ways. It will try to correct small variants, local missamblies, collapsed repeats etc. depending on the settings. This Pilon pipeline uses BWA MEM to map the Illumina reads to the assembly, and utilizes the read file to find and correct mismatches between the two sets of data.

This polisher should also be used to be used iteratively to progressively improve the quality. The reason being, a corrected assembly will attract more mapped reads which will allow progressively better assemblies. This fact is used by the **Unicycler** polish module.

This polish is a great choice to progressively increase the quality of your assemblies as a last step, but since it was developed for bacterial genomes specifically and is uses very sensitive parameters, it is not recommended for the first polish of an assembly. However, it should be fine to use for polishing of the small variants of the genome.

Shell scripts:
{{ :meteora_bwash.doc |}}
{{ :pilonsh.doc |}}

Usage:

First, make a BWA index of the assembly you wish to map onto by using the following command:
<code>
bwa index assembly_to_polish.fasta
</code>
Next, use the meteora_bwa.sh script to map the short reads onto your assembly. This will create a sorted.bam file. In this example, two paired-end read files will be mapped:
<code>
bwa mem -t 16 assembly_to_polish.fasta /scratch2/user/path/to/trimmedreads_1_PairNtrim.fastq.gz /scratch2/user/path/to/trimmedreads_2_PairNtrim.fastq.gz | samtools view -Sb - | samtools sort >  piloninput.sorted.bam
</code>
Once this is finished, use Pilon.sh to make changes in the assembly and generate a new consensus sequence. Pilon.sh can be formatted like so:
<code>
java -Xmx16G -jar /scratch2/software/pilon/pilon-1.22.jar --genome assembly_to_polish.fasta --frags piloninput.sorted.bam --output P2x --outdir Pilon2x --threads 16
</code>
You may run into an error where Pilon does not recognize the bam file created from the previous step as being indexed. To fix this, run:
<code>
samtools index /path/to_bam_file
</code>
This will return a .bam.bai file. This file needs to be in the same folder as Pilon.sh, but does not need to be placed in the script.

**Alternatively:** you can use **Unicycler** to do iterative steps of Pilon. This saves you the trouble of having to edit the shell script with the names of new reads. Unicycler can be used to do a huge variety of tasks. Here, we only use it to utilize Pilon. Unicycler will use several different Pilon settings to polish the assembly, and requires paired end trimmed Illumina reads, the trimmed long reads, and the assembly to be fixed. [[https://github.com/rrwick/Unicycler]]

Shell script:
{{ :unicyclersh.docx |}}

Formatting:
<code>
/scratch2/software/Python-3.6.0/bin/unicycler_polish -1 /scratch2/user/path_to/Illumina_shortreads_1_PairNtrim.fq -2 /scratch2/user/path_to/Illumina_shortreads_2_PairNtrim.fq --long_reads longreads_trimmed.fastq.gz -a assembly_to_clean.fasta --pilon=/scratch2/software/pilon/pilon-1.22.jar --samtools=/opt/perun/bin/samtools --threads 16
</code>

*Unicycler does not like colons ":" in headers. Replace them with a sed command such as:
<code>
sed -i 's/:/_/g' filename
</code>

There is also a handy perl script "translate\_fasta\_headers" that can be used to generate clean headers that are accepted by most programs (https://github.com/nylander/translate_fasta_headers). This porgram keeps the old headers in a tabulated file so that they can be restored at the users request.

## USAGE

./translate_fasta_headers.pl [--tabfile=<tabfile.tab>] [--in=<in.fas>] [--out=<out.fas>] <in.fas>

#### From long to short labels:

./translate_fasta_headers.pl --out=out.fas in.fas

#### An back, using a translation table:

./translate_fasta_headers.pl --tabfile=out.fas.translation.tab out.fas

----

**Nanopolish** [[https://github.com/jts/nanopolish]]

Nanopolish uses the raw event data that is extracted from FAST5 files. The reads are aligned using bwa mem and the consensus algorithm polishes the assembly in pieces. Finally the pieces are merged to create a final assembly. The final assembly appears to be more accurate but it is very computer intensive. The resulting polished assemblies appear to reach about 99.5% or higher but are definitely not perfect. People have suggested that running Racon 2 rounds before doing nanopolish will speed up this. Since Racon is comparatively quick it is strongly recommend that you try this.

Please be aware! Headers of fasta file should not contain any commas (:) or pipes (|), or it will fail in the nanopolish_consensus step!!.  So any assembly out of Abruijn will need to have commas removed.

If illumina reads are available it might be possible to skip nanopolish altogether and go directly to Pilon polishing after Racon. This has been exemplified in the Solanum penellii pre-print where nanopolish simply was not feasible.

Location: /scratch2/software/nanopolish

Scripts:

{{ :nanopolish_align_busseltonsh.doc |}}
{{ :nanopolish_consensus_busselton_jun6sh.doc |}}
{{ :nanopolish_extractsh.doc |}}
{{ :nanopolish_merge_jun5sh.doc |}}

A quick overview of the order to use each script, and what they do:

nanopolish_extract - extracts data from a directory of FAST5 files. Do not do this if you used albacore 2.0 or higher to basecall.

nanopolish_align - uses bwa mem -x ont2d mappings.

nanopolish_consensus - computes the new consensus in batches of 50 kb.

nanopolish merge - merges the pieces into new a new consensus.

**Updated Nanopolish protocol (as of Dec 1st 2017):**

First, index your unchopped, raw reads file:
<code>
nanopolish index -d /path/to/raw_fast5s/ albacore_output.fastq
</code>

This will create four files with the following extensions: .index, .index.gzi, .index.fai, .index.readdb .

Next, index your draft genome (which has already been polished with Racon once or twice). You do not need to do this if you're using **minimap2** to create the alignment .bam file. I recommend splitting your draft genome into 5Mb pieces before doing this. **pyfasta split** is a handy tool for splitting. Either through the command line or a shell script, do:
<code>bwa index piece_of_draft_genome.fa </code>
Then alignment file .bam must be made:
<code>bwa mem -x ont2d -t 8 draft.fa reads.fastq | samtools sort -o reads.sorted.bam -T reads.tmp
samtools index reads.sorted.bam </code> "reads" must correspond to the name of the raw fastq and its nanoplish index file names. The "draft.fa" must also correspond to the split piece of the draft genome you bwa indexed.

Finally, build the actual consensus sequence:
<code>
python nanopolish_makerange.py draft.fa | parallel --results nanopolish.results -P 8 \
    nanopolish variants --consensus polished.{1}.fa -w {1} -r reads.fa -b reads.sorted.bam -g draft.fa -t 4 --min-candidate-frequency 0.1
</code>
This will create new pieces of consensus sequences in the current working directory (not in the nanopolish.results folder indicated in the script above)  named polished.*.fa. Merge them by:
<code>
python nanopolish_merge.py polished.*.fa > polished_genome.fa </code>

Since you'll have several pieces of polished_genome.fa because you split the original draft genome, just **cat** them all into one file in the end.

Now, move on to Pilon or Unicycle!

Please fill all the letters into the box to prove you're human. N S D U O Please keep this field empty:

Edit summary

Note: By editing this page you agree to license your content under the following license: CC Attribution-Noncommercial-Share Alike 4.0 International