Estimating multiple infections from variant data obtained using massively parallel sequencing experiments can be done in several ways using moimix.

The core function in moimix is binomix which is an interface to a binomial mixture model estimation procedure using the R package flexmix. The input to this function is a matrix of read counts supporting the reference and alternate alleles for all the variants in data.

In this vignette, we will introduce some key functions in moimix and how to set up your variant data to estimate some measures of multiplicity of infection.

Key functions in `moimix`

Function	Description
binommix	estimate MOI using binomial mixture model
alleleCounts	extract a matrix of read counts
getBAF	estimate B-allele frequencies from read counts
plotBAF	plot the B-allele frequency spectra genome-wide
getMAF	esimate minor allele frequencies from read counts
getFws	estimate within isolate diversity according to Manske et. al, 2012
callMajor	call major haplotypes using B-allele frequencies
extractCOIL	Output a COIL barcode file
extractPED	Create haploid PLINK PED file

Data Formats

moimix assumes that variants have been called either from whole-genome/exome or RNA sequencing and the variant calling algorithm used outputs a variant call format (VCF) version 4.1 file. We also assume that each sample in the VCF file has been annotated with a AD tag at the variant sites for both the alternate and reference alleles. Currently, we have tested the functions using output from GATK UnifiedGenotyper and HaplotypeCaller algorithms.

moimix requires that the VCF file is then converted to the Genomic Data Storage (GDS) format. The GDS format provides a fast and convienent way of manipulating large VCF files in R. We make use of the SeqArray package for both VCF conversion and manipulating GDS files.

If you are just interested in running the binomix procedure, then you do not need to have a GDS file. For each sample construct a matrix with 2 columns correpsonding the read counts for the alternate and referene alleles. Then pass this matrix to the binomix function.

Converting a VCF to GDS format

To convert a tabixed gzipped VCF file to the GDS format use the SeqArray package.

library(SeqArray)
seqVCF2GDS("/path/to/my/file.vcf.gz", "/path/to/my/file.gds")

To read the GDS file into R, use the seqOpen function and to check that its output matches the orginal VCF file use seqSummary.

my_vcf <-seqOpen("/path/to/my/file.gds")
seqSummary(my_vcf)

Case study: Plasmodium falciparum isolates from The Gambia.

Here we present an example analysis from data obtained from the Malaria Genomics Consoritum community project. The data consists of 3 Plasmodium falciparum isolates from The Gambia with approximately 60,000 high-quality biallelic SNVs on 14 nuclear chromosomes and the apicoplast.

For details on how the data was processed please see the README file located at ftp://ngs.sanger.ac.uk/production/pf3k/release_4/

First, we use SeqArray to load our GDS file and look at it.

library(SeqArray)

## Loading required package: gdsfmt

gds_file <- system.file("extdata", "gambian_isolates.gds", package = "moimix")
isolates <- seqOpen(gds_file)



seqSummary(isolates)

