Background

Viral quasispecies (VQS) refers to a genetically closely related population of viral variants. Measurement of VQS diversity, especially in highly mutating RNA viruses, is vital for understanding their micro-evolution. Long-read sequencing facilitates retrieving full-length sequences of specific genes or genomes from each variant, yet challenges like sequencing errors and low coverage depth hinder diversity metric computation and comparison. Here, we present methodologies to overcome these challenges, enabling comparisons of VQS diversity.

Individual sample preparation

Sequencing error/noise minimization and down-sampling

Sequencing errors create plenty of artificial single nucleotide variants (SNVs), leading to an overestimated number of singleton haplotypes (haplotype = group of identical reads, singleton haplotype = haplotype with only a single read). Our package aims to minimize potential sequencing errors and potential noises from rare mutations along long reads defined as SNV(s) or nucleotide base(s) with a frequency less than n % at each nucleotide position. Such problematic base(s) will be replaced with the consensus or majority base of that particular position.

First, load the “longreadvqs” package and a FASTA file of example read alignment (sample1.fasta). To choose the % cut-off for noise minimization, use the “pctopt” function that shows % singleton haplotypes in the read alignment across different % cut-offs. Then, use the “vqssub” function for noise minimization, in this case, we choose a 10% cut-off that decreases singleton haplotypes to less than 10% of all reads. This function will give a summary of VQS diversity metrics calculated by the “QSutils” package (Gregori, et al., 2016) after noise minimization.

library(longreadvqs)
library(ggplot2)

#Load sample 1.
sample1 <- system.file("extdata", "sample1.fasta", package = "longreadvqs")
#Check which % cut-off could effectively minimize noises by assessing % singleton haplotypes.
x <- pctopt(sample1, pctsing = 0, label = "sample1")
ggplot(x, aes(x=pct, y=pctsingleton)) + geom_line() + geom_point()

#VQS diversity metrics of noise minimized (10% cut-off) read alignment
vqssub(sample1, pct = 10, label = "sample1")
#>     label  method samplingwhen pct fulldepth depth haplotypes nsingleton
#> 1 sample1 conbase        after  10       311   311         47         17
#>   pctsingleton polymorph mutations  shannon norm_shannon gini_simpson      FAD
#> 1     5.466238         6       116 2.917347    0.7577235    0.8819625 6.475025
#>           Mfe         Pie         Mfm         Pim
#> 1 0.002580953 0.002994924 6.79165e-09 0.001803544

To compare VQS diversity across samples, it’s important to standardize the read alignment depth. Since sequencing results can vary, we perform down-sampling after minimizing noises with the “vqssub” function. In our case, for sample1, sample2, and sample3 with original depths of 311, 316, and 655, respectively, we randomly down-sample them to a common depth of 300.

#Load samples 2 and 3.
sample2 <- system.file("extdata", "sample2.fasta", package = "longreadvqs")
sample3 <- system.file("extdata", "sample3.fasta", package = "longreadvqs")

#Noise minimization (10% cut-off) and down-sampling (depth of 300)
a <- vqssub(sample1, pct = 10, samsize = 300, label = "sample1")
b <- vqssub(sample2, pct = 10, samsize = 300, label = "sample2")
c <- vqssub(sample3, pct = 10, samsize = 300, label = "sample3")

#Compare VQS diversity across three samples after Noise minimization and down-sampling.
rbind(a, b, c)
#>     label  method samplingwhen pct fulldepth depth haplotypes nsingleton
#> 1 sample1 conbase        after  10       311   300         47         17
#> 2 sample2 conbase        after  10       316   300         57         22
#> 3 sample3 conbase        after  10       655   300         32          8
#>   pctsingleton polymorph mutations  shannon norm_shannon gini_simpson       FAD
#> 1     5.666667         6       116 2.942932    0.7643685    0.8844816  6.475025
#> 2     7.333333         7       138 3.188717    0.7886908    0.9053735 10.214067
#> 3     2.666667         5        67 2.606215    0.7519947    0.8698997  2.530071
#>           Mfe         Pie          Mfm         Pim
#> 1 0.002580953 0.002994924 7.190260e-09 0.001836544
#> 2 0.002467944 0.003199896 6.170471e-09 0.001951459
#> 3 0.002325433 0.002550475 6.186794e-09 0.001566417

We can also repeatedly down-sample the alignments and check the consistency of each VQS diversity metric at the different depths like the following example.

