[BioC] single color using Limma : very low p-value "1.40456052121752e-60"

[Alex Gutteridge]

 Agreed. Specifying the second coefficient (which excludes the intercept 
 term) in topTable should do what you want:

 result=topTable(ebayes,coef=2,number=4608, p=0.05,adjust = "BH")

 Alex

 On Wed, 10 Aug 2011 13:21:00 -0400, Sean Davis wrote:
> I didn't go through the code in great detail, but since these are
> 1-color data, make sure that you are not testing that the intercept 
> is
> different than 0 (and it looks like you probably are).
>
> Sean
>
>
> On Wed, Aug 10, 2011 at 5:11 AM, bigoun <bigounn at yahoo.fr> wrote:
>> Dear all,
>>
>> I'm analysing data single color data (only intensities at F532) 
>> using limma.
>>
>> I have 14 arrays in one group and 5 in the other group. I want to 
>> compare them
>> so my targets file is as follows
>> FileName Mal Sain
>> M2.gpr 1 0
>> M3.gpr 1 0
>> M4.gpr 1 0
>> M5.gpr 1 0
>> M6.gpr 1 0
>> M7.gpr 1 0
>> M8.gpr 1 0
>> M9.gpr 1 0
>> M10.gpr 1 0
>> M11.gpr 1 0
>> M12.gpr 1 0
>> M13.gpr 1 0
>> M14.gpr 1 0
>> M15.gpr 1 0
>> S1.gpr 0 1
>> S2.gpr 0 1
>> S3.gpr 0 1
>> S4.gpr 0 1
>> S5.gpr 0 1
>>
>> I use some correction and normalisation methods as recommanded in 
>> limma user
>> guide. Afterthat, I make my design and as I have 2 technical 
>> replicate (2 spot
>> for the same miRNA) for each array, I wanted to compute 
>> duplicatecorrelation
>> coefficient
>>
>>
>> .GAL file as follows :
>> ATF 1
>> 52 5
>> Type=GenePix ArrayList V1.0
>> BlockCount=48
>> BlockType=0
>> URL=http://genome-www4.stanford.edu/cgi-bin/SGD/locus.pl?locus=[ID]
>> "Block1=2645, 2295, 100, 12, 300, 8, 300"
>> "Block2=7145, 2295, 100, 12, 300, 8, 300"
>> "Block3=11555, 2235, 100, 12, 300, 8, 300"
>> "Block4=16055, 2235, 100, 12, 300, 8, 300"
>>
>>
>> The  corfit$consensus is not good (0,36). Do I still use this 
>> parameters for
>> lmfit???
>>
>> Anyway, the more strange is with or without this parameters I have 
>> p-adjust
>> value very very low.
>>
>> first line of my result :
>> Block Row Column ID Name X.Intercept. factor.pData.population.sain 
>> AveExpr F P.Value adj.P.Val
>>
>> 39 2 10 42554 hsa-miR-923 13.4125293049973 -0.0517387609947947 
>> 13.3989138415776 25264.4487822867 1.20563134868457e-63 
>> 1.40456052121752e-60
>>
>>
>> Is anybody help me and what is wrong.
>>
>> As you can see in UNFILTERED section, I have less than 50% genes for 
>> the
>> analysis. Could it be the explanation of these strange result??
>>
>> Thanks in advance
>>
>> My srcipt:
>> library(limma)
>> myFilter = function(X) {
>> H_threshold=2
>> okFLAG = X$Flags > -49;
>> okm1 =abs(X[,"F532 Median"]-X[,"F532 Mean"])
>> okm2 = 0.5*(X[,"F532 Median"]+X[,"F532 Mean"])
>> okH = ((okm1)/okm2) <H_threshold
>> as.numeric(okFLAG & okH)}
>>
>> targets=readTargets(file="targets.txt", path=NULL, sep="\t")
>>
>> E 
>> =read.maimages(targets,source="genepix",wt.fun=myFilter,columns=list(E="F532
>> Mean",Eb="B532 Mean",Names="Name"))
>>
>> #Taux de filtrage
>> unFiltered = apply(E$weights,MARGIN=2,FUN=mean)
>> round(unFiltered,2)
>>
>> #  M2   M3   M4   M5   M6   M7   M8   M9  M10  M11  M12  M13  M14 
>>  M15   S1   S2
>>
>> #0.41 0.43 0.34 0.29 0.51 0.35 0.48 0.35 0.29 0.56 0.38 0.56 0.61 
>> 0.52 0.48 0.34
>>
>> #  S3   S4   S5
>> #0.36 0.46 0.47
>>
>> #Boite a moustache avant correction
>> boxplot(as.data.frame(E$E),main="Mean intensities")
>>
>> #Correction background - normexp
>> Enorm <- backgroundCorrect(E, method="normexp",offset=1)
>> boxplot(as.data.frame(Enorm$E), main="Mean intensities - normexp")
>> #Normalisation quantile
>> NormE <- normalizeBetweenArrays(Enorm,method="quantile")
>> boxplot(as.data.frame(NormE$E), main="Normalized intensities")
>>
>> NormE$E <- log2(NormE$E)
>>
>> pData <- data.frame(population = c('mal','mal', 'mal',
>> 
>> 'mal','mal','mal','mal','mal','mal','mal','mal','mal','mal','mal','sain','sain','sain','sain','sain'))
>>
>>
>> # create design matrix
>> design <- model.matrix(~factor(pData$population))
>>
>> # In my .gpr files all miRNAs contain the string "miR" in their name
>> # so the grep function can be used to extract all of these, removing
>> # all control signals and printing buffers etc.
>> # You need to check your .gpr files to find which signals you should 
>> extract.
>> miRs <- c(grep("-miR-", NormE$genes$Name), grep("-let-", 
>> NormE$genes$Name))
>> NormE.final <- NormE[miRs, ]
>> NormE.final <-NormE.final[order(NormE.final$genes$Name), ]
>> NormE.final$genes$Name
>>
>> # calculate duplicate correlation between the 2 replicates on each 
>> array
>> corfit <- duplicateCorrelation(NormE.final, design, ndups=2)
>> corfit$consensus
>> # show a boxplot of the correlations
>> boxplot(tanh(corfit$atanh.correlations))
>>
>> # fit the linear model, including info on duplicates and correlation
>> fit <- lmFit(NormE.final, design, ndups=2, 
>>  correlation=corfit$consensus)
>> #fit <- lmFit(NormE.final, design)
>> # calculate values using ebayes
>> ebayes <- eBayes(fit)
>>
>> # output a list of top differnetially expressed genes...
>> result=topTable(ebayes,number=4608, p=0.05,adjust = "BH")
>> genes=result$Status=="genes"
>> write.table(result,"gene.txt",quote=FALSE,sep="\t",row.names=FALSE,
>> col.names=TRUE)
>>
>>
>> mikel
>>        [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> Bioconductor mailing list
>> Bioconductor at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/bioconductor
>> Search the archives: 
>> http://news.gmane.org/gmane.science.biology.informatics.conductor
>>
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives:
> http://news.gmane.org/gmane.science.biology.informatics.conductor

-- 
 Alex Gutteridge