I am uncertain. Below is my r code:
> library(GEOquery)
> my_id <- "GSE30929"
> gse <- getGEO(my_id)
> length(gse)
[1] 1
> gse <-gse[[1]]
> gse
gs0 <- pData(gse)$`subtype:ch1`
> table(gs0)
gs0
dedifferentiated myxoid myxoid/round cell
40 17 11
pleomorphic well-differentiated
20 52
>
> keep <- c("dedifferentiated", "well-differentiated")
> sample_idx <- which(gs0 %in% keep)
>
> gse <- gse[, sample_idx]
>
> gs <- pData(gse)$`subtype:ch1`
> table(gs)
gs
dedifferentiated well-differentiated
40 52
> pData(gse)$data_processing[1]
summary(exprs(gse))
library(dplyr)
> sampleInfo <- pData(gse)
> head(sampleInfo)
table(sampleInfo$characteristics_ch1.1)
sampleInfo <- select(sampleInfo, characteristics_ch1.1)
> sampleInfo <- rename(sampleInfo, sample = characteristics_ch1.1)
>
> head(sampleInfo)
dim(sampleInfo)
[1] 92 1
> sampleInfo$sample
sampleInfo <- sampleInfo[,-1]
> sampleInfo
library(limma)
> design <- model.matrix(~0 + sampleInfo$group)
design
cutoff <- median(exprs(gse))
>
> ## TRUE or FALSE for whether each gene is "expressed" in each sample
> is_expressed <- exprs(gse) > cutoff
>
> ## Identify genes expressed in more than 2 samples
>
> keep <- rowSums(is_expressed) > 3
>
> ## check how many genes are removed / retained.
> table(keep)
## subset to just those expressed genes
> gse <- gse[keep,]
> aw <- arrayWeights(exprs(gse),design)
> aw
fit <- lmFit(exprs(gse), design,
+ weights = aw)
head(fit$coefficients)
fit2 <- contrasts.fit(fit, contrasts)
fit2 <- eBayes(fit2)
topTable(fit2)
top_genes <- topTable(fit2, coef = 1, n = 2000)
write.csv(top2000, "dge_data.csv", row.names = FALSE)
Hi,
https://www.ncbi.nlm.nih.gov/geo/geo2r/?acc=GSE30929
I have performed DGE analysis between WDLPS and DDLPS tumour samples in R and have the following results table. I ended up with the same results as the GEO2R.
logFC AveExpr t P.Value adj.P.Val B
207175_at -5.7 8.8 -14 3.1e-25 7.0e-21 47
213706_at -3.5 8.9 -14 1.0e-24 1.1e-20 46
203548_s_at -4.7 10.1 -14 2.8e-24 2.1e-20 45
203549_s_at -4.3 10.1 -14 9.0e-24 4.4e-20 43
205913_at -4.5 9.9 -14 9.9e-24 4.4e-20 43
218736_s_at -2.5 8.7 -13 4.3e-23 1.6e-19 42
203571_s_at -4.0 9.3 -13 2.0e-22 6.3e-19 40
49452_at -3.0 9.1 -13 4.3e-22 1.2e-18 40
43427_at -2.8 8.9 -13 5.3e-22 1.3e-18 39
219398_at -3.3 9.6 -13 9.3e-22 2.1e-18 39
However, does the logFC of -5.7 in the first row of the table indicate that the gene is unregulated or down regulated in WDLPS or DDLPS?
Thanks!
1 answer
It will depend upon how did you defined two groups. Would you please tell which is your control and which is test?
I am not sure, I think DDLPS might be the test and WDLPS the control. Is this correct? I just performed a DGE analyses between the two tumour subtypes.
Please show me your design table. It looks like you are comparing between subtype: dedifferentiated vs subtype: well-differentiated.
head(design)
sample id: ALT476 sample id: ALT479 sample id: ALT486 sample id: ALT497
1 0 0 0 0
2 0 0 0 0
3 0 0 0 0
4 0 0 0 0
5 0 0 0 0
6 0 0 0 0
sample id: ALT500 sample id: DD0728 sample id: DD103 sample id: DD112
1 0 1 0 0
2 0 0 1 0
3 0 0 0 1
4 0 0 0 0
5 0 0 0 0
6 0 0 0 0
sample id: DD172 sample id: DD176 sample id: DD188 sample id: DD216
1 0 0 0 0
2 0 0 0 0
3 0 0 0 0
4 1 0 0 0
5 0 1 0 0
6 0 0 1 0
sample id: DD226 sample id: DD330 sample id: DD344 sample id: DD3516
1 0 0 0 0
2 0 0 0 0
3 0 0 0 0
4 0 0 0 0
5 0 0 0 0
6 0 0 0 0
sample id: DD382 sample id: DD391 sample id: DD442 sample id: DD454
1 0 0 0 0
2 0 0 0 0
3 0 0 0 0
4 0 0 0 0
5 0 0 0 0
6 0 0 0 0
sample id: DD456 sample id: DD474 sample id: DD506 sample id: DD513
1 0 0 0 0
2 0 0 0 0
3 0 0 0 0
4 0 0 0 0
5 0 0 0 0
6 0 0 0 0
sample id: DD595 sample id: DD597 sample id: DD610 sample id: DD611
1 0 0 0 0
2 0 0 0 0
3 0 0 0 0
4 0 0 0 0
5 0 0 0 0
6 0 0 0 0
sample id: DD6213 sample id: DD639 sample id: DD668 sample id: DD670
1 0 0 0 0
2 0 0 0 0
3 0 0 0 0
4 0 0 0 0
5 0 0 0 0
6 0 0 0 0
sample id: DD672 sample id: DD69 sample id: DD712 sample id: DD714
1 0 0 0 0
2 0 0 0 0
3 0 0 0 0
4 0 0 0 0
5 0 0 0 0
6 0 0 0 0
sample id: DD715 sample id: DD718 sample id: DD719 sample id: DD7533
1 0 0 0 0
2 0 0 0 0
3 0 0 0 0
4 0 0 0 0
5 0 0 0 0
6 0 0 0 0
sample id: DD8379 sample id: DD88 sample id: ML138 sample id: ML161
1 0 0 0 0
2 0 0 0 0
3 0 0 0 0
4 0 0 0 0
5 0 0 0 0
6 0 0 0 0
sample id: ML163 sample id: ML165 sample id: ML169 sample id: ML253
1 0 0 0 0
2 0 0 0 0
3 0 0 0 0
4 0 0 0 0
5 0 0 0 0
6 0 0 0 0
sample id: ML257 sample id: ML283 sample id: ML292 sample id: ML298
1 0 0 0 0
2 0 0 0 0
3 0 0 0 0
4 0 0 0 0
5 0 0 0 0
6 0 0 0 0
sample id: ML304 sample id: ML414 sample id: ML593 sample id: ML618
1 0 0 0 0
2 0 0 0 0
3 0 0 0 0
4 0 0 0 0
5 0 0 0 0
6 0 0 0 0
attr(,"assign")
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$`factor(phenoData$characteristics_ch1.2)`
[1] "contr.treatment"
Your design in not correct. I suggest you to invest your time to learn this. Please follow the tutorial I have linked here which will help you understand how it is done-
DDLPS WDLPS
GSM766533 1 0
GSM766534 1 0
GSM766535 1 0
GSM766536 1 0
GSM766537 1 0
GSM766538 1 0
GSM766539 1 0
GSM766540 1 0
GSM766541 1 0
GSM766542 1 0
GSM766543 1 0
GSM766544 1 0
GSM766545 1 0
GSM766546 1 0
GSM766547 1 0
GSM766548 1 0
GSM766549 1 0
GSM766550 1 0
GSM766551 1 0
GSM766552 1 0
GSM766553 1 0
GSM766554 1 0
GSM766555 1 0
GSM766556 1 0
GSM766557 1 0
GSM766558 1 0
GSM766559 1 0
GSM766560 1 0
GSM766561 1 0
GSM766562 1 0
GSM766563 1 0
GSM766564 1 0
GSM766565 1 0
GSM766566 1 0
GSM766567 1 0
GSM766568 1 0
GSM766569 1 0
GSM766570 1 0
GSM766571 1 0
GSM766572 1 0
GSM766621 0 1
GSM766622 0 1
GSM766623 0 1
GSM766624 0 1
GSM766625 0 1
GSM766626 0 1
GSM766627 0 1
GSM766628 0 1
GSM766629 0 1
GSM766630 0 1
GSM766631 0 1
GSM766632 0 1
GSM766633 0 1
GSM766634 0 1
GSM766635 0 1
GSM766636 0 1
GSM766637 0 1
GSM766638 0 1
GSM766639 0 1
GSM766640 0 1
GSM766641 0 1
GSM766642 0 1
GSM766643 0 1
GSM766644 0 1
GSM766645 0 1
GSM766646 0 1
GSM766647 0 1
GSM766648 0 1
GSM766649 0 1
GSM766650 0 1
GSM766651 0 1
GSM766652 0 1
GSM766653 0 1
GSM766654 0 1
GSM766655 0 1
GSM766656 0 1
GSM766657 0 1
GSM766658 0 1
GSM766659 0 1
GSM766660 0 1
GSM766661 0 1
GSM766662 0 1
GSM766663 0 1
GSM766664 0 1
GSM766665 0 1
GSM766666 0 1
GSM766667 0 1
GSM766668 0 1
GSM766669 0 1
GSM766670 0 1
GSM766671 0 1
GSM766672 0 1
attr(,"assign")
[1] 1 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"
Is this design better? Thanks!
It looks better than what you had before.
Log in to answer this question.