*2.6. Bioinformatic Analysis*

Reads obtained from Illumina were de-multiplexed and converted to fastq files. All fastq files were quality and adaptor trimmed using Trimmomatic (v0.39) [20]. Only highquality paired reads (-phred 33 -leading 3 -trailing 3- slidingWindow: 4:15) with 150 bp were used for further analysis. FASTQC tools were further used to assess whether any adaptors remained [21]. High-quality reads were then mapped to the HPV 16 reference genome from the Papillomavirus Episteme (PaVE) [22] using bwa (v0.7.17) [23], to create a sam file. Due to the circular HPV genome, the reference genome was modified by adding the 258 nucleotides from the beginning to the end of the genome sequence to not lose coverage of amplicons 46 and 47. SAMtools (v1.14) [24] was then used to convert files from sam to bam and to curate files for the variant calling. BCFtools (v1.14), mpileup and consensus tools were used for the variant calling and for the generation of a consensus sequence [25], using default parameters. Positions not covered were annotated as Ns.

New consensus files were aligned using MAFFT (v7.490) with default parameters [26]. A manual edit was performed when required. Maximum likelihood trees were inferred using RaxML (v2.0.8) [27] with the GTR substitution model (ML + transfer bootstrap expectation + consensus, 1 run, 100 reps). Visualization of the trees generated by RaxML was performed using Figtree (v1.4.4). Each sample was assigned with a sub-lineage corresponding to the nearest neighbor.

Sub-lineages references were obtained from the PAVE for each of the HPV 16 sublineages: A1 (K02718.1), A2 (AF536179.1), A3 (HQ644236.1), A4 (AF534061.1), B1 (AF536180.1), B2 (HQ644298.1), B3 (HQ644298.1), B4 (KU053914.1), C1 (AF472509.1), C2 (HQ644244.1), C3 (KU053920.1), C4 (KU053925.1), D1 (HQ644257.1), D2 (AY686579.1), D3 (AF402678.1) and D4 (AF402678.1) A sub-lineage assignment was performed for all specimens excluding those with <100× median depth or low genome coverage (<80% genome coverage).
