Validating Whole Genome Nanopore Sequencing, using Usutu Virus as an Example

NOTE: List of software tools to be used: usearch v11.0.667; muscle v3.8.1551; porechop 0.2.4; cutadapt 2.5; minimap2 2.16-r922; samtools 1.9; trimmomatic 0.39; bbmap 38.33; spades v3.13.1; kma-1.2.8

1. Primer design

1. Start with downloading or retrieving a set of relevant reference whole genome sequences from public or private data collections. For instance, retrieve all full length USUV genomes (taxid64286) from the NCBI database²². USUV encodes a genome of around 11,000 nucleotides so only retrieve the sequences with a sequence length of 8,000-12,000 nucleotides. Do this using the following search entry:
   – taxid64286[Organism:noexp] AND 8000[SLEN]:12000[SLEN].
   1. Click on Send to | Complete Record | File; use Format = FASTA and create the File.
2. To downsize the set of reference sequences, remove duplicate sequences or sequences with over 99% nucleotide identity from the dataset. Do this using the cluster fast option from usearch²³. On the command line enter:
   – usearch -cluster_fast All_USUV.fasta -id 0.99 -centroids All_USUV_dedup.fasta
3. To generate the primers, sequences need to be aligned. This is done using MUSCLE²⁴. On the command line enter:
   – muscle -in All_USUV_dedup.fasta -out All_USUV_dedup_aligned.fasta -log log_muscle.txt
   NOTE: It is essential to manually inspect the alignment to check for discrepancies. These can be manually corrected if needed and the ends can be trimmed according to the length of most whole genome sequences.
