cgeb2001's DokuWiki!

User Tools

Site Tools

Trace: • gene_prediction_just_augustus • multi-gene_phylogeny_pipeline
nanopore_tools_for_polishing

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Both sides previous revisionPrevious revision
Next revision		Previous revision
nanopore_tools_for_polishing [2018/07/17 12:54] – 129.173.91.14nanopore_tools_for_polishing [2024/08/07 13:01] (current) – 134.190.232.164
Line 1: Line 1:
 ====== Polishing your MinION assembly ====== ====== Polishing your MinION assembly ======
-Documentation by Jon Jerlstrom Hultqvist and Shelby Williams+Documentation by Jon Jerlström Hultqvist and Shelby Williams 
 +(updates by Joran Martijn) 
 + 
 +**Be aware that some scripts and commands might not be working any longer on Perun due to the switch to the new conda-environment system. Sections will be progressively updated to reflect this.**
  
 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. 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.
Line 17: Line 20:
 {{ :bman_racon_1sh.doc |}} {{ :bman_racon_1sh.doc |}}
  
-Formatting for this script:+The following is a script using default settings:
 <code> <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 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> </code>
-Note: for each polish performed, the script must be changed to reflect the new file names.+**Note 1**: for each polish performed, the script must be changed to reflect the new file names.
  
 +For genomes with lots of stubborn indels that persist after Nanopolishing and Unicycling, play around with the flags -e (error threshold), -w (window size), and -q (coverage threshold). I found that increasing window size and decreasing coverage threshold works well. The following is a script that contains one round of minimap2 mapping (read **Note 2** on interleaved short reads), Racon polishing, then Bowtie2 mapping (for viewing bam files):
 +<code>
 +#!/bin/bash
 +#$ -S /bin/bash
 +. /etc/profile
 +#$ -cwd
 +#$ -pe threaded 8
 +
 +cd /your_working_directory/
 +input=your_input_genome.fasta
 +output=your_output_racon_round1.fasta
 +minimap2 -t 8 $input interleavedshortreads.fq > temporary.paf
 +echo "minimap2 done"
 +racon -u -e 0.1 -w 500 -q 1 -t 8 interleavedshortreads.fq temporary.paf $input >$output
 +echo "racon done"
 +rm temporary.paf
 +source activate bowtie2
 +bowtie2-build $output $output
 +bowtie2 -k 2 --threads 8 -x $output --very-sensitive \
 +-1 your_forward_short_reads.fq \
 +-2 your_reverse_short_reads.fq \
 +-S $output.sam
 +echo "Bowtie done"
 +conda deactivate
 +samtools view -F 4 -bS $output.sam |samtools sort > $output.sorted.bam
 +samtools index $output.sorted.bam > $output.sorted.bam.bai
 +rm $output.sam
 +rm $output.*.bt2
 +</code>
 +You can modify the above script to loop for multiple rounds, but if you're looking for the best settings for correcting errors always check the bam files by eye after each round. 
 +
 +**Note 2**: The "**interleavedshortreads.fq**" must not contain blanks/white spaces in headers. You can simply concatenate forward and reverse reads into one read file, but making them interleaved is useful for other read correction tools. To obtain interleaved reads:
 +<code>
 +source activate fastx_toolkit
 +fastq-interleave read1.fq read2.fq > interleaved.fq
 +</code>
  
 +**WARNING**: Keep an eye on genome/contig size difference after each round of polishing. Using window size (-w) above 5000 may remove repetitive regions or regions with no short-read coverage. 
 ---- ----
  
Line 42: Line 82:
  
 First, make a BWA index of the assembly you wish to map onto by using the following command: First, make a BWA index of the assembly you wish to map onto by using the following command:
 +
 <code> <code>
 bwa index assembly_to_polish.fasta bwa index assembly_to_polish.fasta
 </code> </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: 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> <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+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 sort --threads 16 -o piloninput.sorted.bam
 </code> </code>
 +
 +UPDATE: You can now run bwa-mem2, which is an optimized version of bwa mem. It generates the exact same output, but is 2-4x faster:
 +
 +<code>
 +bwa-mem2 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 sort --threads 16 -o 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: 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> <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+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> </code>
 +
 +UPDATE: The --threads option is as of v1.24 no longer maintained. It seems Pilon doesn't use more than 200-300% CPU (i.e. 3 threads) at most, so setting --threads to 4 orso should be sufficient.
 +
 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: 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> <code>
 samtools index /path/to_bam_file samtools index /path/to_bam_file
 </code> </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. 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.
  
