Please see following code.I am wondering DESeq() automatically removes the outliers or should I include or specify the outliers function? When I did the analysis in around 100 samples the results looks fine. But in case of higher numbers of samples, the results looks different. Please advice how can I visualize the normalized data in box plots.

## QNBCvsTNBC ARP
dds <- DESeqDataSetFromMatrix(countData = cts,
                              colData = coldata,
                              design = ~ TNBCstatus)
dds <- DESeq(dds)
##results
colData(dds)[, "TNBCstatus"]
resultsNames(dds)
res <- results(dds, contrast=c("TNBCstatus", "QNBC", "TNBCARP"))
res
mcols(res, use.names = TRUE)
resOrdered <- res[order(res$pvalue),]
write.csv(as.data.frame(resOrdered), 
          file="QNBC_vs_TNBC.csv")
resSig <- subset(resOrdered, padj < 0.05)
resSig
write.csv(resSig,"QNBC_VS_TNBC_significant")
plotMA(res, ylim=c(-2,2))
vsd <- varianceStabilizingTransformation(dds,blind = TRUE)
n <- assay(vsd)
write.csv(n,"normalized data_log2values_QNBC_VS_TNBCdds")
plotMA(resSig, ylim=c(-2,2))

Please take a look at the DESeq2 paper and vignette, which both have sections about how outliers are treated.