## File: /tmp/RtmpabFGac/Rinst15217151ed10/moimix/extdata/gambian_isolates.gds
## Format Version: v1.0
## Reference: file:///lustre/scratch109/malaria/pf3k_methods/resources/Pfalciparum.genome.fasta
## Ploidy: 2
## Number of samples: 3
## Number of variants: 53,379
## Chromosomes:
##     chrPf3D7_01_v3: 1471, chrPf3D7_02_v3: 1958, chrPf3D7_03_v3: 2219, chrPf3D7_04_v3: 3528, chrPf3D7_05_v3: 3142, chrPf3D7_06_v3: 2719
##     chrPf3D7_07_v3: 3751, chrPf3D7_08_v3: 3303, chrPf3D7_09_v3: 3229, chrPf3D7_10_v3: 3915, chrPf3D7_11_v3: 4654, chrPf3D7_12_v3: 4662
##     chrPf3D7_13_v3: 6717, chrPf3D7_14_v3: 8111
## Contigs:
##     Pf3D7_01_v3, 640851
##     Pf3D7_02_v3, 947102
##     Pf3D7_03_v3, 1067971
##     Pf3D7_04_v3, 1200490
##     Pf3D7_05_v3, 1343557
##     Pf3D7_06_v3, 1418242
##     Pf3D7_07_v3, 1445207
##     Pf3D7_08_v3, 1472805
##     Pf3D7_09_v3, 1541735
##     Pf3D7_10_v3, 1687656
##     Pf3D7_11_v3, 2038340
##     Pf3D7_12_v3, 2271494
##     Pf3D7_13_v3, 2925236
##     Pf3D7_14_v3, 3291936
##     Pf3D7_API_v3, 34250
##     Pf_M76611, 5967
## Alleles:
##     NON_REF, Represents any possible alternative allele at this location
##     tabulation: 2, 53379(100.0%)
## Annotation, Quality:
##     Min: 30.2099990844727, 1st Qu: 8685.1298828125, Median: 34331.5, Mean: 620170.906236837, 3rd Qu: 336185, Max: 13424800
## Annotation, FILTER:
##     PASS, All filters passed, 53379(100.0%)
##     Centromere, RegionType == 'Centromere', 0(0.0%)
##     InternalHypervariable, RegionType == 'InternalHypervariable', 0(0.0%)
##     LowQual, Low quality, 0(0.0%)
##     Low_VQSLOD, VQSLOD <= 0.0, 0(0.0%)
##     SubtelomericHypervariable, RegionType == 'SubtelomericHypervariable', 0(0.0%)
##     SubtelomericRepeat, RegionType == 'SubtelomericRepeat', 0(0.0%)
##     VQSRTrancheINDEL99.00to99.90, Truth sensitivity tranche level for INDEL model at VQS Lod: -11.7168 <= x < -0.61, 0(0.0%)
##     VQSRTrancheINDEL99.90to100.00+, Truth sensitivity tranche level for INDEL model at VQS Lod < -24316.1745, 0(0.0%)
##     VQSRTrancheINDEL99.90to100.00, Truth sensitivity tranche level for INDEL model at VQS Lod: -24316.1745 <= x < -11.7168, 0(0.0%)
##     VQSRTrancheSNP99.90to100.00+, Truth sensitivity tranche level for SNP model at VQS Lod < -38542.49, 0(0.0%)
##     VQSRTrancheSNP99.90to100.00, Truth sensitivity tranche level for SNP model at VQS Lod: -38542.49 <= x < -4.5767, 0(0.0%)
## Annotation, INFO variable(s):
##     AF, A, Float, Allele Frequency, for each ALT allele, in the same order as listed
##     FS, 1, Float, Phred-scaled p-value using Fisher's exact test to detect strand bias
##     GC, 1, Float, GC content around the variant (see docs for window size details)
##     HaplotypeScore, 1, Float, Consistency of the site with at most two segregating haplotypes
##     MQ, 1, Float, RMS Mapping Quality
##     MQRankSum, 1, Float, Z-score From Wilcoxon rank sum test of Alt vs. Ref read mapping qualities
##     QD, 1, Float, Variant Confidence/Quality by Depth
##     ReadPosRankSum, 1, Float, Z-score from Wilcoxon rank sum test of Alt vs. Ref read position bias
##     SOR, 1, Float, Symmetric Odds Ratio of 2x2 contingency table to detect strand bias
##     VariantType, 1, String, Variant type description
## Annotation, FORMAT variable(s):
##     GT, 1, String, Genotype
##     AD, ., Integer, Allelic depths for the ref and alt alleles in the order listed
##     DP, 1, Integer, Approximate read depth (reads with MQ=255 or with bad mates are filtered)
##     GQ, 1, Integer, Genotype Quality
## Annotation, sample variable(s):
##     <None>

# save sample identifiers
sample.id <- seqGetData(isolates, "sample.id")

We can look at the number of variants on each of the main nuclear chromosomes.

library(moimix)

## Loading required package: flexmix

# get genomic coordinates of all variants
coords <- getCoordinates(isolates)
head(coords)

##    chromosome position variant.id
## 1 Pf3D7_02_v3   106099          1
## 2 Pf3D7_02_v3   106230          2
## 3 Pf3D7_02_v3   108524          3
## 4 Pf3D7_02_v3   109059          5
## 5 Pf3D7_02_v3   109073          6
## 6 Pf3D7_02_v3   109248          7

table(coords$chromosome)

