#----------------------------------------------PLS-DA分析-----------------------------------------------
library(openxlsx)
#加载R包
library(ropls) # PCA, PLS-DA, OPLS-DA 
data1 = read.xlsx("D:/met.xlsx",1,rowNames = TRUE)  ## Orthogroups.GeneCount.tsv文件
group = read.xlsx("D:/met.xlsx",2,rowNames = TRUE)#读取分组信息
##opls-da分析——opls-da仅适用于二元分类，多元分类适用于pls-da分析
#orthoI:正交分量的数量（仅适用于OPLS）；当设置为0[默认]时，将执行PLS；否则将形成OPLS；当设置为NA时，执行OPLS，并通过使用交叉验证自动计算正交分量的数量（最多9个正交分量）。
#predI:要提取的成分数量（PLS和OPLS情况下的预测成分）；对于OPLS，predI（自动）设置为1；
data1 <- data1[,-15] ##去掉total列
plsda = opls(t(data1), group$group, predI = 2,orthoI=0)
plsda@scoreMN
##提取 VIP 分数，识别重要变量:
vip_scores <- getVipVn(plsda)

# 选择 VIP 分数较高的变量
important_vars <- names(vip_scores[vip_scores > 2]) ##通常以VIP值>1作为筛选标准。阈值选择可灵活变换
#整理绘图数据
plot = data.frame(Sample=rownames(plsda@scoreMN),
                  group=group$group,
                  PC1=plsda@scoreMN[,1],#读取第一列
                  PC2=plsda@scoreMN[,2])#读取第二列
#设置x，y轴标签
xlab = paste0("PC1(",round(plsda@modelDF[1,1]*100,2),"%)")#round()保留小数点几位
ylab = paste0("PC2(",round(plsda@modelDF[2,1]*100,2),"%)")

library(ggplot2)
ggplot(data=plot,#绘图数据
       aes(x=PC1,y=PC2,color=group))+#X，Y轴，以及颜色按照分组分类
  geom_point(aes(shape=group),size=3)+#绘制点，点的形状按照分组分类，大小为3
  geom_vline(xintercept=0,colour="black",linetype="dashed")+#图中横着的虚线
  geom_hline(yintercept=0,colour="black",linetype="dashed")+#图中竖着的虚线
  labs(title="PLS-DA Plot",x=xlab,y=ylab)+#标签
  stat_ellipse(aes(fill=group),alpha=0.2,color=NA,geom="polygon")+#添加置信椭圆
  theme_bw()+
  #主题设定
  theme(panel.grid=element_blank(),#背景空白
        axis.title=element_text(face="bold",colour="black",size=10),#坐标轴标题
        axis.text=element_text(face="bold",colour="black",size=8),#坐标轴文本
        plot.title=element_text(face="bold",colour="black",size=12,hjust=0.5),#图标题设定
        legend.text=element_text(face = "bold",size = 10,color="black"),#图注文本
        legend.title=element_blank())#图注标题

#----------------------------------------------Wilcoxon秩和检验-----------------------------------------------
library(rstatix)
library(tidyverse)
data1 = read.xlsx("D:/as.xlsx",1,rowNames = TRUE) ##计数数据，可以是表达值，也可以是orthofinder结果中的Orthogroups.GeneCount.tsv

res <- list()
##15753为data1的行数
for( i in 1:15753){
  
  fit <- wilcox.test(as.numeric(data1[i,as.numeric(row.names(group[group$group == "LAS", ]))]),as.numeric(data1[i,as.numeric(row.names(group[group$group == "NAS", ]))])) ##1:6列为对照组，7：14列为实验组
  res[[i]] <- data.frame(pvalue = fit$p.value,stats = fit$statistic)
}
#这里用的双尾检验，如果是符合单尾检验的假设，请修改wilcox.test()的参数，具体查看wilcox.test()的说明）

res_df <- do.call(rbind, res)
p.adjust<- p.adjust(res_df[,1],"BH")#对p值进行FDR矫正
res_df$p.adjust=p.adjust

top20 <- read.table("clipboard",header = F,check.names = F,fill=TRUE)

as <- data1[,c( <- ,nas)]
subset(rowSums(as[, nas]) == 0 & rowSums(as[, las]) >= 7)