inst/doc/freqpcr-intro.R

In freqpcr: Estimates Allele Frequency on qPCR DeltaDeltaCq from Bulk Samples

## ---- echo = FALSE------------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ---- eval = FALSE------------------------------------------------------------
#  library(remotes)
#  install_github("sudoms/freqpcr")
#  
#  library(freqpcr)
#  packageVersion("freqpcr")

## ---- eval = FALSE------------------------------------------------------------
#  ** byte-compile and prepare package for lazy loading
#  Error: (converted from warning) package 'cubature' was built under R version 3.6.3

## ---- eval = FALSE------------------------------------------------------------
#  Sys.setenv(R_REMOTES_NO_ERRORS_FROM_WARNINGS="true")
#  install_github("sudoms/freqpcr")

## -----------------------------------------------------------------------------
library(freqpcr)

## ---- results = "hide"--------------------------------------------------------
# Example: R and S are mixed at the exact ratios 1:9, 1:3, 1:1, and 1:0.
# Four dummy bulk samples are generated for each combination. 
# Template DNA amounts follows the gamma distribution.
# K:2, scaleDNA:1e-11, targetScale:1.5, baseChange:0.3, zeroAmount:1e-3,
# sdMeasure:0.3, and EPCR:0.95. 

A <- rep(1, 16)
trueY <- c(rep(0.1, 4), rep(0.25, 4), rep(0.5, 4), rep(1, 4))
housek0 <- c( 19.39, 19.78, 19.28, 19.58,  18.95, 19.91, 19.66, 19.96,
              20.05, 19.86, 19.55, 19.61,  19.86, 19.27, 19.59, 20.21 )
target0 <- c( 19.16, 19.08, 19.28, 19.03,  19.17, 19.67, 18.68, 19.52,
              18.92, 18.79, 18.8, 19.28,   19.57, 19.21, 19.05, 19.15 )
housek1 <- c( 21.61, 21.78, 21.25, 21.07,  22.04, 21.45, 20.72, 21.6,
              21.51, 21.27, 21.08, 21.7,   21.44, 21.46, 21.5, 21.8 )
target1 <- c( 24.3, 24.22, 24.13, 24.13,   22.74, 23.14, 23.02, 23.14,
              21.65, 22.62, 22.28, 21.65,  20.83, 20.82, 20.76, 21.3 )
d.cmp <- data.frame(A, trueY, housek0, target0, housek1, target1)
print(d.cmp)

## -----------------------------------------------------------------------------
# housek0: control sample (intact), housekeeping gene
# target0: control sample (intact), target gene
# housek1: test sample (restriction enzyme), housekeeping gene
# target1: test sample (restriction enzyme), target gene
# trueY  : exact frequency of R allele
# A      : relative concentration of template DNA (not necessary)

result <- knownqpcr(housek0=d.cmp$housek0, target0=d.cmp$target0,
                    housek1=d.cmp$housek1, target1=d.cmp$target1,
                    trueY=d.cmp$trueY, A=d.cmp$A, verbose=FALSE)
print(result)

## -----------------------------------------------------------------------------
# housek1: sample with known ratio, housekeeping gene
# target1: sample with known ratio, target gene, allele specific primer set

result <- knownqpcr(housek1=d.cmp$housek1, target1=d.cmp$target1,
                    trueY=d.cmp$trueY, A=d.cmp$A, verbose=FALSE)
print(result)

## ---- results = "hide"--------------------------------------------------------
A <- rep(1, 16)
trueY <- c(rep(0.1, 4), rep(0.25, 4), rep(0.5, 4), rep(1, 4))
housek0 <- c( 19.39, 19.78, 19.28, 19.58,  18.95, 19.91, 19.66, 19.96,
              20.05, 19.86, 19.55, 19.61,  19.86, 19.27, 19.59, 20.21 )
target0 <- c( 19.16, 19.08, 19.28, 19.03,  19.17, 19.67, 18.68, 19.52,
              18.92, 18.79, 18.8, 19.28,   19.57, 19.21, 19.05, 19.15 )
housek1 <- c( 21.61, 21.78, 21.25, 21.07,  22.04, 21.45, 20.72, 21.6,
              21.51, 21.27, 21.08, 21.7,   21.44, 21.46, 21.5, 21.8 )
target1 <- c( 24.3, 24.22, 24.13, 24.13,   22.74, 23.14, 23.02, 23.14,
              21.65, 22.62, 22.28, 21.65,  20.83, 20.82, 20.76, 21.3 )

# When the combination of (with or without restriction enzyme) is incomplete,
# the R data frame can be prepared in "long" format.
# First, make a complete data frame and then extract a subset.
d.long.all <- data.frame(
    trueY=rep(trueY, 4), 
    Digest=c(rep(0, 16 + 16), rep(1, 16 + 16)),
    Gene=c(rep(0, 16), rep(1, 16), rep(0, 16), rep(1, 16)),
    A=rep(1, 16*4), 
    Cq=c(housek0, target0, housek1, target1) )

# For example, samples without restriction enzyme (Digest == 0) are only available if trueY == 1.
d.long <- d.long.all[d.long.all$Digest == 1 | d.long.all$trueY == 1, ]
print(d.long)

result <- knownqpcr_unpaired(   Digest=d.long$Digest, Gene=d.long$Gene,
                                trueY=d.long$trueY, Cq=d.long$Cq, A=d.long$A   )

