### 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)
网友评论