[BioC] Another limma factorial that needs controlling and pairing
Sabrina AT Case
hxs61 at case.edu
Mon Aug 16 15:33:22 CEST 2010
Hi, Karl:
Thank you for the detailed explanation. I need to understand the
Treatment-contrast parametrization since I always used the design
matrix as what you did the first time. I need to digest more about
it.
Just a quick comment. I noticed that you used
decideTests(data, methods="global",adj="BH"...)
In Limma's user's guide $10.3, it mentioned that there is no theory
prove about the combination of BH and global will correctly control
the FDR for negatively correlated contrasts (though simulation is ok),
I wonder if you have tried other ways, such as separate topTable etc.
Does it make difference? Thanks!
For some reason, I posted the same question as you suggested to the
mailinglist, it never got posted, I hope this time it works,,,
Sabrina
On Wed, Aug 11, 2010 at 9:59 AM, Karl Brand <k.brand at erasmusmc.nl> wrote:
> Hi Sabrina,
>
> Yes, i believe i did figure out a contrast matrix which allowed me to use
> limma to indentify DE genes between two factors of interest ('Tissue' and
> 'Pperiod' in my case) whilst controlling for a third factor ('Time' in my
> case where the samples were collected at different times of the day as part
> of another question incorporated in the same experiment but not discussed
> here). So to be clear, my goal was to find DE genes between all pairwise
> comparisons of 6 data sets, namely, R.S, R.E, R.L, C.S, C.E and C.L and also
> for any genes significant for interaction between 'Tissue' and 'Pperiod'.
>
> In my original post "[BioC] Another limma factorial that needs controlling
> and pairing" i'd outlined an approach based closly on the first example
> discussed in Gorodn Smyth's chapter 23*, pg412. In the end this approach
> didn't work for me. Instead i succeeded using the "treatment-contrast"
> parametrization (pg414) example. You'll find a representation of the code i
> used for my analysis below. I haven't run it, but for the purpose of
> illustration, together with my targets file should suffice in explaining
> exactly how i used limma. The most diffciult part for me was understanding
> the way R sets up this "treatment-contrast" parametrization. To be sure
> about how i was setting up the contrasts, i manually calculated the Log2FC
> for each contrast and compared this with the limma output. This confirmation
> was critical for a dilatante like myself to be sure my costrasts were indeed
> what i thought they were. I would advise similar such corss checking of your
> contrasts.
>
> As you'll see from my code, i used 'decide tests' since overall significance
> of genes across all contrasts was of greater help to me in understanding the
> biology. As to which correction method you use, there is plenty of online
> info (BioC archives, limma user guide) to clarify this. In short, "global"
> seemed the best balanced method for my question(s).
>
> Final note on controlling for the paring of tissues (being obtained from the
> same animal): i omitted this for simplicty in the code below. In fact Prof.
> Smyth clarified out to deal with this using duplicateCorrelation in a later
> BioC post if this is relevant for you; and i can also repost if needed.
>
> hth,
>
> Karl
>
> *http://www.statsci.org/smyth/pubs/limma-biocbook-reprint.pdf
>
> ###not run
> library(limma)
> targets <- read.delim("RNA_Targets.txt")
> Tissue <- factor(targets$Tissue, levels = c("R", "C"))
> Pperiod <- factor(targets$Pperiod, levels = c("E", "L", "S"))
> Time <- factor(targets$Time)
>
> design <- model.matrix(~Tissue * Pperiod + Time)
> colnames(design)
>
> fit <- lmFit(rma.pp, design)
> cont.matrix <- cbind(
>
> R.E_R.S = c(0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
> R.E_R.L = c(0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
> R.S_R.L = c(0,0,-1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
>
> C.E_C.S = c(0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1),
> C.E_C.L = c(0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0),
> C.S_C.L = c(0,0,-1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,1),
>
> R.S_C.S = c(0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1),
> R.E_C.E = c(0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
> R.L_C.L = c(0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0),
>
> R.E_R.S__C.E_C.S = c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1),#int
> R.E_R.L__C.E_C.L = c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0))#int.2
>
> fit2 <- contrasts.fit(fit, cont.matrix)
> fit3 <- eBayes(fit2)
>
> #write resutls of contrasts for 0.05
> results <- (decideTests(fit3, method="global", adjust.method="BH",
> p.value=0.05))
> write.fit(fit3, results, file="w.int.glob.BH.0.05.csv", sep=",")
>
> #write resutls of contrasts for 0.01
> results <- (decideTests(fit3, method="global", adjust.method="BH",
> p.value=0.01))
> write.fit(fit3, results, digits=10, file="w.int.glob.BH.0.01.csv", sep=",")
>
> ###end
>
> "RNA_Targets.txt":
>
> FileName Tissue Pperiod Time Animal
> 01file.CEL R S 1 1
> 02file.CEL C S 1 1
> 03file.CEL R S 2 2
> 04file.CEL C S 2 2
> 05file.CEL R S 3 3
> 06file.CEL C S 3 3
> 07file.CEL R S 4 4
> 08file.CEL C S 4 4
> 09file.CEL R S 5 5
> 10file.CEL C S 5 5
> 11file.CEL R S 6 6
> 12file.CEL C S 6 6
> 13file.CEL R S 7 7
> 14file.CEL C S 7 7
> 15file.CEL R S 8 8
> 16file.CEL C S 8 8
> 17file.CEL R S 9 9
> 18file.CEL C S 9 9
> 19file.CEL R S 10 10
> 20file.CEL C S 10 10
> 21file.CEL R S 11 11
> 22file.CEL C S 11 11
> 23file.CEL R S 12 12
> 24file.CEL C S 12 12
> 25file.CEL R S 13 13
> 26file.CEL C S 13 13
> 27file.CEL R S 14 14
> 28file.CEL C S 14 14
> 29file.CEL R S 15 15
> 30file.CEL C S 15 15
> 31file.CEL R S 16 16
> 32file.CEL C S 16 16
> 33file.CEL R E 1 17
> 34file.CEL C E 1 17
> 35file.CEL R E 2 18
> 36file.CEL C E 2 18
> 37file.CEL R E 3 19
> 38file.CEL C E 3 19
> 39file.CEL R E 4 20
> 40file.CEL C E 4 20
> 41file.CEL R E 5 21
> 42file.CEL C E 5 21
> 43file.CEL R E 6 22
> 44file.CEL C E 6 22
> 45file.CEL R E 7 23
> 46file.CEL C E 7 23
> 47file.CEL R E 8 24
> 48file.CEL C E 8 24
> 49file.CEL R E 9 25
> 50file.CEL C E 9 25
> 51file.CEL R E 10 26
> 52file.CEL C E 10 26
> 53file.CEL R E 11 27
> 54file.CEL C E 11 27
> 55file.CEL R E 12 28
> 56file.CEL C E 12 28
> 57file.CEL R E 13 29
> 58file.CEL C E 13 29
> 59file.CEL R E 14 30
> 60file.CEL C E 14 30
> 61file.CEL R E 15 31
> 62file.CEL C E 15 31
> 63file.CEL R E 16 32
> 64file.CEL C E 16 32
> 65file.CEL R L 1 33
> 66file.CEL C L 1 33
> 67file.CEL R L 2 34
> 68file.CEL C L 2 34
> 69file.CEL R L 3 35
> 70file.CEL C L 3 35
> 71file.CEL R L 4 36
> 72file.CEL C L 4 36
> 73file.CEL R L 5 37
> 74file.CEL C L 5 37
> 75file.CEL R L 6 38
> 76file.CEL C L 6 38
> 77file.CEL R L 7 39
> 78file.CEL C L 7 39
> 79file.CEL R L 8 40
> 80file.CEL C L 8 40
> 81file.CEL R L 9 41
> 82file.CEL C L 9 41
> 83file.CEL R L 10 42
> 84file.CEL C L 10 42
> 85file.CEL R L 11 43
> 86file.CEL C L 11 43
> 87file.CEL R L 12 44
> 88file.CEL C L 12 44
> 89file.CEL R L 13 45
> 90file.CEL C L 13 45
> 91file.CEL R L 14 46
> 92file.CEL C L 14 46
> 93file.CEL R L 15 47
> 94file.CEL C L 15 47
> 95file.CEL R L 16 48
> 96file.CEL C L 16 48
>
>
>
> On 8/10/2010 4:21 PM, sabrina wrote:
>>
>>
>>
>> Hi, Karl:
>> I have a similar question as yours, did you figure it out since your last
>> post?
>>
>> Other question I have is about using topTable. In your previous
>> design, did you use topTable for each contrast, meaning for each
>> particular
>> contrast, how you find the top list of genes? If that is the case, what is
>> the
>> best way to sort the genes? Did you use B value or adj.p? Or you just used
>> decideTests?
>>
>> Thanks
>>
>> Sabrina
>>
>> _______________________________________________
>> Bioconductor mailing list
>> Bioconductor at stat.math.ethz.ch
>> https://stat.ethz.ch/mailman/listinfo/bioconductor
>> Search the archives:
>> http://news.gmane.org/gmane.science.biology.informatics.conductor
>
> --
> Karl Brand <k.brand at erasmusmc.nl>
> Department of Genetics
> Erasmus MC
> Dr Molewaterplein 50
> 3015 GE Rotterdam
> P +31 (0)10 704 3409 | F +31 (0)10 704 4743 | M +31 (0)642 777 268
>
More information about the Bioconductor
mailing list