rm(list=ls())

library(pacman)

pacman::p_load(tidyverse,psych,reshape,ggtree,aplot)

#对genus进行读取，转置，转为数据框

table1 <- read.delim("genus.xls",header =T,sep="\t",

           row.names = 1,check.names = F) %>% 

  t() %>% as.data.frame()

#对环境因子进行读取

table2 <- read.delim("env.xls",header =T,sep="\t",

                     row.names = 1,check.names = F)

#对2个数据对象，相同的行名，进行计算相关性及p数值与校正。

pp <- corr.test(table1,table2,method="pearson",adjust = "fdr")

#索取r数值与p数值

cor <- pp$r 

pvalue <- pp$p

#定义计算显著性表示的函数

myfun <- function(pval) {

  stars = ""

  if(pval <= 0.001)

    stars = "***"

  if(pval > 0.001 & pval <= 0.01)

    stars = "**"

  if(pval > 0.01 & pval <= 0.05)

    stars = "*"

  if(pval > 0.05 & pval <= 0.1)

    stars = ""

  stars

}

#melt函数将宽表转为长表，rename设置列名，mutate函数选取整列，计算，添加列。

heatmap <- melt(cor) %>% rename(replace=c("X1"="sample","X2"="gene",

                                          "value"="cor")) %>%

  mutate(pvalue=melt(pvalue)[,3]) %>%

  mutate(signif = sapply(pvalue, function(x) myfun(x)))



write.table (heatmap,file ="heatmap.xls", sep ="\t", row.names = F)

#对cor的距离矩阵进行聚类

phr <- hclust(dist(cor)) %>% 

  ggtree(layout="rectangular", branch.length="none")

#对列cor的距离矩阵进行聚类

phc <- hclust(dist(t(cor))) %>% ggtree() + layout_dendrogram()

#对cor的数值进行画热图

pp <- ggplot(heatmap,aes(gene,sample,fill=cor)) + 

  geom_tile()+theme_minimal()+

  scale_fill_viridis_c()+

  geom_text(aes(label=signif),size=5,color="white",hjust=0.5,vjust=0.5)+

  scale_y_discrete(position="right")+xlab(NULL) + ylab(NULL)+

  theme(axis.text.x=element_text(angle =90,hjust=1,vjust=0.5,

                                 family = "Times",face = "italic",colour = "black",size=12),

        axis.text.y=element_text(family= "Times",

                                 face = "plain",colour = "black",size=12),

        legend.text=element_text(face="plain",family = "Times",

                                 colour = "black",size = 12))+

  guides(fill = guide_colorbar(direction = "vertical",reverse = F, 

                               barwidth = unit(.5, "cm"),

                               barheight = unit(8, "cm")))+labs(fill= "")

#热图与聚类图结合在一起

pp %>% insert_left(phr, width=.2) %>%

  insert_top(phc, height=.1)