Calculating Contrast Matrix in Limma for Differential Expression
0
0
Entering edit mode
4.1 years ago
sp29 ▴ 50

About the sample: There are 40 samples in which 10 samples are healthy controls and 30 samples are diseased (10 of type-1 + 10 of type-2 + 10 of type-3). I wish to perform differential expression analysis between the 30 diseased and 10 controls, i.e. 30_diseased_vs_10_control. The accession id of the sample is GSE65127, through which I have retrieved the expression-set (using exprs(x)), pheno-data (using pData(x)) and feature-data (using fData(x)) respectively.

Reproducible Code for differential expression that I am using

e_data <- read.csv("GSE65127_expression_set.csv", row.names = 1)

p_data <- read.csv("GSE65127_pheno_set.csv",row.names = 1)

f_data <- read.csv("GSE65127_feature_set.csv",row.names = 1)

design_labels <- p_data$source_name_ch1

design_factors <- factor(design_labels)

design <- model.matrix(~0+design_factors)

colnames(design) <- levels(design_factors)

fit <- lmFit(e_data, design)

cont.matrix <- makeContrasts(disease_vs_healthy = "(Lesional+Non.Lesional+Peri-Lesional)-(Healthy.Volunteer)",levels=design)

Error: While making the cont.matrix, I am getting the following message

Error in makeContrasts(disease_vs_healthy = "(Lesional+Non.Lesional+Peri-Lesional)-(Healthy.Volunteer)", :

The levels must by syntactically valid names in R, see help(make.names). Non-valid names: Healthy Volunteer,Non Lesional,Peri-Lesional

I am using this workflow: https://www.bioconductor.org/help/course-materials/2005/BioC2005/labs/lab01/estrogen/

limma R Micro-array Contrast-Matrix DEG • 2.0k views
ADD COMMENT
1
Entering edit mode

Rewrite your design_factors; eg, use design_factors <- factor(make_names(design_labels)). The column names of your design matrix are not the same as the levels used in your makeContrasts string.

ADD REPLY
0
Entering edit mode

Hi! thank you, it worked. However, after running the testing with the following code cont.matrix <- makeContrasts(disease_vs_healthy = (Lesional+Non.Lesional+Peri.Lesional)-(Healthy.Volunteer),levels=design)

fit2 <- contrasts.fit(fit, cont.matrix)

fit2 <- eBayes(fit2)

results <- decideTests(fit2)

summary(results)

All my genes are coming to be upregulated in summary()

disease_vs_healthy

Down 0

NotSig 0

Up 54675

When I do the same testing with Geo2R, I get significant results.

ADD REPLY
0
Entering edit mode

Is that surprising though? Your code works, but it is wrong. Let's suppose I'm a gene that is expressed at the same level in your four groups. Say my fitted value for treatment-group G, is x + (noise_G). With your contrast as currently defined, the value of your contrast would be (x + noise_lesional) + (x + noise_non_lesional) + (x + noise_peri_lesional) - (x + noise_healthy) = 2x + some noise. If x is appreciably bigger than the inter-group noise, this contrast will be significant. Fix your contrast.

ADD REPLY
0
Entering edit mode

Remove the quotation marks within the contrast function and please show the content of design.

ADD REPLY
0
Entering edit mode

Hi! Thanks for the quick reply, I have tried by removing the "" but still got the same error

Here is the design

> design

Healthy Volunteer Lesional Non Lesional Peri-Lesional

1 1 0 0 0

2 1 0 0 0

3 1 0 0 0

4 1 0 0 0

5 1 0 0 0

6 1 0 0 0

7 1 0 0 0

8 1 0 0 0

9 1 0 0 0

10 1 0 0 0

11 0 1 0 0

12 0 1 0 0

13 0 1 0 0

14 0 1 0 0

15 0 1 0 0

16 0 1 0 0

17 0 1 0 0

18 0 1 0 0

19 0 1 0 0

20 0 1 0 0

21 0 0 1 0

22 0 0 1 0

23 0 0 1 0

24 0 0 1 0

25 0 0 1 0

26 0 0 1 0

27 0 0 1 0

28 0 0 1 0

29 0 0 1 0

30 0 0 1 0

31 0 0 0 1

32 0 0 0 1

33 0 0 0 1

34 0 0 0 1

35 0 0 0 1

36 0 0 0 1

37 0 0 0 1

38 0 0 0 1

39 0 0 0 1

40 0 0 0 1

attr(,"assign")

[1] 1 1 1 1

attr(,"contrasts")

attr(,"contrasts")$design_factors

[1] "contr.treatment"

I am new to Biostars. Hope this is the correct format for pasting matrix results

ADD REPLY

Login before adding your answer.

Traffic: 2289 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