[BioC] Limma: What to choose?

James W. MacDonald jmacdon at med.umich.edu
Fri Sep 5 15:17:14 CEST 2008


Hi Ingrid,

I am a bit confused by what you have done here. I understand the first 
analysis (basic), and given your data I think that is probably not the 
way you want to go.

For the second analysis it appears that either you have re-ordered your 
data, or the design matrix is completely wrong. The design matrix has to 
mimic your data regardless of the analysis, and if the data are still 
NG, NG, NG, HG, HG, etc then this design matrix isn't going to make the 
correct comparisons. So to do this correctly you likely need to use the 
design matrix from the basic analysis.

For the third analysis I think you got confused. You set up the 
treatment and target factors, but then you don't ever use the targets 
factor. In addition, you don't have a contrasts matrix so you are simply 
testing to see if the average expression of a particular group is 
different from zero rather than comparing different groups. To do the 
paired analysis correctly you will need to do something like

designMA <- model.matrix(~ 0 + Treatment + Target)

You will then have a six-column design matrix looking something like this:

 > designMA
    Treatment1 Treatment2 Treatment3 Treatment4 Target2 Target3
1           1          0          0          0       0       0
2           1          0          0          0       0       1
3           1          0          0          0       1       0
4           0          1          0          0       1       0
5           0          1          0          0       0       0
6           0          1          0          0       0       1
7           0          0          1          0       0       0
8           0          0          1          0       0       1
9           0          0          1          0       1       0
10          0          0          0          1       1       0
11          0          0          0          1       0       0
12          0          0          0          1       0       1
attr(,"assign")
[1] 1 1 1 1 2 2
attr(,"contrasts")
attr(,"contrasts")$Treatment
[1] "contr.treatment"

attr(,"contrasts")$Target
[1] "contr.treatment"

the Target2 and Target3 coefficients capture the average differences 
between patient 29 and 30 and 29 and 34, respectively, which 
algebraically is the same as a paired analysis.

You then want a contrast matrix that takes these extra coefficients into 
account:

 > contrast <- matrix(c(1,0,-1,0,0,0,0,1,0,-1,0,0,1,-1,0,0,0,0),
ncol = 3, dimnames = list(c("NG","HG","NG_ben","HG_ben","P1","P2"), 
c("NG-NG_ben","HG-HG_ben","NG-HG")))
 > contrast
        NG-NG_ben HG-HG_ben NG-HG
NG             1         0     1
HG             0         1    -1
NG_ben        -1         0     0
HG_ben         0        -1     0
P1             0         0     0
P2             0         0     0

And the usual analysis from there on out.

Best,

Jim



