


> Date: Thu, 20 Oct 2011 11:32:30 -0400
> From: jmacdon@med.umich.edu
> To: bahhab@hotmail.com
> CC: bioconductor@r-project.org
> Subject: Re: [BioC] removeBatchEffect problem
> 
> Hi John,
> 
> I don't know what to tell you.
> 
>  > x
>           [,1]     [,2]     [,3]     [,4]     [,5]     [,6]
> [1,] 21.34228 18.43027 21.29862 15.48104 23.16786 20.83750
> [2,] 27.00744 23.49814 27.41335 23.03560 24.70404 22.56948
> [3,] 27.02660 24.11282 26.70383 23.81685 22.84880 23.33938
> [4,] 21.71527 21.39641 21.07932 20.88411 20.31454 20.50233
> [5,] 32.57084 30.33254 31.80184 30.28329 31.11213 30.74922
> [6,] 22.51718 22.04644 22.17471 21.84274 20.94986 21.76465
> 
>  > design <- model.matrix(~factor(c(1,2,1,2,3,4)))
>  > removeBatchEffect(x, rep(c(1,2,1), each=2), design=design)
>           [,1]     [,2]     [,3]     [,4]     [,5]     [,6]
> [1,] 20.59406 17.68205 22.04684 16.22926 23.16786 20.83750
> [2,] 26.99328 23.48398 27.42751 23.04976 24.70404 22.56948
> [3,] 26.87191 23.95813 26.85852 23.97153 22.84880 23.33938
> [4,] 21.42821 21.10935 21.36638 21.17117 20.31454 20.50233
> [5,] 32.36628 30.12798 32.00640 30.48785 31.11213 30.74922
> [6,] 22.38064 21.90990 22.31125 21.97928 20.94986 21.76465
> 
> Works for me with the small subset of data you supply.

Thanks Jim,

My issue was that if I then ran removeBatchEffect with a different design element 
(just to prove to myself I was doing this right) I got the same results; i.e. the design 
element appeared to be having no effect.  But that didn't happen with just six cols 
of data (following on from what you have above)...

> design <- model.matrix(~factor(c(1,1,1,2,2,2)))
> removeBatchEffect(x, rep(c(1,2,1), each=2), design=design)
           V1       V2       V3       V4       V5       V6
[1,] 20.49073 17.57872 23.00172 17.18414 22.31631 19.98595
[2,] 27.26734 23.75804 26.89355 22.51580 24.96394 22.82938
[3,] 27.33608 24.42230 26.08487 23.19789 23.15828 23.64886
[4,] 21.71513 21.39627 21.07960 20.88439 20.31440 20.50219
[5,] 32.52130 30.28300 31.90092 30.38237 31.06259 30.69968
[6,] 22.58024 22.10950 22.04858 21.71661 21.01292 21.82771

so I tried again with 16 columns and there are differences between
    design <- model.matrix(~factor(c(1,2,1,2,3,4,3,4,5,6,5,6,7,8,7,8)))
and  

    design<- model.matrix(~factor(c(1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2) )) 
it's just they are so small they get lost in the rounding.

So I guess it is working as it should but my method of convincing myself that it works
was seriously flawed.

Thanks for your time.

> x<-head(e.qn) #the first 6 rows of my data; see below
> design<- model.matrix(~factor(c(1,2,1,2,3,4,3,4,5,6,5,6,7,8,7,8) ))
> y<-removeBatchEffect(x, rep(c(1,1,2,2), 4), design=design)
> design<- model.matrix(~factor(c(1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2) ))
> z<-removeBatchEffect(x, rep(c(1,1,2,2), 4), design=design)
> max(y-z)
[1] 7.105427e-15
> format(z[1,1], digits=20)
[1] "21.342279433232679509"
> format(y[1,1], digits=20)
[1] "21.342279433232683061"
> x
         [,1]     [,2]     [,3]     [,4]     [,5]     [,6]