#Randomly down-sample from depth of 655 to 300, and 100 (10 iterations for each sample size).
set.seed(123)
c_full <- vqssub(sample3, pct = 10, label = "full_depth") #non-sampled alignment
c_300 <- vqsresub(sample3, iter = 10, pct = 10, 
                  samsize = 300, label = "depth_300") #down-sample to depth of 300
c_100 <- vqsresub(sample3, iter = 10, pct = 10,
                  samsize = 100, label = "depth_100") #down-sample to depth of 100

all_c <- rbind(c_full, c_300, c_100) #combine data

#Load packages for visualization.
library(cowplot)

all_c$depth <- as.character(all_c$depth)
depthorder <- c("655", "300", "100")

#Plot box-plots of three VQS metrics variation through different sampling depths.
haplotypes <- ggplot(all_c, aes(x = factor(depth, level=depthorder), y = haplotypes, 
    group = interaction(depth, label), fill = label)) + geom_boxplot() + 
  xlab("depth") + theme(legend.position = "none")
shannon <- ggplot(all_c, aes(x = factor(depth, level=depthorder), y = shannon, 
    group = interaction(depth, label), fill = label)) + geom_boxplot() + 
  xlab("depth") + theme(legend.position = "none")
gini_simpson <- ggplot(all_c, aes(x = factor(depth, level=depthorder), y = gini_simpson, 
    group = interaction(depth, label), fill = label)) + geom_boxplot() + 
  xlab("depth") + theme(legend.position = "none")
plot_grid(haplotypes, shannon, gini_simpson, nrow = 1)

From this plot, the number of haplotypes obviously decreases as the depth decreases, while Shannon entropy and the Gini-Simpson index remain relatively stable with high variation at a depth of 100.

VQS data preparation for in-depth comparison

The previous “vqssub” function just summarizes key VQS diversity metrics after noise minimization and down-sampling. However, these metrics do not provide us with comprehensive data for further comparative analyses, such as SNV positions and full read alignment after down-sampling. The “vqsassess” function will generate a complete VQS profile for each sample, and the outputs can be used for other functions.

In this example, we repeat noise minimization and down-sampling for samples 1-3 but apply the “vqsassess” function instead of the “vqssub” function. We then compare the SNV profiles between these samples using the “snvcompare” function.

#Generate complete VQS data for further comparison for each sample.
s1 <- vqsassess(sample1, pct = 10, samsize = 300, label = "sample1")
s2 <- vqsassess(sample2, pct = 10, samsize = 300, label = "sample2")
s3 <- vqsassess(sample3, pct = 10, samsize = 300, label = "sample3")
#Compare SNV profile between three samples.
snvcompare(samplelist = list(s1, s2, s3), ncol = 1)

This plot shows that SNV positions are well-aligned between three samples. In some cases, SNV plot can highlight sequencing error concentrating at a particular range in read alignment like the additional sample added to the following example.

#Load sample 4.
sample4 <- system.file("extdata", "mock.fasta", package = "longreadvqs")
s4 <- vqsassess(sample4, pct = 10, samsize = 300, label = "sample4")
#SNV profile of sample 4
s4$snv

The plot above captures atypical SNVs towards the 3’ end of the sample 4 alignment, which is potentially caused by sequencing error. We can address this issue using the “vqscustompct” function, which allows us to adjust the % cut-off for noise minimization at specific range(s) in the alignment. In the case of sample 4, we increase the % cut-off from position 972 onward, raising it from 10% to 30%.

#Use "vqscustompct" function to increase % cut-off after position 972 to 30%.
s4_fix <- vqscustompct(sample4, pct = 10, brkpos = c("1:971","972:982"), lspct = c(10,30), 
                       samsize = 300, label = "sample4")
#SNV profile of sample 4 after % cut-off adjustment
s4_fix$snv

Comparing VQS diversity across multiple samples

Once we have all samples prepared, we can combine them to not only compare within-sample diversity but also explore between-sample relationships. Using the “vqscompare” function, noise-minimized down-sampled reads from all samples are pooled to identify 1) identical haplotypes and 2) genetically closely related haplotypes between samples. In this example, we combine the previously prepared outputs of samples 1 to 4 into the “vqscompare” function, allowing it to generate a summary comparative plot.

#List outputs from "vqsassess" or "vqscustompct" that we want to compare into "vqscompare" function.
#Set the number of new OTU groups based on k-means clustering to 10 groups (kmeans.n = 10).
set.seed(1234)
comp <- vqscompare(samplelist = list(s1, s2, s3, s4_fix),
                   lab_name = "Sample", kmeans.n = 10)
#The most important output of the "vqscompare" function is the summary plot.
comp$summaryplot