Hi Aaron :)

Thank you so much for your response and help. Could I please ask another question? I have noticed that in my top differentially expressed genes (long non-coding RNA genes), that the batch effect adjustment using the removeBatcheffect function (edgeR) significantly alters the expression values for some samples. This was more pronounced for samples in batch 4 where there is an uneven sample distribution.

For example, in this boxplot, there are 3 samples in group 2 where the expression increases significantly when removeBatcheffect was used however not with ComBat. This was the case for most of the top differentially expressed genes found by including the batch factor in the design matrix with edgeR. 

Could you please shed some light on this matter? Thank you! 

These are the codes that were used:

# s2c is my sample to conditions table
# Generating the batch effect adjusted matrix using removeBatchEffect
status <- s2c$condition2 %>% factor
design <- model.matrix(~0 + status)
colnames(design) <- gsub("status", "", colnames(design))
TPMrbe <- removeBatchEffect(TPMmatrix , batch = s2c$batch, design = design)

# Generating the batch effect adjusted matrix using ComBat
batch = factor(s2c$batch)
modcombat = model.matrix(~1, data = s2c)
modPDL1 = model.matrix(~condition1, data = s2c)
TPMcombat <- ComBat(TPMmatrix, batch = batch , mod = modcombat, par.prior=TRUE, prior.plots=FALSE)

# Performing differential expression analysis using edgeR.
status <- factor(s2c$condition2)
batch <- factor(s2c$batch)
design <- model.matrix(~0 + status + batch)
colnames(design) <- gsub("status", "", colnames(design))
xE <- DGEList(counts = countsMatrix, group = status)

xE <- calcNormFactors(xE, method = "TMM")
xE <- estimateDisp(xE, design, robust = TRUE) 
fit <- glmQLFit(xE, design, robust = TRUE)
contr.matrix <- makeContrasts( group = group2 - group1, levels = design)
resE <- glmTreat(fit, contrast = contr.matrix, lfc=1)
ttE <- topTags(resE , n=Inf, p.value=0.05, adjust.method="BH")