=== Please help ======== This is was a code used to analyze DEG of an microarray data All steps are working fine But error in the final step while trying to write the result to .txt file

ERROR

> write.fit(fit2, file = "de-results.txt",adjust="BH")
Error in (function (..., row.names = NULL, check.rows = FALSE, check.names = TRUE,  : 
  arguments imply differing number of rows: 21801, 27416

library(GEOquery)
library(limma)
filenm <- "GSE39956_series_matrix.txt.gz"
gse <- getGEO(filename=filenm)
gse
head(exprs(gse))
boxplot(exprs(gse),outline=FALSE)
pd <- pData(gse)
library(genefilter)
design <- matrix(nrow=21,ncol=2)
design[,1] <- c(1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
design[,2] <- c(0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
colnames(design) <- c("Control","Salt300")
gse.expFilt <- varFilter(gse)
fit <- lmFit(exprs(gse.expFilt), design)
names(fit)
dim(fit$coefficients)
head(fit$coefficients)
contrasts <- makeContrasts(logFC = Salt300 - Control, levels=design)
fit2 <- contrasts.fit(fit, contrasts)
head(fit2$coeff)
fit2 <- eBayes(fit2)
fit2
topTable(fit2)
library(org.At.tair.db)
anno <- data.frame(SYMBOL=sapply(contents(org.At.tairSYMBOL), paste, collapse=", "), DESC=sapply(contents(org.At.tairGENENAME), paste, collapse=", "))
fit2$genes <- anno
topTable(fit2)
decideTests(fit2)
table(decideTests(fit2))
sum(abs(decideTests(fit2))==1)
volcanoplot(fit2)
volcanoplot(fit2,highlight=10,names = fit2$genes$"Symbol")
write.fit(fit2, file = "control_salt300.txt",adjust="BH")

To add to what Gordon Smyth has already said, the way you are annotating your data is woefully naive. You are assuming that you can just dump the contents from the OrgDb into a data.frame and assume that the ordering of your annotation will magically be the same as the ordering of your genes. This is not a good assumption! You are in essence assuming that the OrgDb has all the same data as this Agilent array, and in the same order that you have, after filtering the data! I don't know why you would assume that, but it's not a reasonable assumption at all.

Also, do note that the data you got from GEO are conveniently already annotated in such a way that limma will happily use the data. All you really need to do is remove some of the extraneous columns.

> gse <- getGEO("GSE39956")[[1]]
> fData(gse) <- fData(gse)[,c("REFSEQ","GENE_SYMBOL","GENE_NAME")]
> head(fData(gse))
               REFSEQ GENE_SYMBOL
A_84_P10000                      
A_84_P10001 NM_102184   AT1G23350
A_84_P10002 NM_116719       MSRB4
A_84_P10003 NM_116748   AT4G05090
A_84_P10004 NM_116765   AT4G05260
A_84_P10005 NM_116785   AT4G05460
                                                                            GENE_NAME
A_84_P10000                                                                          
A_84_P10001 plant invertase/pectin methylesterase inhibitor domain-containing protein
A_84_P10002                                 peptide methionine sulfoxide reductase B4
A_84_P10003                                         putative PAP-specific phosphatase
A_84_P10004                                                  ubiquitin family protein
A_84_P10005                                                      F-box protein SKIP19

In addition, do note that using varFilter prior to eBayes is not a good idea. The idea behind eBayes is that you can get a better estimate of the within-sample variance by using all the genes. If you filter out the low variance genes, then you are likely to be artificially increasing the variance estimate that eBayes is computing, which may have the effect of causing your tests to be too conservative.

The data.frame anno has more rows than fit2 The rows of anno should match the genes in your dataset.