Thanks but I wouldn't ask if it was that simple, well I am not interested in the slight difference between treatments, my interest is in the general picture of the cell control vs treated(any treatment method)

Let suppose I didn't mention the treatment_methods and only add the cell_lines information to the second design would it be more accurate in this case compare to the first one even if I am only using the cell_class binary factor ( control/treated) as coef in topTable ?     

Thank you for taking the time Aaron, you resumed pretty well the situation in your code, I just simplified and I changed it a little to reflect the fact that I have a least 2 control per cell line, some cell lines are more represented than others and a few cell line have more than one treatment method.

cell_lines <- factor(c(rep(LETTERS[1:4], each=4),rep('D',4))) # Cell line A to D
treatment_method <- factor(c(0,0,1,1,0,0,1,1,0,0,2,2,0,0,2,2,0,0,1,1)) # Most cell lines contain 1 treatment method, control is 0
cell_class <- factor(treatment_method == 0) # the contrast I am mostly interest in

To answer your first question; the coefficients in the design matrix from the additive model:

design <- model.matrix(~cell_class + cell_lines + treatment_method)

... are not interpretable due to linear dependencies between cell_classTRUE, treatment_method1 and treatment_method2. If you want to interpret each coefficient, you have to discard one of them. limma chooses (more-or-less arbitrarily, in this case) to discard treatment_method2 for you, as lmFit doesn't try to estimate it.

Now, consider samples for cell line D. (I've chosen this cell line because it has both treatment conditions, which makes my point easy to demonstrate.) We describe the systematic component for each type of sample (i.e., control/treatment 1/treatment 2), based on whether the coefficient has a value of 1 or 0 in the design matrix for the row corresponding to that type of sample:

Control = Intercept + cell_classTRUE + cell_linesD
Treatment_1 = Intercept + treatment_method1 + cell_linesD
Treatment_2 = Intercept + cell_linesD

It should be fairly obvious that cell_classTRUE represents the difference (i.e., the log-fold change) between control and treatment 2, while treatment_method1 represents the difference between treatment 1 and treatment 2.

To answer your second question; probably yes, because your old model (using only cell_class) is likely to be flawed. I don't know what treatment methods you're using, but I find it hard to believe that different cell lines would not have systematic differences in behaviour. These differences will inflate the variance estimates in your old model.

To answer your third question; try to figure out what the values of the contrast vector mean. I get:

    A0     A1     B0     B1     C0     C2     D0     D1     D2
-0.250  0.250 -0.250  0.250 -0.250  0.250 -0.250  0.125  0.125

... which can also be generated from this:

con <- makeContrasts((A1 + B1 + C2 + (D1+D2)/2)/4 - (A0 + B0 + C0 + D0)/4, 
    levels=design)

In other words, this contrast compares the average of all treatment samples to the average of all control samples. I've averaged the two different treatment types in cell line D so that each cell line contributes the same number of coefficients to the contrast.