(I asked this as a side question before, but I though it deserves it's own post)

I have a bulk-RNA experiment with 3 tissue types: Lesional, Non-Lesional and Controls, whereas Lesional and Non-Lesional are taken from the same patients. Most of the Controls were processed in batch#2. To batch-correct, two technical replicates of lesional samples and two technical replicates of non-lesional samples were also included in batch #2. All samples have Patient ID and Sample ID, like so:

.

I've read some comments regarding batch correction in the forum, and a similar problem was raised here with a more thorough explanation here, but I couldn't find an answer for a case of batch correction followed by repeated measures. I can use sva and ComBat_seq (or RUVSeq), but the comments recommend utilizing limma's duplicateCorrelation. This is a valid option, but I don't know which blocking parameter to use - because of the repeated-measures nature of the experiment, I need to use PatientID, but this might not suit the batch correction, which should match SampleID.

My current approach is something like this, which I'm not sure is the "right" way.

exprs.mat.batchcorrected = ComBat_seq(counts=exprs.mat, batch=phenodata$batch, group = phenodata$SampleID, full_mod = T)
DGE.cpm<-DGEList(counts = exprs.mat.batchcorrected, samples = phenodata)
...  # filtering genes etc ... #
design<-model.matrix(~0 + Sample.Type + batch, data=DGE.cpm$samples)
voom.result <-voomWithQualityWeights(DGE.cpm, design, plot=T, normalize="quantile")
corfit.rna<-duplicateCorrelation(voom.result, design, block=DGE.cpm$samples$PatientID)
voom.result.randeffect <-voomWithQualityWeights(DGE.cpm, design, plot=TRUE, block=DGE.cpm$samples$PatientID, correlation=corfit.rna$consensus, normalize="quantile")
corfit.rna2<-duplicateCorrelation(voom.result.randeffect, design, block=DGE.cpm$samples$PatientID)
fit2.rna<-lmFit(voom.result.randeffect, design, block=voom.result.randeffect$targets$PatientID, correlation=corfit.rna2$consensus)

Thank you so much for helping!

You should just do a standard limma analysis as recommended in the previous posts.

design <- model.matrix(~0 + Sample.Type + batch)
i <- filterByExpr(exprs.mat, group=Sample.Type)
fit <- edgeR::voomLmFit(exprs.mat[i,], design, block=PatientID, sample.weights=TRUE)

Here I am assuming that exprs.mat is your count matrix (just counts, not cpm or anything else).

That's all you need. You should not use ComBat to pre batch-correct.

A caveat is that this analysis treats the four technical replicates as if they were biological replicates, but that can't be helped, and it will be relatively harmless given the number of samples you have. I don't know any practical way to do better.