[BioC] HTqPCR question

[Heidi Dvinge]

Hello Simon,

your questions seems to have made it to the Bioconductor list with a few
days delay, so I'm continuing our off-list discussion from last week here
(everything copied in below).

Regarding the normalisation, it sounds like you might want to try either
the "norm.rankinvariant" or "scale.rankinvariant" methods in
normalizeCtData. Both of them identifies a set of rank-invariant genes,
i.e. don't rely on e.g. actin. "scale.rankinvariant" simply scales each
qPCR card by a fixed factor identified based on the rank-invariant genes.
"norm.rankinvariant" fits a smoothing line through the rank-invariant
genes, and then applies this smoothing to the rest of the genes.

The latter works best if there is more than a few rank-invariant genes.
For a comparison of the data before/after normalisation, you can do a
simple scatterplot, similar to figure 7 in the vignette (also included in
the examples in ?normalizeCtData). This gives you an idea of how each
method performs, and what the differences between them are.

Regarding how replicates are treated, then this depends on the function.
Functions such as setCategory() for marking some genes as potentially
unreliable can account for technical replicates within each qPCR plate by
setting replicates=TRUE. However, replicates from different qPCR plates
are treated identically, whether they come from biological or technical
sample replicates.

When testing for differential expression, technical vs biological
replicates can be accounted for using the function limmaCtData, which
incorporates funcitonality from the "limma" package. There are some
examples of this in the limma user manual, accessible via
limmaUsersGuide() in R, especially in section 8.2.

In your case that would involve calling duplicateCorrelation() first,
followed by limmaCtData(). However, if you give some more details about
which samples you have + what type of replicates they are, we can probably
help you construct you design matrix if you haven't used limma before.

HTH
\Heidi

======================
Off-list conversation copied in below, for the archives
======================

Hi Simon,

> Thanks Heidi,
> I posted my questions (plus a few others), on the forum as well.

hm, I seemed to have overlooked those in my Bioconductor email box. Do you
know which date you sent your questions? I that case I can probably chase
them up.

> So for
> the flags, HTqPCR needs them to be either "Passed" or "Flagged"? Then you
> just substitute using rep as you suggest below?

Actually, the flags can be anything. For an example, I just used the
default qPCRset object that comes with HTqPCR. You can read in data
without specifying flags, and alter that later as you wish.

In the example I email you, the flag should be read in automatically:

test <- readCtData(files="1131135160.csv", # file name
+     n.features=2304, # if you use the 48.48 assay
+     flag=9, # the column containing flags
+     feature=5, # column containing gene names
+     type=6, # column containing probe type
+     Ct=7, # column containing the actual values
+     skip=24, # additional argument to read.table; the number of lines to
skip before hte data starts. Probably different in your file!
+     sep=",") # additional argument to read.table; the data is "," not
"\t" separated.
Warning message:
In readCtData(files = "1131135160.csv", n.features = 2304, flag = 9,  :
  2304 gene names (rows) expected, got 2314
show(test)
An object of class "qPCRset"
Size:  2304 features, 1 samples
Feature types:		 NRC, Test
Feature names:		 T09F3_3_HK T09F3_3_HK F33H1_2_HK ...
Feature classes:
Feature categories:	 OK
Sample names:		 1131135160 NA NA ...
head(flag(test))
  1131135160
1       Pass
2       Pass
3       Pass
4       Pass
5       Pass
6       Pass
table(flag(test))

Fail Pass
 835 1469

> Also, if you could respond to the post on the forum (with potential
> answers to the questions about technical versus biological, and potential
> normalizers, that would be great).

I'm not exactly sure what your question here is. You say that you're not
clear about technical versus biological reps. Do you mean how many you
need, whether to normalise one kind different from another or something
else?

Cheers
\Heidi

>> Hi Simon,
>>
>> thanks for your interest in HTqPCR. In principle, there's no reason why it
>> shouldn't be applicable to Fluidigm data. I used it briefly at some point
>> for one of the 96.96 assays, and I've also had other questions from people
>> on the Bioconductor list using HTqPCR for Fluidigm. I've been meaning to
>> add a section about that to the HTqPCR vignette, but didn't have time for
>> the before the previous bioconductor release.
>>
>> For specifying the flags, you can use the flag(<your_object>) <-
>> <your_flags> syntax. The flags have to be in a data.frame, but apart from
>> that they can be anything. Although I must admit that I never really
>> tested this... E.g. for changing just the first column, if flags are
>> already present:
>>
>> data(qPCRraw)
>> head(flag(qPCRraw))
>>  sample1 sample2 sample3 sample4 sample5 sample6
>> 1  Passed  Passed  Passed Flagged  Passed  Passed
>> 2  Passed  Passed  Passed Flagged  Passed  Passed
>> 3  Passed  Passed  Passed  Passed  Passed  Passed
>> 4  Passed  Passed  Passed  Passed  Passed  Passed
>> 5  Passed  Passed  Passed  Passed  Passed  Passed
>> 6 Flagged Flagged Flagged Flagged Flagged Flagged
>> new.flags <- data.frame(rep("OK", 384))
>> flag(qPCRraw)[,1] <- new.flags
>> head(flag(qPCRraw))
>>  sample1 sample2 sample3 sample4 sample5 sample6
>> 1      OK  Passed  Passed Flagged  Passed  Passed
>> 2      OK  Passed  Passed Flagged  Passed  Passed
>> 3      OK  Passed  Passed  Passed  Passed  Passed
>> 4      OK  Passed  Passed  Passed  Passed  Passed
>> 5      OK  Passed  Passed  Passed  Passed  Passed
>> 6      OK Flagged Flagged Flagged Flagged Flagged
>>
>> Alternatively, it should be possible to read them into a qPCRset object at
>> the same time as your data. My mac has done strange things to the line
>> ending in your file, so this might not work directly for you, but
>> something along the lines of:
>>
>> test <- readCtData(files="1131135160.csv", # file name
>>    n.features=2304, # if you use the 48.48 assay
>>    flag=9, # the column containing flags
>>    feature=5, # column containing gene names
>>    type=6, # column containing probe type
>>    Ct=7, # column containing the actual values
>>    skip=24, # additional argument to read.table; the number of lines to
>> skip before hte data starts. Probably different in your file!
>>    sep=",") # additional argument to read.table; the data is "," not "\t"
>> separated.
>>
>> HTH
>> \Heidi

>>> I'm trying to use the package developed by Heidi Dvinge (HTqPCR), and
I've
>>> a couple of basic questions. The package seems very nice, however, I'm
>>> using a high throughput real-time PCR system, the Biomark from Fluidigm,
>>> which can carry out upto 40,000 PCR's in a single run (approximately).
>>> Most of my runs use about 10,000 PCRs, with multiple technical
replicates,
>>> and occasionally I use several plates if I'm assying many samples, and
>>> dozens of genes.
>>>
>>> Typically, I've been using the geNorm approach for identifying potential
>>> normalizers in a suite of genes in any given biological paradigm, as a
>>> priori its always best to identify a gene which doesnt change in the
>>> question of interest (as opposed to using something like actin, and
>>> assuming it doesnt change). Is there a way of identifying the "best"
>>> normalizer in HTqPCR? Also, I'm having trouble specifying flags. I
>>> understand how to designate the appropriate column in the file as the
flag
>>> column, however, I dont understand how to specify the pass/fail criteria.
>>> In the Biomarks output, each PCR reaction is specified as "Pass" or
>>> "Fail". I've been over the vignette for HTqPCR, and cant see how to
define
>>> the "Faliure" criteria.
>>>
>>> In addition, I'm a bit fuzzy on how HtQPCR handles technical reps as
>>> opposed to biological reps. There is some discussion of this in the
>>> vignette, but I'm having trouble following it.
>>>
>>> Here are a few lines of our typical output from our Biomark machine in
>>> case this helps.
>>>
>>>
>>> Chip Run Filename
>>>
>>> X:\Real-time studies\Biomark\HK_Mar2010\ChipRun.bml
>>>
>>>
>>>
>>>
>>>
>