
> Hi Mark,
> 
> I was able to run RP from the M values but have couple of doubts. Thank you for suggesting that to me. I really appreciate if you could help me out with it. First question - is it possible to use only significantly differentially expressed genes obtained from the linear model (LIMMA ) for RP analysis. I was unable to select the subset of genes from topTable command.  Second question is something very simple but I am unable to get it to work. I want to see the gene ID/gene names in the RP output tables. For some reason it is unable to recognize the gene names column from the file. Session info is pasted below.
> 
> Thank you, Mark. 
> 
> 
> Sincerely, Neel
> 
> 
> Session info
> 
> > library(limma)
> > getwd()
> [1] "/Users/Neel"
> > setwd("/Users/Neel/agilent/limma")
> > getwd()
> [1] "/Users/Neel/agilent/limma"
> >  targets = readTargets()
> > targets = readTargets("Targets.txt", row.names = "Name")
> > spottypes = readSpotTypes()
> > myfun = function(x) as.numeric(x$ControlType > -50.5)
> > RG = read.maimages(targets, source="agilent", wt.fun=myfun)
> Read conta.txt 
> Read contb.txt 
> Read contc.txt 
> Read contd.txt 
> Read pcba.txt 
> Read pcbb.txt 
> Read pcbc.txt 
> Read pcbd.txt 
> > spottypes = readSpotTypes()
> > RG$genes$Status = controlStatus(spottypes, RG)
> Matching patterns for: ProbeName GeneName 
> Found 42990 gene 
> Found 14 BLANK 
> Found 604 Blank 
> Found 320 positive 
> Found 153 negative 
> Found 130 flag1 
> Found 120 flag2 
> Setting attributes: values Color 
> > RG.b = backgroundCorrect(RG, method="normexp", offset=50)
> Green channel
> Corrected array 1 
> Corrected array 2 
> Corrected array 3 
> Corrected array 4 
> Corrected array 5 
> Corrected array 6 
> Corrected array 7 
> Corrected array 8 
> Red channel
> Corrected array 1 
> Corrected array 2 
> Corrected array 3 
> Corrected array 4 
> Corrected array 5 
> Corrected array 6 
> Corrected array 7 
> Corrected array 8 
> > MA.p = normalizeWithinArrays(RG.b, method="loess")
> > MA.pAq = normalizeBetweenArrays(MA.p, method = "Aquantile")
> > design = cbind(CONT=c(0,0,0,0,1,1,1,1), PCB=c(1,1,1,1,0,0,0,0))
> > isGene = RG$genes$Status=="gene"
> > fit = lmFit(MA.p[isGene,], design)
> > cont.matrix = makeContrasts(PCBvsCONT= CONT-PCB, levels=design)
> > fit2 = contrasts.fit(fit, cont.matrix)
> > fit2 = eBayes(fit2)
> > top1 = topTableF(fit2, number=300, genelist=fit$genes, adjust.method="BH", sort.by="F", p.value=1)
> > write.table(top1, file= "Fstatistic.txt", quote=FALSE, sep="\t", row.names=FALSE, col.names=TRUE)
> > cl = c(rep(0,4), rep(1,4))
> > RP.out = RP(MA.p, cl, logged=TRUE, rand=123)
> Rank Product analysis for two-class case 
>  
> Starting 100 permutations... 
> Computing pfp .. 
> Outputing the results .. 
> > topGene (RP.out, cutoff=0.05, logged = TRUE)
> Table1: Genes called significant under class1 < class2 
> 
> Table2: Genes called significant under class1 > class2 
> 
> $Table1
>      gene.index  RP/Rsum FC:(class1/class2)    pfp P.value
> [1,]      25521   1.1388             0.0296 0.0000       0
> [2,]      26068   1.8476             0.0335 0.0000       0
> [3,]      21312 153.5698             0.3200 0.0000       0
> [4,]       2577 303.9498             0.3469 0.0200       0
> [5,]       6418 330.7369             0.3833 0.0200       0
> [6,]       7471 338.5279             0.3641 0.0167       0
> [7,]      31524 415.9073             0.3912 0.0314       0
> 
> $Table2
>      gene.index  RP/Rsum FC:(class1/class2)    pfp P.value
> [1,]       4970 138.9352             3.5370 0.0000       0
> [2,]       6578 152.4629             3.5384 0.0000       0
> [3,]       9065 177.9296             3.3864 0.0000       0
> [4,]       8198 336.6878             2.9670 0.0200       0
> [5,]       5300 369.8948             2.8839 0.0280       0
> [6,]      15935 393.9698             2.7581 0.0350       0
> [7,]       7262 445.7020             2.5756 0.0486       0
> [8,]       3089 471.5363             2.6675 0.0488       0
> > names(MA.p)
> [1] "weights" "targets" "genes"   "source"  "M"       "A"      
> > topGene (RP.out, num.gene = 10, logged = TRUE, gene.names = MA.p$genes)
> Warning: gene.names should have the same length as the gene vector. 
> No gene.names are assigned 
> 
> 
> On Dec 17, 2009, at 7:09 PM, Mark Cowley wrote:
> 
>> no worries Neel,
>> you should be able to run RP using the M values, ie the R-G log2 ratios that you fitted the linear models to.
>> mark
>> 
>> On 17/12/2009, at 2:07 PM, Neel Aluru wrote:
>> 
>>> Hi Mark,
>>> 
>>> Thanks for your mail. I was able to figure out how to get the replicate values. However, I haven't tried working on rank product analysis yet. Will get to it in a couple of days. Thanks a lot for enquiring about my progress. The mailing list is really helpful to learn R packages.
>>> 
>>> Sincerely, Neel
>>> 
>>> At 07:57 PM 12/16/2009, you wrote:
>>>> hi Neel,
>>>> how did you go with this problem? let me know if you still need it  
>>>> answering
>>>> mark
>>>> On 11/12/2009, at 1:45 PM, Neel Aluru wrote:
>>>> 
>>>>> Hello,
>>>>> 
>>>>> I am hoping someone will be able to help me with this. I am  
>>>>> analyzing the two color common reference design agilent arrays using  
>>>>> Limma. I have four replicates and two treatments(control and PCB). I  
>>>>> am trying to extract fold change values (PCB/control) for the four  
>>>>> replicates separately. lmfit() command after normalizing between  
>>>>> arrays is pooling all the samples and I get one fold change value.  
>>>>> Is there anyway I can get the values for individual replicates  
>>>>> separtely? I want to use the replicate values in Rank Product  
>>>>> analysis. Any help will be really helpful.
>>>>> 
>>>>> Sorry for bothering everyone with my questions. I have made a lot of  
>>>>> progress since I started using this package and it is mainly due to  
>>>>> this mailing list.
>>>>> Sincerely, Neel
>>>>> 
>>>>> 
>>>>> Neel Aluru
>>>>> Postdoctoral Scholar
>>>>> Biology Department
>>>>> Woods Hole Oceanographic Institution
>>>>> Woods Hole, MA 02543
>>>>> USA
>>>>> 508-289-3607
>>>>> 
>>>>> _______________________________________________
>>>>> Bioconductor mailing list
>>>>> Bioconductor@stat.math.ethz.ch
>>>>> https://stat.ethz.ch/mailman/listinfo/bioconductor
>>>>> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor
>>>> 
>>> Neelakanteswar Aluru Ph.D.
>>> Post doctoral Scholar
>>> Biology Department
>>> Redfield 304 (MS#32)
>>> Woods Hole Oceanographic Institution
>>> Woods Hole MA 02543 USA
>>> Phone: (508) 289-3607 [Office]
>>> 774-392-3727 [Cell] 
>>> RID: A-7237-2009
>>> 
>> 
> 
> Neel Aluru
> Postdoctoral Scholar
> Biology Department
> Woods Hole Oceanographic Institution
> Woods Hole, MA 02543
> USA
> 508-289-3607
> 
> 
> 

Neel Aluru
Postdoctoral Scholar
Biology Department
Woods Hole Oceanographic Institution
Woods Hole, MA 02543
USA
508-289-3607




	[[alternative HTML version deleted]]

