hi ! I have a two factor RNA seq ,each factors have multiple levels(4 developmental stages and 16 different tissues), here is a part of my data,

There are 10 of these samples having replicates.(total 58 combinations,not 64 combinations since there are some tissues are not developed in certain developmental stage) 
Here are my assumptions: 

1. There are common DE geneto guide development regardless tissues in order to response same developmental signals, 
2. Different tissues have different DE genes to guide differentiations even if these different tissues are at the same developmental stage.

I would like to ask How to design a formalu(   moderl.matrix(~time+tissue?) ) to answer my previous questions?:
1. To find differential expressed genes regardless tissues among four periods?? 
In my thoughts, I should perform DE analysis in each single tissue among 4 periods, afterwards I just pick the co-occurrence genes in top (maybe 50?) DEed genes ? what would be the difference between code below:
d1 <- model.matrix(~time+tissue)
fit <- glmFit(y, d1)
lrt <- glmLRT(fit, coef=2:4)

I am new to edgeR and biostatics, thanks in advance.

The simplest approach would be to define your design matrix with a one-way layout, as below:

grouping <- factor(paste0(time, ".", tissue))
design <- model.matrix(~0 + grouping)
colnames(design) <- levels(grouping)

This is more appropriate than the additive model (i.e., model.matrix(~time+tissue)) if you're interested finding tissue-specific DE genes between timepoints. More specifically, if you use your additive model, you'll be making some assumptions about the consistency of the time effect across different tissues. This would defeat the purpose of looking for time-induced changes that are different between tissues. Now, if y is your DGEList object, you can then run:

y <- estimateDisp(y, design)
fit <- glmFit(y, design)

The contrasts will depend on what you're interested in. If you want to find DE genes that change between time points in a similar manner across tissues, you can do something like:

con.brain <- makeContrasts(E13.brain - P00.brain, levels=design)
lrt.brain <- glmLRT(fit, contrast=con.brain)
con.heart <- makeContrasts(E13.heart - P00.heart, levels=design)
lrt.heart <- glmLRT(fit, contrast=con.heart)

Then, you intersect the set of DE genes in lrt.heart and lrt.brain to get common changes in both tissue types (note that you'll need to manually check that the sign of the log-fold change is the same in both comparisons). Obviously, you can include more contrasts in the intersection, e.g., for kidney or lung, to get common changes across more tissue types. However, be aware that this'll be more conservative, so you'll get fewer DE genes in the final set.

Your other aim refers to tissue-specific DE genes across developmental stages. This can be accommodated by running:

con.brain2heart <- makeContrasts(E13.brain - P00.brain - (E13.heart - P00.heart), levels=design)
lrt.heart <- glmLRT(fit, contrast=con.brain2heart)

This'll find the set of genes where the DE between timepoints in brain is different to that in the heart. If you want to include more tissues, then you can do an ANOVA-style analysis:

con.tissue <- makeContrasts(E13.brain - P00.brain - (E13.heart - P00.heart),
    E13.brain - P00.brain - (E13.lung - P00.lung),
    E13.brain - P00.brain - (E13.kidney - P00.kidney),
    E13.brain - P00.brain - (E13.liver - P00.liver),
    levels=design)
lrt.tissue <- glmLRT(fit, contrast=con.tissue)

This'll find genes that respond differently between these two timepoints in any tissue.

Thanks, Aaron. Here is part of my definition of "time" and "tissue" vectors, is it like this? In each "time", There are at least two replicates of different time~tissue combinations.

Thanks, Aaron! 

"This is more appropriate than the additive model (i.e., model.matrix(~time+tissue)) if you're interested finding tissue-specific DE genes between timepoints. "

In my opinion, "between timepoints" here means a comparison between only two time points, 

eg. con.brain <- makeContrasts(E13.brain - P00.brain, levels=design)

so, how can I find tissue-specific DE genes among 4 time points? Will this kind of genes be important since they always respond to devlopmental signal ?