diff --git a/script/run_deseq.R b/script/run_deseq.R
index f55bf9a..fd2e68b 100644
--- a/script/run_deseq.R
+++ b/script/run_deseq.R
@@ -10,6 +10,7 @@ sizeFactorsFile  <- Args[7]
 qcPdfFile        <- Args[8]
 filterPercentile <- as.numeric( Args[9] )
 
+
 # Get data and samples
 countData     <- read.table(   countFile, header=TRUE, row.names=1 )
 samples       <- read.table( samplesFile, header=TRUE, row.names=1 )
@@ -42,7 +43,7 @@ if (filterPercentile) {
 # Create DESeqDataSet (with design according to number of factors)
 dds <- DESeqDataSetFromMatrix(countData, samples, design = ~ condition)
 if (numFactors == 2) {
-    design(dds) <- formula(~ group + condition)
+    design(dds) <- formula(~ group * condition)
 }
 
 # Ensure control level (usually "sibling") is first level (i.e. before "mutant")
@@ -51,14 +52,34 @@ colData(dds)$condition <- factor(colData(dds)$condition,
 
 # Differential expression analysis
 dds <- DESeq(dds)
-res <- results(dds)
+
+# Find the condition term (does not contain '.' and constains 'condition')
+results_names=resultsNames(dds)
+non_interaction_terms = which(grepl("\\.",results_names,perl =T) == FALSE)
+all_condition_terms_index = which(grepl("condition",results_names) == TRUE)
+condition_terms_index = intersect(non_interaction_terms,all_condition_terms_index)
+
+# create matrix with pvalues of all terms in the model
+pvalues_matrix = matrix(ncol=0,nrow=nrow(dds))
+col_names<-c();
+for (i  in 1:length(results_names) )
+{
+	equation_term = results_names[i]; 
+	results = results(dds, name= equation_term)
+	palues = as.matrix(results$pvalue)
+	padj =   as.matrix(results$padj)
+	pvalues_matrix = cbind(pvalues_matrix,palues)
+	pvalues_matrix = cbind(pvalues_matrix,padj)
+	col_names<-c(col_names,paste(equation_term,"_pvalue",sep=""),paste(equation_term,"_pvalue_adjusted",sep=""))
+}
+colnames(pvalues_matrix)=col_names;
+
+
 
 # Write output
-out <- data.frame(pvalue=res$pvalue, padj=res$padj, row.names=rownames(res))
-write.table( out, file=outputFile, col.names=FALSE, row.names=TRUE,
-    quote=FALSE, sep="\t" )
-write.table( sizeFactors( dds ), file=sizeFactorsFile, col.names=FALSE,
-    row.names=FALSE, quote=FALSE, sep="\t" )
+out <- data.frame(pvalues_matrix, row.names=rownames(res))
+write.table( out, file=outputFile, col.names=TRUE, row.names=TRUE, quote=FALSE, sep="\t" )
+write.table( sizeFactors( dds ), file=sizeFactorsFile, col.names=FALSE,row.names=FALSE, quote=FALSE, sep="\t" )
 
 # Data transformations for QC
 rld <- rlogTransformation(dds, blind=TRUE)