此篇简书是线性相关分析（秩相关），可视化是热图；
还有其他相关分析方法（如基于峰值模型和降维手段的CCA分析）及对应的图，应该适合不同生物类群
install.packages(pheatmap)
library("pheatmap")
install.packages("psych")
library("psych")
species_table <-read.table("otu_table.Genus.relative.txt",comment.char = "",check.names=F,stringsAsFactors=F,header=T,sep="\t")
#comment.char设置注释字符，以这个字符开头的行将被忽略；设置为空时，就不会忽略数据中的信息
#check.names=F,不会自动修改行名和列名；
#stringsASFactor=F；不会自动把strings设置成因子
#header=true设置第一个有效行为列名
#sep="\t",设置制表符为分隔符；"\"表示制表符；
species_table$Taxonomy[duplicated(species_table$Taxonomy)]<-paste(species_table$Taxonomy[duplicated(species_table$Taxonomy)],"_1",sep="")
#此命令用来区分重复的属名，由于数据库不同科下的属会出现同名不同类的情况，因此需要加以区分
#"[]"表示索引；索引true项
#duplicated函数表示判断是否重复；是重复的‘TRUE’；非重复‘FALSE’
species_table<-species_table[-nrow(species_table),]
#表示去除最后一行；
#表格索引‘[a,b]’,a代表表格第a行，b代表表格第b列；
#‘-’表示去除这一行或这一列；
#nrow表示计算这个表格总行数；
#ncol表示计算这个表格的总列数；
rownames(species_table)<-species_table$Taxonomy
#将第一列作为行名；
species_table<-species_table[,-c(1,ncol(species_table))]
#去掉第一列和最后一列，不需要的分类信息
map<-read.table("mapping_file.txt",sep = "\t",na.strings = "",header = T,row.names = 1,comment.char = "",check.names = F,stringsAsFactors = F)
#含有准备好含有环境变量因子的map文件
#na.strings=""设置缺失值的字符
map<-map[,10:14]
#提取需要分析的环境因子
species_table<-t(species_table)
#转置
sum_of_species<-colSums(species_table)
#计算每个样品中每个物种的总丰度，用于筛选
species_table<-species_table[match(rownames(map),rownames(species_table)),]
#重排species_table行的顺序，使其与map保持一致；
merged_table<-data.frame(map,species_table,check.names=F,check.rows = T )
#合并环境因子map表格和otu表格
correlation_results<-corr.test(merged_table,method = "spearman",adjust = "fdr")
#计算spearman秩相关；spearman也可以换成pearson
correlation_results$ci
#'$'在表格后面表示返回某一列值;
#查看ci置信区间，如果置信区间太大（上限和下限），则相关系数不太可信；
r<-correlation_results$r
#提取相关矩阵
p<-correlation_results$p
#提取p值矩阵
r<-r[-c(1:5),-c(6:192)]
p<-p[-c(1:5),-c(6:192)]
#剔除不关心的相关系数，环境与环境之间的相关系数；生物与生物之间的相关系数；
selected_position_of_species<-head(order(colSums(t(p)<0.05),decreasing = T),20)
#colSums(t(p)<0.05)，每列所得出的逻辑值的和，也可换成rowSums(p<0.05)
#order对值进行排序（升序或降序），从大到小或从小到大得出值的位置，并返回值所在位置；方便索引
#sort也是排序，是返回从大到小或从小到大的值
#head是取前多少的值（20为例）
rowSums(p<0.05)
#每行所得出的逻辑值相加
r<-r[selected_position_of_species,]
p<-p[selected_position_of_species,]
#得到筛选后的相关系数矩阵和p值矩阵；
#“[selected_position_of_species,]”索引排名前20的物种的行（行名是物种名）的位置
sig_label<-function(x){ifelse(x<0.001,"***",ifelse(x<0.01,"**",ifelse(x<0.05,"*","")))}
#自定义显著性标记函数
#ifelse函数；（test，a，b），test为真返回a，为假返回b
sig_matrix<-sig_label(p)
pheatmap(r,fontsize = 15,border_color="black",
display_numbers = sig_matrix,fontsize_row = 15,fontsize_col = 15,
fontsize_number = 22,
#显著性标记的大小
cluster_rows=T,clustering_distance_rows = "correlation",
#指明行聚类、聚类依据的距离；
cluster_cols=T,clustering_distance_cols = "euclidean",
#指明列聚类、聚类依据的距离；
clustering_method = "centroid")
#指明聚类方法
#可视化
