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]]

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