rm(list = ls()) ## 魔幻操作，一鍵清空~

options(stringsAsFactors = F)

#######進行細菌和病毒相關性分析

###

library(ggplot2)

library(ggpubr)

library(reshape2)

library(vegan)

library(ggrepel)

library(aplot)

library(gridExtra)

library(readxl)

library(corrplot)

library(psych)

library(tidyverse)

library(psych)

library(pheatmap)

library(magrittr)

library(scico)

############行為樣本

bacteria <- read.csv("細菌.csv",sep = ",",header = TRUE)

row.names(bacteria) <- bacteria$sample

bacteria <- bacteria[,-1]

virus <- read.csv("病毒.csv",sep = ",",header = TRUE)

row.names(virus) <- virus$sample

virus <- virus[,-1]

bacteria <- apply(bacteria, 2, as.numeric)

virus <- apply(virus, 2, as.numeric)

#############用spearman/person

# corr <- cor (bacteria, virus,method="pearson")

pp <- corr.test(bacteria, virus, method = "spearman", adjust = "fdr")

cor <- pp$r # 獲取相關系數矩陣

pvalue <- pp$p # 獲取p-value矩陣

df <- melt(cor) %>%

? mutate(pvalue = melt(pvalue)[, 3],

? ? ? ? p_signif = symnum(pvalue, corr = FALSE, na = FALSE,?

? ? ? ? ? ? ? ? ? ? ? ? ? cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),

? ? ? ? ? ? ? ? ? ? ? ? ? symbols = c("***", "**", "*", "", " "))) %>%

? set_colnames(c("bacteria", "virus", "r", "p", "p_signif"))

####存儲結果

write.csv(df,"./bacteria_virus_corrplot.csv",row.names = FALSE)

#將相關系數矩陣轉換為寬格式，行名為環(huán)境變量，列名為物種，值為相關系數

rvalue <- df %>%

? select(1, 2, 3) %>%

? pivot_wider(names_from = "bacteria", values_from = r) %>%

? column_to_rownames(var = "virus")

# 將顯著性符號矩陣轉換為寬格式，行名為環(huán)境變量，列名為物種，值為顯著性符號

pvalue <- df %>%

? select(1, 2, 5) %>%

? pivot_wider(names_from = "bacteria", values_from = p_signif) %>%

? column_to_rownames(var = "virus")

mycol <- scico(100, palette = "vik")

p <- pheatmap(rvalue, scale = "none", cluster_row = TRUE, cluster_col = TRUE, border = NA,

? ? ? ? ? ? ? display_numbers = pvalue, fontsize_number = 12, number_color = "white",

? ? ? ? ? ? ? main = " ",

? ? ? ? ? ? ? cellwidth = 21, cellheight = 20, color = mycol)