专栏首页优雅R如何在相关矩阵图上添加p-value

如何在相关矩阵图上添加p-value

这篇文章可以说是《怎么分析和展示RNAseq基因表达数据中基因的相关性》的延续。上次绘制了下图:

可以发现只有两个基因的表达表现出了较强的相关(ETV3-ELK4)。 一般教材描述相关性大小为:

相关系数r 是否是:·|r|>0.95 存在显著性相关; ·|r|≥0.8 高度相关; ·0.5≤|r|<0.8 中度相关;·0.3≤|r|<0.5 低度相关;·|r|<0.3 关系极弱,认为不相关 计算公式为:

公式

可见这两个基因属于中度相关。

但是我们知道皮尔逊相关系数表示的是两组数据线性相关的程度,但是如果两者在统计学上不存在相关性呢?那这个指标还有什么意义?因此,我们在评判相关的时候需要同时考量p值和r相关系数大小。

一个网络博主是这样认为的:

看两者是否算相关要看两方面:显著水平以及相关系数 (1)显著水平,就是P值,这是首要的,因为如果不显著,相关系数再高也没用,可能只是因为偶然因素引起的,那么多少才算显著,一般p值小于0.05就是显著了;如果小于0.01就更显著;例如p值=0.001,就是很高的显著水平了,只要显著,就可以下结论说:拒绝原假设无关,两组数据显著相关也说两者间确实有明显关系.通常需要p值小于0.1,最好小于0.05设甚至0.01,才可得出结论:两组数据有明显关系,如果p=0.5,远大于0.1,只能说明相关程度不明显甚至不相关.起码不是线性相关. (2)相关系数,也就是Pearson Correlation(皮尔逊相关系数),通常也称为R值,在确认上面指标显著情况下,再来看这个指标,一般相关系数越高表明两者间关系越密切.

从上面看来,在进行相关分析考量相关系数r(或者R2)前,先考量显著性是有必要的。不过,如果你已经能看到两变量有很明显的线性关系了,你不看p值也无伤大雅,因为那个时候p值肯定少于0.05。

参考R包corrplot文档对上次的函数进行优化,代码如下:

gene_exp.corr <- function(gene.list, project_code, project.clinical, project.exp, outdir, ID_transform=TRUE, conf.level=0.95){
    # Arguments:
    # gene.list: a list of gene you want to analyze their expression correlation
    # project_code: data project name or name you want to specify this analysis
    # project.clinical: clinical information about samples, data.frame format
    # project.exp: normalized gene expression (RNA seq) about samples, data.frame format
    # ID_transform: sometimes clinical information use "-" as separate symbol for sample ID,
    #               we need it to be the same as it in project.exp data
    # one sample ID example: in clinical information, one sample may be marked by "TCGA-3N-A9WB-06",
    #                        in RNA seq data.set, this sample is "TCGA.3N.A9WB.06". If it is not, set ID_transform=FALSE. 

    # note: you need to install "corrgram" package before use this function

    gene_exp.list <- subset(project.exp, sample%in%gene.list)
    rownames(gene_exp.list) <- gene_exp.list[,1]
    gene_exp.list <- gene_exp.list[,-1]
    gene_exp.list <- t(gene_exp.list)
    # gene_exp.list <- gene_exp.list[,c(5,1,2,3,4,6,7,8,9,10)]
    library(corrplot)
    # combine with significance test
    cor.mtest <- function(mat, conf.level = 0.95){
        mat <- as.matrix(mat)
        n <- ncol(mat)
        p.mat <- lowCI.mat <- uppCI.mat <- matrix(NA, n, n)
        diag(p.mat) <- 0
        diag(lowCI.mat) <- diag(uppCI.mat) <- 1
        for(i in 1:(n-1)){
            for(j in (i+1):n){
                tmp <- cor.test(mat[,i], mat[,j], conf.level = conf.level)
                p.mat[i,j] <- p.mat[j,i] <- tmp$p.value
                lowCI.mat[i,j] <- lowCI.mat[j,i] <- tmp$conf.int[1]
                uppCI.mat[i,j] <- uppCI.mat[j,i] <- tmp$conf.int[2]
            }
        }
        return(list(p.mat, lowCI.mat, uppCI.mat))
    }

    if(ID_transform){
        project.clinical$sampleID = gsub("-",".",project.clinical$sampleID, fixed = TRUE)
    }
    n.gene <- ncol(gene_exp.list)
    # all samples
    M1 <- cor(gene_exp.list)
    res1 <- cor.mtest(gene_exp.list, conf.level)

    pdf(paste(outdir,project_code,"_all_sample_genelist_expression_corrgram.pdf", sep=""))
    corrplot(M1, order = "AOE", tl.pos = "d", p.mat = res1[[1]], insig = "p-value")
    title(paste("Corrgram of ", n.gene," Genes Expression in ", project_code, sep = ""))
    dev.off()

    # choose tumor sample
    # table(project.clinical$sample_type)
    primary_tumor <- "Primary Tumor"
    Metast_tumor  <- "Metastatic"
    primary_tumor.id <- project.clinical[project.clinical$sample_type==primary_tumor,]$sampleID
    Metast_tumor.id  <- project.clinical[project.clinical$sample_type==Metast_tumor,]$sampleID

    if(length(primary_tumor.id)<2 & length(Metast_tumor.id<2)){
        stop("Maybe your data have something wrong. Please check it!")
    }else{
        if(length(primary_tumor.id)<2){
            stop("I don't think it's reasonable that there are less than 2 primary tumor samples.")}

        gene_exp.list.primary <- subset(gene_exp.list, rownames(gene_exp.list)%in%primary_tumor.id)
        M2 <- cor(gene_exp.list.primary)
        res2 <- cor.mtest(gene_exp.list.primary, conf.level)

        pdf(paste(outdir,project_code,"_primary_tumor_sample_genelist_expression_corrgram.pdf", sep=""))
        corrplot(M2, order = "AOE", tl.pos = "d", p.mat = res2[[1]], insig = "p-value")
        title(paste("Corrgram of ", n.gene," Genes Expression in ", project_code, sep = ""))
        dev.off()

        if(length(Metast_tumor.id)<2){
            cat("It seems has no Metastatic sample in this analysis. \n")
            return(0)
        }

        gene_exp.list.Metast  <- subset(gene_exp.list, rownames(gene_exp.list)%in%Metast_tumor.id)
        M3 <- cor(gene_exp.list.Metast)
        res3 <- cor.mtest(gene_exp.list.Metast, conf.level)
        pdf(paste(outdir,project_code,"_Metastatic_sample_genelist_expression_corrgram.pdf", sep=""))
        corrplot(M3, order = "AOE", tl.pos = "d", p.mat = res3[[1]], insig = "p-value")
        title(paste("Corrgram of ", n.gene," Genes Expression in ", project_code, sep = ""))
        dev.off()}
}