4. Primal is used to make a draft selection of the primers which can be used for full length amplicon sequencing¹⁹. Upload the alignment to the primal website (http://primal.zibraproject.org/) and select the preferred amplicon length and overlap length between the different amplicons. Go to primal.zibraproject.org, fill in the Scheme name, upload the aligned fasta file, select the amplicon length, overlap size, and generate the scheme.
5. Align the complete set of available complete USUV sequences (not the downsized or deduplicated set). On the command line enter:
   – muscle -in All_USUV.fasta -out All_USUV_aligned.fasta -log log_muscle.txt
   NOTE: Map the generated primers against the complete alignment (do not use the deduplicated alignment), manually correct errors and include a maximum of 5 degenerative primer positions.

2. Multiplex PCR

1. Perform the multiplex PCR using the designed primers and nanopore and Illumina sequencing. The multiplex PCR for USUV was performed as previous described^19,21.
2. Perform basecalling with flip-flop version 3.0.6.6+9999d81.

3. Data analysis to generate consensus sequences from nanopore data

1. Several samples can be multiplexed on a single nanopore sequencing run. After performing the sequence run, demultiplex the nanopore data. Use Porechop²⁵ for this. To prevent contamination and enhance accuracy, use the require_two_barcodes flag. On the command line enter:
   – porechop -i Run_USUV.fastq -o Run_USUV_demultiplex –require_two_barcodes
2. After demultiplexing, remove primer sequences (indicated in the file Primers_Usutu.fasta in both orientations) using cutadapt²⁶. In addition, remove sequences with a length shorter than 75 nucleotides. The primers have to be removed since they can introduce artificial biases in the consensus sequence. On the command line enter:
   – cutadapt -b file:Primers_USUV.fasta -o BC01_trimmed.fastq BC01.fastq -m 75
3. Demultiplexed sequence reads can be mapped against a panel of distinct reference strains using minimap2²⁷ and a consensus sequence can be generated using samtools²⁸. Follow the example below which shows the procedure of a reference-based alignment and the consensus sequence generation of one sample: BC01. On the command line enter:
   – minimap2 -ax map-ont Random_Refs_USUV.fasta BC01_trimmed.fastq > BC01.bam
   – samtools sort BC01.bam > BC01_sorted.bam
   – bcftools mpileup -Ou -f Random_Refs_USUV.fasta BC01_sorted.bam | bcftools call -mv -Oz -o BC01.vcf.gz
   – bcftools index BC01.vcf.gz
   – cat Random_Refs_USUV.fasta | bcftools consensus BC01.vcf.gz > BC01_consensus.fasta
4. For reference-based alignments it is essential that a closely related reference sequence is used. Therefore, perform a BlastN search with the generated consensus sequence to identify the closest reference strain. After that, repeat the reference-based alignment with the closest reference strain as reference (step 3.3 and 3.4). On the command line enter:
   – minimap2 -ax map-ont Ref_USUV_BC01.fasta BC01_trimmed.fastq > BC01_ref.bam
   – samtools sort BC01_ref.bam > BC01_sorted_ref.bam
   – bcftools mpileup -Ou -f Ref_USUV_BC01.fasta BC01_sorted_ref.bam | bcftools call -mv -Oz -o BC01_ref.vcf.gz
   – bcftools index BC01_ref.vcf.gz
   – cat Ref_USUV_BC01.fasta | bcftools consensus BC01_ref.vcf.gz > BC01_ref_consensus.fasta

4. Analysis of the Illumina data

1. These sequences are automatically demultiplexed after sequencing. Reads can be quality controlled using trimmomatic²⁹. For paired-end Illumina sequences, use the commonly used cut-off median PHRED score of 33 and a minimal read length of 75 to get accurate, high quality reads. On the command line enter:
   – trimmomatic PE -phred33 9_S9_L001_R1_001.fastq.gz 9_S9_L001_R2_001.fastq.gz 9_1P.fastq 9_1U.fastq 9_2P.fastq 9_2U.fastq LEADING:3 TRAILING:3 SLIDINGWINDOW:3:15 MINLEN:75
2. Remove primers (indicated in the file Primers_Usutu.fasta in both orientations), since they can introduce artificial biases, using cutadapt²⁶. In addition, remove sequences with a length shorter than 75 nucleotides using the commands below. On the command line enter:
   – cutadapt -b o 9_1P_trimmed.fastq -p 9_2P_trimmed.fastq 9_1P.fastq 9_2P.fastq -m 75
3. Before de novo assembly, the sequence reads can be normalized for an even coverage across the genome. This is essential since de novo assemblers like SPAdes take the read coverage into account when assembling sequence reads. Normalize reads to a read coverage of 50 using BBNorm from the BBMap package³⁰. On the command line enter:
   – bbmap/bbnorm.sh target=50 in=9_1P_trimmed.fastq in2=9_2P_trimmed.fastq out=Sample9_FW_norm.fastq out2=Sample9_RE_norm.fastq
4. The normalized reads are de novo assembled using SPAdes³¹. Default settings are used for the assembly using all different kmers (21, 33, 55, 77, 99 and 127). On the command line enter:
   – spades.py -k 21,33,55,77,99,127 -o Sample9 -1 Sample9.qc.f.fq -2 Sample9.qc.r.fq
5. Map the QC reads against the obtained consensus sequence using minimap2 and programs like Geneious, Bioedit or Ugene to curate the alignment. It is important to check the beginning and the end of the contig.
   1. Align the QC reads against the obtained consensus sequencing using minimap2.
   2. Import the alignment in Geneious/Bioedit/UGene.
   3. Manually inspect, correct and curate especially the beginning and the end of the genome.

5. Determining the required read coverage to compensate for the error profile in nanopore sequencing using Illumina data as gold standard

1. Select sequence reads mapping to one amplicon, in this case amplicon 26. Subsequently, map the nanopore reads against this amplicon using minimap2. Use Samtools to select only the reads mapping to amplicon 26 and to convert the bam file into fastq. On the command line enter:
   – minimap2 -ax map-ont -m 150 Amplicon26.fasta BC01_trimmed.fastq > BC01.bam
   – samtools view -b -F 4 BC01.bam > BC01_mapped.bam
   – samtools bam2fq BC01_mapped.bam | seqtk seq – -> BC01_mapped.fastq
2. Randomly select subsets of for instance 200 sequence reads one thousand times. For example, changing it to 10 will result in the random selection of one thousand times a subset of 10 sequence reads. The script is provided as Supplementary File 1. On the command line enter:
   – python Random_selection.py
3. All randomly selected sequence reads are aligned to amplicon 26. Use KMA³² to map the sequence reads and to immediately generate a consensus sequence. Use optimized settings for nanopore sequencing, indicated by the -bcNano flag. On the command line enter:
   – kma index -i Amplicon26.fasta
   – for file in random_sample*; do
   – sampleID=${file%.fastq}
   – kma -i ${sampleID}.fastq -o ${sampleID} -t_db Amplicon26.fasta -mem_mode -mp 5 -mrs 0.0 -bcNano
   – done
4. Inspect the generated consensus sequences on the command line using:
   – cat *.fsa > All_genomes.fsa
   – minimap2 -ax map-ont Amplicon26.fasta All_genomes.fsa > All_genomes.bam
   – samtools sort All_genomes.bam > All_genomes_sorted.bam
   – samtools stats All_genomes_sorted.bam > stats.txt
   1. The error rate is displayed in the stats.txt under the heading error rate #mismatches / bases mapped. Display it on the screen with the following command:
      – grep ^SN stats.txt | cut -f 2-
   2. The amount of indels is displayed under the heading #Indels per cycle. Display it on the screen with the following command:
      – grep ^IC stats.txt | cut -f 2-