[1,] 21.38993 18.47792 21.25097 15.43339 23.21551 20.88515
[2,] 27.19205 23.68276 27.22873 22.85098 24.88866 22.75410
[3,] 26.83332 23.91955 26.89710 24.01012 22.65553 23.14611
[4,] 21.46395 21.14509 21.33064 21.13543 20.06322 20.25101
[5,] 32.35783 30.11952 32.01486 30.49631 30.89911 30.53621
[6,] 21.98646 21.51572 22.70543 22.37346 20.41914 21.23393
         [,7]     [,8]     [,9]    [,10]    [,11]    [,12]
[1,] 22.96322 19.39987 23.04609 20.16191 23.91040 20.03716
[2,] 24.79286 22.31044 27.18536 24.02368 27.96298 23.38279
[3,] 24.08869 24.10216 24.31715 23.68728 24.82202 24.48001
[4,] 20.86551 21.10980 20.98385 20.38379 22.12131 21.91251
[5,] 30.27100 31.07816 30.83535 30.68225 32.35783 30.89911
[6,] 21.75568 22.39880 21.24631 20.87767 22.90156 22.74874
        [,13]    [,14]    [,15]    [,16]
[1,] 20.96133 19.92078 24.62166 19.67955
[2,] 27.51171 23.99884 26.76161 22.99287
[3,] 25.50943 24.19461 24.78192 24.17332
[4,] 21.55765 21.09746 21.33941 21.15254
[5,] 31.54015 29.93768 32.12096 31.07816
[6,] 21.13245 21.39392 21.68297 21.73048
> y
         [,1]     [,2]     [,3]     [,4]     [,5]     [,6]
[1,] 21.34228 18.43027 21.29862 15.48104 23.16786 20.83750
[2,] 27.00744 23.49814 27.41335 23.03560 24.70404 22.56948
[3,] 27.02660 24.11282 26.70383 23.81685 22.84880 23.33938
[4,] 21.71527 21.39641 21.07932 20.88411 20.31454 20.50233
[5,] 32.57084 30.33254 31.80184 30.28329 31.11213 30.74922
[6,] 22.51718 22.04644 22.17471 21.84274 20.94986 21.76465
         [,7]     [,8]     [,9]    [,10]    [,11]    [,12]
[1,] 23.01087 19.44752 22.99844 20.11426 23.95805 20.08481
[2,] 24.97748 22.49506 27.00074 23.83906 28.14760 23.56740
[3,] 23.89542 23.90888 24.51042 23.88055 24.62875 24.28674
[4,] 20.61419 20.85848 21.23517 20.63511 21.86999 21.66119
[5,] 30.05798 30.86514 31.04836 30.89527 32.14481 30.68609
[6,] 21.22496 21.86808 21.77704 21.40839 22.37084 22.21802
        [,13]    [,14]    [,15]    [,16]
[1,] 20.91368 19.87313 24.66931 19.72720
[2,] 27.32709 23.81422 26.94623 23.17748
[3,] 25.70270 24.38789 24.58865 23.98005
[4,] 21.80897 21.34879 21.08809 20.90122
[5,] 31.75316 30.15070 31.90794 30.86514
[6,] 21.66317 21.92464 21.15225 21.19976
> z
         [,1]     [,2]     [,3]     [,4]     [,5]     [,6]
[1,] 21.34228 18.43027 21.29862 15.48104 23.16786 20.83750
[2,] 27.00744 23.49814 27.41335 23.03560 24.70404 22.56948
[3,] 27.02660 24.11282 26.70383 23.81685 22.84880 23.33938
[4,] 21.71527 21.39641 21.07932 20.88411 20.31454 20.50233
[5,] 32.57084 30.33254 31.80184 30.28329 31.11213 30.74922
[6,] 22.51718 22.04644 22.17471 21.84274 20.94986 21.76465
         [,7]     [,8]     [,9]    [,10]    [,11]    [,12]
