[R] create pcrfit model with qpcR package

Luigi Marongiu marongiu.luigi at gmail.com
Sun Mar 15 12:05:55 CET 2015


Dear all,
I am trying to set a model using the pcrfit() function of the qpcR
package. I have defined a dataframe with 45 cycles (rows) and 4
readings (columns) derived from a TaqMan qPCR with relative
quantification. The data is based on the normalized fluorescent
results (main fluorescence minus ROX).
The example reported herein has a couple of readings with abnormal
amplification (which would be the usual kind of results for the PCR I
am carrying out) and two with a proper amplification; the plot section
of the example is meant to describe the data itself.
For each implementation, the error returned is:

Error in model.frame.default(formula = ~Fluo + Cycles, data = DATA,
weights = WEIGHTS,  :
  variable lengths differ (found for '(do.optim)')

However each column has 45 rows.

Any tips on how to problem-solve this issue?
Many thanks,
Luigi

PS:
The following example writes a file in the current folder, uses it for
the pcrfit modelling and then removes it.
>>>
# define working directory
setwd(getwd())

# qpcR
if(!is.element('qpcR', installed.packages()[,1]))
{
  install.packages('qpcR')
} else {
  library(qpcR)
}

# create dataframe
col.0<-1:45
col.a<-c(-0.033,  -0.027,    -0.015,    -0.007,    -0.001,    0.001,
 0.003,    0.004,    0.005,    0.003,    0.004,    0.001,    0.002,
0.002,    0.001,    0.002,    0.002,    0.001,    0.001,    0.001,
0.004,    0.002,    -0.001,    0.001,    0.001,    -0.001,    -0.002,
  -0.001,    0.001,    -0.001,    -0.002,    -0.001,    -0.002,
-0.001,    -0.003,    -0.004,    -0.002,    -0.001,    -0.002,
-0.001,    0.001,    0.001,    0.002,    0.003,    0.006)
col.b<-c(-0.128,    -0.116,    -0.089,    -0.068,    -0.048,
-0.032,    -0.018,    -0.006,    0.004,    0.009,    0.015,    0.02,
 0.024,    0.025,    0.028,    0.029,    0.028,    0.029,    0.029,
0.026,    0.025,    0.024,    0.023,    0.019,    0.018,    0.018,
0.015,    0.012,    0.01,    0.009,    0.006,    0.003,    -0.001,
-0.003,    -0.006,    -0.01,    -0.012,    -0.014,    -0.017,
-0.02,    -0.023,    -0.024,    -0.027,    -0.029,    -0.029)
col.c<-c(-0.045,    -0.048,    -0.043,    -0.036,    -0.023,    -0.01,
   0.006,    0.02,    0.029,    0.034,    0.037,    0.032,    0.026,
 0.017,    0.01,    -0.002,    -0.013,    -0.021,    -0.029,
-0.034,    -0.025,    0.006,    0.065,    0.163,    0.301,    0.477,
 0.68,    0.89,    1.092,    1.254,    1.41,    1.546,    1.664,
1.772,    1.865,    1.929,    1.983,    2.032,    2.066,    2.09,
2.108,    2.126,    2.129,    2.13,    2.142)
col.d<-c(-0.049,    -0.04,    -0.026,    -0.013,    -0.006,    0.002,
  0.007,    0.009,    0.01,    0.01,    0.009,    0.007,    0.007,
0.007,    0.005,    0.003,    0,    -0.002,    -0.005,    -0.006,
-0.006,    -0.008,    -0.005,    0,    0.01,    0.032,    0.072,
0.14,    0.237,    0.361,    0.51,    0.677,    0.847,    1.016,
1.181,    1.336,    1.481,    1.619,    1.742,    1.858,    1.962,
2.056,    2.136,    2.212,    2.27)
x.file<-cbind(col.0, col.a, col.b, col.c, col.d)

# write data working file
write.table(
  x.file,
  file="x.file.txt",
  sep="\t",
  eol = "\n",
  na = "NA",
  dec = ".",
  quote = FALSE,
  col.names = TRUE,
  row.names = FALSE,
  append = FALSE
)

# create objects
obj<-pcrimport2(
  file="x.file.txt",
  sep="\t",
  dec=".",
  header=TRUE,
  colClasses="numeric",
  quote=""
)

# calculate Cy0
Cy<-1:4
model<-1:4
model[1]<-pcrfit(obj, cyc=1, 2, model=l4, do.optim=TRUE, robust=FALSE)
model[2]<-pcrfit(obj, cyc=1, 3, model=l4, do.optim=TRUE, robust=FALSE)
model[3]<-pcrfit(obj, cyc=1, 4, model=l4, do.optim=TRUE, robust=FALSE)
model[4]<-pcrfit(obj, cyc=1, 5, model=l4, do.optim=TRUE, robust=FALSE)
Cy[1]<-Cy0(model[1], plot=FALSE)
Cy[2]<-Cy0(model[2], plot=FALSE)
Cy[3]<-Cy0(model[3], plot=FALSE)
Cy[4]<-Cy0(model[4], plot=FALSE)

# plot results as control
plot(col.a ~ col.0)
plot(col.b ~ col.0)
plot(col.c ~ col.0)
plot(col.d ~ col.0)

# remove working file
unlink("x.file.txt",
       recursive = FALSE,
       force = FALSE
)



More information about the R-help mailing list