I am analysing an Affymetrix microarray dataset using Affy and limma but have encountered a problem with my code.
The experiment has two genotypes (RAB and N50) and two treatments (C and D), with three reps of each combination. I would like to know 1) which genes are differentially expressed due to treatment, 2) which genes are differentially expressed due to genotype and 3) which genes are differentially expressed due to a genotype*treatment interaction.
Firstly, I have normalised and corrected for background etc using gcrma and written a mydata file which is a matrix of expression data for all genes for each of the 12 samples.
I have also made a targets file, which specifies the file name meanings as below:
FileName Genotype Treatment N50C2 N50 C N50C4 N50 C N50C7 N50 C N50D1 N50 D N50D4 N50 D N50D8 N50 D RABC3 RAB C RABC8 RAB C RABC12 RAB C RABD2 RAB D RABD9 RAB D RABD10 RAB D > targets <- readTargets("targets.txt") > GT <- paste(targets$Genotype, targets$Treatment, sep=".") > design <- model.matrix(~Genotype*Treatment) > Genotype <- factor(targets$Genotype, levels=c("N50", "RAB")) > Treatment <- factor(targets$Treatment, levels=c("C", "D")) > design <- model.matrix(~Genotype*Treatment) > fit <- lmFit(eset, design) > colnames(design)  "(Intercept)" "Genotype1" "Treatment1"  "Genotype1:Treatment1"
My question is how to extract the contrasts which I am interested in at this point? From reading the manual, I understand I need to make a contrast matrix but I am unsure how to do this. Once I have this line of code, I think the next steps are:
> fit2 <- contrasts.fit(fit, cont.matrix) > fit2 <- eBayes(fit2) > topTable (fit2, coef=1) # and so on for other coefficients
Is this correct? Any help you could give me would be very much appreciated. If I need to give any more information then please do let me know.