| Title: | Estimating mulitplicity of infection from high-throughput sequencing data |
|---|---|
| Description: | Evaluate multiplicity of infection (MOI) of organisms using binomial mixture models and other population genetics methods. moimix can be used to simulate single nucleotide variant (SNV) data from samples with known MOI. |
| Authors: | "Stuart Lee" [aut, cre] |
| Maintainer: | "Stuart Lee" <[email protected]> |
| License: | GPL-3 + file LICENSE |
| Version: | 0.0.2.9001 |
| Built: | 2026-10-09 10:21:49 UTC |
| Source: | https://github.com/bahlolab/moimix |
Return allele counts from a genofile based on read counts
alleleCounts(gdsfile)alleleCounts(gdsfile)
gdsfile |
a |
If ref.allele = NULL, the function returns a list
of matrices corresponding to the read counts of the reference and alternate
alleles respectively, while if ref.allele = 0L or ref.allele = 1L
the function returns a matrix corresponding the reference or alternate
read counts. Currently, the AD tag formats from varscan2re and GATK
are supported.
an alleleCounts matrix object, which is a list containg the
Compute B-allele frequency spectrum genome-wide
bafMatrix(gdsfile)bafMatrix(gdsfile)
gdsfile |
a |
Constructs a bafMatrix object which is a list consisting of a matrix of estimated B-allele frequencies and
a bafMatrix object
Generate GenomicRanges object from GDS file
barcodes(publication)barcodes(publication)
publication |
a character string indicating malaria publication (default Volkmann2008) |
Return a GenomicRanges object containing Pf barcode information from Volkmann, 2008.
Fit binomial mixture model on coverage data
binommix(counts_matrix, sample.id, k, coverage_threshold = 0L, niter = 1000)binommix(counts_matrix, sample.id, k, coverage_threshold = 0L, niter = 1000)
counts_matrix |
an |
sample.id |
character sample.id to fit model on. |
k |
vector of mixture components to fit |
coverage_threshold |
exclude sites with total coverage < coverage_threshold |
niter |
number of iterations to run |
CDF for binomial mixture model
binommixCDF(model, x, N)binommixCDF(model, x, N)
model |
flexmix object |
x |
vector of read counts supporting each SNV |
N |
vector of coverage at SNV |
Vector of theoretical quantiles
Extract major alleles from a GDS file using the non-reference allele frequency (NRAF) estimated from read counts directly. In the case of RAF = 0.5 we select the major allele randomly from Bernoulli distribution if use.hets is FALSE. Other wise hets are coded as 3. If get.nucleotides is TRUE a character matrix is retruned with the alleles coded, otherwise if use.hets the allele is coded as N. Missing genotype calls are coded as X.
callMajor(gdsfile, get.nucleotides = FALSE, use.hets = FALSE)callMajor(gdsfile, get.nucleotides = FALSE, use.hets = FALSE)
gdsfile |
a |
get.nucleotides |
FALSE logical indicating whether to replace allele with nucleotide bases |
use.hets |
FALSE logical indicating whether to code hets sepeartely |
a binary matrix or character matrix containing reference and alternate alleles
Generate Confidence Intervals for parameter estimates using estimated Fisher Information Matrix
ciMM(y, model, alpha = 0.05)ciMM(y, model, alpha = 0.05)
y |
filtered counts matrix |
model |
flexmix object |
alpha |
confidence level (default = 0.05) |
a data.frame containing parameter estimates and lower and upper confidence limits
Compute the proportion of samples within each SNP that are covered to at least threshold.
coverageBySample(coverage_list, threshold)coverageBySample(coverage_list, threshold)
coverage_list |
an |
threshold |
coverage threshold |
a vector of n.snps length
Compute the proportion of SNPs within each sample that are covered to at least threshold.
coverageBySNP(coverage_list, threshold)coverageBySNP(coverage_list, threshold)
coverage_list |
an |
threshold |
coverage threshold |
a vector of n.samples length
Density, distribution function, quantile function and random generation for the truncated Beta distribution with parameters shape1 and shape2 (and optional truncation bounds)
dtbeta(x, shape1, shape2, lower = 0, upper = 1) ptbeta(x, shape1, shape2, lower = 0, upper = 1) qtbeta(p, shape1, shape2, lower = 0, upper = 1) rtbeta(n, shape1, shape2, lower = 0, upper = 1)dtbeta(x, shape1, shape2, lower = 0, upper = 1) ptbeta(x, shape1, shape2, lower = 0, upper = 1) qtbeta(p, shape1, shape2, lower = 0, upper = 1) rtbeta(n, shape1, shape2, lower = 0, upper = 1)
x |
vector of quantiles |
shape1, shape2
|
non-negative parameters of Beta distributions |
lower, upper
|
bounds for truncation |
p |
vector of probabilities |
n |
number of observations |
ptbeta: Cumulative distribution function for truncated beta
qtbeta: Quantile function
rtbeta: Random number generation
Density function for random generation for truncated exponential distribution
dtexp(x, theta, lower, upper) lltexp(x, theta, lower, upper) mleTexp(x, lower, upper)dtexp(x, theta, lower, upper) lltexp(x, theta, lower, upper) mleTexp(x, lower, upper)
x |
numeric vector |
theta |
rate parameter for exponential |
lower, upper
|
boundary points for truncation |
lltexp: Truncated exponential distrubition
mleTexp: Truncated exponential distrubition
Extract barcode for use in COIL program
extractBarcode(gdsfile, variant.id, barcode.file)extractBarcode(gdsfile, variant.id, barcode.file)
gdsfile |
a |
variant.id |
a numeric vector of variant identifiers |
barcode.file |
character name of output file |
Generate SNP barcode using supplied variant ids. This is constructed by collapsing alleles over the variant ids, forming a barcode for each sample. Then a barcode file containing sample ids + barcode is generated and can be supplied to the COIL web-app or to the command line COIL program.
data.frame containing barcode for each sample
Need to look at how COIL constructed it's orginal barcode. Esp. het calls.
Galinsky, Kevin, et al. "COIL: a methodology for evaluating malarial complexity of infection using likelihood from single nucleotide polymorphism data." Malaria journal 14.1 (2015): 4.
Extract PED files from gds
extractPED(gdsfile, moi.estimates = NULL, use.hets = FALSE, outfile = NULL)extractPED(gdsfile, moi.estimates = NULL, use.hets = FALSE, outfile = NULL)
gdsfile |
a |
moi.estimates |
a vector of MOI values obtained by |
use.hets |
FALSE include heterozygote genotypes |
outfile |
prefix of PLINK files for output (default NULL) |
This function writes a plink .ped and .map file for a given gdsfile. If moi.estimates is set then use.hets is redundant. It will set the sex in each the ped file to 2 if MOI > 1 and set heterozygote genotypes to missing for MOI = 1 calls. This is for use in isoRelate. If the use.hets option is true the genotypes are used as is, other wise heterzygotes are set to missing. This function is slow if there are a large number of variants in the GDS file.
Variant calls of polymorphic variants from 65 P.falciparum isolates in The Gambia, from Amambua-Ngwa A, Tetteh KK et al., 2012 http://www.ncbi.nlm.nih.gov/pubmed/23133397. The data was extracted from VCF files and converted to GDS format from Malaria Genomics Consortium Pf3k version 4. See the README file at ftp://ngs.sanger.ac.uk/production/pf3k/release_4/ for details. The GDS object was obtained by filtering the MalariaGen VCF files and keeping only biallelic SNPs that were polymorphic in The Gambia and passed all filters.
A SeqVarGDSClass object.
ftp://ngs.sanger.ac.uk/production/pf3k/release_4/
Construct non-overlapping Genomic Windows by Chromosome
generateWindows(gdsfile, window_size)generateWindows(gdsfile, window_size)
gdsfile |
a |
window_size |
size of overlapping window in base-pairs |
Helper functions for extracting variant IDs corresponding to previously published barcodes.
getBarcodeVariants(gdsfile, publication = "Volkman2008")getBarcodeVariants(gdsfile, publication = "Volkman2008")
gdsfile |
a |
publication |
a character string representing a publication |
Helper function for producing data frame of genomic coordinates
getCoordinates(gdsfile)getCoordinates(gdsfile)
gdsfile |
a |
Helper function for producing matrix of allele dosages
getDosage(gdsfile)getDosage(gdsfile)
gdsfile |
a |
Compute within-host diversity statistic
getFws(gdsfile)getFws(gdsfile)
gdsfile |
a |
Compute the within host diversity statistic according to the
method devised in Manske et.al, 2012. Briefly, within sample heterozygosity
and within population heterozygosity are computed and assigned to ten equal sized
MAF bins [0.0.05]...[0.45,0.5]. For each bin the mean within sample and population
heterozygosity is computed. A regression line of these values through the orgin
is computed for each sample. The is then .
Currently only works on GATK derived gdsfiles. Needs to be updated to define populations.
Manske, Magnus, et al. "Analysis of Plasmodium falciparum diversity in natural infections by deep sequencing." Nature 487.7407 (2012): 375-379.
getHeterozygosity, getHeterozygosityBySample
Obtain GATK info tags relevant to best practices.
getGATKInfo(gdsfile)getGATKInfo(gdsfile)
gdsfile |
a |
a data.frame with the following columns: chromosome position variant.id qd - (average quality by depth) mq - RMS mapping quality fs - Fisher strand bias mqrank - MQ rank haplotypescore readposranksum
http://gatkforums.broadinstitute.org/gatk/discussion/2806/howto-apply-hard-filters-to-a-call-set
Compute population level heterozygosity
getHeterozygosity(gdsfile)getHeterozygosity(gdsfile)
gdsfile |
a |
Need to update so population ids can be submitted.
Compute heterozygosity by sample
getHeterozygosityBySample(gdsfile)getHeterozygosityBySample(gdsfile)
gdsfile |
a |
Follows the method outlined in Manske et al. 2012. Briefly, for each sample, the reference and and alternate allele frequencies are computed as the proportion of reads covering each allele. Then heterozygosity at a SNP is 1 - (raf^2 + aaf^2)
a numeric matrix of size l by n where l is the number of samples and n is the number of SNPs.
Compute minor allele frequency
getMAF(gdsfile)getMAF(gdsfile)
gdsfile |
a |
Coverage based estimation of minor allele frequences. The MAF is computed by computing the minimum of the proportion of reads covering the reference and alternate allele over all samples.
Currently only supports gds files obtained from GATK callers
Model assessment for fitted model using MSE
getMSE(fitted_models)getMSE(fitted_models)
fitted_models |
a moimix object |
Return estimated model parameters
getTheta(model, k = NULL, criterion = NULL)getTheta(model, k = NULL, criterion = NULL)
model |
stepFlexmix or flexmix object |
k |
if stepFlexmix choose model with k components |
criterion |
if stepFlexmix choose model according to information criterion |
Estimate the Fisher information matrix for a bionimal mixture model
infomat(y, model)infomat(y, model)
y |
a two column matrix of read counts |
model |
flexmix object |
A 2k by 2k Fisher information matrix
A dataset containing information about exonic SNPs in Plasmodium falciparum from Manske et al., 2012 http://www.ncbi.nlm.nih.gov/pubmed/22722859.
A data frame with 29621 rows and 10 variables:
http://www.malariagen.net/data
moimix enables the estimation of MOI from SNV data and to evaluate the effect of MOI on SNV data.
Compute sum of squares for mixture components based on known truth.
mse(model, theta)mse(model, theta)
model |
estimated mixture model |
theta |
true mixture-model parameters of form (pi_1,..,pi_k, mu_1,...mu_k) |
Compute per site per sample coverage
perSiteAlleleCoverage(gdsfile)perSiteAlleleCoverage(gdsfile)
gdsfile |
a |
a matrix of allele counts
Compute per variant depth accross all samples.
perSiteCoverage(gdsfile)perSiteCoverage(gdsfile)
gdsfile |
a |
an integer vector
Compute the proportion of samples covered up any number of bases
perSiteCoverageBySample(gdsfile, threshold)perSiteCoverageBySample(gdsfile, threshold)
gdsfile |
a |
threshold |
integer lower bound for coverage |
a vector of length equal to number of variants with the proportion of samples covered up to threshold
Compute proportion of SNPs covered up to a lower bound of bases within each sample.
perSiteCoverageBySNP(gdsfile, threshold)perSiteCoverageBySNP(gdsfile, threshold)
gdsfile |
a |
threshold |
integer lower bound for coverage |
a vector of length equal to number of samples with the proportion of SNPs covered up to threshold
Plot bafMatrix object
## S3 method for class 'bafMatrix' plot( x, sample.id, assignments = NULL, y = NULL, xlab = "", ylab = "SNV frequency", ylim = c(0, 1), pch = 16, ... )## S3 method for class 'bafMatrix' plot( x, sample.id, assignments = NULL, y = NULL, xlab = "", ylab = "SNV frequency", ylim = c(0, 1), pch = 16, ... )
x |
a bafMatrix object |
sample.id |
character name of sample to plot |
assignments |
integer vector of cluster memberships (NULL) |
y |
unnused argument for plot.bafMatrix |
xlab |
Axis label along chromosome regions |
ylab |
Axis label for SNV frequency |
ylim |
Limits for y axis default is between 0 and 1 |
pch |
Plot symbol default is 16 |
... |
other parameters to pass to |
Plots the genome-wide signal of MOI within an isolate from B-allele frequencies.
Parse estMOI log files in a directory into tidy data frame
readEstMOI(path)readEstMOI(path)
path |
path containing .txt files containing estMOI estimates |
a data.frame
Assumes that outPrefix set as option in estMOI is set to sample identifer.
Simulate random numbers from a truncated exponential distribution
rtexp(n, theta, t)rtexp(n, theta, t)
n |
number of random values to return |
theta |
rate parameter for exponential |
t |
cut-off threshold for truncation |
Simulate SNP alternate allele frequencies from a truncated exponential distribution Code is modified from that given at https://stat.ethz.ch/pipermail/r-help/2010-June/242615.html
n random values from TEXP
Generate iid binomial mixture random variables for testing EM implementation
sampleMM(n, N, k, mu, pi)sampleMM(n, N, k, mu, pi)
n |
number of realisations |
N |
number of trials |
k |
number of components |
mu |
true mixture components |
pi |
true mixture weights |
sum number of samples/snps meeting coverage threshold scaled by total number of SNPs/proportions
scaledProportion(x, threshold, scale)scaledProportion(x, threshold, scale)
x |
integer vector of read counts |
threshold |
integer lower bound for coverage |
scale |
integer to scale results. |
se is computed as the inverse of the Fisher information matrix at the EM estimates
seMM(y, model)seMM(y, model)
y |
a two column matrix of read counts |
model |
flexmix object |
Generate read count data in support of SNVs partioned by mixtures of clones within a host. This is done by setting true clonal mixture proportions and the underlying sampling probabilities for each clone. Then a random sample of underlying SNV probabilities is sampled from either a rbeta or rtexp distrubition with given parameters and read counts are generated by bionimial conditonal on the underlying SNV probability with the clone sampling probabilities.
simulateMOI( n.samples, n.snps, moi, mean_coverage, error_coverage, pi.true = NULL, mu.true = NULL )simulateMOI( n.samples, n.snps, moi, mean_coverage, error_coverage, pi.true = NULL, mu.true = NULL )
n.samples |
number of infected hosts in population |
n.snps |
number of SNPs observed |
moi |
number of clonal infections |
mean_coverage |
average coverage of reads |
error_coverage |
deviation in coverage of reads |
pi.true |
optional matrix of true mixture proportions |
mu.true |
optional matrix of true mixture components |
A list containing the following elements pi.true an n.samples by moi matrix containing true clonal proportions mu.true an n.samplps by moi matrix containing true genotype proportions aaf an n.snps vector of under SNV probabilities read.counts an n.samples by n.snps matrix containing read counts supporting each SNV error.counts an n.samples by n.snps matrix containing number of error reads