Line 63: Line 134:
 Shell script: Shell script:
 {{ :unicyclersh.docx |}} {{ :unicyclersh.docx |}}
 +
 +<code>
 +#!/bin/bash
 +#$ -S /bin/bash
 +. /etc/profile
 +#$ -cwd
 +#$ -pe threaded 16
 +
 +#cd /scratch2/jon/MinION/BMAN/assemblies/Unicycler_polish/
 +
 +echo "Starting"
 +
 +unset PYTHONPATH
 +export PATH=/scratch2/software/gcc-6.3.0/bin:/scratch2/software/Python-3.6.0/bin:$PATH
 +export LD_LIBRARY_PATH=/scratch2/software/gcc-6.3.0/lib64:/scratch2/software/Python-3.6.0/lib:$LD_LIBRARY_PATH
 +
 +/scratch2/software/Python-3.6.0/bin/unicycler_polish -1 /scratch2/shelbyw/RCL_Unicycler/RCL_1_PairNtrim.fq -2 /scratch2/shelbyw/RCL_Unicycler/RCL_2_PairNtrim.fq --long_reads RCL_MinION.CutAdapt75.3000.chop.fastq.gz -a RCL_unclean_AB_assembly_fix_Racon2_Pilon3.fasta --pilon=/scratch2/software/pilon/pilon-1.22.jar --samtools=/opt/perun/bin/samtools --threads 16
 +
 +
 +echo "Done!"
 +
 +</code>
  
 Formatting: Formatting:
Line 100: Line 193:
 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. 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+Location: /scratch2/software/anaconda/envs/nanopolish-0.12/bin
  
 Scripts: Scripts:
Line 119: Line 212:
 nanopolish merge - merges the pieces into new a new consensus. nanopolish merge - merges the pieces into new a new consensus.
  
-**Updated Nanopolish protocol (as of Dec 1st 2017):** +**Updated Nanopolish protocol (as of March 8 2020):** 
  
 First, index your unchopped, raw reads file.  First, index your unchopped, raw reads file. 
 Use the sequencing_summary.txt produced by albacore during basecalling to speed up this step significantly. If you have several sequencing_summary.txt files these can be placed in a fof-file with the  path to the txt-file and called by -f. This also works in case of a single-file.: Use the sequencing_summary.txt produced by albacore during basecalling to speed up this step significantly. If you have several sequencing_summary.txt files these can be placed in a fof-file with the  path to the txt-file and called by -f. This also works in case of a single-file.:
 +for the following step **DO NOT** **use more than 1** thread because the program is not threaded!
 <code> <code>
 #!/bin/bash #!/bin/bash
Line 128: Line 222:
 . /etc/profile . /etc/profile
 #$ -cwd #$ -cwd
-#$ -pe threaded 4 +#$ -pe threaded 1
 cd $PWD cd $PWD
  
-export PATH=/scratch2/software/anaconda/bin:$PATH +fast5path=/scratch2/path2/fast5/ 
-source activate nanopolish-python3 +fastq=/path2fastqlongreads.fastq 
- +seqsummary=/path2tosequencing_summary.txt 
-/scratch2/software/anaconda/envs/nanopolish-python3/bin/nanopolish index \ +source activate nanopolish-0.13.2 
--d /path/to/fast5/directory/ \ +export PATH=/scratch2/software/anaconda/envs/nanopolish-0.13.2/bin:$PATH 
--f summary_files.fof \ +nanopolish index -d $fast5path -s $seqsummary $fastq 
-/path/to/reads.fastq+conda deactivate
  
-source deactivate 
  
 </code> </code>
Line 165: Line 257:
 samtools index reads.sorted.bam samtools index reads.sorted.bam
  
-source deactivate+conda deactivate
  
 </code> </code>
Line 185: Line 277:
 cd $PWD cd $PWD
  
-export PATH=/scratch2/software/anaconda/bin:$PATH +export PATH=/scratch2/software/anaconda/envs/nanopolish-0.12/bin:$PATH 
-source activate nanopolish-python3+ 
 +source activate nanopolish-0.12
  
-python /scratch2/software/anaconda/envs/nanopolish-python3/bin/nanopolish_makerange.py \ +nanopolish_makerange.py reference.fasta | parallel --results nanopolish.results -P 20 \ 
-reference.fasta | parallel --results nanopolish.results -P 20 \ +nanopolish variants --consensus -o polished.{1}.fa -w {1} \ 
-/scratch2/software/anaconda/envs/nanopolish-python3/bin/nanopolish variants --consensus polished.{1}.fa -w {1} \ +-r /path/to/reads.fastq -b reads.sorted.bam -g reference.fasta -t 10 --min-candidate-frequency 0.1
--r /path/to/reads.fastq \ +
--b reads.sorted.bam -g reference.fasta -t 10 --min-candidate-frequency 0.1+
  
-source deactivate+conda deactivate
  
 </code> </code>
nanopore_tools_for_polishing.1531842878.txt.gz · Last modified: by 129.173.91.14