Hi I have a quick question about a paired samples and batch effect correction. Especially I am seeking a solution for edgeR or limma/voom.

A few days ago, I fount the following post.

correction for batch effects and paired-design of RNA-seq samples

In his case, subject and batch are identical, so "There is no need to do any batch correction. The baseline differences between subjects, including any batch effect, is already accounted for in the paired analysis."

group <- factor(c("wt","wt","wt","ko","ko","ko"))
subject <- as.factor(c(1,2,3,1,2,3))
batch <- as.factor(c(1,2,3,1,2,3))

But how about the following example? In this case, a given subject's wt/ko samples are from different batches. 

group <- factor(c("wt","wt","wt","ko","ko","ko"))
subject <- as.factor(c(1,2,3,1,2,3))
batch <- as.factor(c(1,1,1,2,2,1))

What is the best model design to correct batch effect and to predict DEG between WT and KO groups? Could you suggest the model fomula? 

This design is quite poor, it's one sample away from being broken (batch is almost confounded with group). As it is, you can squeeze out 1 residual degree of freedom:

design <- model.matrix(~group + subject + batch)

... which should allow you to compare between WT and KO. This will not give you a lot of power - the dispersion will not be stably estimated, and a lot of information will be spent estimating the subject/batch blocking factors rather than being used to compare between genotypes.

Indeed, a careful analysis of the design indicates that the third WT sample is effectively useless; it will not contribute to dispersion estimation or DE testing, and will only be used by edgeR when computing the average. The two samples in batch 2 will mostly be uninvolved in DE testing, meaning that most of your analysis will rely on a 2 WT vs 1 KO comparison - not great.