Question: Limma- deleting control probes from differential expression analysis
0
gravatar for Neel Aluru
9.9 years ago by
Neel Aluru450
United States
Neel Aluru450 wrote:
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
transcription assign • 658 views
ADD COMMENTlink written 9.9 years ago by Neel Aluru450
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 122 users visited in the last hour