Ingrid H. G. Østensen wrote:
> Hi
> 
> I have some data that I am going to analyze, but I have some problems regarding what to choose/do.
> 
> The dataset consists of 4 treatments groups: NG, NG+benfo, HG and HG+benfo and it has been used samples from 3 donors in each group: donor 29, 30 and 34. This means that it is 4 groups consisting of 3 samples in each. The platform that has been used is Illumina Human WG6, and the comparisons that are wanted is: NG - NG_benfo, HG - HG_benfo, NG – HG.
> 
> The data was read into R using lumiR, quality controlled, normalized, log2 transformed and the expression values was extracted resulting in a data matrix called dataSet_Norm_exp_log2.
> 
> First I run a quality control which showed me that the distribution is fine, but the samples were divided into donors not treatment in the PCA plot and hierarchical clustering plot. For me this indicates the differences between the donors are bigger than the differences between the treatments. And because of this I did not expect a “basic” analysis in limma to give me any good results, see ########### Basic ############
> 
> Since the donor effect was so visible I thought I should try to block the donor effect, but that actually gave me worse adjusted p-value. Is it wrong to do this, am I using it incorrectly or can the higher p-values being explained?
> Example in ############ Block #################
> 
> The three donors are in each group and the scientists also want a paired test: Donor 29 NG vs. donor 29 NG_benfo etc. I looked in the limma user guide and found an example of paired t-test, but I am not sure I am doing it correct and how to interpret the results. Example in ####### Pair ###########
> 
> In this situation would a paired test or a “group” test be the best? Is there a way to “see” what analysis is the best before one start the job? I do not like to try too many analyses because eventually the data will give the desired answer; I just have to torture it enough. And when can one say that there are no significant findings in the data, when to stop the torture? 
> 
> Regards,
> Ingrid
> 
> 
>> ############################## Basic ##################################
>>
>> sampleType <- c('NG','NG','NG','HG','HG','HG','NG_benfo','NG_benfo','NG_benfo','HG_benfo','HG_benfo','HG_benfo')     
>>    
>> designMa <- model.matrix(~-1+factor(sampleType))
>>
>> rownames(designMa) <- sampleType
>>       
>> colnames(designMa) <- c("HG_s","HG_benfo_s","NG_s","NG_benfo_s")     
>>
>> fitDesMa <- lmFit(dataSet_Norm_exp_log2, designMa)
>>
>> designMa
>          HG_s HG_benfo_s NG_s NG_benfo_s
> NG          0          0    1          0
> NG          0          0    1          0
> NG          0          0    1          0
> HG          1          0    0          0
> HG          1          0    0          0
> HG          1          0    0          0
> NG_benfo    0          0    0          1
> NG_benfo    0          0    0          1
> NG_benfo    0          0    0          1
> HG_benfo    0          1    0          0
> HG_benfo    0          1    0          0
> HG_benfo    0          1    0          0
> attr(,"assign")
> [1] 1 1 1 1
> attr(,"contrasts")
> attr(,"contrasts")$`factor(sampleType)`
> [1] "contr.treatment"
> 
>>                                                                     
>> contrast.matrix <- makeContrasts(NG_s - NG_benfo_s,HG_s - HG_benfo_s, NG_s-HG_s, levels = designMa)                               
>>
>> fitContr <- contrasts.fit(fitDesMa, contrast.matrix)
>>
>> fitContr2 <- eBayes(fitContr)
>>
>> top1 <- topTable(fitContr2, coef= 1, number = 20, adjust= "BH", sort.by="p")
>> top1
>          ID      logFC   AveExpr         t      P.Value adj.P.Val         B
> 2426   2426 -1.0843868 10.716173 -8.938389 5.482310e-06 0.2675587 1.8529767
> 17320 17320 -0.6227589  6.281775 -8.116261 1.261991e-05 0.2784645 1.5080171
> 36670 36670  1.3093614  7.984147  7.600320 2.202967e-05 0.2784645 1.2574042
> 2752   2752 -0.4334247  7.923580 -7.518626 2.412326e-05 0.2784645 1.2150191
> 25602 25602  0.6119058  9.056356  7.316350 3.029915e-05 0.2784645 1.1066829
> 41399 41399  0.5812922  8.270799  7.209620 3.423464e-05 0.2784645 1.0475136
> 14334 14334  0.5852656  5.977673  6.889422 4.977374e-05 0.3036840 0.8612701
> 2427   2427 -1.0252133 10.001603 -6.784891 5.639001e-05 0.3036840 0.7975167
> 7100   7100  0.7322154  7.625884  6.688233 6.336271e-05 0.3036840 0.7372208
> 47668 47668  0.6990338  8.337938  6.545533 7.542180e-05 0.3036840 0.6457823
> 2592   2592 -0.9705914  6.310054 -6.527260 7.713743e-05 0.3036840 0.6338612
> 31680 31680 -1.0568705  6.573984 -6.504740 7.931020e-05 0.3036840 0.6191020
> 1859   1859  1.0862890  6.535754  6.488754 8.089278e-05 0.3036840 0.6085795
> 7059   7059  0.6369413  5.951195  6.423283 8.774005e-05 0.3045383 0.5650911
> 42018 42018 -0.8482179 11.479166 -6.350866 9.605183e-05 0.3045383 0.5162425
> 41398 41398  0.6482108  8.106491  6.285794 1.042501e-04 0.3045383 0.4716697
> 35368 35368  1.6933744  9.167829  6.272022 1.060805e-04 0.3045383 0.4621527
> 4815   4815  0.6222491 11.210936  6.120113 1.287430e-04 0.3322279 0.3552235
> 3965   3965 -0.8698427  8.916196 -6.116510 1.293404e-04 0.3322279 0.3526427
> 155     155 -0.8420104  8.091452 -5.983258 1.536784e-04 0.3750060 0.2557538
> 
>> top2 <- topTable(fitContr2, coef= 2, number = 20, adjust= "BH", sort.by="p")
>> top2
>          ID      logFC   AveExpr         t      P.Value adj.P.Val          B
> 25602 25602  0.6544021  9.056356  7.824463 1.723540e-05 0.5557462 -0.9515403
> 36670 36670  1.2095773  7.984147  7.021113 4.261105e-05 0.5557462 -1.1474102
> 2426   2426 -0.8373037 10.716173 -6.901732 4.905181e-05 0.5557462 -1.1804654
> 43939 43939  1.5657869  9.799353  6.858864 5.161622e-05 0.5557462 -1.1926108
> 25601 25601  0.6822865  7.984298  6.767874 5.755464e-05 0.5557462 -1.2188870
> 38431 38431  0.3340007  5.317240  6.626234 6.832385e-05 0.5557462 -1.2611723
> 35368 35368  1.7394872  9.167829  6.442817 8.563363e-05 0.5970377 -1.3185448
> 35394 35394  0.4021078  6.140240  6.222546 1.129480e-04 0.6492061 -1.3915750
> 43459 43459  0.8113102  9.200194  6.057577 1.395481e-04 0.6492061 -1.4494001
> 35519 35519 -0.3788567  7.472821 -5.964382 1.575082e-04 0.6492061 -1.4833109
> 3896   3896  0.3074361  5.432543  5.940904 1.624160e-04 0.6492061 -1.4919995
> 36692 36692  0.8720572  7.116580  5.899135 1.715600e-04 0.6492061 -1.5076044
> 42124 42124 -0.6213031  7.597989 -5.893084 1.729301e-04 0.6492061 -1.5098808
> 27084 27084  0.5015841  5.690324  5.752551 2.083097e-04 0.7149992 -1.5638914
> 5019   5019  0.4185694  5.712520  5.691197 2.261369e-04 0.7149992 -1.5881724
> 42018 42018 -0.7514740 11.479166 -5.626515 2.467244e-04 0.7149992 -1.6142433
> 41399 41399  0.4513182  8.270799  5.597586 2.565770e-04 0.7149992 -1.6260629
> 31680 31680 -0.9018399  6.573984 -5.550571 2.735040e-04 0.7149992 -1.6454852
> 3019   3019 -0.5514822 11.274571 -5.533892 2.797945e-04 0.7149992 -1.6524395
> 35615 35615  0.8136739  6.947605  5.475631 3.030174e-04 0.7149992 -1.6769968
> 
>> top3 <- topTable(fitContr2, coef= 3, number = 20, adjust= "BH", sort.by="p")
>> top3
>          ID      logFC  AveExpr         t      P.Value adj.P.Val         B
> 39854 39854 -0.2387620 5.392370 -4.855615 0.0007294172 0.9999495 -1.982072
> 24894 24894  0.3724066 6.205436  4.835843 0.0007508215 0.9999495 -1.992056
> 14122 14122 -0.2238631 5.325398 -4.710873 0.0009025967 0.9999495 -2.056463
> 32475 32475 -0.2275235 5.393603 -4.687553 0.0009343768 0.9999495 -2.068734
> 37028 37028 -0.2340294 6.164496 -4.671185 0.0009573938 0.9999495 -2.077396
> 12902 12902  0.2264364 5.285410  4.604682 0.0010573110 0.9999495 -2.112996
> 32023 32023 -0.2136657 5.541409 -4.549203 0.0011491544 0.9999495 -2.143203
> 12908 12908  0.1898669 5.417419  4.518870 0.0012029254 0.9999495 -2.159916
> 7786   7786 -0.1972366 5.218460 -4.481958 0.0012719856 0.9999495 -2.180443
> 26520 26520 -0.2213202 5.278165 -4.452604 0.0013299104 0.9999495 -2.196916
> 4727   4727 -0.2274668 5.314431 -4.436143 0.0013636118 0.9999495 -2.206211
> 44298 44298  0.2378503 6.214914  4.394183 0.0014536867 0.9999495 -2.230096
> 18945 18945 -0.2495693 5.582078 -4.347342 0.0015617426 0.9999495 -2.257083
> 1900   1900 -0.2107123 5.414634 -4.244593 0.0018297003 0.9999495 -2.317484
> 46220 46220 -0.1773682 5.446607 -4.234698 0.0018579597 0.9999495 -2.323388
> 29257 29257 -0.2125592 5.599029 -4.229925 0.0018717555 0.9999495 -2.326242
> 4590   4590  0.1954713 5.340609  4.175245 0.0020377623 0.9999495 -2.359194
> 19285 19285  0.1960713 5.378438  4.139800 0.0021536291 0.9999495 -2.380809
> 15867 15867  0.1582910 5.298160  4.127889 0.0021941097 0.9999495 -2.388117
> 33071 33071 -0.2201014 5.464051 -4.103212 0.0022805505 0.9999495 -2.403331
>> ########################################################################
>  
>> ######################## Block #######################################
>> SS <- read.table(file = "GEXSampleSheet.csv",sep = ",",colClasses = "character",header = T,skip = 8)
>> SS
>    Sample_Name Sample_Well Sample_Plate Sample_Group Sentrix_ID Pool_ID Sentrix_Position X X.1 X.2
> 1         1_NG                                    NG 4254964042    LD29                A          
> 2         7_HG                                    HG 4254964042    LD30                B          
> 3   2_NG_benfo                              NG_benfo 4254964042    LD29                C          
> 4   8_HG_benfo                              HG_benfo 4254964042    LD30                D          
> 5         3_HG                                    HG 4254964042    LD29                E          
> 6         9_NG                                    NG 4254964042    LD34                F          
> 7   4_HG_benfo                              HG_benfo 4254964032    LD29                A          
> 8  10_NG_benfo                              NG_benfo 4254964032    LD34                B          
> 9         5_NG                                    NG 4254964032    LD30                C          
> 10       11_HG                                    HG 4254964032    LD34                D          
> 11  6_NG_benfo                              NG_benfo 4254964032    LD30                E          
> 12 12_HG_benfo                              HG_benfo 4254964032    LD34                F          
>> dataSet_Norm_exp_log2_ordnet <- dataSet_Norm_exp_log2[,SS[,1]]
>>
>> S_gr <- SS[,4]
>> s_gr
> Error: object "s_gr" not found
>> S_gr_faktor <- as.factor(S_gr)
>> S_gr_faktor
>  [1] NG       HG       NG_benfo HG_benfo HG       NG       HG_benfo NG_benfo NG       HG       NG_benfo HG_benfo
> Levels: HG HG_benfo NG NG_benfo
>> designMa <- model.matrix(~-1+S_gr_faktor)
>> colnames(designMa) <- levels(S_gr_faktor)
>>
>> rownames(designMa) <- S_gr_faktor
>>
>> blokk <- SS[,6]
>> blokk
>  [1] "LD29" "LD30" "LD29" "LD30" "LD29" "LD34" "LD29" "LD34" "LD30" "LD34" "LD30" "LD34"
>> corfit <- duplicateCorrelation(dataSet_Norm_exp_log2_ordnet, design = designMa, ndups = 1, block = as.factor(blokk))
>>
>> fitDesMa <- lmFit(dataSet_Norm_exp_log2, design = designMa,block = as.factor(blokk),cor = corfit$consensus)
>>
>> designMa
>          HG HG_benfo NG NG_benfo
> NG        0        0  1        0
> HG        1        0  0        0
> NG_benfo  0        0  0        1
> HG_benfo  0        1  0        0
> HG        1        0  0        0
> NG        0        0  1        0
> HG_benfo  0        1  0        0
> NG_benfo  0        0  0        1
> NG        0        0  1        0
> HG        1        0  0        0
> NG_benfo  0        0  0        1
> HG_benfo  0        1  0        0
> attr(,"assign")
> [1] 1 1 1 1
> attr(,"contrasts")
> attr(,"contrasts")$S_gr_faktor
> [1] "contr.treatment"
> 
>> contrast.matrix <- makeContrasts(NG - NG_benfo,HG - HG_benfo, NG-HG, levels = designMa)                               
>>
>> fitContr <- contrasts.fit(fitDesMa, contrast.matrix)
>>
>> fitContr2 <- eBayes(fitContr)
>>
>> top1 <- topTable(fitContr2, coef= 1, number = 20, adjust= "BH", sort.by="p")
>> top1
>          ID      logFC  AveExpr         t      P.Value adj.P.Val         B
> 7611   7611 -0.2527342 5.594157 -5.439053 0.0003198243 0.9999325 -1.521966
> 32680 32680  0.3197339 5.304647  5.213128 0.0004380912 0.9999325 -1.630612
> 13175 13175  0.2373210 5.436111  4.912886 0.0006730955 0.9999325 -1.786076
> 12484 12484 -0.2244960 5.545403 -4.887328 0.0006985769 0.9999325 -1.799920
> 48359 48359  0.2393309 5.316981  4.863403 0.0007233641 0.9999325 -1.812967
> 14456 14456  0.2406121 5.618726  4.858584 0.0007284695 0.9999325 -1.815605
> 21982 21982 -0.2415139 5.354844 -4.799576 0.0007942028 0.9999325 -1.848198
> 4668   4668  0.2189805 5.294603  4.752691 0.0008509389 0.9999325 -1.874474
> 24165 24165 -0.2138513 5.438563 -4.723256 0.0008887563 0.9999325 -1.891143
> 8228   8228 -0.3325391 6.138562 -4.615691 0.0010429202 0.9999325 -1.953203
> 23640 23640  0.2089381 5.580925  4.596901 0.0010726550 0.9999325 -1.964230
> 12265 12265  0.1636480 5.471247  4.590664 0.0010827233 0.9999325 -1.967903
> 16029 16029  0.2862802 6.143381  4.584987 0.0010919753 0.9999325 -1.971252
> 26564 26564  0.2162782 5.452064  4.572604 0.0011124478 0.9999325 -1.978572
> 25896 25896  0.2122054 5.509182  4.568471 0.0011193727 0.9999325 -1.981022
> 12563 12563  0.2087823 5.555015  4.517547 0.0012085463 0.9999325 -2.011420
> 29041 29041 -0.1691662 5.467616 -4.501427 0.0012383247 0.9999325 -2.021129
> 8086   8086 -0.2228617 6.070106 -4.489329 0.0012611850 0.9999325 -2.028443
> 286     286  0.1975515 6.102632  4.480072 0.0012789776 0.9999325 -2.034055
> 47658 47658  0.2078376 5.293210  4.459124 0.0013202373 0.9999325 -2.046807
> 
>> top2 <- topTable(fitContr2, coef= 2, number = 20, adjust= "BH", sort.by="p")
>> top2
>          ID      logFC  AveExpr         t      P.Value adj.P.Val         B
> 23091 23091  0.2354301 5.244294  5.224066 0.0004313957   0.99981 -1.625191
> 32259 32259  0.2518479 5.332262  5.165712 0.0004684309   0.99981 -1.654305
> 31838 31838 -0.2355327 5.340717 -5.126989 0.0004948792   0.99981 -1.673889
> 44090 44090 -0.2251488 5.383418 -5.040692 0.0005597503   0.99981 -1.718304
> 45163 45163  0.2894191 8.376503  4.892174 0.0006936669   0.99981 -1.797287
> 33415 33415  0.2636604 5.306506  4.855887 0.0007313444   0.99981 -1.817084
> 43956 43956  0.2483797 5.346141  4.831985 0.0007573518   0.99981 -1.830232
> 48396 48396  0.1977520 5.400935  4.678863 0.0009492160   0.99981 -1.916536
> 35709 35709 -0.5484775 9.632759 -4.643015 0.0010012341   0.99981 -1.937266
> 23341 23341  0.1741108 5.346627  4.563827 0.0011272060   0.99981 -1.983776
> 41288 41288  0.3055193 5.338962  4.538377 0.0011711893   0.99981 -1.998935
> 24890 24890 -0.2025487 5.238161 -4.496772 0.0012470673   0.99981 -2.023940
> 47443 47443  0.1997024 5.298141  4.381631 0.0014855767   0.99981 -2.094596
> 13408 13408 -0.2320642 5.335356 -4.380533 0.0014880710   0.99981 -2.095280
> 19419 19419  0.2029474 5.785537  4.364322 0.0015254257   0.99981 -2.105404
> 33589 33589  0.1680122 5.321035  4.329689 0.0016086005   0.99981 -2.127177
> 21838 21838  0.2162749 6.009970  4.308605 0.0016615763   0.99981 -2.140528
> 28720 28720  0.2007788 5.305068  4.288266 0.0017144334   0.99981 -2.153476
> 32627 32627 -0.1571718 5.386997 -4.279245 0.0017384451   0.99981 -2.159241
> 36979 36979 -0.1851754 5.272498 -4.249410 0.0018204265   0.99981 -2.178404
> 
>> top3 <- topTable(fitContr2, coef= 3, number = 20, adjust= "BH", sort.by="p")
>> top3
>          ID      logFC  AveExpr         t      P.Value adj.P.Val         B
> 20549 20549 -0.2562633 5.297783 -5.787760 0.0001995769 0.9998617 -1.367343
> 46835 46835 -0.2730167 5.579581 -5.741662 0.0002122093 0.9998617 -1.386916
> 18794 18794  0.2400783 5.352057  5.411186 0.0003323510 0.9998617 -1.534995
> 15282 15282 -0.3525527 5.574132 -5.284227 0.0003964783 0.9998617 -1.595671
> 19419 19419 -0.2395820 5.785537 -5.152140 0.0004775244 0.9998617 -1.661145
> 31093 31093  0.2084498 5.679371  5.033290 0.0005657241 0.9998617 -1.722164
> 14231 14231  0.1895001 5.504983  4.938160 0.0006488719 0.9998617 -1.772483
> 13615 13615  0.2638662 5.282304  4.835343 0.0007536390 0.9998617 -1.828380
> 24296 24296  0.1750666 5.267973  4.828252 0.0007615022 0.9998617 -1.832294
> 14926 14926 -0.2954851 5.358565 -4.734460 0.0008741543 0.9998617 -1.884782
> 29229 29229  0.2305711 5.673516  4.731604 0.0008778519 0.9998617 -1.886402
> 31572 31572 -0.1835065 5.353449 -4.668348 0.0009641692 0.9998617 -1.922596
> 10319 10319  0.2386456 6.346361  4.649946 0.0009909460 0.9998617 -1.933243
> 32211 32211  0.2222265 5.483684  4.619019 0.0010377456 0.9998617 -1.951256
> 20957 20957  0.1950235 5.440708  4.599657 0.0010682380 0.9998617 -1.962610
> 31581 31581 -0.2053997 5.378321 -4.592595 0.0010795955 0.9998617 -1.966765
> 31817 31817  0.1761750 5.382717  4.575112 0.0011082689 0.9998617 -1.977088
> 12484 12484 -0.2067044 5.545403 -4.500002 0.0012409931 0.9998617 -2.021989
> 29205 29205  0.2347584 5.334616  4.457761 0.0013229690 0.9998617 -2.047639
> 32048 32048  0.1927972 5.312840  4.436783 0.0013658003 0.9998617 -2.060484
> ######################################################################## 
> 
>> ############################## Pair ####################################
>> targets <- readTargets("targets2.txt")
>> targets
>    FileName Donor Treatment
> 1     File1    29        NG
> 2     File2    34        NG
> 3     File3    30        NG
> 4     File4    30        HG
> 5     File5    29        HG
> 6     File6    34        HG
> 7     File7    29   NGbenfo
> 8     File8    34   NGbenfo
> 9     File9    30   NGbenfo
> 10   File10    30   HGbenfo
> 11   File11    29   HGbenfo
> 12   File12    34   HGbenfo
>> Target <- factor(targets$Donor)
>> Target
>  [1] 29 34 30 30 29 34 29 34 30 30 29 34
> Levels: 29 30 34
>> Treat <- factor(targets$Treatment, levels=c("NG","HG","NGbenfo","HGbenfo"))
>> Treat
>  [1] NG      NG      NG      HG      HG      HG      NGbenfo NGbenfo NGbenfo HGbenfo HGbenfo HGbenfo
> Levels: NG HG NGbenfo HGbenfo
>> design <- model.matrix(~-1+Treat)
>> design
>    TreatNG TreatHG TreatNGbenfo TreatHGbenfo
> 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")$Treat
> [1] "contr.treatment"
> 
>> fit <- lmFit(dataSet_Norm_exp_log2, design)
>> fit2 <- eBayes(fit)
>>
>> top1 <- topTable(fit2, coef= 1, number = 20, adjust= "BH", sort.by="p")
>> top1
>          ID    logFC  AveExpr        t      P.Value    adj.P.Val        B
> 29377 29377 14.91596 14.92944 516.8806 7.735991e-23 3.558893e-18 35.95591
> 7792   7792 14.77356 14.80328 467.7334 2.036389e-22 3.558893e-18 35.77056
> 46748 46748 14.54269 14.56708 453.1413 2.768223e-22 3.558893e-18 35.70514
> 7793   7793 14.90254 14.90698 441.9898 3.523945e-22 3.558893e-18 35.65132
> 28882 28882 14.28753 14.29856 435.5766 4.059943e-22 3.558893e-18 35.61875
> 33208 33208 14.14750 14.13732 426.6177 4.965377e-22 3.558893e-18 35.57112
> 28580 28580 14.87370 14.83600 424.1510 5.252285e-22 3.558893e-18 35.55755
> 41917 41917 14.03380 14.05301 419.5783 5.833773e-22 3.558893e-18 35.53186
> 4769   4769 14.41497 14.44560 410.1179 7.276004e-22 3.945535e-18 35.47637
> 41907 41907 13.77121 13.74080 400.9056 9.067137e-22 4.313991e-18 35.41912
> 26904 26904 14.29154 14.33758 398.0242 9.723362e-22 4.313991e-18 35.40053
> 46749 46749 14.50236 14.47209 393.3830 1.089334e-21 4.430320e-18 35.36986
> 29730 29730 14.34559 14.32901 383.3841 1.397896e-21 5.020324e-18 35.30060
> 42046 42046 14.48283 14.51144 380.4832 1.504625e-21 5.020324e-18 35.27966
> 43841 43841 13.71735 13.74382 378.2740 1.591940e-21 5.020324e-18 35.26344
> 41965 41965 14.45261 14.38145 376.9752 1.645873e-21 5.020324e-18 35.25380
> 42047 42047 14.37644 14.35028 372.5960 1.843121e-21 5.216996e-18 35.22067
> 42019 42019 14.27887 14.32325 370.1575 1.964159e-21 5.216996e-18 35.20180
> 32161 32161 13.95997 14.02817 368.8802 2.031041e-21 5.216996e-18 35.19179
> 42021 42021 14.23735 14.15566 365.6109 2.213963e-21 5.402512e-18 35.16579
> 
>> top2 <- topTable(fit2, coef= 2, number = 20, adjust= "BH", sort.by="p")
>> top2
>          ID    logFC  AveExpr        t      P.Value    adj.P.Val        B
> 29377 29377 14.94293 14.92944 517.8155 7.601752e-23 3.468187e-18 35.95899
> 7792   7792 14.81788 14.80328 469.1366 1.978147e-22 3.468187e-18 35.77657
> 46748 46748 14.56020 14.56708 453.6869 2.736145e-22 3.468187e-18 35.70769
> 7793   7793 14.90254 14.90698 441.9898 3.523945e-22 3.468187e-18 35.65132
> 28882 28882 14.31171 14.29856 436.3137 3.993984e-22 3.468187e-18 35.62256
> 33208 33208 14.14569 14.13732 426.5633 4.971519e-22 3.468187e-18 35.57082
> 28580 28580 14.78715 14.83600 421.6828 5.557769e-22 3.468187e-18 35.54377
> 41917 41917 14.07125 14.05301 420.6980 5.685086e-22 3.468187e-18 35.53822
> 4769   4769 14.39271 14.44560 409.4848 7.385709e-22 4.005024e-18 35.47254
> 41907 41907 13.70521 13.74080 398.9841 9.499104e-22 4.251956e-18 35.40676
> 26904 26904 14.31293 14.33758 398.6198 9.583542e-22 4.251956e-18 35.40440
> 46749 46749 14.44334 14.47209 391.7820 1.133228e-21 4.608840e-18 35.35906
> 29730 29730 14.32854 14.32901 382.9285 1.414093e-21 5.099289e-18 35.29734
> 42046 42046 14.51584 14.51144 381.3506 1.471801e-21 5.099289e-18 35.28596
> 43841 43841 13.70964 13.74382 378.0612 1.600642e-21 5.099289e-18 35.26187
> 41965 41965 14.39025 14.38145 375.3486 1.716279e-21 5.099289e-18 35.24160
> 42019 42019 14.38992 14.32325 373.0365 1.822145e-21 5.099289e-18 35.22404
> 42047 42047 14.28467 14.35028 370.2178 1.961065e-21 5.099289e-18 35.20227
> 32161 32161 13.99289 14.02817 369.7503 1.985216e-21 5.099289e-18 35.19862
> 42021 42021 14.15355 14.15566 363.4590 2.344258e-21 5.720459e-18 35.14836
> 
>> top3 <- topTable(fit2, coef= 3, number = 20, adjust= "BH", sort.by="p")
>> top3
>          ID    logFC  AveExpr        t      P.Value    adj.P.Val        B
> 29377 29377 14.96991 14.92944 518.7503 7.470078e-23 3.536707e-18 35.96205
> 7792   7792 14.83129 14.80328 469.5612 1.960886e-22 3.536707e-18 35.77838
> 46748 46748 14.60366 14.56708 455.0410 2.658281e-22 3.536707e-18 35.71397
> 7793   7793 14.85291 14.90698 440.5176 3.639697e-22 3.536707e-18 35.64396
> 28882 28882 14.30565 14.29856 436.1290 4.010404e-22 3.536707e-18 35.62160
> 33208 33208 14.18397 14.13732 427.7176 4.843054e-22 3.536707e-18 35.57711
> 28580 28580 14.84366 14.83600 423.2943 5.356163e-22 3.536707e-18 35.55279
> 41917 41917 14.04286 14.05301 419.8492 5.797406e-22 3.536707e-18 35.53340
> 4769   4769 14.48384 14.44560 412.0774 6.947685e-22 3.767498e-18 35.48813
> 41907 41907 13.69640 13.74080 398.7278 9.558428e-22 4.265616e-18 35.40510
> 26904 26904 14.30819 14.33758 398.4878 9.614331e-22 4.265616e-18 35.40354
> 46749 46749 14.47947 14.47209 392.7620 1.106132e-21 4.498641e-18 35.36568
> 29730 29730 14.33028 14.32901 382.9748 1.412436e-21 5.031474e-18 35.29767
> 42046 42046 14.54513 14.51144 382.1202 1.443337e-21 5.031474e-18 35.29153
> 43841 43841 13.75475 13.74382 379.3054 1.550501e-21 5.042971e-18 35.27104
> 41965 41965 14.35647 14.38145 374.4674 1.755804e-21 5.042971e-18 35.23494
> 32161 32161 14.08381 14.02817 372.1527 1.864498e-21 5.042971e-18 35.21726
> 42047 42047 14.34979 14.35028 371.9053 1.876547e-21 5.042971e-18 35.21535
> 42019 42019 14.27952 14.32325 370.1744 1.963291e-21 5.042971e-18 35.20193
> 42021 42021 14.15562 14.15566 363.5122 2.340936e-21 5.712351e-18 35.14880
> ########################################################################
> 
>> sessionInfo()
> R version 2.7.1 (2008-06-23) 
> 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] splines   tools     stats     graphics  grDevices utils     datasets  methods   base     
> 
> other attached packages:
>  [1] illuminaHumanv3ProbeID.db_1.1.1 statmod_1.3.6                   GOstats_2.6.0                   Category_2.6.0                  genefilter_1.20.0               survival_2.34-1               
>  [7] RBGL_1.16.0                     graph_1.18.1                    annaffy_1.12.1                  KEGG.db_2.2.0                   GO.db_2.2.0                     annotate_1.18.0               
> [13] AnnotationDbi_1.2.2             RSQLite_0.6-9                   DBI_0.2-4                       xtable_1.5-2                    RColorBrewer_1.0-2              limma_2.14.5                   
> [19] lumi_1.6.2                      mgcv_1.4-0                      affy_1.18.2                     preprocessCore_1.2.0            affyio_1.8.0                    Biobase_2.0.1                  
> 
> loaded via a namespace (and not attached):
> [1] cluster_1.11.11
> 
> 
> 
> 	[[alternative HTML version deleted]]
> 
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at stat.math.ethz.ch
> 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
Hildebrandt Lab
8220D MSRB III
1150 W. Medical Center Drive
Ann Arbor MI 48109-0646
734-936-8662



More information about the Bioconductor mailing list