Dear Wei,

Sorry for replying late. I was on vacation and went back to home. thanks for your suggestion.
I have tried reading some previous mails from bioc forum to find a way to add batch factor
in the linear model. But did not clearly understand how to do that.

In my linear model I am comparing three conditions. So I have created a contrast matrix:

    fit1<- lmfit(normalized_matrix, design, array_weights)
    contrast<- make.contrast(group2-gropu1,group3-group2,group3-group1,levels= design)
    confit<- contrasts.fit(fit1 ,contrast)
    efit<- eBayes()
       

I have tried the same R script  with a new batch of 96 samples from the same 3 groups.
These 96 sample does not really have that much variation among them . Yet again the top
list of genes have a very high adjusted p.value.

Q1. Does this mean that there is not not enough evidence of differential expression.  ??
Q2. Is it a very good idea to analyze 96 samples in one go, or should I subset the file( ExpressionSet)
     to work with smaller set of data. Cause I have got this impression that for this large data
     set variations in one of the slides (illumina HT12 containing 12 arrays) do a lot of damage
     when the whole data set is normalized. ...?????



(I have done this part to check if my algorithm is ok or not  ..)
I have tested the R script for contrasting among male and female group. and found the
evidence of differential expression among the genes found in Y chromosome.
This finding gave me a idea that my R script is ok (anyway this is very straight forward)
****Please any suggestion about this assumption.

  	contrast<- make.contrast(male-female, levels= design)

Can you please suggest me few more gene filtering methods. I will start with the gene filter package.

Thanks in advance.

regards,
Mamun



