ViralFlow: An Automated Tool for Assembling COVID-19 Genomes, Identifying Variants, and Tracking Virus Changes Within Patients

The workflow was developed within an Ubuntu 20.04.2 LTS Docker environment (https://hub.docker.com/_/ubuntu↗, accessed on 16 November 2021) and is composed of six steps used to analyze SARS-CoV-2 Illumina paired-end amplicon sequencing data (Figure 1A): reference genome indexing, quality control, consensus generation, intrahost variant analysis, virus lineage assignment and mutation analysis, and assembly metrics analytics. This workflow can be used in different computational environments (Figure 1B).

First, the reference genome is processed and indexes are obtained using the BWA index [23]; we recommend the SARS-CoV-2 reference genome Wuhan-Hu-1, NCBI refseq NC_045512.2 code. In this step, the pangolin tool is updated. The quality control step for the raw sequencing reads is performed with the fastp v.0.20.1 tool [24], where the paired-end reads are trimmed using a minimum read quality threshold (Phred score = 20). The adapters or the PCR primers and the minimum length threshold for the trimmed reads should be defined by the user. In addition to the paired-end treated data, the fastp tool generates an html file information from pre- and post-treatment steps with associated statistics (Figure 1C).

The generation of the consensus genome is performed using a reference-guided assembly strategy. In this step, the paired-end libraries are mapped against a reference genome with BWA v.0.7.17 [23]. Following the alignment step, the Samtools [25] sort and index parameters are used to sort and index the BAM files. Next, minor variant analyses are performed, the Samtools v.1.9 and iVar v.1.3.1 [26] tools are used for the correct recovery of SNPs and indels, and two consensus are generated: one with the majority allele present in every nucleotide position along the genome (iVar consensus -t 0) and another version with ambiguous nucleotide characters, in cases where the majority intrahost Single Nucleotide Variants (iSNVs) encompass up to 60% of allele frequencies (iVar consensus -t 0.6). Only mapped bases with quality equal to or greater than 30 (-q 30) were used in iVar counts. The minimum depth threshold to consider a position with supported intrahost variants can be defined by the user (default equals 100x).

An extra step for intrahost variant calls is necessary, considering that iVar does not provide an option to generate consensus harboring all iSNVs with two or more alleles found in low frequency (≤49% of the reads). Given that the consensus genome with minor iSNVs is essential to understanding the effect of intrahost variants, we developed an in-house python script (intrahost_script.py) that uses the allele frequencies per position output of bamreadcount v.0.8.0 [27] to detect only positions with two or more alleles and to generate a consensus harboring all minor supported alleles. To avoid the recognition of sequencing artifacts as intrahost variants, genomic positions were selected:(i)The minor allele frequency represented at least 5% of the total allele depth;(ii)The minor alleles had at least 100 reads of depth (default depth);(iii)The minor allele nucleotides were supported by reads of both senses (at least 5% of depth should come from each read sense).

Combining the first two requirements stated above, a sequenced depth of 2000 reads is required to detect iSNVs present at the minimum frequency of 5%. However, the user can set a different minimum depth threshold, if necessary.

The virus lineages signature is performed with Pangolin. Pangolin and all information about current and new lineages are updated at the moment of the Docker or Singularity image creation to avoid using outdated data and software versions for the analysis. When using interactive containers, the command "pangolin--update" is strongly recommended. The consensus quality and set of mutations are evaluated using nextclade v.0.14.2 [28]. If the analyzed sample shows intrahost signals, Pangolin and nextclade analyses are performed for both consensus versions (with major and minor allele frequencies), while it will run only for a single consensus genome in the absence of iSNV sites. In the last step, the assembly metrics, such as depth and coverage, are extracted with bamdst v.1.0.6 [29].

Infrastructure and computational experience are heterogeneous in the different research groups working with SARS-CoV-2 genomic data; therefore, we evaluated our workflow in two use case scenarios: Case I: Using an average personal computer to install all dependencies or using Docker or Singularity container services. Case II: Using a multi-user computational server to install all dependencies or using Docker or Singularity container services.

Case I was run on a personal laptop with the following configurations: Ubuntu 20.04.2 LTS, 02 × RAM 8 Gb DDR4 2667 MHz and CPU AMD^® Ryzen™ 7-3750H 2.88 GHz. For Case II, we used a computational server with the following configurations: a node with 191 Gb of RAM DRAM 2933 MHz and two CPUs Intel(R) Xeon(R) Gold 5220R CPU @ 2.20 GHz totaling 96 threads. The scalability of ViralFlow according to the number of threads provided was evaluated on both computational resource scenarios using a Singularity container. The performance of the workflow was accessed using two datasets. The first one is a public dataset of 86 Brazilian SARS-CoV-2 Illumina paired-end libraries generated by the amplicon sequencing method using the Illumina COVIDSeq protocol, available under the EMBL-EBI study accession PRJEB47823. It was used to evaluate the computational resources and the run time of the workflow. The second one is an artificial dataset (Supplementary File S1) created with the ART [30] of five paired-end libraries simulating a simultaneous infection (coinfection/codetection) of different SARS-CoV-2 lineages in a single sample (sample information in Table S1), used to evaluate the capacity of the workflow to detect intrahost variants.

To compare the genome assembly coverage breadth, depth and lineage assignment of SARS-CoV-2 lineages, we assembled the same 86 samples with HAVoC [31], a published workflow that performs similar steps of read processing (fastp) and mapping (BWA). We performed two tests with HAVoC:

Case I: fastp with -q (qualified_quality_phred) equal to 20 and parsing the adapters sequence file, to perform the same treatment of ViralFlow

Case II: fastp with -q equal to 15 and without an adapters sequence file, to run the native HAVoC.

Both tests were performed considering a minimum of 5× coverage depth to generate consensus and 75 as the minimum read length threshold in the fastp step.

A set of 86 samples was generated in an amplicon-sequencing strategy using the Illumina COVIDSeq Test kit, generating paired-end reads of 150 nucleotides that encompass PCR positive samples for SARS-CoV-2 from the Pernambuco state in Brazil from August 2020 to May 2021 (Table S2). It was assembled and tested in both environments (see the Materials and Methods section, Case I and Case II). The benchmarks for the 86 samples show that it is possible to process this dataset in ~20 min using ≥6 threads (Figure 2A) on a personal computer. Such time can be reduced to ~2 min or less on a dedicated server using ≥50 threads (Figure 2B). Benchmarks also revealed that using more than one thread per sample generally decreases performance. Most of the ViralFlow steps are handled by a single thread, and additional threads imply extra waiting time for those processes to finish. As a general rule of thumb, one thread per sample should scale better in almost every scenario. Therefore, scaling up using a computer/server with a large number of cores/threads is ideal for speeding up ViralFlow results. The total RAM (Random Access Memory) used was ~0.70 Gb and ~0.80 Gb in Case I and Case II, respectively. This number was similar and did not change significantly depending on the number of threads. We observed that the tools used in our workflow relied more on CPU usage than RAM. ViralFlow scales well with the number of threads available and is able to process life-size sequence batches in a few minutes, even using a modest computational infrastructure.

We detected a low number of iSNVs, from zero to two, with a mean of zero (stdev = 0.43; for details, see Table S2), among the 86 "non-artificial" samples, supporting published estimates of low intrahost variant variability of SARS-CoV-2 [19,32]. To evaluate the performance of the workflow on samples with a large number of iSNVs, we generated five artificial coinfection/codetection samples (see the Section 2). The workflow was able to consistently detect a large amount of well-supported iSNVs (47 iSNVs per sample) (Figure 3A, Table S3). These results show the capacity of the workflow to rapidly detect and generate a range of useful information that is important to generate new insights, such as single consensus and coinfection of different SARS-CoV-2 lineages in a single sample. Moreover, the intrahost multi-allele frequencies can also be used to detect sample contamination in a scenario where most samples show the same intrahost pattern found in a negative control sample.

Two key pieces of information are necessary for the deployment of SARS-CoV-2 outbreak control strategies: the virus lineage assignment and the mutation characterization [9,33]. The workflow generates two tabular files for each sample containing this key information: '.pango.csv' and '.nextclade.csv', which provide the virus lineage and the mutations found, respectively. The lineages identified in the 86 non-artificial samples (Figure 4) correspond to the set of mutations and to the expected lineages circulating at the sampling date in each location (Table S4). For the five artificial samples, the virus lineages and set of mutations correspond to the allele frequencies present in iSNV multi-allele frequencies (Figure 3B, Table S3). To show the precision of our workflow in detecting indel regions, five random non-artificial samples were assigned as P.1 with the deletion of 11288–11297 in ORF1a and an insertion into the intergenic region at 28,262 positions, which were manually investigated using Integrative Genomics Viewer [34] (Supplementary File S2).

In addition to the intrahost, lineage and mutations information generated by different tools present in our workflow, the '.fastp.html' generated by fastp and the 'coverage.report' and 'chromosomes.report' files generated by bamdst can be used to assess the quality of mapping and assembly steps and for feedback to wet lab staff. This information can also be crossed with the information of the 'qc.overallScore' and 'qc.overallStatus' columns present in the '.nextclade.csv' table.

The comparison with HAVoC showed a similar genome coverage () between ViralFlow (coverage mean 99.70 with stdev equal to 3.24), HAVoC Test I (coverage mean 99.71 with stdev equal to 2.9) and HAVoC Test II (coverage mean 99.71 with stdev equal to 2.71). The mean depth between ViralFlow and HAVoC Test I are similar (342.2 (stdev equal to 66.47) and 341.54 (stdev equal to 67.42), respectively), and HAVoC Test II showed a high mean depth (404.19 with stdev equal to 80.76). The higher depth of HAVoC Test II can be explained by the lower quality threshold than those used in ViralFlow and HAVoC Test I. Table S5

When investigating the consensus genomes, we noticed that HAVoC fails to assemble a deletion region of nine nucleotides in ORF1ab of Gamma (P.1) genomes. In the 34 samples with this deletion, the HAVoC Test I fails to correctly assemble all 34 samples bearing this deletion (), and the HAVoC Test II correctly assembled only 12 of the 34 samples (). Moreover, we found 23 samples with pangolin lineage incongruencies between ViralFlow and HAVoC Test I (), where HAVoC fails to recover a mutation in position 22,812 of the SARS-CoV-2 genome, probably owing to the presence of minor iSNVs in these genomic loci that support the reference base. In these cases, we performed a manual curation to replace the iSNV with a minor frequency to an iSNV with a major frequency. After the curation, the pangolin signed the correct lineage (). Supplementary File S3A Supplementary File S3B Table S5 Supplementary File S3C and Table S5