Help with limma contrast/design matrices
Entering edit mode
mat149 ▴ 70
Last seen 13 months ago
United States


I was hoping someone could give this script a look-through and determine if I have made an error at some point.  This regards the analysis of Affymetrix 1.1 st zebrafish gene arrays.

I have three groups of samples and I am detecting differentially expressed genes between all three contrasts (1: morphant-control 2: rescue-control 3: morphant-rescue).  I have 8x control samples, 4x treatment#1 (morphant) samples, and 4x treatment#2 (rescue) samples.  I am using limma's lmFit and eBayes for this.  When I get my topTable results back, the fold change signs are inverted. The genes that should be downregulated are returned as positive logFC integers and vice versa . I know from analyzing this data previously that, for instance, KEGG pathway enrichment for phototransduction should be downregulated but in this case it will be returned as upregulated.

It seems like a silly inquiry, but I cant figure it out.



#### read samples in

CELlist <- list.celfiles("C:\\Users\\mat149\\Desktop\\GG\\CEL", full.names=TRUE, pattern=NULL, all.file=FALSE)
CELdat <- read.celfiles(filenames = CELlist,experimentData=TRUE,phenoData=pdat,verbose=TRUE)

### RMA

eset<-rma(CELdat, background=TRUE, normalize=TRUE, subset=NULL, target="core")     
eset <- annotateEset(eset, annotation(eset))

fd <- fData(eset)
fd$ENTREZID <- mapIds(, as.character(fd$SYMBOL), "ENTREZID","SYMBOL")
fData(eset) <- fd




ph = CELdat@phenoData 
ph@data[ ,1] = c("WT1","WT2","WT3","WT4","WT5","WT6","WT7","WT8","MO1","MO2","MO3","MO4","RS1","RS2","RS3","RS4")
ph@data[ ,2] = c("control","control","control","control","control","control","control","control","morphant","morphant","morphant","morphant","rescue","rescue","rescue","rescue")
groups = ph@data$Treatment
f = factor(groups,levels=c("control","morphant","rescue"))


#### LIMMA

design = model.matrix(~ 0 + f)
design = lmFit(eset,design)

contrast.matrix = makeContrasts(morphant-control,rescue-control,morphant-rescue,levels=design) =,contrast.matrix) = eBayes(





> sessionInfo()
R version 3.3.1 (2016-06-21)
Platform: i386-w64-mingw32/i386 (32-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1

[1] LC_COLLATE=English_United States.1252 
[2] LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

attached base packages:
[1] stats4    parallel  stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] limma_3.28.21         AnnotationDbi_1.34.4    
 [4] affycoretools_1.44.3     pd.zebgene.1.1.st_3.12.0 RSQLite_1.0.0           
 [7] DBI_0.5-1                oligo_1.36.1             Biostrings_2.40.2       
[10] XVector_0.12.1           IRanges_2.6.1            S4Vectors_0.10.3        
[13] Biobase_2.32.0           oligoClasses_1.34.0      BiocGenerics_0.18.0     

loaded via a namespace (and not attached):
 [1] httr_1.2.1                    munsell_0.4.3                
 [3] latticeExtra_0.6-28           BSgenome_1.40.1              
 [5] dichromat_2.0-0               R.utils_2.5.0                
 [7] PFAM.db_3.3.0                 httpuv_1.3.3                 
 [9] R6_2.2.0                      ensembldb_1.4.7              
[11] graph_1.50.0                  affxparser_1.44.0            
[13] BiocInstaller_1.22.3          data.table_1.9.6             
[15] reshape_0.8.6                 annotate_1.50.1              
[17] xtable_1.8-2                  gdata_2.17.0                 
[19] tools_3.3.1                   stringr_1.1.0                
[21] rtracklayer_1.32.2            mime_0.5                     
[23] GSEABase_1.34.1               shiny_0.14.2                 
[25] chron_2.3-47                  R.oo_1.21.0                  
[27] GOstats_2.38.1                foreach_1.4.3                
[29] digest_0.6.10                 KernSmooth_2.23-15           
[31] GO.db_3.3.0                   GenomeInfoDb_1.8.7           
[33] codetools_0.2-15              GGally_1.2.0                 
[35] ff_2.2-13                     GenomicAlignments_1.8.4      
[37] gplots_3.0.1                  genefilter_1.54.2            
[39] scales_0.4.1                  stringi_1.1.2                
[41] locfit_1.5-9.1                R.methodsS3_1.7.1            
[43] gcrma_2.44.0                  lattice_0.20-33              
[45] AnnotationForge_1.14.2        interactiveDisplayBase_1.10.3
[47] biovizBase_1.20.0             Rcpp_0.12.7                  
[49] OrganismDbi_1.14.1            caTools_1.17.1               
[51] Hmisc_4.0-0                   Formula_1.2-1                
[53] ggplot2_2.1.0                 htmlTable_1.7                
[55] Category_2.38.0               grid_3.3.1                   
[57] ReportingTools_2.12.2         GenomicRanges_1.24.3         
[59] preprocessCore_1.34.0         plyr_1.8.4                   
[61] RBGL_1.48.1                   survival_2.39-4              
[63] edgeR_3.14.0                  acepack_1.4.1                
[65] affy_1.50.0                   rpart_4.1-10                 
[67] magrittr_1.5                  SummarizedExperiment_1.2.3   
[69] VariantAnnotation_1.18.7      gridExtra_2.2.1              
[71] affyio_1.42.0                 biomaRt_2.28.0               
[73] htmltools_0.3.5               ggbio_1.20.2                 
[75] knitr_1.15                    nnet_7.3-12                  
[77] gtable_0.2.0                  zlibbioc_1.18.0              
[79] colorspace_1.2-7              geneplotter_1.50.0           
[81] cluster_2.0.4                 gtools_3.5.0                 
[83] RCurl_1.95-4.8                DESeq2_1.12.4                
[85] bitops_1.0-6                  RColorBrewer_1.1-2           
[87] Matrix_1.2-6                  foreign_0.8-66               
[89] bit_1.1-12                    GenomicFeatures_1.24.5       
[91] hwriter_1.3.2                 reshape2_1.4.2               
[93] XML_3.98-1.4                  AnnotationHub_2.4.2          
[95] iterators_1.0.8               splines_3.3.1                
[97] BiocParallel_1.6.6            Rsamtools_1.24.0

limma contrast matrix • 588 views
Entering edit mode

Several points:

  1. Make a limma tag, otherwise the maintainers don't get notified.
  2. Post a specific example of a gene that's behaving in the unexpected manner that you describe, i.e., its expression values across all samples and the log-fold changes that are reported by limma. The code looks (mostly) fine so, right now, there's only your word that something's amiss.
Entering edit mode

In addition: take a gene that "you know" the sign is inverted for and plot its expression from your data among the two groups implicated in your contrast.

I agree with Aaron in that your code looks mostly correct, so apart from perhaps a mislabelling of the conditions to samples, positive logFCs in your coef=1 result will give you genes higher expressed in "morphant" compared to their expression levels in "control"


Login before adding your answer.

Traffic: 436 users visited in the last hour
Help About
Access RSS

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6