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)
#将相关系数矩阵转换为宽格式，行名为环境变量，列名为物种，值为相关系数
rvalue <- df %>%
select(1, 2, 3) %>%
pivot_wider(names_from = "bacteria", values_from = r) %>%
column_to_rownames(var = "virus")
# 将显著性符号矩阵转换为宽格式，行名为环境变量，列名为物种，值为显著性符号
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)
