limma package for DEGs analysis
1
0
Entering edit mode
@lawardeankita1-14844
Last seen 10 months ago
Estonia

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

 

limma • 1.3k views
ADD COMMENT
0
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

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.

 

ADD REPLY
0
Entering edit mode
@gordon-smyth
Last seen 58 minutes ago
WEHI, Melbourne, Australia

There are actually no errors here from the limma package. You are just making some programming errors in R after you've run limma. It's not my responsibility to debug your code, but here is a hint. If a matrix X has only one column, then trying to access the second column by X[,2] will give an error. That is why fit3$p.value[,2] gives an error.

Note also that you never need to call p.adjust() yourself, because limma does this for you with topTable() etc. Anyway, I don't think you are running p.adjust() correctly -- you are asking it do perform Bonferroni adjustment, which I don't think is really what you want.

ADD COMMENT
0
Entering edit mode

Hello Sir,

Thank you sir, the explanation was very helpful.

ADD REPLY

Login before adding your answer.

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