Another limma factorial that needs controlling and pairing
0
0
Entering edit mode
k. brand ▴ 420
@k-brand-1874
Last seen 9.6 years ago
Dear BioC's, I'd very much appreciate input on my design to achieve a sound analysis with limma. Initially i thought i could get away with a two-factor design: -treatments, 'S', 'E' & 'L' -tissues, 'R' & 'C' (targets file at bottom) With 6 groups i adhered to the example in the limma user guide sec.8.7, thus- targets <- readTargets("targets_no_time_no_animal.txt") TP <- paste(targets$Tissue, targets$Pperiod, sep = ".") TP <- factor(TP, levels = c("R.S", "C.S", "R.E", "C.E", "R.L", "C.L")) design <- model.matrix(~0 + TP) colnames(design) <- levels(TP) fit <- lmFit(rma.pp, design) cont.matrix <- makeContrasts( R.12t12vsR.6t18 = R.E - R.S, R.12t12vsR.18t6 = R.E - R.L, R.6t18vsR.18t6 = R.S - R.L, C.12t12vsC.6t18 = C.E - C.S, C.12t12vsC.18t6 = C.E - C.L, C.6t18vsC.18t6 = C.S - C.L, R.6t18vsC.6t18 = R.S - C.S, R.12t12vsC.12t12 = R.E - C.E, long.int.tis = (R.E - R.S) - (C.E - C.S), short.int.tis = (R.E - R.L) - (C.E - C.L), levels = design) fit2 <- contrasts.fit(fit, cont.matrix) fit3 <- eBayes(fit2) results <- (decideTests(fit3, method="global", adjust.method="BH", p.value=0.05)) etc. This worked fine for my contrasts of interest until discussing with a statistician (thank you Renee) just what i was trying to get away with- within each of the 6 groups are 16 samples collected over a 24 hour period. These group time course's focus on another aspect of this experiment(*fyi see at bottom). What i'm trying to acheive now with limma is not concerned with 'time', merely identifying differentially expressed genes (thus not changing in time) between the 6 groups. However i understand this factor needs to be controlled for, not ignored as i did above. And this is where my R & limma skills prove inadequate- controlling for the factor, Time. Scouring the user guide and BioC list for examples of what i was attempting, i tried- targets <- read.delim("RNA_Targets.txt") Tissue <- factor(targets$Tissue) Pperiod <- factor(targets$Pperiod) Time <- factor(targets$Time) design <- model.matrix(~ Tissue + Pperiod + Time, data=targets) colnames(design) fit <- lmFit(rma.pp, design) fit2 <- eBayes(fit) topTable(fit2, coef=2) _But_ i've failed to extract my contrasts of interest as detailed in the first design. My primary appeal for help is then- how to obtain my contrast's of interest? Either from the second design, and/or how to setup an appropriate design which incorporates factor=time so that i can control for this? Additionally, in attempting to really do the 'right thing' setting up this design, there is the 4th factor- pairing of the tissue's. As is evident form the targets "RNA_Targets.txt", each 'R' & 'C' pair come from the same animal. I expect incorporating this pairing into the model would add power. Any help/suggestions here would also be greatly appreciated. Limma was built to easily handle all these parameters, but the limma/R expertise in this corner of the world are just too thin for me to take advantage of it w/o further help. I look forward to any suggestions with thanks in advance, cheers, Karl "targets_no_time_no_animal.txt"- same as "RNA_Targets.txt" w/o columns 'Time' & 'Animal' "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 *fyi- The group time course's, consisisting of 16 samples collected every 1.5 hours over 24 hours, attempt to identify 'circadian expression profiles' by F-testing the 16 samples for goodness-of-fit to a sinus curve. There is only one replicate per time point, relying instead on the density of samples, as opposed to replication. Clearly it would be ideal to incorporate this model into limma, as i understand is possible. But in failing to employ limma for 'routine' DE identification can't even dream of attempting such in the time available. Anyone who could do this would be a co-author on the paper. And *VERY* popular in the circadian field. -- Karl Brand <k.brand-ampersand-erasmusmc.nl> Department of Genetics Erasmus MC Dr Molewaterplein 50 3015 GE Rotterdam lab +31 (0)10 704 3409 fax +31 (0)10 704 4743 mob +31 (0)642 777 268
limma limma • 813 views
ADD COMMENT

Login before adding your answer.

Traffic: 905 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6