I prefer to perform PCA outside of DESeq2 (but inspired by the code): Assume we start from a count matrix counts that should be on log scale. That could be your matrix of batch-corrected counts. I did not read the entire question because it was too extensive so tell me if I got you wrong.

# Find most variable genes, say top 500:
rv <- rowVars(as.matrix(counts))
select <- order(rv, decreasing = TRUE)[seq_len(min(500, length(rv)))]

# perform PCA on these genes:
pca <- prcomp(t(counts[select, ]))

# calculate the explained variation per principal component:
percentVar <- pca$sdev^2/sum(pca$sdev^2)

# feed into ggplot:
d <- data.frame(PC1 = pca$x[, 1], PC2 = pca$x[, 2], stringsAsFactors = FALSE)

# Most simplistic plot:
ggplot(data=d, aes(x = PC1, y = PC2)) + geom_point()

You can add legends etc. to indicate which data point is which sample etc.