Bayesian Analysis of qRT- PCR Data

12 Pages·2013·1.31 MB·English

No Control Genes Required: Bayesian Analysis of qRT-
PCR Data
Mikhail V. Matz1*, Rachel M. Wright2, James G. Scott3
1DepartmentofIntegrativeBiology,UniversityofTexasatAustin,Austin,Texas,UnitedStatesofAmerica,2InstituteforCellandMolecularBiology,UniversityofTexasat
Austin,Austin,Texas,UnitedStatesofAmerica,3McCombsSchoolofBusinessandDivisionofStatisticsandScientificComputing,UniversityofTexasatAustin,Austin,
Texas,UnitedStatesofAmerica
Abstract
Background:Model-based analysis of data from quantitative reverse-transcription PCR (qRT-PCR) is potentially more
powerfulandversatilethantraditionalmethods.Yetexistingmodel-basedapproachescannotproperlydealwiththehigher
sampling variances associated with low-abundant targets, nor dothey provide a natural way to incorporate assumptions
aboutthe stability ofcontrol genesdirectly intothe model-fitting process.
Results:Inourmethod,rawqPCRdataarerepresentedasmoleculecounts,anddescribedusinggeneralizedlinearmixed
models under Poisson-lognormal error. A Markov Chain Monte Carlo (MCMC) algorithm is used to sample from the joint
posteriordistributionoverallmodelparameters,therebyestimatingtheeffectsofallexperimentalfactorsontheexpression
ofeverygene.ThePoisson-basedmodelallowsforthecorrectspecificationofthemean-variancerelationshipofthePCR
amplificationprocess,andcanalsogleaninformationfrominstancesofnoamplification(zerocounts).Ourmethodisvery
flexible with respect to control genes: any prior knowledge about the expected degree of their stability can be directly
incorporatedintothemodel.Yetthemethodprovidessensibleanswerswithoutsuchassumptions,oreveninthecomplete
absenceofcontrolgenes.WealsopresentanaturalBayesiananalogueofthe‘‘classic’’analysis,whichusesstandarddata
pre-processingsteps(logarithmictransformationandmulti-genenormalization)butestimatesallgeneexpressionchanges
jointlywithinasinglemodel.Thenewmethodsareconsiderablymoreflexibleandpowerfulthanthestandarddelta-delta
Ct analysisbased onpairwise t-tests.
Conclusions:Ourmethodologyexpandstheapplicabilityoftherelative-quantificationanalysisprotocolallthewaytothe
lowest-abundance targets, and provides a novel opportunity to analyze qRT-PCR data without making any assumptions
concerning target stability. Theseprocedures have beenimplemented asthe MCMC.qpcr package inR.
Citation: Matz MV, Wright RM, Scott JG (2013) No Control Genes Required: Bayesian Analysis of qRT-PCR Data. PLoS ONE 8(8): e71448. doi:10.1371/
journal.pone.0071448
Editor:AamirAhmed,UniversityCollegeLondon,UnitedKingdom
ReceivedDecember13,2012;AcceptedJuly2,2013;PublishedAugust19,2013
Copyright: (cid:2) 2013 Matz et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits
unrestricteduse,distribution,andreproductioninanymedium,providedtheoriginalauthorandsourcearecredited.
Funding:ThisworkwassupportedbytheNationalScienceFoundationgrantDEB-1022196.Thefundershadnoroleinstudydesign,datacollectionandanalysis,
decisiontopublish,orpreparationofthemanuscript.
CompetingInterests:AuthorM.MatzisaPLOSONEEditorialBoardmember.Thisdoesnotaltertheauthors’adherencetoallthePLOSONEpoliciesonsharing
dataandmaterials.
*E-mail:[email protected]
Introduction of which is assumed to be constant. This approach is attractive
becauseofitsmathematicaleleganceandtheabilitytocorrectfor
Real-time quantitative PCR [1] is a gold standard for
unequal amounts of biological material (‘‘template loading’’)
quantifying the abundances of nucleic acid targets (DNA or
betweensamplesbyusingthecontrolgeneasaninternalstandard.
RNAmoleculesofaparticularsequence).Oneofitsmostcommon
As a disadvantage, however, it relies on pairwise comparisons of
implementations, quantitative reverse transcription-PCR (qRT-
samplesandthereforemakesitdifficulttohandlemoreelaborate
PCR), is a gene expression quantification method extensively
experimental designs, particularly involving interactions between
applied to test specific hypotheses suggested by genome-scale
factors[3].Furthermore,Pfaffletal[4]pointedoutanadditional
approaches (microarrays or RNA-seq), as well as for analyzing
complication:theneedtoaccountforthedifferenceinefficiencyof
diagnostic gene expression signatures in studies ranging from
amplification of the control and target gene. This problem was
medicinetoecology.DespitethewidespreaduseofqRT-PCR,the addressedbyreplacingtheoriginal22DDCTequationwithafour-
dataprocessingandstatisticalanalysisproceduresarestillinflux,
storyformulaincorporatingtheefficienciesofPCRforcontroland
withmanyalternative approaches coexistingintheliterature and
targetgene[4].Theseefficiencies(amplificationfactorspersingle
new methodologies continuously developed. One of the earliest
thermal cycle) are typically determined by qPCR analysis of
qRT-PCR analysis methods still very commonly used is relative dilution series [4–6], although other methods, based on the
quantification using the 22DDCT (‘‘delta-delta Ct’’) method [2]. It analysis of individual product accumulation curves, have been
compares pairs of samples to see whether a target gene became suggested [7–10].
moreorlessabundantrelativetothecontrolgene,theexpression
PLOSONE | www.plosone.org 1 August2013 | Volume 8 | Issue 8 | e71448
BayesianAnalysisofqRT-PCRData
The need to base qPCR analysis on dilution series data made the sample failed to amplify simply because it contained zero
relativequantificationpracticallyequivalenttotheotherflavorof target molecules. Finally, no solution was provided to directly
qRT-PCR analysis, absolute quantification [11]. Under this include information about control genes into the model-fitting
approach, the raw Cq (‘‘cycle of quantification’’) values are process.
transformed into concentrations of the target per qRT-PCR Here,themixed-modelingapproachofSteibeletalisextended
reaction,basedonthecalibrationcurveconstructedacrossaseries to account for all these issues. First, a generalized linear mixed
of known target concentrations. The efficiency of amplification is model based on the Poisson-lognormal distribution replaces the
implicitlytakenintoaccountduringthisconversion.SinceinqRT- original Gaussian model. This properly handles zero counts, as
PCRtheknowledgeoftheabsolutetargetamountisusuallynotas well as the shot-noise variance associated with low-abundant
important as the knowledge of its variation across samples, so- targets. Second, the model fitting process involves a Bayesian
called relative calibration curves, created by diluting an arbitrary MCMC sampling scheme, and can directly incorporate informa-
amountoftarget,areoftenused[12,13].Theserelativecalibration tionaboutcontrolgenesintheformofpriors.Ourimplementation
curvesprovideessentiallythesameinformationasthecalibration ofthemethodleveragestheMCMCglmmpackageinR[29]and
curves for determining PCR efficiency [4]. To account for ispresentedintheformofaspecializedRpackage,MCMC.qpcr.
variationintemplateloadingacrosssamples,theprocedurecalled
normalizationisperformedwheretheinferredtargetamountsare Results
divided by the abundance of a control gene. In this way, all the
targetabundancesbecomeexpressedasfolddifferencesrelativeto Motivating example
theabundanceofthecontrolgene[14].However,hardlyanyone Thedatasetthatwaschosenforre-analysisaddressedtheeffects
gene remains perfectly stable [15]. Vandesompele etal [16] of heat-light stress and recovery in a reef-building coral Porites
proposedthatmoreaccuratenormalizationcouldbeachievedby astreoides [23]. Briefly, eight individual colonies of the coral were
using multiple nearly-stable control genes: in this case, the target fragmented into 4 pieces each and allowed to acclimate in
abundances are divided by the geometric average of the control common benign conditions for four days. On the fifth morning,
gene abundances. The same paper also introduced a non- two fragments of each colony were placed into a stressful
parametric method, geNorm, for identification of the most stable environment(elevatedheatandlight).Atmidday,whenthestress
control genes based on covariance across samples, which quickly intensity was the highest, one stressed fragment and one control
became a standard in the qRT-PCR field [17]. Another fragment (remaining in the benign conditions) were sampled,
commonly used method for identifying stable genes is the representing the first sampling timepoint. In the end of the same
parametric NormFinder algorithm [18], which makes use of the day, the second stressed fragment was put back into the benign
factthatthelog-transformedqRT-PCRdatasatisfythenormality environment for recovery. This fragment, along with the
criterion[19]andusesmomentsequationstocalculatethestability remaining control fragment, was sampled at midday on the
of each gene independently of other genes. Following normaliza- followingday(thesecondtimepoint).Expressionof15genes,5of
tion, several authors used further parametric approaches such as whichwere putative control genes, wasassayed byqRT-PCR on
ANOVA [6,19] and linear mixed models [20–23], applied on a LightCycler 480 (Roche) with SYBR-based detection, corrected
gene-by-gene basis, to achieve maximum versatility in analysis of for amplification efficiencies, normalized by the 3 genes that
complicated designs. The workflow involving correction for proved to be most stable according to geNorm test [16], and
amplification efficiency followed by multi-gene normalization analyzed using linear mixed models applied on a gene-by-gene
and gene-by-gene analysis with t-tests, ANOVA, or linear basis[23].
modeling represents the current consensus of qRT-PCR data This dataset is interesting from the analytical standpoint
processing [24]. because of three reasons. First, oneof the maineffects of interest
Wewoulddrawananalogywiththeliteratureondataanalysis istheinteractionterm,Condition:Timepoint,describingthegene
forDNAmicroarrays.Here,ithasbeenrepeatedlyarguedthatthe regulation in coral fragments that were first stressed and then
joint analysis of the whole dataset is more appropriate. Such an allowedtorecover.Evaluationoftheinteractiontermnecessitates
approachcanborrowinformationacrossmultiplegenes,improve the use of linear models or ANOVA rather than non-parametric
the precision of gene-specific estimates, and properly account for methods or pairwise t-tests [3]. Second, the dataset includes an
complex experiment designs, as well as both biological and important random factor: the identity of the coral colony from
technical replication [25]. In the qRT-PCR field, this approach which the experimental fragments were obtained. This factor
has been explored by Steibel etal [3], who developed a accounts for variation in the baseline levels of gene expression
parameter-rich linear mixed model to jointly estimate all the between individualcorals,andprompts theuseof alinear mixed
gene-specific effects from non-normalized qRT-PCR data. The model rather than a simple linear model. Finally, several genes
most notable feature of this model is its inclusion of unobserved were so low-abundant under some conditions that they became
random effects, commontoall genes ina sample.Theserandom undetectable in a considerable number of trials, precluding the
effects account for unequal template loading between samples, straightforward use of log-transformation typical of qRT-PCR
thereby achievinga functional equivalent of normalization. analysis[16,19].
Despite the attractiveness of this approach, three major issues
remained unresolved. First, the approach of Steibel etal. Poisson-lognormal mixed models
disregardsheteroscedasticity,ortheincreaseinsamplingvariance ThestandardpracticeforqRT-PCRanalysisistoanalyzeeach
atthelowerendoftargetabundances.Inthisrespectitissimilarto gene individually, using control genes to estimate the required
essentiallyallexistingqPCRanalysispipelines:thestatisticalmodel normalization factors. In contrast, we build a hierarchical model
ignores the discrete nature of the amplification process. This that can be used to jointly estimate the effects of experimental
heteroscedasticity arises because qRT-PCR is fully capable of treatments on the expression of all genes. Under such an
amplifyingjustafewtargetmoleculeswithineachtrial[26,27].It approach, control genes can sharpen estimates of model param-
thusbecomespronetoPoisson-like‘‘shotnoise’’[28].Second,the eters, but are not strictly necessary, as all normalization happens
methodcannoteasilyderiveinformationfromPCRtrialsinwhich withinthemodel.
PLOSONE | www.plosone.org 2 August2013 | Volume 8 | Issue 8 | e71448
BayesianAnalysisofqRT-PCRData
In this respect, our approach is similar to that proposed by A notable feature of the Poisson-lognormal model is its
Steibel etal [3], who model the cycle of quantification (Cq) using overdispersion relative tothePoisson:ify*PLN(cid:4)l,s2(cid:5),then
linear mixed models (LMMs) with Gaussian errors. Since C is
q
proportionaltothenegativelogarithmofagene’sinitialtranscript
E(y)~exp(lzs2=2)
copy number, this assumption implies that copy number (an
integerquantity)isbeingmodeledwithalog-normaldistribution. var(y)~exp(lzs2=2)zexp(2l)fexp(2s2){exp(s2)g§E(y):
Our approachdiffersinthatwedirectlymodeltheinitialcopy
number using generalized linear mixed models (GLMMs) with
Poisson-lognormalerrors[30].Thisisaverynaturalwaytomodel This is important for adequately describing the technical
multiplicative fold-changes. It is also more appropriate for count variability of qPCR measurements, which need not match the
data than a log-normal model, and can flexibly accommodate a strict mean-variance relationship implied by the Poisson distribu-
wide range of mean-variance relationships. Moreover, unlike the tion. Observe that in the limit as the log-variance goes to 0, the
Gaussian model for Cq, it also gleans information from samples model becomes Poisson.
thatfailtoamplify(Cq=‘),asthesecorrespondtocountsofzero. For the purpose of model-fitting, we appeal to (1) and re-write
Athirdadvantageofourmodelisthatitnaturallyaccountsfor the original model in an equivalent hierarchical form, which
the so-called ‘‘shot noise’’ that arises from the discrete nature of facilitates computationvia Markov-chain Monte Carlo:
PCR amplification. Shot noise refers to Poisson-like fluctuations
that become discernible when the number of target molecules is (cid:4)y Dl ,e (cid:5)*Pois(cid:4)l eegijkr(cid:5)
smallenoughsothatsuchfluctuationsarethedominantsourceof gijkr gijk gijkr gijk
variability after signal amplification. To see why this occurs for l ~exp(y )
gijk gijk
weaksignals,observethatifthetruenumberoftargetmoleculesin
a sample is Poisson distributed, then the absolute magnitude of egijk*N(0,s2g):
shot-noise variation grows like the square root of the expected
number of molecules. This is much slower than linear growth,
meaningthattherelativecontributionofshotnoisedecreasesand The e terms, which are mutually independent, appear as
thesignal-to-shot-noise ratioincreasesas theexpectednumber of random effects in the hierarchical specification. This gives the
countsgetslarger.Thisexplainswhyshotnoiseismorefrequently appearanceofasaturatedorevennon-identifiedmodel.Butthey
observedwhenamplifyingsampleswithveryfewtargetmolecules, are best thought of as merely data-augmentation variables that
as weillustrate experimentally below. yield a computationally efficient way to recover the original
In describing our model, we use the following subscript Poisson-lognormal specification (which is neither saturated nor
conventions: unidentified).
– g:gene Structure of the regression model
– i:level ofa treatment condition Our model forlog-rate termygijk~log(lgijk)takes theform
– j:levelofagroupingvariable,e.g.block,plot,geneticline,etc.
– k:biological replicate (a singleRNA sample) y ~I zB zt za zs : ð2Þ
gijk g ig k jg kg
– r: technicalreplicate.
Thusy isthecount(initialtranscriptcopynumber)forgeneg Wedescribeeachcomponentofthemodelinmoredetail,along
gijkr
undertreatmenti,groupj,samplek,andtechnicalreplicater.Our with thepriorsused fortherandom effects.
model assumes that y arises from a Poisson-lognormal
gijkr
distribution: – Igrepresentsagene-specificintercept.Thisisalwayspartofthe
model;itrepresentsthebaselinelevelofeachgene’sexpression
(cid:2) (cid:3) underthe‘‘reference’’orbaselinecombinationofexperimental
ygijkr*PLN lgijk,s2g treatments (fixed factors), to which all other treatment effects
will be compared. In our coral example, a natural choice for
where PLN(m,v) denotes the Poisson-lognormal distribution with the reference combination is the control condition at the first
rate parameter m and log-variance v. In our model, the log- timepoint.
varianceisgene-specific,andtheratetermsinvolveregressionson – B isthefixedeffectoftreatmentiongeneg.Thistermderives
ig
both fixed and random effects, as detailed in the subsequent fromtheexperimentaldesignmatrix,andconsistsofaseriesof
section. gene-specific fixed effects that we are primarily interested in.
There is no closed-form expression for the density of the We denote this term B to signify that it captures the main
Poisson-lognormal, but it may be interpreted as a mixture of biological relevance of the study. In our coral example, this
Poissons.Specifically, supposethat corresponds to a series of interaction terms: gene/condition,
gene/timepoint, andgene/timepoint/condition.
y*PoisðleeÞ, e*N(cid:4)0,s2(cid:5): ð1Þ – t isarandomeffectmeanttocaptureunequaltemplateloading
k
(hencet)inbiologicalreplicatek.Intuitively,thisaccountsfor
the reality that, all else being equal, biological replicates may
Then the marginal distribution of y is Poisson-lognormal with stilldiffer systematically inthetranscriptcopynumbersacross
parameters l and s2. Intuitively, the PLN is similar to the allgenesduetovariationinthegrossamountand/orqualityof
negative-binomial distribution, which can also be expressed as a RNA among samples. These effects are modeled with a
mixtureofPoissons,butwhich(unlikethePLN)hasaclosed-form Gaussianpriorwherethegene-specificvarianceiscommonto
density. all replicates.
PLOSONE | www.plosone.org 3 August2013 | Volume 8 | Issue 8 | e71448
BayesianAnalysisofqRT-PCRData
– a isthegene-specificrandomeffectassociatedwiththejthlevel estimate. In the simplified notation that we use throughout this
jg
of some grouping variable, such as block, plot, litter as in section,wewilllistthenamesofrandomfactorsinsquarebrackets,
Steibel etal [3], or coral colony in our motivating example. to discriminate them from the factors of primary interest (‘‘fixed
Depending on the experimental design, there may be more factors’’) that wediscussed before:
than one grouping variable, or none at all. These effects are
modeled with a Gaussian prior whose variance may, in ln(rate)*genezgene:Conditionzgene:Timepointz
principle,be gene-specific.
– s captures residual variation across different biological gene:Timepoint:Conditionz½sample(cid:2):
kg
samples, assuming that some genes might vary more than
others.TheseeffectsaremodeledwithaGaussianpriorwhose
Notethat,sincethevariationincDNAqualityand/orquantity
variance mayalsobegene-specific.
affectsallgenesinasampleinthesameway,thisrandomfactoris
notgene-specific.Theintroductionofthisrandomfactorintothe
Explaining the model using coral example qRT-PCR model was perhaps the most important innovation in
To make sure we are understood not only by statisticians but themodel of Steibeletal [3].
also by qRT-PCR practitioners, below we explain in more The experimental design might have involved additional
colloquial terms how the model (2) is constructed for our ‘‘groupingfactors’’thatarenotdirectlyrelatedtotheexperimental
motivating example. treatments being studied but still might be responsible for a
Themodelhasasingleresponsevariable,thetranscriptcount, considerableproportionofvariationandmustbeaccountedforto
whose rate is modeled on a log-linear scale. The most basic achievemoreaccuratepredictions.Thesefactors,ifpresent,would
explanatoryvariableinthemodelis‘gene’,whichcorrespondsto correspond to the term ajg in the formula (2). For example, the
the term I in the formula (2) and accounts for different levels of experiment might have involved repeated measurements of
g
expression between genes.We mayexpress thisinformally as: participatingindividuals,partitioningoftheexperimentalsubjects
between several blocks (plots, tanks) for technical reasons, or
measurements of all the effects of interest on different genotypes.
ln(rate)*gene,
Thelatteristhecaseinourcoralexample,whereweused8coral
colonieseachsplitintofourclonalfragmentsthatwereexposedto
where ‘‘rate’’ refers to the count rate of the Poisson-lognormal
our experimental treatments. The grouping factors can be
model. Intuitively, this is the likely frequency of transcript counts
specified in the model as additional random factors; however, in
for a particular gene in a particular sample. Such a primitive
contrast to the sample factor, these would be gene-specific since
model,however,wouldbeoflittlevaluesinceinqRT-PCRweare
different genes might be affected by the grouping factors
typicallynotinterestedindifferenceinexpressionbetweengenes,
differently.Inourcase,wewanttoaccountforpossibledifferences
but want to learn how the expression of each gene varies
in baseline level of expression of each of our genes between 8
depending on the experimental treatments. To find this out, we
colonies,and weaugmentour model asfollows:
augmentourmodelwithaseriesoftermsdescribinggene-specific
effectsofexperimentaltreatments,correspondingtothetermB in
ig
theformula(2).Inourcoralexperiment,wehavetwotreatments ln(rate)*genezgene:Conditionzgene:Timepointz
(or, using linear modeling terminology, factors): Condition with gene:Timepoint:Conditionz½sample(cid:2)z
levels‘‘control’’and‘‘heat’’andTimepointwithlevels‘‘one’’and
‘‘two’’, plus their interaction (i.e., we suspect that there might be z½gene:colony(cid:2)
some Timepoint-specific effects of Condition). This experimental
design isincorporated into themodel asfollows:
Once again, the model is flexible in the number of grouping
factors that could beincluded.
The two remaining terms that we still need to add are both
ln(rate)*genezgene:Conditionzgene:Timepointz
error terms, accounting for the residual variation that remained
gene:Timepoint:Condition, unexplained. The first one is specified as a random factor and
reflects the unexplained differences between biological replicates
where the colon indicates interaction, essentially standing for ‘‘- (samples),correspondstotheterms intheformula(2).Itmakes
kg
specificeffectof‘‘.Themodelisfullyflexible,notbeinglimitedto sensetoassumethat thisfactor wouldbegene-specific, i.e.,some
aparticularnumberoffactors,numberoflevelswithineachfactor, genes willvary morethan othersamong samples:
or presence-absence ofinteractions.
Even though this model specification seems to contain all the ln(rate)*genezgene:Conditionz
terms wewanttoestimate, wemust takecare ofother important
sourcesofvariationthat,whilebeingofnorealinteresttous,must gene:Timepointzgene:Timepoint:
be taken into account to ensure that the model is accurate and
Conditionz½sample(cid:2)z
powerful.Themostimportantoftheseistherandomeffectofthe
biological replicate (i.e., an individual RNA sample), accounting z½gene:colony(cid:2)z½gene:sample(cid:2)
forthevariationinqualityand/orquantityofbiologicalmaterial
among samples, which corresponds to the term t in the formula
k
(2). The designation ‘‘random effect’’ implies that we are not Finally, the remaining unexplained variation would be due to
interestedinactualestimatesofeachsample’squalityorquantity, the differences between technical replicates, reflecting the preci-
but simply want to partition out the corresponding variance. sion of the qPCR instrument used. This term corresponds to the
Randomeffectsareimaginedasrandomvariablesdrawnfroman error term e of the general Poisson-lognormal model given by
underlying distribution the variance of which the model will equation(1).WefollowSteibeletal[3],whofoundthatthemodel
PLOSONE | www.plosone.org 4 August2013 | Volume 8 | Issue 8 | e71448
BayesianAnalysisofqRT-PCRData
fitistypicallyimprovedwhenspecifyingthistermasgene-specific: distribution for B provides a direct answer to questions about
ig
probable effect size – or example, via a 95% posterior credible
ln(rate)*genezgene:Conditionzgene:Timepointz interval for each Big term. Since within a single qRT-PCR
experiment multiple comparisons of this kind are typically
gene:Timepoint:Conditionz½sample(cid:2)z performed, there is a need to correct these results for multiple
z½gene:colony(cid:2)z½gene:sample(cid:2)z½gene:residual(cid:2) testing(butsee[31]foradiscussionofthemultiplicityissueinthis
context). There are also many strategies for full Bayesian model
selectionthatnaturallyaccountformultipletesting,e.g.[32].But
these are computationally intensive, and cannot be straightfor-
It is important to note that separating variances due to
gene:sample from gene:residual is only possible when the dataset wardlyapplied innon-Gaussian settings suchas ours.
contains technical replicates (which, ideally, it should); otherwise AlthoughitislessnaturaltodosoundertheBayesianparadigm,
thevariances collapse intoa singleterm gene:residual. onemayalsouseposteriortailareastoconstructaprocedurethat
behaves very much like a classical significance test. Specifically,
Specification of priors define thetwo-sidedBayestailarea pigastwicethefractionof all
sampledparametervaluesforB thatcrosszerowithrespecttothe
For fixed factors involving genes that are not designated as ig
posterior mean. These values would correspond to p
control genes, a diffuse normal prior is used, with mean =0 and MCMC
very large variance (108) [29]. For control genes, we want to calculated within the MCMCglmm package [29]. We treat the
p ’sasiftheywereclassicalp-valuesandcorrectthemformultiple
specifythattheyshouldbestable(i.e.,havemean =0)withbetter ig
testingusingamethodthatcontrolsthefalse-discoveryrate(FDR).
confidence,whichmeanssomesmallerpriorvariance.Withinthe
For large MCMC samples, our definition of p based on
method implemented in the MCMC.qpcr package, the default ig
posteriortailareasisusuallysufficient.Butthelowestnon-zerop-
setting of the stability parameter for designated control genes
valuethatcanbethusobtainedis2/M,whereMisthesizeofthe
allows them to vary 1.2-fold on average across experimental
MCMC sample. To derive lower p-values based on a limited
conditions, but this value can be increased to relax the stability
assumption, or decreased downto1, whichwould mean that the MCMC run, we approximate the posterior distribution of each
controlgenesareexpectedtobeperfectlystable.Importantly,itis parameter by a normal distribution and calculate a Bayesian z-
also possible to run the analysis in the ‘‘naive’’ form, without score(themeanoftheposteriordividedbyitsstandarddeviation).
specifying any control genes. Thisyields a two-tailed p-valuebased ona standard z-test.
Inadditiontorestrictingthemeanchangeofthecontrolgenes There are at least two reasons why proceeding in this manner
in response to fixed factors, we might also wish to restrict their yields a sensible, though not exact, test. First, the Bernstein–von
variances due to gene-specific random effects, a and s in the Mises theorem implies that, under quite general conditions, the
jg kg
formula (2). MCMC.qpcr package provides an option to fix the jointposteriordistributionbehavesasymptoticallylikeamultivar-
variancecomponentsforthecontrolgenesatsomespecifiedvalue. iate normal distribution centered at the maximum-likelihood
For variance components of all non-control genes, flat non- estimate, and with inverse covariance matrix given by the Fisher
informative priors are used by default, resulting in estimates informationmatrix[33].ThisimpliesthatBayesiantailareasare
approximating the maximum likelihood method [29], with an asymptotically equivalent toclassical p-values (see,e.g.[34]).
option to substitute them for two types of inverse Wishart priors. While this asymptotic guarantee may be cold comfort for
We found that, at least in the experiments described here, the researcherswithmodestsamplesizes,itshouldbeemphasizedthat
inferredeffectsizesandcredibleintervalswerevirtuallyunaffected even purely classical analyses of generalized linear mixed models
bypriorvariancespecificationsandthereforechosenottodiscuss yieldsignificanceteststhatarevalidonlyasymptotically(e.g.[35]
theseoptionshere,eventhoughtheymaycomeusefulinthefuture pg 385). Indeed, the construction of exact significance tests that
intheexperimentsspecificallydesignedtoquantifyvariances(such properly account for the sampling distribution of both random
as,for example,inquantitative genetics). effects and variance components in mixed models is notoriously
challenging, andan openarea of statistical research.
Model estimates and credible intervals Second, the key feature of a p-value is that it has a uniform
One notable advantage of the MCMC-based approach is that distributionunderthenullhypothesis.Weconductedasimulation
pointestimatesandcredibleintervalsforanymodeledeffectscan studytocheckwhetherthisfactholdsforpigcalculatedonthebasis
beeasilycalculatedbasedontheparametervaluessampledbythe of the Bayesian z-score in a realistic scenario. The results were
Markov chain. A credible interval is a Bayesian analogue of the encouraging (see below), seemingly justifying the use of pig as a
confidence interval in frequentist statistics, and is defined as an classical test statistic.
interval that, with a specified posterior probability (e.g. 0.95),
contains the true value of the parameter. Pairwise differences Conversion of qRT-PCR data to counts
between conditions characterized by various factor combinations The central procedure in our method is the transformation of
can also be computed, along with their credible intervals. This is raw Cq values into molecule counts. In principle, this can be
useful for situations when factors have more than two levels. In achieved using absolute quantification curves [11], which would
addition to the fixed effects, credible intervals of variance need to be constructed for all genes. One disadvantage of this
components can be similarly examined; however, it must be approachisthatgenerationofsuchcalibrationcurvesmustrelyon
remembered that interval estimates for variance components are an independent methodof molecule quantification, whichideally
robust only for data sets with many replications at the should be more precise than qPCR. Here we explored an
corresponding level ofthemodel hierarchy. alternative approach, which relies on the knowledge of amplifi-
cation efficiency (E, the factor of amplification per cycle) and the
Testing for statistical significance Cq of a single target molecule (Cq1). The counts can then be
ThequestionofinterestinqPCRanalysisiswhetheraspecific obtained usingtheformula:
treatment had an effect on a specific gene. The posterior
PLOSONE | www.plosone.org 5 August2013 | Volume 8 | Issue 8 | e71448
BayesianAnalysisofqRT-PCRData
Figure1.EffectoftargetconcentrationonCqvarianceandestimationoftheCqforasingletargetmolecule(Cq1).(a–c)Examplesof
amplification of different gene targets from a series of four-fold template dilutions with six-fold technical replication. The red line is the linear
regressionacross6–7mostconcentrateddilutionswheretheincreaseinvariancewasnotyetpronounced.Theamplificationefficiency(E,thefactor
ofamplificationpercycle),calculatedfromtheslopeofthisregression[4],islistedoneachpanelinthelowerleftcorner.Thedottedhorizontalline
marksthevisuallydeterminedCq1.(d)EstimatedCq1valuesforalltestedtargetsplottedagainsttheiramplificationefficiency.Thebluelineislinear
regression,theshadedareacorrespondsto95%confidenceinterval.Theredlineisthecorrelationexpectedfromtheefficiency-basedPCRmodel[4]
andCq1=36foragenewithE=2.
doi:10.1371/journal.pone.0071448.g001
(35).Asexpected,toolargeortoosmallCq1(Fig.2a)affectedpoint-
Count~EðCq1{CqÞ, estimatesconsiderablyonlyfortheleast-abundantgeneswhichwere
zero-bounded(chrom,clect,g3pdh,andhsp16).Forthesegenes,larger
Cq1 resulted in larger inferred fold-change and broader credible
rounded tointeger. (3)
intervals, withtheexception ofhsp16. Forother genes,thepoint-
To directly estimate Cq1 for several gene targets and evaluate estimateswereunaffected,whilethecredibleintervalsexhibitedthe
theextentofitsvariation,weamplifiedseventargetsfromaseries sametendencytoscalewithCq1,buttoamuchsmallerextentthan
of four-fold cDNA dilutions that extended into the range of less forthezero-boundedgenes.AssumingasingleCq1=37forallgenes
thanonetargetmoleculeperamplificationtrialwithsixtechnical had virtually no effect on the inference, compared to the more
replicates per dilution (Fig. 1). As expected, for all targets the accurateformula-approximatedCq1data(Fig.2b).Goodness-of-fit
variancebetweentechnicalreplicatesincreasedastheaverageCq characteristics werevisually indistinguishablebetweenthedataset
per dilution begantoexceed 30(Fig.1 a–c), presumablybecause convertedusingformula-approximatedCq1andthedatasetwiththe
ofmoreandmorepronouncedeffectofshotnoiseduetolowering same Cq1=37 for all genes (Fig. 3). The departures from perfect
number of molecules sampled per trial. We assumed that higher linearity(Fig.3a)andtrendtowardshighervarianceatthelowend
Cq values in dilutions where not all of the six amplification trials ofpredictedvalues(Fig.3b)werenegligible,andthedistributionof
were successful likely corresponded to single-molecule amplifica- lognormal residuals was very close to normal (Fig. 3 c). The
tionevents[26],andestimatedCq1visually(dashedlinesonFig 1 probabilities of Poisson residuals [36] were nearly exactly as
a–c). expected(Fig.3d),supportingthebasicassumptionthatthehigher
The Cq1 values were negatively correlated with the amplifica- varianceat higherCqisdue toPoissonfluctuation innumbers of
tion efficiency estimated from 6–7 higher-concentration dilutions sampledmolecules[28].
(blue line of Fig. 1 d, p=0.005). This correlation was reasonably
closetowhatisexpectedunderthemodelwiththesingleefficiency Effect of control genes and fixation parameters
parameter[4]assumingCq1=36foratargetwithefficiencyE=2 Figure4ashowsthecomparisonbetweenna¨ıvemodelandtwo
(redlineonFig.1d),althoughtheslopeoftherealregressionwas informed models in which either one (nd5) or two (nd5 and rpl11)
notably shallower. Below we show that the results of relative control genes were specified but allowed to change 1.2-fold on
quantification within our method are relatively robust to mis- averageinresponsetothefixedfactors.Inclusionofcontrolgenes
specificationoftheCq1value,sotheempiricalformuladescribing did not have much effect on the point-estimates, but led to
theregressiononFig.1d(Cq1=79–21.5E)providesareasonable considerable narrowing of the credible intervals. The most
approximation if the experimental estimates of gene-specific Cq1 narrowing was seen after adding one control gene; addition of
values are not available. Moreover, as we describe in the the second one had much less effect. Figure 4 b compares the
subsequent section, the results were virtually identical even when results of the na¨ıve and informed model with two control genes
wesimply assumed thesameCq1=37forall genes. (thesameasonFig.4a)withthe‘‘fixed’’modelinwhichthetwo
control genes were required toremain perfectly stable. Complete
Cq1 and goodness of fit fixation of control genes results in very narrow credible intervals;
Themis-specificationofCq1canbeexpectedtoaffectprimarily however, perfect stability of expression of any one gene is
low-abundant genes,whichlieinthe‘‘shotnoisezone’’(Cq.30, considered to be an unrealistic assumption [15], and hence the
Fig.1a–c)andforwhichtherangeofabundancesmightbebounded powergainedinthefixedmodelmaybecomingattheexpenseof
byzero(i.e.,emptyamplificationtrialsinthedata).Fig.2showsthe accuracy.
resultsofana¨ıvemodel(nocontrolgenesspecified)estimatingthe
main effect of heat stress in our coral dataset, fitted to the data Analysis of smaller datasets
convertedeitherwithformula-approximatedCq1(Fig.1d)orthe Goodperformanceofthena¨ıvemodel,notrelyingonanycontrol
sameCq1forallgenes,whichwaseithertoolarge(40),ortoosmall gene information, may not be too surprising for a dataset that
PLOSONE | www.plosone.org 6 August2013 | Volume 8 | Issue 8 | e71448
BayesianAnalysisofqRT-PCRData
Figure 2. Effect ofCq1 setting on the point estimates and credible intervals of the fixed effects of stress in the coral dataset,
accordingtothenaivemodel.Thepointsareposteriormeans,the95%credibleintervalsaredenotedasdashedlinesconnectingupperand
lowerintervallimitsacrossgenes,tobettervisualizechangesintheirwidth.(a)Comparisonoftheresultsbasedonformula-approximatedCq1(‘‘est.’’)
withanalysesassumingthesameinflated(40)ordiminished(35)Cq1forallgenes.Thethreemostaffectedgenesareclect,chrom,g3pdh,andhsp16,
whichweresolow-abundantinatleastoneoftheexperimentalconditionssuchthatmanyoftheirqPCRtrialswereempty.(b)Comparisonbetween
analyseswithformula-basedCq1(‘‘est.’’)anduniformCq1=37forallgenes.
doi:10.1371/journal.pone.0071448.g002
containsmanygenesdemonstratingvariousexpressionpatterns,but fromthecoralstressdatacontainingonlyfourgenes:acontrolgene
willitworkwhengenesarefewandtheirexpressionpatternsare rpl11andthreemosthighlyregulatedgenes,theheatshockproteins
unbalanced?Toexplorethisissue,asmallerdatasetwasextracted hsp16, hsp60 and hsp90 (Fig. 5 a, c). Analysis of this small dataset
recoveredexactlythesameregulationpatternsforthefourgenesas
observedinthewholedataset,irrespectiveofwhetherthecontrol
genewasspecifiedasapriorornot(Fig.5a,c).Moreover,evenwhen
the control gene was removed from the dataset, leaving only the
threegenesthatwerestronglyup-ordown-regulatedinconcert,the
sameresultswererecovered(Fig.5b,d).
The ‘‘classic’’ model
Somedatasetsmaynotconformtotheassumptionofourmodel
(2)thatthevariationintemplateloadingbetweensamples,t,can
k
be modeled as a single Gaussian distribution across all exper-
imental conditions. If for some conditions the RNA samples
systematically show substantially lower concentration and/or
quality, the model will infer down-regulation of all genes under
these conditions. Accounting for such a bias would unavoidably
requirerelianceoncontrolgenesfornormalization.Wetherefore
implementedoursingle-modelMCMC-basedapproachinvolving
the multi-gene normalization procedure [16]. The normalized
data are analyzed using the same model as given by formula (2)
only lacking the t term since this variation is supposed to be
k
subtractedout.Sincenormalizationprocedurewouldprecludethe
useofcountsastheresponsevariable,wehavetofallbacktolog-
transformed expression values and run a lognormal rather than
Poisson-lognormalmodel.Thedataforthisanalysisispreparedby
converting raw Cq values into natural logarithms of relative
abundances(Ra)whilecorrectingfortheefficiencyofamplification
using thetransformationintroduced previously[3,23]:
Figure 3. Goodness-of-fit characteristics of the na¨ıve model
appliedtothecoralstressdataset,fortheCq1=37settingfor
Ra~{Cq:ln(E) ð4Þ
allgenes.(a)Plotoflognormalresidualsagainstpredictedvaluesto
testforlinearity.(b)Scale-locationplottotestforhomoscedasticity.A
goodfitiscorroboratedbythelackofpronouncedmeantrendinthese
two plots (red lines). (c) Plot of quantiles of standardized lognormal We call this model ‘‘classic’’ since it is based on earlier
residuals against theoretical quantiles of the normal distribution. Red developments[22–24]andlacksthemainadvancementsproposed
diagonal corresponds to the exact match. (d) Probabilities of
in this paper, such as the use of generalized linear modeling to
experimental Poisson residuals plotted against their theoretical
account for higher variance of low-abundant targets and the
probabilities for one of the MCMC samples. All MCMC samples show
thesamenearlyperfectfittothePoissonexpectations. possibilitytoanalyzethedatawithoutrelianceuponcontrolgenes.
doi:10.1371/journal.pone.0071448.g003 The one innovation this model offers, however, is a single-model
PLOSONE | www.plosone.org 7 August2013 | Volume 8 | Issue 8 | e71448
BayesianAnalysisofqRT-PCRData
Figure4.Effectofaddingcontrolgenes(highlighted)aspriorsontheestimatesoffixedeffectsofstressinthecoraldataset.The
pointsare posterior means,the 95% credibleintervals are denotedas dashed lines connecting upper and lower intervallimits acrossgenes. (a)
Comparisonbetweenna¨ıvemodel(nocontrolgenesspecified)andtwoinformedmodels,withone(nd5)ortwo(nd5,rpl11)controlgenesspecified
(thecontrolgeneswereallowedtochange1.2-foldonaverageinresponsetofixedfactors).(b)Comparisonofthena¨ıveandtwo-geneinformed
modeltothetwo-genefixedmodel,inwhichthesamecontrolgeneswererequiredtobeabsolutelystable.
doi:10.1371/journal.pone.0071448.g004
ratherthangene-by-gene analysis,whichisexpectedtoboostthe thecredibleintervalstothemostpowerfulmodelofthethree,the
powerconsiderablysincethemodeldrawsevidencefromallgenes fixed model (Fig.6b).
simultaneously [25].
For the coral stress dataset, the ‘‘classic’’ model generated Bayesian p-values
virtually identical point-estimates of fold changes as the full The p-values based on Bayesian z-test agreed well with p
ig
Bayesianmodels(Fig.6a),whilebeingcomparableinthewidthof calculated directly from the posterior tail areas (Fig. 7 a),
supporting the validity of the z-test approach. To explore the
distributionofz-testderivedp-valuesundernullhypothesisunder
three models (na¨ıve, informed, and ‘‘classic’’) we generated three
null datasets based on our coral data by subtracting the main
effectsinferredbyeachofthemodelsfromtheefficiency-corrected
Cqvalues(Cqcorrected = Cq?log (E),the‘‘perfectworld’’Cqvalues
2
thatwouldhavebeenobservedifallgeneswereamplifiedwiththe
efficiencyofexactly2,[23]).Wethenran100re-analysesofthese
datasets while randomly shuffling samples between experimental
conditions and calculated Bayesian p-values for the main effects
andtheirinteractionforeachgene,resultingin2,600p-valuesper
model. The frequency distribution of these p-values was nearly
ideallyuniformunderthena¨ıvemodel(Fig.7b)andonlyslightly
skewed under informed and ‘‘classic’’ models (Fig. 7 c,d). This
result suggests that the test based on Bayesian p-values will be of
correctsizeforthena¨ıvemodel,willbeslightlymoreconservative
than the nominal alpha value under the informed model, while
under ‘‘classic’’ model it might be very slightly less conservative.
These deviations of Bayesian p-values from the uniform distribu-
tion under the null hypothesis are probably negligible for all
practical intents and purposes. In particular, they should not
precludetheapplicabilityofmultipletestingcorrectioncontrolling
forfalse discovery rate[37].
Comparison with other methods
Figure 8 presents the results of MCMC.glmm analysis of the
coralstressdatasetusingthena¨ıvemodel,thistimeshowingboth
theeffectofstressandsubsequentrecovery,whichwerethefocus
Figure 5. Analysis of small unbalanced subsets of the coral of the original analysis based on multi-gene normalization and
stress data. (a, c): Three heat shock protein genes plus one control gene-by-gene linear mixed modeling [23]. Notably, the na¨ıve
gene (rpl11), analyzed using na¨ıve, informed, and fixed model. (b,d): model was able to recapitulate the previous findings without
Analysisofonlythethreeheat-shockproteins.(a,b)–effectsofstress; making any assumptions regarding control genes. As before, the
(c,d)–effectsofrecovery.Thepointsareposteriormeans,thewhiskers
maineffectsofstressincludedverystrongup-regulationofhsp16(a
denote 95% credible intervals. Bayesian analysisinfers the same fold-
smallheat-shockprotein),moremodestup-regulationoflargeheat
changesregardlessofspecificationofthecontrolgeneandeveninits
absence. shockproteinshsp60(Fig.8c)andhsp90,anddown-regulationof
doi:10.1371/journal.pone.0071448.g005 actin(Fig.8b);duringrecoverythispatternwasreversed.Forthese
PLOSONE | www.plosone.org 8 August2013 | Volume 8 | Issue 8 | e71448
BayesianAnalysisofqRT-PCRData
Figure 6. Comparison of the procedure based on multigene normalization followed by single-model MCMC analysis (‘‘classic’’
model)tofullBayesianmodels.(a)Point-estimatesandcredibleintervalsfortheeffectsofstressincoraldatasetinferredbytheinformedmodel
(allowingcontrolgenesrpl11andnd5tochange1.2-foldonaverage)and‘‘classic’’modelbasedonnormalizationusingthesametwocontrolgenes.
Pointsareposteriormeans,thewhiskersdenote95%credibleintervals.(b)Comparisonofthepowerofthe‘‘classic’’model(transformedp-values
alongthehorizontalaxis)tothepoweroffullBayesianmodels(seelegend).Thecoloredlinesarelinearregressionstoillustratetrends(nostatistical
implicationsintended),theblackdottedlineis1:1correspondence.‘‘Classic’’modelgeneratesthethecredibleintervalsthatareasnarrowasunder
thefixedmodel,butreliesonmorerealisticassumptions.
doi:10.1371/journal.pone.0071448.g006
four response genes reported the original paper [23], we plotted factor, treatment, with five levels (control and four different heat
theirinferredfold-changesduringstressandrecoveryunderna¨ıve, stressregimes)ontheexpressionofnineresponsegenesincultured
informed and ‘‘classic’’ models against the originally reported murine fibroblasts. A single control gene (glyceraldehyde-3-
values(Fig.9a),observingaveryclosematch.Plottingthenewp- phosphate dehydrogenase, g3pdh) was used to calculate fold-
values(basedonBayesianz-test)againstthepreviousresults(Fig.9 changes relative to the control condition while correcting for
b)illustratedthatthena¨ıvemodelislesspowerful(mostpointsare amplificationefficiencies[4].Thedatasetincludedthreebiological
below the diagonal), informed model appears about as powerful, replicates per factor level andthree technical replicates per RNA
and ‘‘classic’’ model is considerably more powerful than the sample.Forna¨ıveandinformedmodels,weconvertedtheoriginal
original analysis. The new models, including the na¨ıve one, Ct data to counts assuming Cq1=37 for all genes. The ‘‘classic’’
suggested additional significantly changing genes (such as down- model recapitulated the delta-delta Ct results nearly exactly, as
regulation of glyceraldehyde-3-phosphate dehydrogenase, g3pdh, expected,butitisnotablethatverysimilarfold-changeswerealso
duringrecovery,Fig8a),whichwerenotdetectedpreviouslyand inferred by the na¨ıve model, which did not use any control gene
notincludedintheplotsonFig.9.Themuchhigherpowerofthe information (Fig. 10 a). Moreover, the power of all our models
‘‘classic’’ model compared to the previous gene-by-gene analysis (eventhena¨ıveone)wassubstantiallyhigherthanthepowerofthe
confirmedtheexpectationthatasingle-modelapproach,drawing original analysisbasedon pairwise t-tests (Fig.10b).
evidencefromallgenesatonce,wouldbemorepowerfulevenifit
is based on the same data pre-processing pipeline (log-transfor- Discussion
mation andmultigene normalization).
To see how the new methods compares to the classic delta- qRT-PCR data as counts
delta-Ct procedure [2], we re-analyzed the data comprising ThepurposeofCq-to-countstransformationdescribedhereisto
Figure 6 in[38]. Thisexperiment examinedtheeffect ofa single createadatasetinwhichtheincreaseinvarianceatthelowgene
Figure7.Propertiesofp-valuesbasedonBayesianz-test.(a)Correspondencebetweenp-valuesbasedonposteriortailareas(horizontalaxis)
andz-testbasedp-values(verticalaxis)forthestressandrecoveryeffectsinthecoraldataset.(b–d)Frequencydistributionofz-testbasedp-values
obtainedusingna¨ıve(b),informed(c)and‘‘classic’’(d)modelsfromdatasetssimulatedundernullhypothesis.Thefractionofsimulatedp-valuesthat
arelessthan0.05isgivenaboveeachplot.
doi:10.1371/journal.pone.0071448.g007
PLOSONE | www.plosone.org 9 August2013 | Volume 8 | Issue 8 | e71448
BayesianAnalysisofqRT-PCRData
Figure 9. Comparison of the Bayesian analysis of the coral
stress dataset to the previously published results based on
multigene normalization and gene-by-gene linear mixed
modeling [23]. (a) The match between the fold-changes inferred in
the current reanalysis and previously reported changes, for na¨ıve,
informedand‘‘classic’’models.(b)Correspondencebetweenp-values
underthena¨ıve,informed,and‘‘classic’’models(seelegendonpanela)
and the previously reported p-values. Points above the line indicate
higher power (lower p-values) of the new models. The reanalysis p-
valueswerederivedbytheBayesianz-test.
doi:10.1371/journal.pone.0071448.g009
accommodate, and would handle any number of fixed and
random effects andtheirinteractions.
Figure8.VisualizationoftheresultsofBayesiananalysisofthe The ‘‘classic’’ and lognormal Bayesian models
coral stress dataset. (a) Effects of stress and recovery under na¨ıve The ‘‘classic’’ model represents a viable alternative to the full
model.Itcanbeseenthatrecoverygeneregulationisbasicallyamirror
Bayesiananalysisespeciallyforthecaseswhenthequantityand/or
imageofstressresponse.(b)Transcriptabundancesofselectedgenes
quality of the RNA samples varies systematically across experi-
(see legend) across conditions of interest. The points are posterior
means,thewhiskersdenote95%credibleintervals. mentalconditions,andwhentheexpressionlevelsarenottoolow
doi:10.1371/journal.pone.0071448.g008 (i.e., the majority of Cq values are below 30). The model is
considerably more powerful than previously used procedures
expressionvalues(Fig.1a–c)canbeadequatelyaccountedforby based on the same principles of data processing (Fig. 9 and 10),
the relative quantification model. General reasoning, earlier highlighting the advantage of joint analysis of expression of all
literature [28], and our analysis presented here (Figure 3 d) geneswithinthesamemodel.Itisalsomorepowerfulthatanyof
suggest that this increase in variance is likely due to Poisson- thefull-Bayesianmodelsdescribedhere,withtheexceptionofthe
fixedmodel(Fig.7).Thepowerisexpectedtoimprovewithmore
distributedfluctuationswhenthenumberofsampledmoleculesis
small.Ithasbeenpreviouslyshownthattheamplificationproducts
obtainedfromDNAsamplessohighly-dilutedthatonlyafraction
ofPCRtrialsaresuccessfulindeedcorrespondtoamplificationof
individual molecules [26]. This substantiates our approach of
determiningCq1,thenumberofPCRcyclesrequiredtoamplifya
single target molecule, through analysis ofover-extended dilution
series(Fig.1a–c).Itisnotable,however,thattheCq1valuesthus
determinedshowlesspronouncedcorrelationwiththeefficiencyof
amplification (Fig. 1 d) than expected under the simple qPCR
model with a single efficiency parameter [4]. We must therefore
caution that even though Cq-to-counts data transformation
achieves its primary purpose, enabling quantification of very
low-abundant targets within the general relative quantification
framework, further qPCR model development might be required
toensureaccurateabsolutequantificationofindividualmolecules.
Figure 10. Comparison of the Bayesian method to the delta-
delta-Ct method, reanalyzing data from [38]. (a) Fold-changes
The Poisson-lognormal Bayesian model
inferred by the na¨ıve, informed and ‘‘classic’’ models plotted against
This model is designed to be universally applicable in qRT- fold-changesderivedbythedelta-delta-Ctanalysiswithasinglecontrol
PCR.Itproperlyhandlesthefullrangeoftargetabundancesdown gene.(b)Comparisonofp-valuesderivedbytheBayesianz-testfrom
to the individual molecule level, derives information from no- na¨ıve,informed,and‘‘classic’’models(seelegendonpanela)withp-
amplification trials, and provides a possibility to specify control values obtained by pairwise t-tests within the delta-delta-Ct pipeline.
On both panels, the line denotes 1:1 correspondence. On panel b,
geneswiththedesireddegreeofconfidence.Inthesametime,itis
pointsabovethelineindicatehigherpower(lowerp-values)ofthenew
fully flexible in terms of experimental designs that it can
models.
doi:10.1371/journal.pone.0071448.g010
PLOSONE | www.plosone.org 10 August2013 | Volume 8 | Issue 8 | e71448