rcompanion包中的wilcoxonZ函数提供计算两组wilcoxon检验Z-score计算方法，sign函数可以返回值的正负，-1或1，wilcox.test则进行检验获得p value。

R文档： https://rdrr.io/cran/rcompanion/man/wilcoxonZ.html

wilcoxonZ参数：

wilcox.test参数：

exact: a logical indicating whether an exact p-value should be computed.
correct: a logical indicating whether to apply continuity correction in the normal approximation for the p-value.
两组wilcoxon比较计算z-score和p value

setwd("analysis_ding/06.FunctionPrediction/03_wilcoxon/")
df = read.table("t_input.txt", header=T, sep="\t", row.names=1)
# z-score数据标准化
# df = cbind(data.frame(Group=df$Group), scale(df[, 2:ncol(df)]))
BC = df[which(df$Group=="BC"), c(2:ncol(df))]
JZ = df[which(df$Group=="JZ"), c(2:ncol(df))]
KC = df[which(df$Group=="KC"), c(2:ncol(df))]
SAP = df[which(df$Group=="SAP"), c(2:ncol(df))]
# control
CK = df[which(df$Group=="CK"), c(2:ncol(df))]
# wilcoxon比较
# 函数，计算z-score，wilcoxon p value
library("rcompanion")
fun_zp <- function(v1, v2)
    z = as.numeric(wilcoxonZ(v1, v2, exact=FALSE, correct=FALSE))
    wil = wilcox.test(v1, v2, exact=FALSE, correct=FALSE)
    return(c(z, wil$p.value))
# analysis
length = ncol(BC)
#3.1 BC_CK
#fun_zp(BC[,1], CK[,1])
data = BC
mark = c()
zscore = c()
pvalue = c()
for(i in 1:length)
    res = fun_zp(data[,i], CK[,i])
    mark = c(mark, colnames(CK)[i])
    zscore = c(zscore, res[1])
    pvalue = c(pvalue, res[2])
res_BC = data.frame(mark=mark,zscore=zscore,pvalue=pvalue)
#3.2 JZ_CK
data = JZ
mark = c()
zscore = c()
pvalue = c()
for(i in 1:length)
    res = fun_zp(data[,i], CK[,i])
    mark = c(mark, colnames(CK)[i])
    zscore = c(zscore, res[1])
    pvalue = c(pvalue, res[2])
res_JZ = data.frame(mark=mark,zscore=zscore,pvalue=pvalue)
#3.3 KC_CK
data = KC
mark = c()
zscore = c()
pvalue = c()
for(i in 1:length)
    res = fun_zp(data[,i], CK[,i])
    mark = c(mark, colnames(CK)[i])
    zscore = c(zscore, res[1])
    pvalue = c(pvalue, res[2])
res_KC = data.frame(mark=mark,zscore=zscore,pvalue=pvalue)
#3.4 SAP_CK
data = SAP
mark = c()
zscore = c()
pvalue = c()
for(i in 1:length)
    res = fun_zp(data[,i], CK[,i])
    mark = c(mark, colnames(CK)[i])
    zscore = c(zscore, res[1])
    pvalue = c(pvalue, res[2])
res_SAP = data.frame(mark=mark,zscore=zscore,pvalue=pvalue)
## 整理结果数据
res_BC
res_JZ
res_KC
res_SAP
markdown = c()
signif = res_BC[which(res_BC$pvalue <= 0.05),]
markdown = c(markdown, signif$mark)
signif = res_BC[which(res_JZ$pvalue <= 0.05),]
markdown = c(markdown, signif$mark)
signif = res_BC[which(res_KC$pvalue <= 0.05),]
markdown = c(markdown, signif$mark)
signif = res_BC[which(res_SAP$pvalue <= 0.05),]
markdown = c(markdown, signif$mark)
# 各组比较都显著的target
out = data.frame(table(markdown))
out = out[which(out$Freq==4),]
target = as.character(out$mark)
# 提取target的zscore和pvalue绘图
res_BC
res_JZ
res_KC
res_SAP
input_BC = res_BC[which(res_BC$mark%in%target),]
input_JZ = res_JZ[which(res_JZ$mark%in%target),]
input_KC = res_KC[which(res_KC$mark%in%target),]
input_SAP = res_SAP[which(res_SAP$mark%in%target),]
input = input_BC
input = input_JZ
input = input_KC
input = input_SAP
input$sign = as.character(sign(input$zscore))
#library("ggplot2")
ggplot(input, aes(x=mark, y=zscore, fill=sign)) +
  geom_col(color=NA, width=0.3) +
  scale_fill_manual(
      values = c("-1" = "deepskyblue3", "1" = "indianred3")) +
  labs(x="", y="Z-score") +
  theme_classic() +
  theme(axis.title = element_text(size = 18),
        axis.text = element_text(size = 12),
        axis.line = element_line(size = 1),
        axis.ticks = element_line(size = 1)) +
  coord_flip() +
  scale_y_continuous(limits = c(-8, 8), 
                     expand = c(0, 0), 
                     breaks = c(-4, -2, 0, 2, 4)) +
  geom_hline(yintercept = -1.96, size=0.8, linetype="dashed") +
  geom_hline(yintercept = 1.96, size=0.8, linetype="dashed") +
  geom_hline(yintercept = 0, size=0.8)
ggsave(p, file="wilcoxon_BC.pdf", width=14)
ggsave(p, file="wilcoxon_JZ.pdf", width=14)
ggsave(p, file="wilcoxon_KC.pdf", width=14)
ggsave(p, file="wilcoxon_SAP.pdf", width=14)