[1,] 23.01087 19.44752 22.99844 20.11426 23.95805 20.08481
[2,] 24.97748 22.49506 27.00074 23.83906 28.14760 23.56740
[3,] 23.89542 23.90888 24.51042 23.88055 24.62875 24.28674
[4,] 20.61419 20.85848 21.23517 20.63511 21.86999 21.66119
[5,] 30.05798 30.86514 31.04836 30.89527 32.14481 30.68609
[6,] 21.22496 21.86808 21.77704 21.40839 22.37084 22.21802
        [,13]    [,14]    [,15]    [,16]
[1,] 20.91368 19.87313 24.66931 19.72720
[2,] 27.32709 23.81422 26.94623 23.17748
[3,] 25.70270 24.38789 24.58865 23.98005
[4,] 21.80897 21.34879 21.08809 20.90122
[5,] 31.75316 30.15070 31.90794 30.86514
[6,] 21.66317 21.92464 21.15225 21.19976
> y-z
              [,1]         [,2]          [,3]          [,4]
[1,]  3.552714e-15 0.000000e+00 -3.552714e-15 -1.776357e-15
[2,]  3.552714e-15 0.000000e+00  0.000000e+00  0.000000e+00
[3,]  0.000000e+00 0.000000e+00  0.000000e+00  0.000000e+00
[4,]  0.000000e+00 0.000000e+00  0.000000e+00  0.000000e+00
[5,] -7.105427e-15 0.000000e+00  0.000000e+00 -7.105427e-15
[6,]  7.105427e-15 3.552714e-15  0.000000e+00  0.000000e+00
             [,5]         [,6]          [,7]          [,8]
[1,] 0.000000e+00 0.000000e+00  0.000000e+00  0.000000e+00
[2,] 0.000000e+00 0.000000e+00  0.000000e+00  0.000000e+00
[3,] 0.000000e+00 0.000000e+00  0.000000e+00  0.000000e+00
[4,] 0.000000e+00 0.000000e+00  0.000000e+00  0.000000e+00
[5,] 0.000000e+00 0.000000e+00 -7.105427e-15 -7.105427e-15
[6,] 3.552714e-15 3.552714e-15  0.000000e+00  0.000000e+00
             [,9]        [,10]        [,11]         [,12]
[1,] 3.552714e-15 3.552714e-15 0.000000e+00 -3.552714e-15
[2,] 0.000000e+00 0.000000e+00 0.000000e+00  0.000000e+00
[3,] 0.000000e+00 0.000000e+00 0.000000e+00 -3.552714e-15
[4,] 0.000000e+00 0.000000e+00 0.000000e+00  0.000000e+00
[5,] 0.000000e+00 0.000000e+00 7.105427e-15  0.000000e+00
[6,] 0.000000e+00 3.552714e-15 0.000000e+00  0.000000e+00
             [,13] [,14]         [,15]         [,16]
[1,]  0.000000e+00     0 -3.552714e-15 -3.552714e-15
[2,] -3.552714e-15     0  0.000000e+00  0.000000e+00
[3,] -3.552714e-15     0  0.000000e+00  0.000000e+00
[4,]  0.000000e+00     0  0.000000e+00  0.000000e+00
[5,]  0.000000e+00     0  0.000000e+00 -7.105427e-15
[6,]  0.000000e+00     0 -3.552714e-15 -3.552714e-15
> 


