| popPCR {popPCR} | R Documentation |
EM Mixture Model fitting of dPCR droplet fluorescence
Description
Estimates target concentration by counting positive droplets with Poisson correction. Positive, negative, and rain populations are fitted using EM. Droplets are then classified using Maximum A Posteriori rule
Usage
popPCR(
x,
dist,
volSamp = 20,
volDrp = 0.85,
maxComponents = Inf,
negProbThres = 1e-07,
useOnlyNegProbThres = FALSE
)
Arguments
x |
numeric, vector of droplet fluorescence amplitude |
dist |
character, distribution of the mixture models ("normal", "skewed-normal", "t", "skewed-t") |
volSamp |
numeric, sample volume in microliter |
volDrp |
numeric, droplet (or partition) volume in nanoliter |
maxComponents |
numeric, maximum number of components (e.g. populations) |
negProbThres |
numeric, if only one population was detected, then its assumed as a negative population. Droplets will be classified as positive if its probability given the population < negProbThres. |
useOnlyNegProbThres |
logical, if TRUE, then droplets will be classified as positive if its probability given the leftmost population < negProbThres. Default is FALSE, i.e. classification is done by Maximum A Posteriori rule. |
Value
Returns a result.popPCR S4 class object with attributes
classification - character, vector of droplet classification
dist - character, user-specified parameter for the mixture model
dropletCount - list, droplet classification count
em - list, returned value of EMMIXskew's EmSkew()
estConc - list, estimated target concentration as lambda and sample concentration (with 95% CI)
G - numeric, number of components fitted
memberProb - list, component membership probability of all droplets
Examples
library(popPCR)
# Plot histograms of available data
hist(x_onePop, breaks = 100)
hist(x_twoPop, breaks = 100)
hist(x_multiPop, breaks = 100)
# ---- Mixture model fitting ---- #
# Example 1. One population sample
result <- popPCR(x_onePop, dist = "t")
printSummaryConc(result)
# Output:
# Populations detected : 1
# Total droplets : 8000
# Positive : 1 (0.01%)
# Negative : 7999 (99.99%)
#
# Target copies in sample : 2.9414 ( 95% CI: [ -2.8237 , 8.7064 ] )
# Mean target copies per partition : 1e-04 ( 95% CI: [ -1e-04 , 4e-04 ] )
printSummaryFit(result)
# Output:
# Results of fitting a 1-component t mixture model
#
# Negative Population
# Mix prop. : 1
# Mu : 1024.1614
# Sigma : 35253.1747
# Dof : 2.005
# (Option) increase negProbThres to classify negative droplets more strictly
result <- popPCR(x_onePop, dist = "t", negProbThres = 1e-4)
printSummaryConc(result)
# Output:
# Populations detected : 1
# Total droplets : 8000
# Positive : 691 (8.64%)
# Negative : 7309 (91.36%)
#
# Target copies in sample : 2125.5312 ( 95% CI: [ 1966.9936 , 2284.0688 ] )
# Mean target copies per partition : 0.0903 ( 95% CI: [ 0.0836 , 0.0971 ] )
# Example 2. Two population sample
result <- popPCR(x_twoPop, dist = "t")
printSummaryConc(result)
# Output:
# Populations detected : 2
# Total droplets : 10254
# Positive : 8693 (84.78%)
# Negative : 1561 (15.22%)
#
# Target copies in sample : 44290.3819 ( 95% CI: [ 43215.6408 , 45365.1231 ] )
# Mean target copies per partition : 1.8823 ( 95% CI: [ 1.8367 , 1.928 ] )
printSummaryFit(result)
# Output:
# Results of fitting a 2-component t mixture model
#
# Negative Population
# Mix prop. : 0.1522
# Mu : 2136.7435
# Sigma : 4126.8357
# Dof : 12.3562
#
# Positive Population
# Mix prop. : 0.8478
# Mu : 7580.1275
# Sigma : 42621.1894
# Dof : 2.415
# Example 3. Multiple population sample
result <- popPCR(x_multiPop, dist = "t", maxComponents = 4)
printSummaryConc(result)
# Output:
# Populations detected : 4
# Total droplets : 1814
# Positive : 44 (2.43%)
# Negative : 1252 (69.02%)
# Rain (1) : 258 (14.22%)
# Rain (2) : 260 (14.33%)
#
# Target copies in sample : 8724.5195 ( 95% CI: [ 7999.0578 , 9449.9812 ] )
# Mean target copies per partition : 0.3708 ( 95% CI: [ 0.34 , 0.4016 ] )
# In the output above, we see 2 rain populations! Let's examine its plot.
plot(stats::density(x_multiPop))
# We can see that Rain (1) is very close to the Negative population.
# Let's include droplets in Rain (1) in the negative droplet count.
nNegative <- result@dropletCount$neg + result@dropletCount$rain1
nTotal <- result@dropletCount$total
# Re-estimate concentration as follows
newEstimates <- calculateConc(nNegative, nTotal, volSamp = 20, volDrp = 0.85)
newEstimates
# Output:
# $lambda
# lambda lower upper
# 0.1834247 0.1627763 0.2040731
#
# $conc
# conc lower upper
# 4315.875 3830.031 4801.719