Problem with write.fit function
2
0
Entering edit mode
@tillmannruland-14225
Last seen 5.9 years ago
Germany/Münster/Department of Psychiatry

 

Dear community,

i have an issue with the write.fit function which I use to create a .txt of a limma object because the summary(results) differs from the values in the .txt independent of the 'adjust.method' or the 'method' I choose. (tried different options in decideTests as well as write.fit). write.fit doesn't respond to any change. 

Compared to summary(result) the number of significant genes in contrasts can be reduced by a hundred genes when a file is created with write.fit.  

.....

fit2 <- eBayes(fit2)

results <- decideTests(fit2, method = "global", adjust.method = "BH")

write.fit(fit2, results, "limma_results2.txt", adjust="none", method ="separate")

 

I'm looking forward to your comments!

Tillmann

 

R version 3.5.0 (2018-04-23)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Sierra 10.12.6

Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] de_DE.UTF-8/de_DE.UTF-8/de_DE.UTF-8/C/de_DE.UTF-8/de_DE.UTF-8

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

other attached packages:
[1] data.table_1.11.4    limma_3.36.1         qvalue_2.12.0        BiocInstaller_1.30.0

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.17     grid_3.5.0       plyr_1.8.4       gtable_0.2.0     magrittr_1.5    
 [6] scales_0.5.0     ggplot2_2.2.1    pillar_1.2.3     stringi_1.2.3    rlang_0.2.1     
[11] reshape2_1.4.3   lazyeval_0.2.1   splines_3.5.0    tools_3.5.0      stringr_1.3.1   
[16] munsell_0.5.0    compiler_3.5.0   colorspace_1.3-2 tibble_1.4.2  

 

 

limma write.fit • 1.1k views
ADD COMMENT
0
Entering edit mode
Aaron Lun ★ 28k
@alun
Last seen 9 hours ago
The city by the bay

Works fine for me.

library(limma)
g <- gl(3, 2)
y <- matrix(rnorm(length(g) * 10000), ncol=6)

design <- model.matrix(~g)
fit <- lmFit(y, design)
fit <- eBayes(fit)

results <- decideTests(fit, method = "global", adjust.method = "BH")
write.fit(fit, results, "limma_results1.txt", adjust="none", method = "separate")
write.fit(fit, results, "limma_results2.txt", adjust="BH", method ="separate")

first <- read.delim("limma_results1.txt", header=TRUE, check.names=FALSE)
colnames(first)
##  [1] "A"                   "Coef.(Intercept)"    "Coef.g2"            
##  [4] "Coef.g3"             "t.(Intercept)"       "t.g2"               
##  [7] "t.g3"                "p.value.(Intercept)" "p.value.g2"         
## [10] "p.value.g3"          "F"                   "F.p.value"          
## [13] "Res.(Intercept)"     "Res.g2"              "Res.g3"             
second <- read.delim("limma_results2.txt", header=TRUE, check.names=FALSE)
##  [1] "A"                       "Coef.(Intercept)"       
##  [3] "Coef.g2"                 "Coef.g3"                
##  [5] "t.(Intercept)"           "t.g2"                   
##  [7] "t.g3"                    "p.value.(Intercept)"    
##  [9] "p.value.g2"              "p.value.g3"             
## [11] "p.value.adj.(Intercept)" "p.value.adj.g2"         
## [13] "p.value.adj.g3"          "F"                      
## [15] "F.p.value"               "Res.(Intercept)"        
## [17] "Res.g2"                  "Res.g3"

You can see that if you turn adjust off in write.fit, you don't even get the adjusted p-value columns in first. And obviously, you will get different adjustments between write.fit and decideTests if you change method=.

ADD COMMENT
0
Entering edit mode

Thank you for your answer!

Using the following code which Gordon (Smyth) suggested in a different post, I could get the same results as in summary(results).

tab <- topTable(fit2,coef=n n=Inf,sort="none")

tabUp <- tab[results[,n]<0,]
tabDn <- tab[results[,n]>0,]

For some reason in my write.fit .txt there is an issue with the column names. They appear to be shifted meaning that in my summary(first) there are for example 15 columns but in corresponding write.fit file I have 16 columns and the 16th column containing Res.gx values has no column name.

ADD REPLY
1
Entering edit mode

That's because write.fit is writing out the rownames, and the rownames don't have a column heading. If you don't want to write out rownames, use row.names=FALSE in your call to write.fit.

ADD REPLY
0
Entering edit mode

Thank you for your answer. 

ADD REPLY

Login before adding your answer.

Traffic: 887 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

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

Powered by the version 2.3.6