? ? ? ?本篇文章給大家?guī)頍釄D的實操，大家一定要動手實踐起來哦！學(xué)會了的知識才是自己的！趕緊跟著小易學(xué)習(xí)起來吧！以微生物組與代謝組聯(lián)合分析作為示例，帶大家復(fù)現(xiàn)幾種相關(guān)性熱圖！
1. 準(zhǔn)備輸入數(shù)據(jù)
1) 物種豐度表（*_otu.tsv），格式如下：
說明：第一列為物種名稱，第一行為樣本的名稱，其他的為對應(yīng)樣本的物種豐度，表格用制表格分割。
2）代謝物含量表（metabolites_*.tsv），格式如下：
說明：第一列為代謝物id（或者名稱），第一行為樣本的名稱，其他的為對應(yīng)樣本的代謝物的含量，表格用制表格分割。
2. 腳本及運行
2.1 相關(guān)性聚類熱圖
用法：
Rscript combine_heatmap.R genus_otu.tsv metabolites_differential.tsv prefix
可視化示例圖：
library(psych)
library(vegan)
library(reshape2)
library(ComplexHeatmap)
library(circlize)
options(bitmapType='cairo') #關(guān)閉服務(wù)器與界面的互動響應(yīng)
read_data <- function(data){
? ? data <- read.table(data, header=T, row.names=1, sep="\t", comment.char="", quote="", check.names=FALSE)
? ? #check.names=FALSE 保留原始表頭
? ? data <- data[which(rowSums(data) >0),] #過濾所有豐度都為0得行 ? ?
? ? data <- t(data)
? ? return(data)
}
filter_pvalue_row <- function(data, pmax=0.05){#過濾P值均大于0.05的行
? ? r <- c()
? ? for(i in rownames(data)){
? ? ? ? ? ? if(min(data[i, ]) <= pmax){ #存在P值小于等于0.05的行
? ? ? ? ? ? ? ? ? r <- c(r, i)
? ? ? ? ? ? ? ? ? }else{ ?
? ? ? ? ? ? ? ? ? ? ? ? ? print(i) ? ? ??
? ? ? ? ? ? ? ? ? } ? ?
? ? ? ? ? }
? ? ? ? ? return(r)
}
filter_rvalue_row <- function(data, rmin=0.6){ #過濾R值均小于0.6的行
?? ?r <- c()
? ? for(i in rownames(data)){
? ? ? ? if(max(abs(data[i, ])) >= rmin){ #存在R值大于等于0.6的行
?? ? ? ? ? ?r <- c(r, i)
? ? ? ? }else{?
? ? ? ? ? ?print(i)
?? ? ? ?}
?? ?}
? ? return(r)
}
mark_pvalue <- function(pvalue){#篩選標(biāo)注哪些小于0.01小于0.05并標(biāo)注
? ? if(!is.null(pvalue)){
? ? ? ? site1 <- pvalue<0.01
? ? ? ? pvalue[site1] <-"**"
? ? ? ? site2 <- pvalue >0.01& pvalue <0.05
? ? ? ? pvalue[site2] <-"*"
? ? ? ? pvalue[!site2&!site1]<-""
? ? }else{
?? ? ? ?pvalue <- F
?? ?}
? ? return(pvalue)
}
filter_pvalue <- function(pvalue, rvalue){
? ? indexs <- which(apply(pvalue,1, min, na.rm=TRUE) <0.05)
?? ?pvalue <- pvalue[indexs,]
?? ?rvalue <- rvalue[indexs,]
? ? r <-list(pvalue=pvalue, rvalue=rvalue)
? ? return(r)
}
size_picture <- function(data){
?? ?rows <- nrow(data) ? ? ?#計算行數(shù)
? ? #print(rows)
? ? columns <- length(data[1,]) ? ? #計算列數(shù)
? ? #print(columns)
? ? if(rows<=20) {
? ? ? ? widths <-20
? ? ? ? ratio <-1
? ? }elseif(rows<=50) {
? ? ? ? widths <-25
? ? ? ? ratio <-1.2
? ? }elseif(rows<=80){
? ? ? ? widths <-35
? ? ? ? ratio <-1.6
? ? ? ? }elseif(rows<=100){
? ? ? ? widths <-50
? ? ? ? ratio <-1.6
? ? }else{
? ? ? ? widths <-100
? ? ? ? ratio <-1.8
?? ?}
? ? if(columns<=8) {
? ? ? ? heights <-8
? ? }elseif(columns<=12) {
? ? ? ? heights <-10
? ? }elseif(columns<=20) {
? ? ? ? heights <-15
? ? }else{
? ? ? ? heights <-20
?? ?}
? ? r <-list(widths=widths, heights=heights, ratio=ratio)
}
correlation_analysis <- function(data1, data2, prefix){
?? ?data1 <- read_data(data1)
?? ?data2 <- read_data(data2)
?? ?sample_id <- rownames(data1)
?? ?data2 <- data2[sample_id, ] #匹配樣本
? ? result <- corr.test(data1, data2, method="spearman", adjust="none") #計算相關(guān)性 pearson
?? ?rvalue <- result$r
?? ?pvalue <- result$p
? ? rowpid <- filter_pvalue_row(pvalue,0.05) #獲取有P值小于0.05的行
?? ?pvalue <- t(pvalue[rowpid,])
? ? colpid <- filter_pvalue_row(pvalue,0.05) #獲取有P值小于0.05的列
?? ?rvalue <- rvalue[rowpid, colpid]
? ? rowrid <- filter_rvalue_row(rvalue,0.6) #獲取有R值小于0.6的行
?? ?rvalue <- t(rvalue[rowrid,])
? ? colrid <- filter_rvalue_row(rvalue,0.6) #R值小于0.6的行
?? ?rvalue <- rvalue[colrid, ]
?? ?pvalue <- pvalue[colrid, rowrid]
?? ?r <- filter_pvalue(pvalue, rvalue)
?? ?rvalue <- r$rvalue
?? ?pvalue <- r$pvalue
? ? result <- melt(rvalue, value.name="cor")
? ? result$pvalue <- as.vector(pvalue)
? ? write.table(result, file=paste(prefix,".correlation_pvalue.xls", sep=""), row.names=FALSE, sep="\t", quote=FALSE, na ="")
? ? data.1<-scale(data1)
? ? data.2<-scale(data2)
?? ?rvalue<-t(rvalue)
?? ?pvalue <- mark_pvalue(pvalue)
?? ?pvalue<-t(pvalue)
? ? data.2<- data.2[,match(colnames(rvalue), colnames(data.2))]
? ? data.1<- t(data.1[,match(rownames(rvalue), colnames(data.1))])
?? ?wh <- size_picture(t(data2))
?? ?widths <- wh$widths
?? ?heights <- wh$heights
? ? ratio <- wh$ratio
? ? col_fun = colorRamp2(c(-1,0,1), c("blue","white","red"))
?? ?ha <- rowAnnotation(empty=anno_empty(border=FALSE),
? ? ? ? Microbiome=data.1, col=list(Microbiome=col_fun),
? ? ? ? show_legend=F, annotation_label="Microbiome", foo = anno_text(rownames(data.1)))
? ? p1 <- Heatmap(rvalue, name="Correlation",?
? ? ? ? cluster_columns=T, cluster_rows=T, show_row_names=F,
? ? ? ? show_heatmap_legend=TRUE, width=unit(widths*ratio,"cm"), height=unit(heights*0.55,"cm"),
? ? ? ? column_dend_height=unit(3,"cm"),
?? ? ? ?cell_fun = function(j, i, x, y, width, height, fill) {
? ? ? ? ? ? grid.text(sprintf("%s",pvalue[i, j]), x, y, gp = gpar(fontsize =10))
?? ? ? ?}, right_annotation=ha)
? ? p2 <- Heatmap(data.2, name="Metabolome",?
? ? ? ? column_title="Metabolome",column_title_side="bottom", cluster_columns=F, show_row_names=F, #show_row_names=T 添加樣本名稱?
? ? ? ? cluster_rows=F, width=unit(widths*ratio,"cm"), height=unit(heights*0.55,"cm"))
? ? lgd.list= Legend(col_fun=col_fun, title="Microbiome")
?? ?ht = p1%v%p2
? ? pdf(paste(prefix,".combine_heatmap.pdf", sep=""), width=widths*1.2, height=heights)
? ? ptemp <- dev.cur()
? ? png(paste(prefix,".combine_heatmap.png", sep=""), width=widths*75, height=heights*60)
? ? dev.control("enable")
? ? draw(ht, annotation_legend_list=lgd.list, column_title="", ?#不顯示標(biāo)題column_title="Metabolome"
? ? ? ? column_title_side="bottom", column_title_gp=gpar(fontface='bold',fontsize=20),? ? ? ??
? ? ? ? padding=unit(c(2,2,10,2),"mm"))
? ? decorate_annotation("Microbiome", {
? ? ? ? grid.text("Microbiome", gp=gpar(fontface='bold', fontsize=20) ,
? ? ? ? y=unit(1,"npc") + unit(2,"mm"), just ="bottom")}
?? ?)
?? ?dev.copy(which=ptemp)
?? ?dev.off()
?? ?dev.off()
}
2.2?相關(guān)性層次聚類熱圖
用法：
Rscript correlation_heatmap.R?metabolites_differential.tsv?genus_otu.tsv? prefix
可視化示例圖：
library(psych)
library(pheatmap)
library(reshape2)
options(bitmapType='cairo') #關(guān)閉服務(wù)器與界面的互動響應(yīng)
run_data <- function(data){
? ? data <- read.table(data, header=T, row.names=1, sep="\t", comment.char="",quote="", check.names=FALSE)
? ? #check.names=FALSE 保留原始表頭
? ? data <- data[which(rowSums(data) > 0),] #過濾所有豐度都為0得行
? ? data <- t(data)
? ? return(data)
}
filter_pvalue_row <- function(data, pmax=0.05){
? ? #過濾P值均大于0.05的行
? ? r <- c()
? ? for (i in rownames(data)){
? ? ? ? if(min(data[i, ]) <= pmax){ #存在P值小于等于0.05的行
? ? ? ? ? ? r <- c(r, i)
? ? ? ? }else{
? ? ? ? ? ? print(i)
? ? ? ? }
? ? }
? ? return(r)
}
filter_rvalue_row <- function(data, rmin=0.6){#過濾R值均小于0.6的行
? ? r <- c()
? ? for (i in rownames(data)){
? ? ? ? if(max(abs(data[i, ])) >= rmin){ #存在R值大于等于0.6的行
? ? ? ? ? ? r <- c(r, i)
? ? ? ? }else{
? ? ? ? ? ? print(i)
? ? ? ? }
? ? }
? ? return(r)
}
mark_pvalue <- function(pvalue){
? ? if(!is.null(pvalue)){
? ? ? ? site1 <- pvalue<0.01
? ? ? ? pvalue[site1] <- "**"
? ? ? ? site2 <- pvalue > 0.01& pvalue < 0.05
? ? ? ? pvalue[site2] <- "*"
? ? ? ? pvalue[!site2&!site1]<- ""
? ? } else{
? ? ? ? pvalue <- F
? ? }
? ? return(pvalue)
}
filter_pvalue <- function(pvalue, rvalue){
? ? indexs <- which(apply(pvalue, 1, min, na.rm=TRUE) < 0.05)
? ? pvalue <- pvalue[indexs,]
? ? rvalue <- rvalue[indexs,]
? ? r <- list(pvalue=pvalue, rvalue=rvalue)
? ? return(r)
}
correlation_analysis <- function(data1, data2, prefix){
? ? data1 <- run_data(data1)
? ? data2 <- run_data(data2)
? ? sample_id <- rownames(data1)
? ? data2 <- data2[sample_id, ] #匹配樣本
? ? result <- corr.test(data1, data2, method="spearman", adjust="none") #pearson
? ? rvalue <- result$r
? ? pvalue <- result$p
? ? rowpid <- filter_pvalue_row(pvalue, 0.05) #獲取有P值小于0.05的行
? ? pvalue <- t(pvalue[rowpid,])
? ? colpid <- filter_pvalue_row(pvalue, 0.05) #獲取有P值小于0.05的列
? ? rvalue <- rvalue[rowpid, colpid]
? ? rowrid <- filter_rvalue_row(rvalue, 0.6) #獲取有R值小于0.6的行
? ? rvalue <- t(rvalue[rowrid,])
? ? colrid <- filter_rvalue_row(rvalue, 0.6) #R值小于0.6的行
? ? rvalue <- rvalue[colrid, ]
? ? pvalue <- pvalue[colrid, rowrid]
? ? r <- filter_pvalue(pvalue, rvalue)
? ? rvalue <- r$rvalue
? ? pvalue <- r$pvalue
? ? result <- melt(rvalue, value.name="cor")
? ? result$pvalue <- as.vector(pvalue)
? ? write.table(result, file=paste(prefix, ".correlation.xls", sep=""), row.names=FALSE, sep="\t", quote=FALSE, na ="")
? ? pvalue <- mark_pvalue(pvalue)
? ? mycol <- colorRampPalette(c("blue","white","tomato"))(800)
? ? pheatmap(rvalue, scale="none", cluster_row=T, cluster_col=T,
? ? ? ? ? ? border=NA, display_numbers=pvalue, fontsize_number=12,
? ? ? ? ? ? number_color = "white", cellwidth=12, cellheight=15,
? ? ? ? ? ? color=mycol, filename=paste(prefix, ".correlation_heatmap.pdf", sep=""))
? ? pheatmap(rvalue, scale="none", cluster_row=T, cluster_col=T,
? ? ? ? ? ? border=NA, display_numbers=pvalue, fontsize_number=12,
? ? ? ? ? ? number_color = "white", cellwidth=12, cellheight=15,
? ? ? ? ? ? color=mycol, filename=paste(prefix, ".correlation_heatmap.png", sep=""))
}
2.3 代謝物與微生物的相關(guān)性熱圖
用法：Rscript?ellipse_heatmap.R?metabolites_differential.tsv?genus_otu.tsv? prefix
可視化示例圖：
library(psych)
library(corrplot)
library(reshape2)
options(bitmapType='cairo') #關(guān)閉服務(wù)器與界面的互動響應(yīng)
run_data <- function(data){
? ? data <- read.table(data, sep="\t", header=T, row.names=1, check.names=FALSE)
? ? #check.names=FALSE 保留原始表頭
? ? data <- data[which(rowSums(data) > 0),] #過濾所有豐度都為0得行
? ? data <- t(data)
? ? return(data)
}
filter_pvalue_row <- function(data, pmax=0.05){#過濾P值均大于0.05的行
? ? r <- c()
? ? for (i in rownames(data)){
? ? ? ? if(min(data[i, ]) <= pmax){ #存在P值小于等于0.05的行
? ? ? ? ? ? r <- c(r, i)
? ? ? ? }else{
? ? ? ? ? ? print(i)
? ? ? ? }
? ? }
? ? return(r)
}
filter_rvalue_row <- function(data, rmin=0.6){#過濾R值均小于0.6的行
? ? r <- c()
? ? for (i in rownames(data)){
? ? ? ? if(max(abs(data[i, ])) >= rmin){ #存在R值大于等于0.6的行
? ? ? ? ? ? r <- c(r, i)
? ? ? ? }else{
? ? ? ? ? ? print(i)
? ? ? ? }
? ? }
? ? return(r)
}
max_name <- function(samples){ #獲取樣本的最長名稱長度
? ? maxlen <- 0
? ? for(i in samples){
? ? ? if(maxlen <= nchar(i)){
? ? ? ? ? maxlen <- nchar(i)
? ? ? }
? ? }
? ? return(maxlen)
}
ellipse_heatmap <- function(data1, data2, prefix){
? ? data1 <- run_data(data1)
? ? data2 <- run_data(data2)
? ? sample_id1 <- rownames(data1)? #獲取樣本名稱
? ? sample_id2 <- rownames(data2)? #獲取樣本名稱
? ? sample_id <- intersect(sample_id1, sample_id2)? #提取共有的樣本名稱
? ? print(sample_id)
? ? data1 <- data1[sample_id, ]
? ? data2 <- data2[sample_id, ]
? ? #連續(xù)數(shù)據(jù)，正態(tài)分布，線性關(guān)系，用pearson相關(guān)系數(shù)是最恰當(dāng)
? ? #上述任一條件不滿足，就用spearman相關(guān)系數(shù)，不能用pearson相關(guān)系數(shù)。
? ? result <- corr.test(data1, data2, method="spearman", adjust="none")
? ? pvalue <- result$p
? ? rowpid <- filter_pvalue_row(pvalue, 0.05) #獲取有P值小于0.05的行
? ? pvalue <- t(pvalue[rowpid,])
? ? colpid <- filter_pvalue_row(pvalue, 0.05) #獲取有P值小于0.05的列
? ? rvalue <- result$r
? ? rvalue <- rvalue[rowpid, colpid]
? ? rowrid <- filter_rvalue_row(rvalue, 0.6) #獲取有R值小于0.6的行
? ? rvalue <- t(rvalue[rowrid,])
? ? colrid <- filter_rvalue_row(rvalue, 0.6) #R值小于0.6的行
? ? rvalue <- rvalue[colrid, ]
? ? pvalue <- pvalue[colrid, rowrid]
? ? #rows <- nrow(rvalue) #獲取行數(shù)
? ? rows <-length(colrid)
? ? rowmaxlen <- max_name(colrid) #最長的行名稱
? ? #columns <- length(rvalue[1,])? #獲取列數(shù)目
? ? columns <- length(rowrid)
? ? colmaxlen <- max_name(rowrid) #最長的列名稱
? ? tlcex <- 2
? ? if(rows<=10){
? ? ? ? heights <- 6
? ? }else if(rows<=20){
? ? ? ? heights <- 14
? ? ? ? tlcex <- 1.8
? ? }else if(rows<=50){
? ? ? ? heights <- 18
? ? ? ? tlcex <- 1.6
? ? }else{
? ? ? ? heights <- 22
? ? ? ? tlcex <- 1.4
? ? }
? ? heights <- heights + rowmaxlen*0.15?
? ? if(columns<=4) {
? ? ? ? widths <- 6
? ? }else if(columns<=8) {
? ? ? ? widths <- 10
? ? }else {
? ? ? ? widths <- 10 + (columns-8)*0.2
? ? }
? ? widths <- widths + colmaxlen*0.15
? ? #col1 = colorRampPalette(colors=c("red", "white", "darkgreen"), space="Lab")
? ? col1 = colorRampPalette(c("#7F0000", "red", "#FF7F00", "yellow", "white",
? ? ? ? ? ? ? ? ? ? ? ? ? "cyan", "#007FFF", "blue", "#00007F"))
? ? pdf(paste(prefix, ".ellipse_heatmap.pdf", sep=""), width=widths, height=heights)
? ? a <- dev.cur()
? ? png(paste(prefix, ".ellipse_heatmap.png", sep=""), width=widths*60, height=heights*60)
? ? dev.control("enable")
? ? corrplot(rvalue, p.mat=pvalue, method="ellipse", tl.col="black", tl.cex=tlcex,
? ? ? ? insig="label_sig", sig.level=c(0.001, 0.01, 0.05),
? ? ? ? pch.col="black", pch.cex=0.8) #不顯著相關(guān)性系數(shù)圖塊上有X符號,? col=col1(10)
? ? #corrplot(rvalue, p.mat=pvalue, method="circle", type="lower",
? ? #? ? tl.col="black", tl.cex=tlcex, tl.srt=45,
? ? #? ? insig="label_sig", sig.level=c(0.001, 0.01, 0.05),
? ? #? ? pch.col="black", pch.cex=0.8) #不顯著相關(guān)性系數(shù)圖塊上有X符號,? col=col1(10)
? ? dev.copy(which=a)
? ? dev.off()
? ? dev.off()
}
