Entering edit mode
Hello,
I have two questions with regard to the AgiMicroRna package and they
are FilterMicroRna and SignificantMicroRna functions. I was able to
figure most of the package easily except these two functions. If any
one can give me some hints or suggestions to get these two functions
working, I will really appreciate your help. I think that the
alternative way I used to filter control probes is causing trouble
with the SignficantMicroRna function.
I used this package to analyze Agilent Zebrafish miRNA arrays which
were scanned using GenePix 4000B scanner and data extracted using
Feature Extraction Software 9.5.3. I have posted the session info
below (two functions that need attention are hightlighted).
Session Info
> library("AgiMicroRna")
> targets.micro=readTargets(infile="targets.txt", verbose=TRUE)
Target File
FileName Treatment GErep Subject
36_DMSO_1 36_DMSO_1.txt 36DMSO 1 1
36_DMSO_2 36_DMSO_2.txt 36DMSO 1 2
36_DMSO_3 36_DMSO_3.txt 36DMSO 1 3
36_TCDD_1 36_TCDD_1.txt 36TCDD 2 1
36_TCDD_2 36_TCDD_2.txt 36TCDD 2 2
36_TCDD_3 36_TCDD_3.txt 36TCDD 2 3
60_DMSO_1 60_DMSO_1.txt 60DMSO 3 1
60_DMSO_2 60_DMSO_2.txt 60DMSO 3 2
60_DMSO_3 60_DMSO_3.txt 60DMSO 3 3
60_TCDD_1 60_TCDD_1.txt 60TCDD 4 1
60_TCDD_2 60_TCDD_2.txt 60TCDD 4 2
60_TCDD_3 60_TCDD_3.txt 60TCDD 4 3
> dd.micro=read.maimages(targets.micro$FileName,
columns=list(R="gTotalGeneSignal",G=
"gTotalProbeSignal",Rb="gMeanSignal", Gb="gProcessedSignal"),
annotation=c("ProbeUID","ControlType","ProbeName","GeneName","Systemat
icName",
"sequence", "accessions","probe_mappings",
"gIsGeneDetected","gIsSaturated","gIsFeatNonUnifOL",
"gIsFeatPopnOL","chr_coord","gBGMedianSignal","gBGUsed"))
Read 36_DMSO_1.txt
Read 36_DMSO_2.txt
Read 36_DMSO_3.txt
Read 36_TCDD_1.txt
Read 36_TCDD_2.txt
Read 36_TCDD_3.txt
Read 60_DMSO_1.txt
Read 60_DMSO_2.txt
Read 60_DMSO_3.txt
Read 60_TCDD_1.txt
Read 60_TCDD_2.txt
Read 60_TCDD_3.txt
> cvArray(dd.micro, "MeanSignal", targets.micro, verbose=TRUE)
Foreground: MeanSignal
FILTERING BY ControlType FLAG
RAW DATA: 5335
PROBES without CONTROLS: 4620
----------------------------------
(Non-CTRL) Unique Probe: 490
(Non-CTRL) Unique Genes: 231
----------------------------------
DISTRIBUTION OF REPLICATED NonControl Probes
reps
5 6 7 10
20 18 36 416
------------------------------------------------------
Replication at Probe level- MEDIAN CV
36_DMSO_1 36_DMSO_2 36_DMSO_3 36_TCDD_1 36_TCDD_2 36_TCDD_3 60_DMSO_1
60_DMSO_2 60_DMSO_3
0.078 0.081 0.091 0.081 0.077 0.067
0.076 0.066 0.103
60_TCDD_1 60_TCDD_2 60_TCDD_3
0.073 0.086 0.069
------------------------------------------------------
DISTRIBUTION OF REPLICATED Noncontrol Genes
reps
20
231
------------------------------------------------------
> ddTGS.rma = rmaMicroRna(dd.micro, normalize=TRUE, background=FALSE)
Calculating Expression
> ddPROC = filterMicroRna(ddTGS.rma, dd.micro, control = TRUE,
IsGeneDetected = TRUE, wellaboveNEG = FALSE, limIsGeneDetected = 50,
limNEG = 25, makePLOT = FALSE, targets.micro, verbose = TRUE)
FILTERING PROBES BY FLAGS
FILTERING BY ControlType
Error in matrix(ddFILT$other$gIsGeneDetected, nrow = dim(ddFILT)[1],
ncol = dim(ddFILT)[2]) :
attempt to set an attribute on NULL
> MMM = ddTGS.rma$Rb
> colnames(MMM) = colnames(dd.micro$Rb)
> maintitle='TGS.rma'
> colorfill='blue'
> ddaux=ddTGS.rma
> ddaux$G=MMM
> mvaMicroRna(ddaux, maintitle, verbose=TRUE)
------------------------------------------------------
mvaMicroRna info:
FEATURES : 231
POSITIVE CTRL: 12
NEGATIVE CTRL: 7
STRUCTURAL: 3
> rm(ddaux)
> RleMicroRna(MMM,"RLE TGS.rma", colorfill)
> boxplotMicroRna(MMM, maintitle, colorfill)
> plotDensityMicroRna(MMM, maintitle)
> spottypes = readSpotTypes()
> ddTGS.rma$genes$Status = controlStatus(spottypes, ddTGS.rma)
Matching patterns for: ProbeName GeneName
Found 231 gene
Found 1 BLANK
Found 1 Blank
Found 0 blank
Found 6 positive
Found 0 negative
Found 0 flag1
Found 0 flag2
Found 6 flag3
Found 5 flag4
Found 1 flag5
Setting attributes: values
> i = ddTGS.rma$genes$Status == "gene"
> esetPROC = esetMicroRna(ddTGS.rma[i,], targets.micro,
makePLOT=TRUE, verbose = TRUE)
outPUT DATA: esetPROC
Features Samples
231 12
> design=model.matrix(~-1+treatment)
> print(design)
treatment36DMSO treatment36TCDD treatment60DMSO treatment60TCDD
1 1 0 0 0
2 1 0 0 0
3 1 0 0 0
4 0 1 0 0
5 0 1 0 0
6 0 1 0 0
7 0 0 1 0
8 0 0 1 0
9 0 0 1 0
10 0 0 0 1
11 0 0 0 1
12 0 0 0 1
attr(,"assign")
[1] 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$treatment
[1] "contr.treatment"
> fit=lmFit(esetPROC, design)
> cont.matrix = makeContrasts(treatment36TCDDvstreatment36DMSO =
treatment36TCDD-treatment36DMSO, treatment60TCDDvstreatment60DMSO =
treatment60TCDD-treatment60DMSO,treatment60TCDDvstreatment36TCDD =
treatment60TCDD-treatment36TCDD, treatment60DMSOvstreatment36DMSO =
treatment60DMSO-treatment36DMSO, levels=design)
> print(cont.matrix)
Contrasts
Levels treatment36TCDDvstreatment36DMSO
treatment60TCDDvstreatment60DMSO
treatment36DMSO -1
0
treatment36TCDD 1
0
treatment60DMSO 0
-1
treatment60TCDD 0
1
Contrasts
Levels treatment60TCDDvstreatment36TCDD
treatment60DMSOvstreatment36DMSO
treatment36DMSO 0
-1
treatment36TCDD -1
0
treatment60DMSO 0
1
treatment60TCDD 1
0
> fit2 = contrasts.fit(fit,cont.matrix)
> print(head(fit2$coeff))
Contrasts
treatment36TCDDvstreatment36DMSO
treatment60TCDDvstreatment60DMSO
dre-let-7a 0.038640984
0.013333873
dre-let-7b 0.074038749
-0.031608286
dre-let-7c 0.026244357
-0.005682488
dre-let-7d 0.067340768
0.055567054
dre-let-7e 0.004569306
0.136348664
dre-let-7f 0.042880109
0.085568058
Contrasts
treatment60TCDDvstreatment36TCDD
treatment60DMSOvstreatment36DMSO
dre-let-7a 1.7358343
1.76114142
dre-let-7b 0.1366920
0.24233899
dre-let-7c 0.9920976
1.02402449
dre-let-7d 0.8098432
0.82161694
dre-let-7e 0.1186829
-0.01309647
dre-let-7f 1.1245878
1.08189990
> fit2 = eBayes(fit2)
> fit2 = basicLimma(esetPROC, design, cont.matrix, verbose = TRUE)
DATA
Features Samples
231 12
> DE = getDecideTests(fit2, DEmethod = "separate", MTestmethod =
"BH", PVcut = 0.1, verbose = TRUE)
------------------------------------------------------
Method for Selecting DEGs: separate
Multiple Testing method: BH - pval 0.1
treatment36TCDDvstreatment36DMSO
treatment60TCDDvstreatment60DMSO
UP 0 5
DOWN 0 1
treatment60TCDDvstreatment36TCDD
treatment60DMSOvstreatment36DMSO
UP 56 51
DOWN 80 91
------------------------------------------------------
> pvalHistogram(fit2, DE, PVcut = 0.1, DEmethod ="separate",
MTestmethod="BH",cont.matrix, verbose= TRUE)
> significantMicroRna(esetPROC, ddTGS.rma, targets.micro, fit2,
cont.matrix, DE, DEmethod = "separate", MTestmethod= "BH", PVcut =
0.1, Mcut=0, verbose=TRUE)
------------------------------------------------------
CONTRAST: 1 - treatment36TCDDvstreatment36DMSO
Error in data.frame(PROBE_ID, as.character(GENE_ID),
as.character(chr_coord), :
arguments imply differing number of rows: 231, 0
Thank you very much,
Sincerely, Neel
Neel 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
[[alternative HTML version deleted]]