Лимма для DEG из данных микрочипов
Я использую Limma для анализа DEG моих данных "Данные по экспрессии из опухолевых эндотелиальных клеток гепатоцеллюлярной карциномы", присоединение GEO (GSE51401). Я хочу выяснить разницу в экспрессии опухолевых и неопухолевых клеток из GSE51401, поэтому я использую файлы.CEL для анализа опухолевых и неопухолевых клеток, всего 48 образцов (24 опухолевых и 24 неопухолевых опухоль), я не беру остатки опухоли и неопухолевые образцы.
ниже приведен код R, который я использовал, чтобы найти DEG, используя пакет limma,
library(affy)
library(limma)
library(gplots)
pd <- read.AnnotatedDataFrame("tnt.txt",header = T,sep = "\t")
pData(pd)
my_expressionset <- ReadAffy(filenames =rownames(pData(pd)))
my_expressionset
eset <- rma(my_expressionset)
pData(eset)
# design matrix
design <- model.matrix(~Substance,pData(pd))
design
fit <- lmFit(eset, design)
fit
fit2 <- eBayes(fit)
fit2
dim(fit2)
options(digits = 2)
fit2$p.value
deg <- topTable(fit2, coef = 2, p.value = 0.001, adjust.method = 'fdr',
number = nrow(eset))
dim(deg)[1]
dim(deg)
head(deg)
##heatmap of DEG
idx = rownames(deg)
eset[idx]
eset2 <- exprs(eset[idx])
head(eset2)
hr <- hclust(as.dist(1-cor(t(eset2), method="pearson")),method="complete")
heatmap(eset2, Rowv = as.dendrogram(hr),Colv = NA)
у меня есть вывод моей матрицы дизайна следующим образом,
(Intercept) SubstanceTEC
GSM1244729_BH2010140_5_NT105.CEL 1 1
GSM1244733_BH2010140_5_OT105.CEL 1 1
GSM1244737_BH2010140_5_QT105.CEL 1 1
GSM1244741_BH2010140_5_RT105_2.CEL 1 1
GSM1244746_BH2010338_3_1T105.CEL 1 1
GSM1244749_BH2010338_3_2T105.CEL 1 1
GSM1244754_BH2010338_3_4T105.CEL 1 1
GSM1244758_BH2010338_3_5T105.CEL 1 1
GSM1244730_BH2010140_5_NT31.CEL 1 1
GSM1244734_BH2010140_5_OT31.CEL 1 1
GSM1244738_BH2010140_5_QT31.CEL 1 1
GSM1244742_BH2010140_5_RT31_2.CEL 1 1
GSM1244745_BH2010338_3_1T31.CEL 1 1
GSM1244750_BH2010338_3_2T31.CEL 1 1
GSM1244755_BH2010338_3_4T31.CEL 1 1
GSM1244757_BH2010338_3_5T31.CEL 1 1
GSM1244763_BH2010610-8-8T31+yuan.CEL 1 1
GSM1244765_BH2010610-8-10T31+yuan.CEL 1 1
GSM1244771_BH2010610-8-16T31+yuan.CEL 1 1
GSM1244774_BH2010610-8-18T31+yuan.CEL 1 1
GSM1244779_BH2010610-8-20T31+yuan.CEL 1 1
GSM1244781_BH2010610-8-46T31+.CEL 1 1
GSM1244786_BH2010610-8-58T31+.CEL 1 1
GSM1244789_BH2010610-8-62T31+.CEL 1 1
GSM1244731_BH2010140_5_NP105.CEL 1 0
GSM1244735_BH2010140_5_OP105.CEL 1 0
GSM1244739_BH2010140_5_QP105.CEL 1 0
GSM1244744_BH2010140_5_RP105.CEL 1 0
GSM1244747_BH2010338_3_1P105.CEL 1 0
GSM1244752_BH2010338_3_2P105.CEL 1 0
GSM1244753_BH2010338_3_4P105.CEL 1 0
GSM1244760_BH2010338_3_5P105.CEL 1 0
GSM1244732_BH2010140_5_NP31.CEL 1 0
GSM1244736_BH2010140_5_OP31.CEL 1 0
GSM1244740_BH2010140_5_QP31.CEL 1 0
GSM1244743_BH2010140_5_RP31.CEL 1 0
GSM1244748_BH2010338_3_1P31.CEL 1 0
GSM1244751_BH2010338_3_2P31.CEL 1 0
GSM1244756_BH2010338_3_4P31.CEL 1 0
GSM1244759_BH2010338_3_5P31.CEL 1 0
GSM1244761_BH2010610-8-8P31+yuan.CEL 1 0
GSM1244767_BH2010610-8-10P31+yuan.CEL 1 0
GSM1244772_BH2010610-8-16P31+yuan.CEL 1 0
GSM1244773_BH2010610-8-18P31+yuan.CEL 1 0
GSM1244780_BH2010610-8-20P31+yuan.CEL 1 0
GSM1244784_BH2010610-8-46N31+.CEL 1 0
GSM1244788_BH2010610-8-58N31+yuan.CEL 1 0
GSM1244791_BH2010610-8-62N31+.CEL 1 0
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$Substance
[1] "contr.treatment"
так что мой вопрос: 1. мой код является правильным? 2. моя матрица дизайна верна? 3. Должен ли я сделать контрастную матрицу для анализа? 4. из анализа я получил всего 5201 DEG, так ли это количество DEG правильно? Я в замешательстве, так как получаю столько генов, сколько выражено почтением, так ли это нормально, что такое количество генов выражается почтительно? 5. Затем я должен найти количество генов с повышенной активностью из DEG.
Пожалуйста, направьте меня с этим, Спасибо, я прикрепил ссылку на тепловую карту, которую я сгенерировал, используя R