For a dataset I'm working on, we would like to use the limma package for differential expression analysis. We are using limma::removeBatchEffect (doing it for downstream analysis) before using limma for differential expression analysis. Would it be okay to still include the batch covariate in the model to account for degrees of freedom, even though it's been already corrected for?

Thank you in advance

Using removeBatchEffect and then adding in the covariate to the model is what you would need to do in order to account for the lost degrees of freedom that were used to correct for the batch effect. But that's not really how you should do it. Instead you should just put the batch effect in your model rather than doing this two-step approach.