Hello,
Iam Using affy and lima packages for differential analysis of gene for my data
i have 10 .CEL (Adult_A172.CEL, Adult_LN229.CEL, Adult_SF268.CEL, Adult_U87MG.CEL, Adult_U118MG.CEL, Pediatric_KNS42.CEL, Pediatric_Res186.CEL, Pediatric_Res259.CEL, Pediatric_SF188.CEL, Pediatric_UW479.CEL) files. this data i obtaied from the ArrayExpress website https://www.ebi.ac.uk/arrayexpress/experiments/E-TABM-579/
i have refered to one of the example from bioconductor and the code is given below
but i don't understand why its not working as I am creating the design matrix in similar way. once i fir the data using lmFit till top table its showing no error but i want to select p values less than 0.05 and after doing p value commad its showing error. and i want to crete a heatmap of those Deferentially expressed genes.
i have refered to one of the example from bioconductor and the code is given below
library(affy)
library(limma)
affy_data <- ReadAffy()
eset_rma <- rma(affy_data)
exprSet.nologs <- exprs(eset_rma)
head(exprSet.nologs)
# List the column (chip) names
colnames(exprSet.nologs)
#log transformation
exprSet = log(exprSet.nologs, 2)
head(exprSet)
#To print out our expression matrix
write.table(exprSet, file="MGMT_matrix.txt", quote=F, sep="\t")
XX<- as.matrix(read.table("MGMT_matrix.txt",header = TRUE,sep = "\t"))
head(XX)
minimalSet <- ExpressionSet(assayData = XX)
minimalSet
ex <- exprs(minimalSet)
head(ex)
groups<-as.factor(c(rep("Adult",5),rep("Pediatric",5)))
groups
design_new<-model.matrix(~0+groups)
colnames(design_new)=levels(groups)
design_new
fit<-lmFit(minimalSet, design_new)
head(fit)
fit
cont.matrix<-makeContrasts(Adult-Pediatric, levels=design_new)
cont.matrix
fit2<-contrasts.fit(fit, cont.matrix)
fit2
ebfit<-eBayes(fit2)
head(ebfit)
ebfit
tt_new <-topTable(ebfit, coef=1)
tt_new
select_new <- p.adjust(ebfit$p.value)<0.05 ####error here
select_new
gg <- minimalSet[select_new,]
gg
#tt_row <-nrow(topTable(ebfit, coef=1, lfc=1))
#tt_row
Another way i tried to do the same exercise usig the following code,
data_age <- ReadAffy()
eset <- rma(data_age)
pData(eset)
# matrix design
design.mat <- cbind(c(1,1,1,1,1,0,0,0,0,0),c(0,0,0,0,0,1,1,1,1,1))
colnames(design.mat)<- c("Adult","Pediatric")
design.mat
#design matrix: model.matrix
sample_10 <- factor(rep(c("Adult","Pediatric"),each=5))
design.mat<- model.matrix(~0+sample_10)
colnames(design.mat)<- levels(sample_10)
design.mat
#contrast matrix: makeContrasts()
contrast.mat <- makeContrasts(Diff = Adult-Pediatric,levels = design.mat)
contrast.mat
#fit limma
fit_data <- lmFit(eset, design.mat)
fit_data
fit_data2 <- contrasts.fit(fit_data,contrast.mat)
fit_data2
fit3 <- eBayes(fit_data2)
fit3
dim(fit3)
selected <- p.adjust(fit3$p.value[,2])<0.001 ### error
selected
deg <- topTable(fit3,coef = 'Diff',p.value=0.05, adjust.method='fdr',lfc = log2(1.5),number = nrow(eset))
dim(deg)[1]
#length(deg)
#fit4 <- exprs(fit3)
it again is not woring
Please be more specific by providing the output of the function that apparently is not doing what you expect it do do, because "it is not working" is not really informative... and people don't have the time to reproduce your lines of code without having access to the raw data.
In addition, please realize that the RMA algorithm returns normalized data that is already log2 transformed. So there is no need to do this again.
Hello Sir,
I'm so sorry for the incomplete message, i couldn't find any attachment option to attach the raw data. i want to find the diferentially expressed genes of my data.
today i tried to do do this DEGs analysis using the new_pheno. txt file such that
new_pheno.txt file format is as follows
sample Substance
A172 sensitive
LN229 sensitive
SF268 sensitive
Res259 sensitive
U87MG resistant
U118MG resistant
KNS42 resistant
i have given output of all the commands.
using the following code,
library(affy)
library(limma)
cels<-ReadAffy()
cels
pheno<-read.AnnotatedDataFrame("new_pheno.txt",header = TRUE,sep = "\t")
pheno
eset<-justRMA(phenoData = pheno)
eset
ExpressionSet (storageMode: lockedEnvironment)
assayData: 54675 features, 10 samples
element names: exprs, se.exprs
protocolData
sampleNames: A172.CEL KNS42.CEL ... UW479.CEL (10 total)
varLabels: ScanDate
varMetadata: labelDescription
phenoData
sampleNames: A172.CEL KNS42.CEL ... UW479.CEL (10 total)
varLabels: Substance
varMetadata: labelDescription
design<-model.matrix(~Substance,pData(eset))
design
(Intercept) Substancesensitive
A172.CEL 1 1
KNS42.CEL 1 1
LN229.CEL 1 1
Res186.CEL 1 1
Res259.CEL 1 0
SF188.CEL 1 0
SF268.CEL 1 0
U118MG.CEL 1 0
U87MG.CEL 1 0
UW479.CEL 1 0
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$Substance
[1] "contr.treatment"
fit<-lmFit(eset,design)
fit
ebase<-eBayes(fit)
ebase
ebase$p.value
(Intercept) Substancesensitive
1007_s_at 1.793398e-10 4.408522e-01
1053_at 3.946614e-11 6.189166e-01
select<-p.adjust(ebase$p.value[,1]) < 0.001
select
select_2 <- p.adjust(ebase$p.value[,2]) < 0.001
head(select_2) # this gives false as there are no values less than 0.001 in the second column
esetsel<-eset[select,]
esetsel
ExpressionSet (storageMode: lockedEnvironment)
assayData: 50413 features, 10 samples
element names: exprs, se.exprs
protocolData
sampleNames: A172.CEL KNS42.CEL ... UW479.CEL (10 total)
varLabels: ScanDate
varMetadata: labelDescription
phenoData
sampleNames: A172.CEL KNS42.CEL ... UW479.CEL (10 total)
de<-exprs(esetsel)
> dim(de)
[1] 50413 10
tt <- topTable(ebase,coef=2, number = nrow(eset))
tt
library(gplots)
a <- de[1:100,1:10]
heatmap.2(a)
i hope this will help to understand the problem.
thank you so much for considering to help me earlier.