> 
> Best,
> 
> Jim
> 
> 
> 
> On 10/20/2011 11:02 AM, John Coulthard wrote:
> > Thanks Jim, but still regardless of the design matrix I give it the output is the same?!  I must be missing something else.
> >
> >> design<- model.matrix(~factor(c(1,2,1,2,3,4,3,4,5,6,5,6,7,8,7,8) ))
> >> head(removeBatchEffect(e.qn,c(1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2),design=design)[,1:6])
> >           [,1]     [,2]     [,3]     [,4]     [,5]     [,6]
> > [1,] 21.34228 18.43027 21.29862 15.48104 23.16786 20.83750
> > [2,] 27.00744 23.49814 27.41335 23.03560 24.70404 22.56948
> > [3,] 27.02660 24.11282 26.70383 23.81685 22.84880 23.33938
> > [4,] 21.71527 21.39641 21.07932 20.88411 20.31454 20.50233
> > [5,] 32.57084 30.33254 31.80184 30.28329 31.11213 30.74922
> > [6,] 22.51718 22.04644 22.17471 21.84274 20.94986 21.76465
> >> head(removeBatchEffect(e.qn,c(1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2))[,1:6])
> >           [,1]     [,2]     [,3]     [,4]     [,5]     [,6]
> > [1,] 21.34228 18.43027 21.29862 15.48104 23.16786 20.83750
> > [2,] 27.00744 23.49814 27.41335 23.03560 24.70404 22.56948
> > [3,] 27.02660 24.11282 26.70383 23.81685 22.84880 23.33938
> > [4,] 21.71527 21.39641 21.07932 20.88411 20.31454 20.50233
> > [5,] 32.57084 30.33254 31.80184 30.28329 31.11213 30.74922
> > [6,] 22.51718 22.04644 22.17471 21.84274 20.94986 21.76465
> >> design<- model.matrix(~factor(c(1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2) ))
> >> head(removeBatchEffect(e.qn,c(1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2),design=design)[,1:6])
> >           [,1]     [,2]     [,3]     [,4]     [,5]     [,6]
> > [1,] 21.34228 18.43027 21.29862 15.48104 23.16786 20.83750
> > [2,] 27.00744 23.49814 27.41335 23.03560 24.70404 22.56948
> > [3,] 27.02660 24.11282 26.70383 23.81685 22.84880 23.33938
> > [4,] 21.71527 21.39641 21.07932 20.88411 20.31454 20.50233
> > [5,] 32.57084 30.33254 31.80184 30.28329 31.11213 30.74922
> > [6,] 22.51718 22.04644 22.17471 21.84274 20.94986 21.76465
> >
> >
> >
> >> Date: Thu, 20 Oct 2011 09:33:01 -0400
> >> From: jmacdon@med.umich.edu
> >> To: bahhab@hotmail.com
> >> CC: bioconductor@r-project.org
> >> Subject: Re: [BioC] removeBatchEffect problem
> >>
> >> Hi John,
> >>
> >> On 10/20/2011 6:42 AM, John Coulthard wrote:
> >>> Hi
> >>>
> >>> I'm having troubles with the removeBatchEffect function.  It seems that whatever I put in as the 'design' I get the same result out of removeBatchEffect(); Perhaps I'm not using the right syntax for the 'design'?
> >>>
> >>> Below is each line of the function run individually and the output from the line using 'design' with various inputs; Always the same output?!
> >>>
> >>> My design is 4x2x2; 4 tissues(T) x 2 subgroups(S) x 2 days(D) and the col order in e.qn is...
> >>> T1S1D1, T1S2D1, T1S1D2, T1S2D2, T2S1D1, T2S2D1, T2S1D2, T2S2D2, T3S1D1, T3S2D1, T3S1D2, T3S2D2, T4S1D1, T4S2D1, T4S1D2, T4S2D2,
> >>>
> >>> And it's the day(D) effect I'm trying to remove.
> >>>
> >>>
> >>> Can anyone point me to the right way to run this.
> >>>
> >>> Many thanks
> >>> John
> >>>
> >>> the command I started with...
> >>> e.qn_batch_removed<-removeBatchEffect(e.qn,c(1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2),design=matrix(c(1,2,1,2,3,4,3,4,5,6,5,6,7,8,7,8),16,1))
> >> That isn't the correct design matrix. What you most likely want is
> >>
> >> design<- model.matrix(~factor(c(1,2,1,2,3,4,3,4,5,6,5,6,7,8,7,8) ))
> >>
> >> Best,
> >>
> >> Jim
> >>
> >>
> >>> then line by line...
> >>>> removeBatchEffect
> >>> function (x, batch, design = matrix(1, ncol(x), 1))
> >>> {
> >>>       x<- as.matrix(x)
> >>>       batch<- as.factor(batch)
> >>>       contrasts(batch)<- contr.sum(levels(batch))
> >>>       X<- model.matrix(~batch)[, -1, drop = FALSE]
> >>>       X<- qr.resid(qr(design), X)
> >>>       qrX<- qr(X)
> >>>       t(qr.resid(qrX, t(x)))
> >>> }
> >>> <environment: namespace:limma>
> >>>
> >>>
> >>>> x<- as.matrix(e.qn)
> >>>> batch<- as.factor(c(1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2))
> >>>> contrasts(batch)<- contr.sum(levels(batch))
> >>>> X<- model.matrix(~batch)[, -1, drop = FALSE]
> >>>> qr.resid(qr(matrix(c(1,2,1,2,3,4,3,4,5,6,5,6,7,8,7,8),16,1)), X)
> >>>      batch1
> >>> 1       1
> >>> 2       1
> >>> 3      -1
> >>> 4      -1
> >>> 5       1
> >>> 6       1
> >>> 7      -1
> >>> 8      -1
> >>> 9       1
> >>> 10      1
> >>> 11     -1
> >>> 12     -1
> >>> 13      1
> >>> 14      1
> >>> 15     -1
> >>> 16     -1
> >>>> qr.resid(qr(matrix(c(1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2),16,1)), X)
> >>>      batch1
> >>> 1       1
> >>> 2       1
> >>> 3      -1
> >>> 4      -1
> >>> 5       1
> >>> 6       1
> >>> 7      -1
> >>> 8      -1
> >>> 9       1
> >>> 10      1
> >>> 11     -1
> >>> 12     -1
> >>> 13      1
> >>> 14      1
> >>> 15     -1
> >>> 16     -1
> >>>> qr.resid(qr(matrix(1,16,1)), X)
> >>>      batch1
> >>> 1       1
> >>> 2       1
> >>> 3      -1
> >>> 4      -1
> >>> 5       1
> >>> 6       1
> >>> 7      -1
> >>> 8      -1
> >>> 9       1
> >>> 10      1
> >>> 11     -1
> >>> 12     -1
> >>> 13      1
> >>> 14      1
> >>> 15     -1
> >>> 16     -1
> >>>
> >>>    		 	   		
> >>> 	[[alternative HTML version deleted]]
> >>>
> >>> _______________________________________________
> >>> Bioconductor mailing list
> >>> Bioconductor@r-project.org
> >>> https://stat.ethz.ch/mailman/listinfo/bioconductor
> >>> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor
> >> -- 
> >> James W. MacDonald, M.S.
> >> Biostatistician
> >> Douglas Lab
> >> University of Michigan
> >> Department of Human Genetics
> >> 5912 Buhl
> >> 1241 E. Catherine St.
> >> Ann Arbor MI 48109-5618
> >> 734-615-7826
> >>
> >> **********************************************************
> >> Electronic Mail is not secure, may not be read every day, and should not be used for urgent or sensitive issues
> >>
> >   		 	   		
> > 	[[alternative HTML version deleted]]
> >
> > _______________________________________________
> > Bioconductor mailing list
> > Bioconductor@r-project.org
> > https://stat.ethz.ch/mailman/listinfo/bioconductor
> > Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor
> 
> -- 
> James W. MacDonald, M.S.
> Biostatistician
> Douglas Lab
> University of Michigan
> Department of Human Genetics
> 5912 Buhl
> 1241 E. Catherine St.
> Ann Arbor MI 48109-5618
> 734-615-7826
> 
> **********************************************************
> Electronic Mail is not secure, may not be read every day, and should not be used for urgent or sensitive issues 
> 
 		 	   		  
	[[alternative HTML version deleted]]