一方面增加了检验部分,另一方面修改了画图函数。如果你想用这个函数绘制更多自定义的图,可以参考R包文档进行修改。

corrplot_demo.png

这跟上次的图类似。偏蓝色代表正相关,偏红色代表负相关。用圆圈大小和颜色鲜艳程度辅助color legend可以很好的区分和找出相关性明显的。corrgram的左下角和右上角是对称的,标有数字的显示的是p值,我这里默认设定0.05为阈值,大于0.05都会显示出来,这些值说明对应的两个基因在统计学上是没有相关性的。

本文分享自微信公众号 - 优雅R(elegant-r),作者:王诗翔

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-09-26

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 怎么分析和展示RNAseq基因表达数据中基因的相关性

    TCGA是癌症基因组分析中相当流行的数据库,针对里面数据的挖掘结果、软件工具发表了许多CNS文章,不过现在已经被整合进GDC数据平台了。虽然现在测序技术发展的很...

    王诗翔呀
  • 「R」UCSCXenaTools 文档(英文)

    UCSCXenaTools is an R package for downloading and exploring data from UCSC Xena ...

    王诗翔呀
  • TidyFriday 快为你的 R 语言基础绘图系统设定绘图主题吧!

    basetheme 的魔力就在于,配置好基础绘图系统的主题之后你可能会再次爱上 R 语言的基础绘图系统!

    王诗翔呀
  • Qt编写地图综合应用13-获取边界点

    获取边界点一般和行政区划搭配起来使用,比如用户输入一个省市的名称,然后自动定位到该省市,然后对该轮廓获取所有边界点集合输出到js文件,最后供离线使用,获取边界点...

    feiyangqingyun
  • Qt编写地图综合应用9-行政区划

    行政区划在地图应用中非常有用,行政区划是行政区域划分的简称,是国家为了进行分级管理而实行的区域划分,百度地图提供的内置的函数类支持传入行政区划的名称来获取对应的...

    feiyangqingyun
  • Python学习中的点点滴滴

    明天依旧可好
  • 爬虫案例|抓取拉勾网职位需求关键词,并生成统计图

    知己知彼,方可百战不殆。在学习技术的时候我们往往面临太多选择而不知所措,可能是各个方面都有涉猎,对某个领域没有深入研究,看似什么都会,真要让你做个什么东西的时候...

    统计学家
  • 怎么分析和展示RNAseq基因表达数据中基因的相关性

    TCGA是癌症基因组分析中相当流行的数据库,针对里面数据的挖掘结果、软件工具发表了许多CNS文章,不过现在已经被整合进GDC数据平台了。虽然现在测序技术发展的很...

    王诗翔呀
  • Android 模拟登陆网站实现移动客户端

    叶应是叶
  • 超几何分布与二项分布及其期望

    惊奇的发现选修2-3上有期望的介绍,不过我没有课本啊qwq。只能去网上找资料了。。

    attack

扫码关注云+社区

领取腾讯云代金券