## 
## Pf3D7_01_v3 Pf3D7_02_v3 Pf3D7_03_v3 Pf3D7_04_v3 Pf3D7_05_v3 Pf3D7_06_v3 
##        1471        1958        2219        3528        3142        2719 
## Pf3D7_07_v3 Pf3D7_08_v3 Pf3D7_09_v3 Pf3D7_10_v3 Pf3D7_11_v3 Pf3D7_12_v3 
##        3751        3303        3229        3915        4654        4662 
## Pf3D7_13_v3 Pf3D7_14_v3 
##        6717        8111

Furthermore, we can extract high quality (‘true positive’) SNPs using the seqGetData function and obtaining filters based on the GATK best practices. Here we look at the distributions of the QD and SOR tags in the original VCF file.

# extract tags from GATK Haplotype Caller using convienence funciton
qd <- seqGetData(isolates, "annotation/info/QD")
sor <- seqGetData(isolates, "annotation/info/SOR") 

hist(qd)

hist(sor)

Estimating the BAF matrix

The first step in our clonality analysis is to estiamte the B-allele frequencies for each isolate directly from the sequencing depth. We construct a bafMatrix object which stores the estimates of the B-allele frequencies in a matrix over the entire cohort and the genomic coordinates of the variants.

isolate_baf <- bafMatrix(isolates)

## # of selected variants: 53,379
## # of selected samples: 3
## # of selected variants: 53,379

class(isolate_baf)

## [1] "bafMatrix"

str(isolate_baf)

## List of 3
##  $ baf_matrix: num [1:3, 1:53379] 0 0 0 0 1 0 0 0 1 0 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ sample : chr [1:3] "PA0008-C" "PA0021-C" "PA0022-C"
##   .. ..$ variant: chr [1:53379] "1" "2" "3" "5" ...
##  $ baf_site  : num [1:53379] 0 0.9308 0.0128 0 0.6336 ...
##  $ coords    :'data.frame':  53379 obs. of  3 variables:
##   ..$ chromosome: chr [1:53379] "Pf3D7_02_v3" "Pf3D7_02_v3" "Pf3D7_02_v3" "Pf3D7_02_v3" ...
##   ..$ position  : int [1:53379] 106099 106230 108524 109059 109073 109248 109343 109601 109623 109634 ...
##   ..$ variant.id: int [1:53379] 1 2 3 5 6 7 8 9 10 11 ...
##  - attr(*, "class")= chr "bafMatrix"

You can plot the genome-wide B-allele frequencies by passing the bafMatrix object and a sample identifer to plot.

Here’s an example for an isolate that is from a single-clone infection.

plot(isolate_baf, sample.id = "PA0022-C")

And an example of multiple-clone infection.

plot(isolate_baf,  sample.id = "PA0021-C")

And a slighly more complicated example.

plot(isolate_baf,  sample.id = "PA0008-C")

## Estimating MOI with binommix

To estimate MOI we use the binomix function. MOI is estimated on a per-isolate basis using a binomial mixture model implemented in the R package flexmix. In order to fit a model for an isolate, we need a matrix of read counts and to set the value of \(k=2c\) where \(c\) is the number clonal infections.

We fit a mixture model assuming there are 2 infections on “PA0021-C” and report the results.

set.seed(2002)
# first generate matrix of read counts for ref and alt alleles
counts <- alleleCounts(isolates)

# fit a model on sample PA0021-C, multiple clone infection
m1 <- binommix(counts, sample.id = "PA0021-C", k = 2)

## 2 : * * * * *

# to inspect the model and view parameter estimates
param_estimates <- getTheta(m1$fits)

The estimated clonal proportions under the binomial mixture model are:

param_estimates

##      pi.hat    mu.hat
## 1 0.4411062 0.2615370
## 2 0.5588938 0.7393493

We can visualise the model output as follows:

plot(isolate_baf, "PA0021-C")
abline(h = param_estimates$mu.hat, col = "blue", lty = 2)

Or by plotting a histogram.

bf <- isolate_baf$baf_matrix["PA0021-C",]
hist(bf, breaks = 20, main = "MOI = 2 infection", xlab = "SNV frequency")
abline(v = param_estimates$mu.hat, col = "blue", lty = 2)

We can also examine the more complicated mixed infection, which appears to be a mixture of 3 or 4 clones. We can use the Bayesian Information Criterion to select the number of components.

m2 <- binommix(counts, sample.id = "PA0008-C", k = 2:4)

