[BioC] Limma- deleting control probes from differential expression analysis
Neel Aluru
naluru at whoi.edu
Thu Dec 10 22:37:08 CET 2009
Hello,
I am doing analysis on the two color data set with common reference sample (Agilent). I am able to generate list of differentially expressed genes but having trouble to remove control probes from the list. I am having difficulty in understanding "asMatrixWeights()" command. I have column in my data titled "ControlType" with 1 or -1 corresponding to control probes and I want to eliminate them from my analysis. "ControlType =0" are the actual probes. Is there anyway to write "isGene" command so that only ControlType=0 can be considered for fitting subsequent models.
I also have a text file with the list of various control probes. Is there any way I can pull the info from there and assign them as not relevant for analysis.
> library(limma)
> getwd()
[1] "/Users/Neel/agilent/limma"
> setwd("/Users/Neel/agilent/limma")
> getwd()
[1] "/Users/Neel/agilent/limma"
> targets = readTargets()
> targets = readTargets("Targets.txt", row.names = "Name")
> RG = read.maimages(targets, source ="agilent")
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 44407 gene
Found 14 BLANK
Found 604 Blank
Found 320 positive
Found 153 negative
Found 130 flag1
Found 120 flag2
Setting attributes: values Color
> 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
> MA = normalizeWithinArrays(RG, method="loess", weights=NULL)
> 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))
> fit = lmFit(MA, design)
Warning message:
Partial NA coefficients for 1269 probe(s)
> cont.matrix = makeContrasts(PCBvsCONT=PCB-CONT, levels=design)
> fit2 = contrasts.fit(fit, cont.matrix)
> fit2 = eBayes(fit2)
> topTable(fit2, number = 10, adjust="BH")
Row Col Start Sequence ProbeUID ControlType ProbeName GeneName
25521 306 66 0 GCTAACATGGTCTCCAATGGCAGTTTGTTGTTGCATAGCATATTTATGATGTTGGCAAAA 17654 0 A_15_P100578 cyp1a
26068 313 23 0 GCTAACATGGTCTCCAATGGCAGTTTGTTGTTGCATAGCATATTTATGATGTTGGCAAAA 17654 0 A_15_P100578 cyp1a
22261 267 65 0 GTCTGTGTTTTTGGGCTAGAGATGTACAGAGATGATTCTTATGCTTAAAAGGAGGATTTT 16235 0 A_15_P116040 si:dkey-13a21.15
25361 304 72 0 ACCAAAAGAGTGTGAAGGTAAAAGATTGAAGACTGTGTTGGGTTTTTATTGGATTCGTGA 17593 0 A_15_P119653 drd3
37362 448 41 0 CTGAACTAAGAAGTGGCTACCAGCACTGTCAAATAAGCATATACAGGAAAAACTAAATAA 5780 0 A_15_P102530 ugt1aa
32555 391 7 0 GCTTCATCTCAGATTGCAGCAGTGTTATGTGTTCTATGCATTGAAATAAAATAATGCACA 20037 0 A_15_P118356 gabpb2
28154 338 38 0 ACCAGTAACCACAGGATTAAGTTCATTACCATCACCGTATCTTGCTCTCCTGAAGTCCAT 18652 0 A_15_P105185 elp4
3413 41 65 0 GGTTTATTGTTGATTATGTAGATCCCCATTAGCCACTATCAACATAGTGGTTACTCGTCC 3040 0 A_15_P107325 hsf1
34336 412 26 0 ATAAATCCTGTGAATCTCAATACAGTGTCTGAAAACATGTGTCCATTTTTACTCTGTTGC 20460 0 A_15_P102917 wdr12
34569 415 6 0 CCCCCCAAAATGCTGTAAATTACTATATTTTCATTTAAAATTTGGAAAGAATTGGTATCT 12610 0 A_15_P108015 mapk14b
SystematicName
25521 NM_131879
26068 NM_131879
22261 NM_001083843
25361 BC076352
37362 BC100055
32555 NM_212712
28154 NM_001017638
3413 NM_131600
34336 NM_199557
34569 NM_178223
Description logFC
25521 "ref|Danio rerio cytochrome P450, family 1, subfamily A (cyp1a), mRNA [NM_131879]" -5.231141
26068 "ref|Danio rerio cytochrome P450, family 1, subfamily A (cyp1a), mRNA [NM_131879]" -5.094103
22261 "ref|Danio rerio si:dkey-13a21.15 (si:dkey-13a21.15), mRNA [NM_001083843]" 2.994005
25361 "gb|Danio rerio dopamine receptor D3, mRNA (cDNA clone MGC:92896 IMAGE:7118973), complete cds. [BC076352]" -2.144771
37362 "gb|Danio rerio UDP glucuronosyltransferase 1 family a, a, mRNA (cDNA clone IMAGE:7284571), partial cds. [BC100055]" -2.557859
32555 "ref|Danio rerio GA repeat binding protein, beta 2 (gabpb2), mRNA [NM_212712]" 2.211286
28154 "ref|Danio rerio elongation protein 4 homolog (S. cerevisiae) (elp4), mRNA [NM_001017638]" 3.755883
3413 "ref|Danio rerio heat shock transcription factor 1 (hsf1), mRNA [NM_131600]" -3.005454
34336 "ref|Danio rerio WD repeat domain 12 (wdr12), mRNA [NM_199557]" 1.511836
34569 "ref|Danio rerio mitogen-activated protein kinase 14b (mapk14b), mRNA [NM_178223]" -2.155142
AveExpr t P.Value adj.P.Val B
25521 12.597193 -30.025678 3.565209e-09 0.0001534394 2.82418014
26068 12.481079 -26.545758 9.021446e-09 0.0001941325 2.76354911
22261 5.474685 10.205127 2.777888e-05 0.3985158471 0.52258102
25361 5.733082 -6.881963 1.611209e-04 0.6678374033 0.30028005
37362 6.635666 -7.874792 1.369250e-04 0.6678374033 -0.01742959
32555 5.982746 6.187841 3.238173e-04 0.7066158126 -0.08376019
28154 5.480283 7.643937 1.638003e-04 0.6678374033 -0.08963515
3413 4.849755 -7.322943 2.117418e-04 0.6951770618 -0.19736474
34336 7.124021 5.990976 3.987536e-04 0.7066158126 -0.20474318
34569 4.560762 -7.065892 2.618034e-04 0.7066158126 -0.29029278
Once I got this, now I am trying to remove control spots... Any ideas will be appreciated.
Thank you,
Neel
Neel Aluru
Postdoctoral Scholar
Biology Department
Woods Hole Oceanographic Institution
Woods Hole, MA 02543
USA
508-289-3607
More information about the Bioconductor
mailing list