Yes, I think I start to understand howthe correlation works.

 

One last thing I am not sure I understand:

Previously, by doing analysis without separating male-female but treating
them together as 6-plicate ( the differences are very few anyway):

 

MA<-normalizeWithinArrays(RG,RG$printer,method="control",controlspots=CONTRO
Lpos2)

MA2<-normalizeBetweenArrays(MA,method="quantile")

 

fit<-lmFit(MA2)

fit<-eBayes(fit)

Final<-topTable(fit, number=389307,sort.by="M")

Final[1:5,]

                        ProbeID        M        A        t      P.Value
adj.P.Val        B

103220 CHR5:100995274-100995773 5.485304 12.64761 22.74918 1.380686e-10
1.126337e-06 14.71703

103219 CHR5:100995274-100995773 5.363803 12.20548 20.44722 4.343994e-10
1.126337e-06 13.68467

103209 CHR5:100993710-100995269 5.344564 12.01568 19.62738 6.732815e-10
1.126337e-06 13.28067

103214 CHR5:100993710-100995269 5.213038 12.00956 17.93719 1.760080e-09
1.236496e-06 12.37856

103211 CHR5:100993710-100995269 5.187434 11.67345 20.33639 4.604543e-10
1.126337e-06 13.63125

 

I was getting ca 31% of the spots statistically significant
(adj.p.value<=0.05)

 

Now, having specified the male/female design, the spots I get with
adj.p.value<=0.05 are “only” 23%

I get exactly the same numbers If I use a contrasts matrix(B) or not(A) as
you suggested:

 

design<-read.table("design.txt",sep="\t",header=TRUE,row.names=1)

design

# MALE FEMALE

#A_M    1      0

#B_M    1      0

#F_M    1      0

#C_F    0      1

#D_F    0      1

#E_F    0      1

 

 

#A

fit2 <- lmFit(MA2, design)

fit3<-eBayes(fit2)

FinalA<-topTable(fit2, number=389307,sort.by="M")

FinalA

                   ProbeID     MALE   FEMALE        A        F      P.Value
adj.P.Val

103220 CHR5:100995274-100995773 5.991492 4.979116 12.64761 392.7107
2.267845e-10 1.539717e-05

103222 CHR5:100995274-100995773 4.945359 3.989958 13.60286 364.9386
3.279371e-10 1.539717e-05

103211 CHR5:100993710-100995269 5.719710 4.655158 11.67345 313.0060
7.087562e-10 1.539717e-05

103219 CHR5:100995274-100995773 5.899835 4.827771 12.20548 308.5935
7.610504e-10 1.539717e-05

103209 CHR5:100993710-100995269 5.932802 4.756326 12.01568 307.1157
7.795880e-10 1.539717e-05

 

#B

cont.matrix <- makeContrasts(overallEnrich = (MALE + FEMALE)/2,FemalevMale =
FEMALE - MALE,levels=design)

fit4 <- contrasts.fit(fit2,cont.matrix)

fit5 <- eBayes(fit4)

FinalB<-topTable(fit5, number=389307,sort.by="M")

FinalB[1:5,]

                        ProbeID overallEnrich FemalevMale        A        F
P.Value    adj.P.Val

103220 CHR5:100995274-100995773      5.485304  -1.0123756 12.64761 392.7107
2.267845e-10 1.539717e-05

103222 CHR5:100995274-100995773      4.467658  -0.9554016 13.60286 364.9386
3.279371e-10 1.539717e-05

103211 CHR5:100993710-100995269      5.187434  -1.0645527 11.67345 313.0060
7.087562e-10 1.539717e-05

103219 CHR5:100995274-100995773      5.363803  -1.0720642 12.20548 308.5935
7.610504e-10 1.539717e-05

103209 CHR5:100993710-100995269      5.344564  -1.1764765 12.01568 307.1157
7.795880e-10 1.539717e-05

 

I think that perhaps the two sets of adj.p.values might represent different
statistical tests so they are both usable but for answering different
questions.

Can this be the case or what might be causing this “worsening” of my
results?

 