## ---- eval = FALSE------------------------------------------------------------
#  > print(d.long)
#     trueY Digest Gene A    Cq
#  13  1.00      0    0 1 19.86
#  14  1.00      0    0 1 19.27
#  15  1.00      0    0 1 19.59
#  16  1.00      0    0 1 20.21
#  29  1.00      0    1 1 19.57
#  30  1.00      0    1 1 19.21
#  31  1.00      0    1 1 19.05
#  32  1.00      0    1 1 19.15
#  33  0.10      1    0 1 21.61
#  34  0.10      1    0 1 21.78
#  35  0.10      1    0 1 21.25
#  36  0.10      1    0 1 21.07
#  37  0.25      1    0 1 22.04
#  38  0.25      1    0 1 21.45
#  39  0.25      1    0 1 20.72
#  40  0.25      1    0 1 21.60
#  41  0.50      1    0 1 21.51
#  ...
#  all following rows are Digest = 1, i.e., qPCR with restricted enzyme.

## -----------------------------------------------------------------------------
targetScale <- 1.2
sdMeasure <- 0.2
scaleDNA <- 1e-06
baseChange <- 0.2

# For the following parameters, we input arbitrary value
P <- 0.15
K <- 4
EPCR <- 0.97
zeroAmount <- 0.0016

## -----------------------------------------------------------------------------
dmy_cq <- make_dummy(   rand.seed=71, P=P, K=K, ntrap=4, npertrap=8,
                        scaleDNA=scaleDNA, 
                        targetScale=targetScale, 
                        baseChange=baseChange,
                        EPCR=EPCR, 
                        zeroAmount=zeroAmount,
                        sdMeasure=sdMeasure, 
                        diploid=FALSE   )
print(dmy_cq)

## ---- results = "hide"--------------------------------------------------------
# Access the keys of an S4 object with @ instead of $. They are extracted as vector.
N <- dmy_cq@N
housek0 <- dmy_cq@housek0
target0 <- dmy_cq@target0
housek1 <- dmy_cq@housek1
target1 <- dmy_cq@target1

# freqpcr() with beta assumption, assuming haploidy, print every optimization step
result <- freqpcr(  N=N, housek0=housek0, target0=target0,
                    housek1=housek1, target1=target1,
                    EPCR=EPCR, zeroAmount=zeroAmount, 
                    beta=TRUE, diploid=FALSE, print.level=2  )
print(result)

## ---- eval = FALSE------------------------------------------------------------
#  > print(result)
#  An object of class "CqFreq"
#  Slot "report":
#                                      Estimate Fixed  (scaled) (scaled.SE)       2.5%       97.5%
#  P (R-allele frequency)            0.09773189     0 -2.222684  0.60914923 0.03178086   0.2633220
#  K (gamma shape parameter)        20.92728471     0  3.041054  1.83522515 0.57354355 763.5884712
#  targetScale (rel. DNA quantity)   1.11922896     0  0.112640  0.08911953 0.93985370   1.3328388
#  sdMeasure (Cq measurement error)  0.20973065     0 -1.561931  0.32845070 0.11017528   0.3992451
#  EPCR (amplification per cycle)    0.97000000     1        NA          NA         NA          NA
#  
#  Slot "obj":
#  $minimum
#  [1] 6.094915
#  
#  $estimate
#  [1] -2.222684  3.041054  0.112640 -1.561931
#  
#  $gradient
#  [1] -3.400973e-05 -8.275889e-05 -5.170087e-05  8.878422e-05
#  
#  $hessian
#              [,1]        [,2]       [,3]       [,4]
#  [1,]  2.71023719  0.05094362   1.168535 -0.1766756
#  [2,]  0.05094362  0.37630056   2.045198 -0.6539723
#  [3,]  1.16853469  2.04519835 147.389558  6.4578190
#  [4,] -0.17667556 -0.65397228   6.457819 11.1504636
#  
#  $code
#  [1] 1
#  
#  $iterations
#  [1] 12
#  
#  
#  Slot "cal.time":
#     user  system elapsed
#     0.72    0.15    0.88

## ---- eval = FALSE------------------------------------------------------------
#  result <- freqpcr(  N=N, housek0=housek0, target0=target0,
#                      housek1=housek1, target1=target1,
#                      K = 1,
#                      EPCR=EPCR, zeroAmount=zeroAmount, beta=TRUE, print.level=1  )

## ---- eval = FALSE------------------------------------------------------------
#  result <- freqpcr(  N=N, housek1=housek1, target1=target1,
#                      targetScale=1.2, sdMeasure=0.2,
#                      EPCR=EPCR, zeroAmount=zeroAmount, beta=TRUE, print.level=1  )

## -----------------------------------------------------------------------------
# Insert the expected size for target1 
x <- 16.2
N <- c(20, 20, 20)
housek0 <- c(20.5, 25.0, 25.6)
target0 <- c(20.5, 25.5, 26.5)
housek1 <- c(25.0, 26.0, 25.5)
target1 <- c(35.3, 26.0 + x, 34.9)
result <- freqpcr(  N=N, housek0=housek0, target0=target0,
                    housek1=housek1, target1=target1,
                    EPCR=0.9, zeroAmount=0.00005, targetScale=0.6,
                    sdMeasure=0.4, beta=TRUE, diploid=TRUE, print.level=0  )

# Try the expected size + 10
x <- 16.2 + 10
target1 <- c(35.3, 26.0 + x, 34.9)
result <- freqpcr(  N=N, housek0=housek0, target0=target0,
                    housek1=housek1, target1=target1,
                    EPCR=0.9, zeroAmount=0.00005, targetScale=0.6, 
                    sdMeasure=0.4, beta=TRUE, diploid=TRUE, print.level=0  )

## -----------------------------------------------------------------------------
citation("freqpcr")