導讀

相關pheatmap圖，cytoscape網絡圖霹菊，一文打盡杠河。

一、模擬輸入文件

ko_abun = as.data.frame(matrix(abs(round(rnorm(200, 100, 10))), 10, 20))
colnames(ko_abun) = paste("KO", 1:20, sep="_")
rownames(ko_abun) = paste("sample", 1:10, sep="_")
ko_abun

二浇辜、組內相關分析

1 寫函數(shù)

library(psych)
library(stringr)

correlate = function(other, metabo, route)
{
    #  讀取方式：check.name=F, row.names=1, header=T
    # 計算相關性:
    #other = data
    #metabo = env
    #route="gut"
    result=data.frame(print(corr.test(other, metabo, use="pairwise", method="spearman", adjust="fdr", alpha=.05, ci=TRUE, minlength=100), short=FALSE, digits=5))
    # FDR矯正
    result_raw=data.frame(print(corr.test(other, metabo, use="pairwise", method="spearman", adjust="none", alpha=.05, ci=TRUE, minlength=100), short=FALSE, digits=5))
    # 原始P value

    # 整理結果
    pair=rownames(result)  # 行名
    result2=data.frame(pair, result[, c(2, 4)])  # 提取信息

    # P值排序
    # result3=data.frame(result2[order(result2[,"raw.p"], decreasing=F),])

    # 格式化結果【將細菌代謝物拆成兩列】
    result4=data.frame(str_split_fixed(result2$pair, "-", 2), result2[, c(2, 3)], p_value=result_raw[, 4])
    colnames(result4)=c("feature_1", "feature_2", "r_value", "fdr_p_value", "raw_p_value")

    # 保存提取的結果
    write.table(result4, file=paste(route, "Correlation_result.txt", sep="/"), sep="\t", row.names=F, quote=F)
}

2 相關分析

dir.create("Result")  # 創(chuàng)建結果目錄
correlate(ko_abun, ko_abun, "Result")

三、pheatmap熱圖

1 寫函數(shù)

library(reshape2)
library(pheatmap)

correlate_pheatmap = function(infile, route)
{
    data=read.table(paste(route, infile, sep="/"), sep="\t", header=T)

    data_r=dcast(data, feature_1 ~ feature_2, value.var="r_value")
    data_p=dcast(data, feature_1 ~ feature_2, value.var="raw_p_value")
    rownames(data_r)=data_r[,1]
    data_r=data_r[,-1]
    rownames(data_p)=data_p[,1]
    data_p=data_p[,-1]
    
    # 剔除不顯著的行
    del_row = c()
    for(i in 1:length(data_p[, 1]))
    {
        if(all(data_p[i, ] > 0.05))
        {
            del_row = c(del_row, i)
        }
    }

    # 剔除不顯著的列
    del_col = c()
    for(j in 1:length(data_p[1, ]))
    {
        if(all(data_p[, j] > 0.05))
        {
            del_col = c(del_col, j)
        }
    }
    
    # null值處理
    if(is.null(del_row) && !(is.null(del_col)))
    {
        data_p = data_p[, -del_col]
        data_r = data_r[, -del_col]
    }else if(is.null(del_col) && !(is.null(del_row)))
    {
        data_p = data_p[-del_row,]
        data_r = data_r[-del_row,]
    }else if(is.null(del_row) && is.null(del_col))
    {
        print("delete none")
    }else if(!(is.null(del_row)) && !(is.null(del_col)))
    {
        data_p = data_p[-del_row, -del_col]
        data_r = data_r[-del_row, -del_col]
    }
    
    # data_p = data_p[-del_row, -del_col]
    # data_r = data_r[-del_row, -del_col]
    write.csv(data_p, file=paste(route, "data_p.csv", sep="/"))
    write.csv(data_r, file=paste(route, "data_r.csv", sep="/"))

    # 用"*"代替<=0.05的p值唾戚，用""代替>0.05的相對豐度
    data_mark=data_p
    for(i in 1:length(data_p[,1])){
        for(j in 1:length(data_p[1,])){
            #data_mark[i,j]=ifelse(data_p[i,j] <= 0.05, "*", "")
            if(data_p[i,j] <= 0.001)
            {
                data_mark[i,j]="***"
            }
            else if(data_p[i,j] <= 0.01 && data_p[i,j] > 0.001)
            {
                data_mark[i,j]="**"
            }
            else if(data_p[i,j] <= 0.05 && data_p[i,j] > 0.01)
            {
                data_mark[i,j]="*"
            }
            else
            {
                data_mark[i,j]=""
            }
        }
    }
    write.csv(data_mark, file=paste(route, "data_mark.csv", sep="/"))

    pheatmap(data_r, display_numbers=data_mark, cellwidth=20, cellheight=20, fontsize_number=18, filename=paste(route, "Correlation_result.pdf", sep="/"))
    pheatmap(data_r, display_numbers=data_mark, cellwidth=20, cellheight=20, fontsize_number=18, filename=paste(route, "Correlation_result.png", sep="/"))
}

2 可視化

correlate_pheatmap("Correlation_result.txt", "Result")

3 結果

三個繪圖文件和兩幅圖

四柳洋、cytoscape網絡圖

1 文件準備

先刪除自我相關，也就是r值為1的相關叹坦。
然后熊镣，更新行號。
接著，建feature1-feature2的box绪囱。
遍歷所有feature2-feature1與box[-1]匹配测蹲，如果匹配成功則是重復，記錄行號鬼吵，最后刪除扣甲。

data = read.table("Result/Correlation_result.txt", sep="\t", header=T)
# 刪除自我相關
data = data[data$r_value != 1,]
# 刪除重復相關
rownames(data) = 1:nrow(data)
box = paste(data$feature_1, data$feature_2, sep="-")
delete = c()
for(i in 1:nrow(data))
{
    tmp = paste(data[i, 2], data[i, 1], sep="-")
    box = box[-1]
    if(tmp %in% box)
    {
        delete = c(delete, i)
    }
}

data = data[-delete,]

# 畫圖文件
data = data[data$raw_p_value <= 0.05,]
r_label = c()
for(i in 1:nrow(data))
{
    if(data[i, 3] < 0)
    {
        r_label = c(r_label, "neg")
    }
    else
    {
        r_label = c(r_label, "pos")
    }
}

data$r_label = r_label
write.table(data, file="Result/input_pre.txt", sep="\t", row.names=F, quote=F)

data$r_value = abs(data$r_value)
input_network = data[, c(1,2,3,6)]
write.table(input_network, file="Result/input_network.txt", sep="\t", row.names=F, quote=F)

2 cytoscape繪圖

繪圖流程見我的另一篇：Cytoscape繪制相關網絡圖
網絡圖結果如下：

圓圈表示KO，相關越多degree齿椅，圓越大
先粗細表示|r|琉挖，紅色負相關，藍色正相關