I apologize for my limited knowledge of statistics and naļve questions and
thank you very much for all the help.

 

Niki

 

 

-------------------------------------------------------

Hi Niki,

I guess it's not exactly clear in the help page for lmFit, but correlation
is only used if ndups > 1 or block is not NULL. Since you shouldn't be using
either of these, correlation won't be used and hence the default value
doesn't matter. This is all explained better in the vignette, under
Technical Replication.

Cheers,
Jenny

At 06:27 AM 3/27/2007, r.athanasiadou wrote:



Hi Jenny,
Thanks for the advice, I tried it and it looks it might be working.
 
I only have one question though:
If I dont calculate the correlation coefficient of my data then the
function lmFit() has a default of cor=0.75.
Isnt there any way I can calculate the average cor. value that is specific
to my arrays?
 
Thanks
Niki Athanasiadou
 
 
 
Hi Niki,

If I'm understanding your experimental design correctly, your 6 arrays are
something like this:

Cy3 vs. Cy5
Male1control vs. Male1chIP
Male2control vs. Male2chIP
...
Female3control vs. Female3chIP

What I don't understand is what you mean by wanting to find sex-related
differences. Do you mean MaleControl vs. FemaleControl? If so, then you will
need to do a separate channel analysis. Instead, if you mean you want to
find genes where the enrichment differs between males and females, as well
as finding those genes that have the same enrichment in males and females,
then I think your analysis below is a lot more complicated than it needs to
be.

First of all, you don't need to do duplicateCorrelation() if you have
biological replicates only. duplicateCorrelation is for technical
replicates. Second of all, I don't understand why you think you need to
separate the data to do male and female comparisons. You should be able to
do it all in one model. I'm not an expert, but off-hand this is the approach
I would take:

# the design matrix you have is good:
> design

    MALE FEMALE

A_M    1      0

B_M    1      0

F_M    1      0

C_F    0      1

D_F    0      1

E_F    0      1

# MALE represents the chIP enrichment in males and FEMALE represents the
chIP enrichment in females.

fit <- lmFit(MA, design)

# You could do eBayes on the fit object to find out which genes are enriched
separately for males and females.
# However, you probably should do contrasts instead:

cont.matrix <- makeContrasts(overallEnrich = (MALE + FEMALE)/2;
                                   MalevFemale = MALE - FEMALE;
                                   levels=design)

# I'm pretty sure that overallEnrich tests for enrichment in chIP samples,
regardless of sex, while  
#MalevFemale tests for the difference in enrichment between males and
females. To finish:

fit2 <- contrasts.fit(fit,cont.matrix)
fit2 <- eBayes(fit2)


Doing a statistical test for differences in enrichment between males and
females is better than just comparing gene lists between the male enrichment
and female enrichment in the fit object. For example, a gene could be
enriched in both sexes, but enriched more in one sex than the other. 

HTH,
Jenny


At 08:30 AM 3/26/2007, r.athanasiadou wrote:

Hi everyone,

I am using limma to analyze my array data. 

I have 6 (Nimblegen) arrays in total. They are two-color chip-on-chip arrays
without dye swaps.

3 arrays are male biological replicates and 3 female biological replicates.
I try to compare them and find sex-related differences (I have a pretty good
idea of what I am expecting to see and about 2% of the spots should be
different) and I also want to analyze them all together and see what is
enriched in both male-female after immunoprecipitations.

 

For the time being I chose to treat them as two-channel arrays (although I
am considering the option of doing single channel analysis).

 

I have two questions regarding the analysis (I apologize for my poor
knowledge of statistics!):

 

1)       Is there a way to use all six arrays for the normalization (I am
doing both within- and between-array normalization) and for the linear fit
model and then separate the male/female again to do eBayes? The reason I
would like to do this is that I observed (as expected) that by using a group
of 6 arrays gives much better statistics at the end than normalizing
separately for the 3 male and 3 female (the differences between them are not
great), but I want to keep the male/female data separate so that I use the
same processed dataset when I am doing the male/female comparison and when I
am searching for consistent ChIP enrichment across the sexes.

 