On 09/30/2009 05:26 AM, Wei Shi wrote:
> Dear Mamun:
>
>     Sorry for replying to you late. Your boxplot of raw data shows 
> that there is a batch effect in your data. You might need to include 
> this batch effect in your linear model fitting.
>
>     Arrays with low weights will not be dropped when you include array 
> weights in your linear model fitting. The idea of giving weights to 
> arrays is to keep all the arrays in the analysis.
>
>     Hope this helps.
>
> Cheers,
> Wei
>
>
> Md.Mamunur Rashid wrote:
>> Dear Wei,
>>
>> Thanks for your reply. I know it's a long email, but I tried to give details about the raw data
>> Please have a look when you have some spare time.
>>
>>
>> 1. About the raw data :
>> -----------------------
>> Yes this arrays are from the same bead chip. The experiments were done by three individuals (each
>> person conducted  experiments with 3 beadchips each consisting 12 arrays, 96 intotal). I think this
>> is the reason for high variance among 3 sets of arrays. 1-32, 33-64, 65-96. But after normalization
>> ("quantile"-lumi package) the data looks quite normal.
>>
>>
>> I have attached two snapshot of boxplot of the raw data and normalized data.
>>
>> 	A. Normalized data -http://www.esnips.com/doc/ad4c8a55-2166-4911-8510-df81eaea63df/data_96_illumina_norm
>>
>> 	B. Raw data -http://www.esnips.com/doc/0a965715-b5ec-44da-b64c-4d025c63827b/96_illumina_raw
>>
>> If you look at the two snapshots above : Expression range for raw data in most cases are below 6. But
>> after normalization it came slightly above 6. Is it possible at all or something is going wrong during
>> the normalization.
>>
>>
>> One question here. Does illumina data require intra-array normalization.?? To the extent of my little
>> study I know that they do not need it as it is already done. Please let me know what do you think.
>>
>>
>> 2. About Arrayweight
>> --------------------
>>
>> As you can see below the the results have not been improved much compare to the old topTable result.
>>
>> 1. Do you think this high variance in arrays are causing this problem.???
>> 2. Does fitting the arrayweight in the linear model along the normalized object drops the downweighted arrays
>>      or that needed to be done manually??
>>
>>
>> 3. Pre-processing with beadarray
>> --------------------------------
>>
>>       I have the done the same analysis using beadArray package(data import and pre processing) and had the similar
>>       result. So there can be two possibilities which can lead to this problem.
>>
>>       A. Problems in the raw data.
>>       B. Or I am doing a wrong linear modeling of the data.(as that is the common part in the both analysis)
>>
>>       I have attached my code fragment here. please have a look:
>>
>>
>> 	*** What are the other bio packages/tests I can use for identifying differentially expressed genes>?
>>
>> Thanks in advance
>>
>> best regards,
>> Md.Mamunur Rashid
>>
>>
>> Code Fragment
>> -----------------------------------
>>
>> About fitting the weights :  I have done the following
>>
>> data_96_nuId_E holds the normalized object
>>
>>
>>
>> data_96_sampleType<- c("I","I","I","I","I","S","S","S","I","C","S","S","C","I","C","I","I","I","I","S","I","S","I","S","C","C","S","I","C","S","I","S","S","C","C","S","S","I","I","I","S","I","S","I","I","C","I","C","S","C","I","I","S","C","S","S","S","C","C","I","S","S","C","I","I","C","C","S","S","S","C","C","S","I","C","C","C","S","C","C","I","I","C","S","C","C","C","S","C","S","C","I","I","S","I","C")
>> data_96_design<- model.matrix(~0+data_96_sampleType)
>> colnames(data_96_design)<- c('C','I','S')
>>
>> ## filtring out the non-expressed genes
>> presentCount<- detectionCall(data_96_raw_nuID)
>> data_96_Matrix<- data_96_Matrix_old[presentCount>   0, ]
>> data_96_probeList<- rownames(data_96_Matrix)
>>
>>
>> arraw_E_96<- arrayWeights(data_96_nuID_E,data_96_design)
>>
>> data_96_contrast<- makeContrasts (C-I,I-S,S-C,levels=data_96_design)
>>
>> fitw_nu_96<- lmFit(data_96_Matrix,data_96_design,weights=arraw_E_96)
>> fitw_nu_con_96<- contrasts.fit(fitw_nu_96,data_96_contrast)
>> fitw_e_nu_96<- eBayes(fitw_nu_con_96)
>> topTable(fitw_e_nu_96,coef=1,adjust="BH")
>>
>>
>>                         ID       logFC  AveExpr         t      P.Value adj.P.Val
>> 3564  oEptU7gl5ULB7ghAf4 -0.11681240 6.341389 -4.749578 7.204301e-06 0.1484950
>> 3900  cUSJ2.kBKqPAk0eZPY  0.31230030 6.881320  4.561405 1.514923e-05 0.1561280
>> 10705 uHk1A0FNxKMV0k7a50 -0.13640964 6.616042 -4.197525 6.070562e-05 0.2895144
>> 14598 WokXUp_QIC4pd_s1.U -0.09998362 6.342947 -4.111694 8.338387e-05 0.2895144
>> 4324  xXLgIpJ_AeI7gRNJ_U  0.11230531 6.478843  4.072272 9.634276e-05 0.2895144
>> 6399  HHteoiUdSJz9wh.e64  0.16605833 7.322017  4.034305 1.106350e-04 0.2895144
>> 13446 x65aXRiHVTIIUQiv00 -0.23601528 6.487643 -3.986742 1.314198e-04 0.2895144
>> 13445 HrlpdGIdVMghRCK.TU -0.20278431 6.470915 -3.986291 1.316337e-04 0.2895144
>> 16912 BgFUhUUgIQu3JDl_sU  0.11902363 6.815447  3.978352 1.354541e-04 0.2895144
>> 9813  NAv6twbtRAh7QuCDro  0.09999394 6.511654  3.958318 1.455733e-04 0.2895144
>>                 B
>> 3564  3.5478599
>> 3900  2.8790555
>> 10705 1.6341835
>> 14598 1.3504633
>> 4324  1.2214803
>> 6399  1.0980586
>> 13446 0.9445695
>> 13445 0.9431198
>> 16912 0.9176269
>> 9813  0.8534456
>>
>>
>>
>>
>>
>> About fitting the arrayweights in the linear model :
>> previously I have produced arrayweights
>>
>> On 09/21/2009 12:59 AM, Wei Shi wrote:
>>    
>>> Dear Mamun:
>>>
>>>      Sorry for replying to you late. I am just back from the holiday.
>>>
>>>      Some of your arrays have pretty low weights. Arrays 37-42 have
>>> weights of around 0.2. Are these arrays from the same bead chip? You
>>> can try put array weights in the linear model fitting  to see if you
>>> can get DE genes.
>>>
>>> Cheers,
>>> Wei
>>>
>>> Md.Mamunur Rashid wrote:
>>>      
>>>> Dear Wei,
>>>>
>>>> Thanks again for your reply. I have measured the array-weights using limma's arrayweights() method. There seem to be
>>>> a good amount of variation among the raw data. But I am not really sure whether it is acceptable or not??.
>>>>
>>>> But the normalized data have an average value of 1.022609
>>>> with min = 0.5641711 max =1.417277.
>>>>
>>>> Also 2 more question :
>>>>
>>>> 1. can you explain a bit more about how can I filter out the genes whose detection p-values
>>>> are larger than cut-off (e.g. 0.01) in all samples??
>>>> 2. how can I check batch effect and chip effect in the raw data.??
>>>>
>>>> Thanks in advance. Any help from anyone also is appreciated.
>>>>
>>>> regards
>>>> Mamun
>>>>
>>>>
>>>> Below is the result from arrayWeigts() method.(Table 1 (weights for raw data) Table 2(weights of normalized data))
>>>> *** Also attached a text file with the weights.
>>>>
>>>> ArrayWeight of Raw Data:
>>>>
>>>> arrayWeights(data_96_raw,data_96_design)
>>>>
>>>>
>>>>           1         2         3         4         5         6         7         8
>>>> 3.2861102 0.6091478 2.6927503 2.2778201 2.4391841 0.7802437 1.1378995 0.5714981
>>>>           9        10        11        12        13        14        15        16
>>>> 0.8091165 0.5260511 0.7565424 0.8354660 0.2929977 0.4950789 0.4512856 0.1283279
>>>>          17        18        19        20        21        22        23        24
>>>> 0.6713553 0.6275333 0.5441217 0.2711778 0.7278274 0.3323277 0.5786381 0.7375983
>>>>          25        26        27        28        29        30        31        32
>>>> 0.5793944 0.4018671 0.6658730 1.6805341 0.7210032 1.2802383 0.8902326 0.6715769
>>>>          33        34        35        36        37        38        39        40
>>>> 0.2959048 0.8854763 0.7587111 1.5067705 0.1766820 0.1639351 0.1739882 0.1833592
>>>>          41        42        43        44        45        46        47        48
>>>> 0.2443871 0.2188281 0.2945594 0.2607380 0.3419603 0.5932775 0.3396003 0.7229979
>>>>          49        50        51        52        53        54        55        56
>>>> 0.8274153 1.7944800 0.8306203 1.0738534 1.5052515 2.9798751 2.0836544 1.5010941
>>>>          57        58        59        60        61        62        63        64
>>>> 2.3827591 2.7575819 1.5967670 3.3397084 2.0576009 2.9660262 3.8417611 1.3469622
>>>>          65        66        67        68        69        70        71        72
>>>> 1.7334743 2.7855044 4.5417142 3.2551286 3.9634064 3.8054377 3.8020479 4.2732205
>>>>          73        74        75        76        77        78        79        80
>>>> 0.2659556 0.3440704 0.6404601 0.8697107 1.0736648 0.8835502 0.9912724 0.4334270
>>>>          81        82        83        84        85        86        87        88
>>>> 0.8801461 1.0309162 2.5352693 1.8366521 2.5501944 4.5388081 4.2603196 1.7665255
>>>>          89        90        91        92        93        94        95        96
>>>> 2.5511084 2.3867529 4.0237262 2.1647555 3.3101943 2.5863558 1.8618118 0.2397050
>>>>
>>>>
>>>>
>>>>
>>>> Array Weights of normalized data:
>>>>
>>>>           1         2         3         4         5         6         7         8
>>>> 0.8908085 1.1978748 1.1454450 1.0215452 1.3605699 0.8141813 1.4172765 1.2958270
>>>>           9        10        11        12        13        14        15        16
>>>> 1.0522608 1.2562065 1.3531495 1.1489293 1.0596229 1.1217643 0.9156484 0.6705476
>>>>          17        18        19        20        21        22        23        24
>>>> 0.9643412 1.2120423 1.0672521 0.9983735 0.8096782 0.9379575 1.0493924 0.7614746
>>>>          25        26        27        28        29        30        31        32
>>>> 1.1573255 1.2735353 1.3795986 0.9499952 1.3602425 1.2549726 1.1772558 1.4158351
>>>>          33        34        35        36        37        38        39        40
>>>> 1.3659316 1.0658774 1.3185647 0.9017821 0.7915704 0.6326567 1.0325512 0.6812818
>>>>          41        42        43        44        45        46        47        48
>>>> 1.0142990 1.1921695 1.1476346 0.8255798 1.2012711 1.1893762 0.9367947 1.1594407
>>>>          49        50        51        52        53        54        55        56
>>>> 1.1072266 1.0561352 0.8488650 1.1756689 1.0554286 0.9326934 1.2836555 0.9665945
>>>>          57        58        59        60        61        62        63        64
>>>> 1.3293081 1.3027377 1.3459009 1.2834152 0.9657114 1.0041629 0.8823282 0.6843356
>>>>          65        66        67        68        69        70        71        72
>>>> 0.8883121 0.9374396 0.9799943 0.9733364 1.2045700 1.0693506 0.7730626 0.9408430
>>>>          73        74        75        76        77        78        79        80
>>>> 0.6685933 1.0492289 0.9952222 0.9690452 1.0150076 1.0148188 0.6687192 0.5641711
>>>>          81        82        83        84        85        86        87        88
>>>> 0.7872238 0.8793660 0.8553373 0.9662603 0.6143880 1.0340769 0.9842437 0.6626941
>>>>          89        90        91        92        93        94        95        96
>>>> 1.0270788 0.8590296 1.0807271 0.8162435 1.0398548 0.8993595 1.2094240 0.5715857
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>        
>>
>>
>> 	[[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
>>    


	[[alternative HTML version deleted]]

