[BioC] HTqPCR

[Heidi Dvinge]

Hello Deborah,

> Good morning Heidi,
>
> thank you for your response.
>
> I have 6 samples with the "LightCycler" format (96 genes one each sample)
> : one control and its duplicate, one "thirty minutes" and its duplicate
> and one "three hours" and its duplicate.
> I did the ttestCtdata, the mannwhitneyCtData and the limmaCtdata on my
> samples for compare all the results.

> I did the commands :
>>PlaqueControle<-readCtData("essai_CT.txt",
>> column.info=list(feature="Name",position="Pos",Ct="Cp"), n.features=96,
>> format="LightCycler")
>>Plaque30min<-readCtData("essai_30min.txt",
>> column.info=list(feature="Name",position="Pos",Ct="Cp"), n.features=96,
>> format="LightCycler")
>>Plaque3h<-readCtData("essai_3h.txt",
>> column.info=list(feature="Name",position="Pos",Ct="Cp"), n.features=96,
>> format="LightCycler")
>>PlaqueControleB<-readCtData("essai_CT_BIS.txt",
>> column.info=list(feature="Name",position="Pos",Ct="Cp"), n.features=96,
>> format="LightCycler")
>>Plaque30minB<-readCtData("essai_30min_BIS.txt",
>> column.info=list(feature="Name",position="Pos",Ct="Cp"), n.features=96,
>> format="LightCycler")
>>Plaque3hB<-readCtData("essai_3h_BIS.txt",
>> column.info=list(feature="Name",position="Pos",Ct="Cp"), n.features=96,
>> format="LightCycler")
>
>>essai<-cbind(PlaqueControle,Plaque30min,Plaque3h,PlaqueControleB,Plaque30minB,Plaque3hB)
>>files_essai<-read.table("files_essai.txt",header=T)

This all looks okay. Incidentally, you can also read in all the files
together, using readCtData(files=c("essai_CT.txt", "essai_30min.txt",
...), n.data=6), in case you don't what to do the cbind() afterwards.

>>essai.cat<- setCategory(essai, groups = files_essai$Treatment,quantile =
>> 0.8)
>>deltaCtnorm <- normalizeCtData(essai.cat, norm = "deltaCt",deltaCt.genes
>> =  c("gene85", "gene86", "gene87", "gene88", "gene89", "gene90",
>> "gene91"))

Just out of curiosity, did you try other normalisation methods as well?
I've seen a few cases where one of hte supposed housekeeper genes has been
really off the chart.
>
> Then I applied the ttestCtData, I compared the "control" samples with the
> "three hours" samples :
>>essai_ttest <- ttestCtData(deltaCtnorm[,c(1,3,6,4)], groups =
>> files_essai$Treatment[1:4],calibrator = "Controle")
>  
> mannwhitneyCtData on the same samples :
>>essaimwtest<-mannwhitneyCtData(deltaCtnorm[,c(1,3,6,4)], groups =
>> files_essai$Treatment[1:4],calibrator = "Controle")
>
Looks okay. I assume the order of your samples is different in your
qPCRset compared to your files_essai$Treatment.
>
> And finally, I did the limmaCtData :
>>essai_design<-model.matrix(~0 + files_essai$Treatment)
>>colnames(essai_design) <- c("Controle", "Trente_m", "Trois_h")
>>print(essai_design)
>
Again, your sample order in files_essai is definitely different from your
object, but you do seem to use c(1,3,6,4,2,5) further down, so I guess it
should work.

>>essai_contrasts <- makeContrasts(Trois_h-Controle,Trois_h - Trente_m,
>> Trente_m - Controle, (Trente_m +Trois_h)/2 - Controle, levels =
>> essai_design)
>>colnames(essai_contrasts) <- c("3h-CT", "3h-30min", "30min-CT",
>> "30min&3h-CT")
>>print(essai_contrasts)
>
>
>>deltaCtnorm2 <- deltaCtnorm[order(featureNames(deltaCtnorm)),]
>
You don't actually have to reorder when ndups=1. Of coruse, it doesn't
hurt though.

