# load  packages
library(ropls)
library(ggplot2)
library(factoextra)
library(ggsci)
library(Cairo)
library(tidyverse)
library(extrafont)
loadfonts()

#load data
metabolites_raw <- read.table(file="metabolites.txt",sep="\t",header=T,check.names=FALSE ,row.names=1)
group <- read.table(file="group.txt",sep="\t",header=T,check.names=FALSE ,row.names=1)
metabolites <- t(metabolites_raw)


#PCA analysis
groupdata <- group[!duplicated(as.character(group$sample)),]
sample <- as.character(unique(groupdata$sample))
group <- as.character(unique(groupdata$group))
groupdata$group <- factor(groupdata$group,levels=group)
data <- subset(metabolites,select=sample)
data <- data[rowSums(data)>1,]
data <- na.omit(data)
data <- log2(data+1)
data <- scale(data)
data <- t(scale(t(data),center=TRUE,scale=F))
pca <- prcomp(t(data),center=FALSE,scale.=FALSE)
pca_mat <- data.frame(pc1=pca$x[,1],pc2=pca$x[,2],sample=groupdata$sample,group=groupdata$group)
p <- ggplot(pca_mat,aes(x=pc1,y=pc2,label=sample,colour=group)) +
  geom_point(size=3) +
  geom_text(colour='black',size=3) +
  labs(x=paste('PC1 ',pc1,sep=''),y=paste('PC2 ',pc2,sep='')) +
  geom_hline(yintercept=0,linetype='dotdash',size=0.8,color='grey') +
  geom_vline(xintercept=0,linetype='dotdash',size=0.8,color='grey') +
  theme_bw() +
  theme(plot.title=element_text(hjust=0.5))
p + stat_ellipse(aes(fill = group), geom = 'polygon', level = 0.95, alpha = 0.1, show.legend = FALSE)




# OPLS-DA analysis
oplsda = opls(df1_oplsda <- opls(metabolites, group$group, predI = 1, orthoI = NA))

# sample scores plot              
data <- as.data.frame(df1_oplsda@scoreMN)
o1 <- df1_oplsda@orthoScoreMN[,1]
data$o1 <- o1
data$group = group$group
data$samples = rownames(data)              
x_lab <- df1_oplsda@modelDF[1, "R2X"] * 100              
col=c("#1597A5","#FFC24B")
p1 <- ggplot(data,aes(x=p1,y=o1,color=group))+#指定数据、X轴、Y轴，颜色
  theme_bw()+#主题设置
  geom_point(size=3)+#绘制点图并设定大小
  theme(panel.grid = element_blank())+
  geom_vline(xintercept = 0,lty="dashed",color="red")+
  geom_hline(yintercept = 0,lty="dashed",color="red")+#图中虚线
  # guides(color=guide_legend(title=NULL))+#去除图例标题
  labs(x=paste0("P1 (",x_lab,"%)"),
       y=paste0("to1"))+#将x、y轴标题改为贡献度
  stat_ellipse(data=data,
               geom = "polygon",level = 0.95,
               linetype = 2,size=0.5,
               aes(fill=group),
               alpha=0.2,
               show.legend = T)+
  scale_color_manual(values = col) +#点的颜色设置
  scale_fill_manual(values = c("#1597A5","#FFC24B"))+
  theme(axis.title.x=element_text(size=12),#修改X轴标题文本
        axis.title.y=element_text(size=12,angle=90),#修改y轴标题文本
        axis.text.y=element_text(size=10),#修改x轴刻度标签文本
        axis.text.x=element_text(size=10),#修改y轴刻度标签文本
        panel.grid=element_blank())#隐藏网格线              

ggsave(p1, filename = 'figures/opls.pdf', 
       width = 5, height = 5, device = cairo_pdf)          


# VIP scores plot
data_VIP <- df1_oplsda@vipVn
data_VIP_select <- data_VIP[data_VIP > 1]           
data_VIP_select <- cbind(otu_raw[names(data_VIP_select), ], data_VIP_select)
names(data_VIP_select)[13] <- "VIP"              
data_VIP_select <- data_VIP_select[order(data_VIP_select$VIP, decreasing = TRUE), ]
data_VIP_select$G = rownames(data_VIP_select)              
p2 <- ggplot(data_VIP_select) +
  geom_segment(aes(x=G, xend=G, y=0, yend=VIP), color="grey",size=1.5,alpha=0.8) +
  geom_point( aes(x=G, y=VIP), size=4,color='red',alpha=0.8) +
  theme_prism(palette = "pearl",
              base_fontface = "plain", 
              base_family = "serif", 
              base_line_size = 0.8,
              axis_text_angle = 45) +
  theme(legend.position = "none") +
  xlab(NULL) +
  ylab(NULL) +
  scale_y_continuous(expand = c(0,0),limits = c(0,2.3))+
  ggtitle("Lollipop Chart",p3<-ggplot(data_VIP_select,aes(x = G, y = VIP),size=2) +
            geom_bar(aes(fill=G),stat = 'identity', width = 0.8)+
            labs(x = NULL, y = "VIP", title = NULL)+
            theme_prism(palette = "candy_bright",
                        base_fontface = "plain", # 字体样式，可选 bold, plain, italic
                        base_family = "serif", # 字体格式，可选 serif, sans, mono, Arial等
                        base_line_size = 0.8, # 坐标轴的粗细
                        axis_text_angle = 45) +# 可选值有 0，45，90，270
            theme(legend.position = "none"))              

p3 <- ggplot(data_VIP_select,aes(x = G, y = VIP),size=2) +
  geom_bar(aes(fill=G),stat = 'identity', width = 0.8)+
  labs(x = NULL, y = "VIP", title = NULL)+
  theme_prism(palette = "candy_bright",
              base_fontface = "plain", # 字体样式，可选 bold, plain, italic
              base_family = "serif", # 字体格式，可选 serif, sans, mono, Arial等
              base_line_size = 0.8, # 坐标轴的粗细
              axis_text_angle = 45) +# 可选值有 0，45，90，270
  theme(legend.position = "none")              

cowplot::plot_grid(p2,p3,ncol = 1)