### Data input
ExM <- as.matrix(log2(tpm + 1))
s1 <- colnames(tpm[,xx:xx]) ### Ctrl groups
s2 <- colnames(tpm[,xx:xx]) ### Treated groups

### Perfom wilcoxon test
cat("wilcox.test\n")
pvalue = log2FoldChange = matrix(0, nrow(ExM), 1)
for(i in 1:nrow(ExM)){
    pvalue[i, 1] = p.value = wilcox.test(ExM[i, s1], ExM[i, s2])$p.value
    log2FoldChange[i, 1] = mean(ExM[i, s1]) - mean(ExM[i, s2])
}
padj = p.adjust(as.vector(pvalue), "fdr", n = length(pvalue))
rTable = data.frame(log2FoldChange, pvalue, padj, row.names = rownames(ExM))
rTable$change <- ifelse(rTable$pvalue < 0.05 & abs(rTable$log2FoldChange) >= 1, 
    ifelse(rTable$log2FoldChange >1, 'Up', 'Down'), 'Stable')

### Draw DEGs
ggplot(rTable, aes(log2FoldChange, -log(pvalue), colour = change)) + 
    geom_point(alpha=0.8, size = 3) + 
    scale_color_manual(values=c("#546de5", "#d2dae2","#ff4757")) + 
    xlab(expression("log"[2]*"FC")) + ylab(expression("-log"[10]*"FDR")) + 
    geom_vline(xintercept=c(-1,1), lty=4, col="black", lwd=0.8) + 
    geom_hline(yintercept = -log10(0.05), lty=4, col="black", lwd=0.8)  +
    theme(legend.position = "none")
table(rTable$change)