>>essai_limma <- limmaCtData(deltaCtnorm2[,c(1,3,6,4,2,5)], design =
>> essai_design,contrasts = essai_contrasts, ndups = 1)
>
>
> I didn't obtain the same conclusions for the same genes in the three
> tests. So I think that I did a lot of error on the commands...
>
>From just glancing over it, it looks okay, although I of course can't
actually see what your initial and resulting objects look like.

> I give you the head screenshots of my results for the three tests in
> attachments.
>
> For example, if I take the "gene10", with limmaCtData, I obtained that
> "gene10" in the "control" samples is significatively different to "gene10"
> in "three hours" samples (p-value = 0.0057). With ttestCtData, "gene10" is
> not significatively different to "gene 10" in the "three hours" samples
> (p-value = 0.063) ; ditto with mannwhitneyCtData (p-value =0.25).
>
> In the part 10.3 of "HTqPCR.pdf", you say about limma that "The result is
> a list with one component per comparison. Each component is similar to the
> result from using ttestCtData." ; so I suppose that my results are not
> consistent.
>
The results from limma are *in principle* similar to a t-test, i.e. each
component is the result from one of the individual tests specified by your
contrasts. However since you're using 3 different types of statistical
tests, it's not surprising that the results you get vary.

Mann-Whitney is mainly used when the data aren't normally distributed.
It's a non-parametric test, so it has less statistical power than the
other two. It's therefore not surprising that the p-values are less
significant.

limma uses a modified, more advanced version of the standard t-test. It
always considers data from all samples, not just the ones being compared
for any given contrast, and it "borrows" data across all features on your
assay to achieve a more robust estimate. It is thus expected to give more
significant values that a standard t-test, which only considers the
samples and genes being compared in each individual test in isolation.

There's no hard'n'fast rule for which of the 3 tests to use, since it also
depends on the (distribution of) your starting data etc. It also depends
on what you're using your data for. Do you want something where you're
100% sure that the genes do indeed differ between your conditions, e.g. as
a validation. Or is this preliminary data, that will be used for further
studies in the lab, so you're more interested in e.g. the top 1-10 hits.
Depending on this, you can always have a look at the actual data for some
of your genes (either the values themselves or some of the plots from
HTqPCR).

The tests mainly differ in how conservative they are. With this in mind,
you can check whether your results are relatively consistent by for
example comparing the order of the resulting p-value. For example, if you
plot the p-values from one test versus another, are they then completely
randomly scattered, or do they follow a general trend.

At the moment it sounds like your results mainly differ in the level of
the p-value, not whether a gene is e.g. up- or down-regulated in a given
comparison. If the latter is the case, then it sounds like there's
definitely something wrong, and I'll have to look into that.

I'm sorry that I can't give you any simple reply, or tell you which of the
3 tests to use. But it really depends on both the purpose of your study,
and your data.

\Heidi

> Furthermore, all my results with the mannwhitneyCtData are not
> significant...  I don't know if it was the good way to use these tests.
>
> Sincerely,
>
> Deborah.
>
>
> ________________________________
>  De : Heidi Dvinge <heidi at ebi.ac.uk>
> À : Deborah Ung <ung.deborah at yahoo.fr>
> Envoyé le : Vendredi 22 juin 2012 11h33
> Objet : Re: HTqPCR
>
> Hello Deborah,
>>
>> Good morning Heidi Dvinge,
>>
>> I am a french student and I am currently a trainee in a biotechnology
>> company. I would like to know if you have other documentations about
>> Limma
>> applied to HTqPCR because I have some problems to analyze my results.
>>
> I'm afraid the only documentation is the examples in the help files for
> limmaCtData, and the examples in the vignette section 10, i.e.
>
> ?limmaCtData
> openVignette("HTqPCR")
>
> If you can email the commands you used and the resulting error messages,
> (or why you think the results don't make sense), we can try to dissect the
> problem.
>
> Best
> \Heidi
>
>> Sincerely,
>>
>> Deborah.
>>