## 2 : * * * * *
## 3 : * * * * *
## 4 : * * * * *
## * * *

plot(isolate_baf, "PA0008-C")
abline(h = getTheta(m2$fits, criterion = "BIC")$mu.hat, col = "blue", lty = 2)

hist(m2$baf, breaks = 20, main = "Complex Infection", xlab = "SNV frequency" )
abline(v = getTheta(m2$fits, criterion = "BIC")$mu.hat, col = "blue", lty = 2)

In this case we end up with 4 components in our model - with estimated clonal proportions as (0.23,0.11,0.38,0.28).

Estimating MOI with \(F_{ws}\)

An alternative way of assessing clonality is to estimte the \(F_{ws}\) statistic. This is useful if you are only interested in identifying clonal infections for filtering in downstream analyses. An \(F_{ws} < 0.99\) is indicative of a clonal infection.

fws_all <- getFws(isolates)

fws_all

##  PA0008-C  PA0021-C  PA0022-C 
## 0.5188443 0.4526363 0.9988770

seqClose(isolates)

References

flexmix package - https://cran.r-project.org/web/packages/flexmix/index.html
SeqArray package - http://bioconductor.org/packages/release/bioc/html/SeqArray.html

Session Info

sessionInfo()

## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 26.04.1 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] moimix_0.0.2.9001 flexmix_2.3-21    SeqArray_1.53.4   gdsfmt_1.49.8    
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.1       farver_2.1.2           dplyr_1.2.1           
##  [4] Biostrings_2.81.9      fastmap_1.2.0          rpart_4.1.27          
##  [7] digest_0.6.39          lifecycle_1.0.5        survival_3.8-12       
## [10] magrittr_2.0.5         compiler_4.6.1         rlang_1.3.0           
## [13] sass_0.4.10            tools_4.6.1            yaml_2.3.12           
## [16] data.table_1.18.6.1    knitr_1.52             RColorBrewer_1.1-3    
## [19] BiocParallel_1.47.0    purrr_1.2.2            GWASExactHW_1.2       
## [22] BiocGenerics_0.59.12   sys_3.4.3              nnet_7.3-21           
## [25] grid_4.6.1             stats4_4.6.1           jomo_2.7-6            
## [28] mice_3.19.0            scales_1.4.0           iterators_1.0.14      
## [31] MASS_7.3-66            mcmc_0.9-8             cli_3.6.6             
## [34] rmarkdown_2.32         crayon_1.5.3           reformulas_0.4.4      
## [37] generics_0.1.4         otel_0.2.0             minqa_1.2.8           
## [40] cachem_1.1.0           modeltools_0.2-25      operator.tools_1.6.3.1
## [43] splines_4.6.1          parallel_4.6.1         XVector_0.53.0        
## [46] vctrs_0.7.3            boot_1.3-32            glmnet_5.1            
## [49] Matrix_1.7-6           jsonlite_2.0.0         SparseM_1.84-2        
## [52] MCMCpack_1.7-1         IRanges_2.47.5         S4Vectors_0.51.10     
## [55] mitml_0.4-5            maketools_1.3.2        foreach_1.5.2         
## [58] jquerylib_0.1.4        tidyr_1.3.2            glue_1.8.1            
## [61] nloptr_2.2.1           pan_2.0                codetools_0.2-20      
## [64] shape_1.4.6.1          GenomicRanges_1.65.4   lme4_2.0-6            
## [67] tibble_3.3.1           pillar_1.11.1          htmltools_0.5.9       
## [70] quantreg_6.1           Seqinfo_1.3.2          R6_2.6.1              
## [73] SeqVarTools_1.51.3     Rdpack_2.6.6           formula.tools_1.7.1   
## [76] evaluate_1.0.5         lattice_0.23-1         Biobase_2.73.2        
## [79] rbibutils_2.4.1        backports_1.5.1        broom_1.0.13          
## [82] bslib_0.12.0           MatrixModels_0.5-4     Rcpp_1.1.2            
## [85] coda_0.19-4.1          nlme_3.1-171           mgcv_1.9-4            
## [88] logistf_1.26.1         xfun_0.61              buildtools_1.0.0      
## [91] pkgconfig_2.0.3

Introduction to moimix package

Key functions in moimix