If you're looking for the genes that are dimorphic based on sex, you've just found them. If you follow through with your "standard" limma dge pipeline, you'll get the genes differentially expressed across the sexes of your samples, ie.
design <- model.matrix(~Sex, data = pData(eset))
fit <- lmFit(eset, design)
fit <- eBayes(fit)
res <- topTable(fit, n = Inf)
The p-values and logFCs reported in the
res data.frame are those for male vs female, the order of which is determined by the
levels(eset$Sex). Just look at the 2nd column name of your
design matrix. Positive logFC's will be those that are more highly expressed in that group vs the intercept (the other sex).
Still, if you really want the residuals matrix, you could do:
resids <- limma::removeBatchEffect(eset, batch = eset$Sex)
Caveat emptor: this should work "just fine" assuming the data in
exprs(eset) is appropriate (ie. on log2 scale and more-or-less normally distributed)