2)       My second question is regarding the calculation of the consensus
correlation to enter in the function  lmFit(). As I already mentioned I have
Nimblegen arrays, thus each gene is printed once but I do have triplicate
biological samples. So I thought that the command:

 

a<-duplicateCorrelation(MA, design=design, block=c(1,1,1,2,2,2)) 

should work, but what I got was:

 

> a$cor

[1] 1

> range(a$atanh.correlations)

[1] Inf Inf 

 

AND

(1)
MA<-normalizeWithinArrays(RG,RG$printer,method="control",controlspots=CONTRO
Lpos2)

(2)> show(RG)

An object of class "RGList"

$R

         A_M     B_M     F_M     C_F     D_F     E_F

[1,] 1603.22  697.67 1079.44  998.67  657.56 1849.00

[2,] 1874.22 1339.56 1331.89 1158.56 1146.00 1749.22

[3,] 1306.67 2033.00 1733.22  691.89 1814.56 1845.33

[4,]  827.00 1417.33 2237.22  694.56 1425.67 1650.56

[5,] 1052.22  565.44 1159.56  933.44  707.56 1091.11

389302 more rows ...

 

$G

         A_M     B_M     F_M     C_F     D_F     E_F

[1,]  521.89  632.33 1402.44  466.56 2093.33 1918.00

[2,]  863.11 1497.22 1733.33 1539.22 2331.56 2510.44

[3,]  730.78 5007.89 1898.56  706.78 4404.00 4502.78

[4,] 1437.44 2538.89 2860.67 1986.11 3349.33 2622.33

[5,]  720.33  779.44 1283.78  830.78 1788.00 1598.00

389302 more rows ...

 

$targets

  ChIP.ID FileName   Cy3 Cy5

1   95159  A_M.txt input ChIP

2  101711  B_M.txt input ChIP

3  101922  F_M.txt input ChIP

4   95854  C_F.txt input ChIP

5  104599  D_F.txt input ChIP

6  101923  E_F.txt input ChIP

 

$source

[1] "generic"

 

$printer

  ngrid.r ngrid.c nspot.r nspot.c

1  389307       1  389307       1

 

$genes

[1] "CHR1:3527111-3527575" "CHR1:3527111-3527575" "CHR1:3527111-3527575"
"CHR1:3527111-3527575" 

 

"CHR1:3527111-3527575"

 

(3)  > design

    MALE FEMALE

A_M    1      0

B_M    1      0

F_M    1      0

C_F    0      1

D_F    0      1

E_F    0      1

 

 


Can what I am getting (ie cor=1, range=inf) be correct? And if not any ideas
of what the problem might be and how I should deal with it?

 

> sessionInfo()

R version 2.4.1 (2006-12-18) 

i386-pc-mingw32 

 

locale:

LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United
Kingdom.1252;LC_MONETARY=English_United
Kingdom.1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252

 

attached base packages:

[1] "stats"     "graphics"  "grDevices" "utils"     "datasets"  "methods"  

[7] "base"     

 

other attached packages:

  limma 

"2.9.1"

 

 

 

Thanks

Niki Athanasiadou

 

ICMB

University of Edinburgh

 

 

 

 


         [[alternative HTML version deleted]]

_______________________________________________
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 

Jenny Drnevich, Ph.D.

Functional Genomics Bioinformatics Specialist
W.M. Keck Center for Comparative and Functional Genomics
Roy J. Carver Biotechnology Center
University of Illinois, Urbana-Champaign

330 ERML
1201 W. Gregory Dr.
Urbana, IL 61801
USA

ph: 217-244-7355
fax: 217-265-5066 
e-mail: drnevich@uiuc.edu

Jenny Drnevich, Ph.D.

Functional Genomics Bioinformatics Specialist
W.M. Keck Center for Comparative and Functional Genomics
Roy J. Carver Biotechnology Center
University of Illinois, Urbana-Champaign

330 ERML
1201 W. Gregory Dr.
Urbana, IL 61801
USA

ph: 217-244-7355
fax: 217-265-5066 
e-mail: drnevich@uiuc.edu


	[[alternative HTML version deleted]]

