Question: Adding Batch Factor in linear model...
0
gravatar for Md.Mamunur Rashid
9.7 years ago by
Md.Mamunur Rashid260 wrote:
-------- Original Message -------- Subject: Re: [BioC] Fitting arrayweights in linear model Date: Tue, 3 Nov 2009 14:23:52 +0000 From: Md.Mamunur Rashid <mamunur.rashid@kcl.ac.uk> To: Wei Shi <shi@wehi.edu.au>, bioconductor <bioconductor@stat.math.ethz.ch> 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-49 11-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.inf ormatics.conductor >> [[alternative HTML version deleted]]
normalization go • 582 views
ADD COMMENTlink written 9.7 years ago by Md.Mamunur Rashid260
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: 142 users visited in the last hour