代码实战系列 | 自动计算数据全特征相关系数的p-value矩阵及其置信矩阵

onlybugs 轻松玩转生信 2023-07-14 20:30 发表于湖北
收录于合集

“ 本文将讲解如何使用R语言基础环境实现相关性矩阵的p-value矩阵以及其置信区间矩阵的计算并可视化。”




本篇文章的主要内容有：

  • 实现相关性矩阵显著性检验计算函数框架。


  • 完善函数并通过操作全局symbol来支持结果自动保存为CSV。


  • 优化代码过程。




关于相关性矩阵可视化的全部内容都展示在上次的文章中：

带显著性标记的组合型相关性矩阵

onlybugs，公众号：轻松玩转生信高级可视化系列 | 带显著性标记的组合型相关性矩阵

部分结果图：

图片


不会绘制相关性矩阵的小伙伴们可以先学习上期推文的内容哦！


接下来，我们就要开始动手实现显著性矩阵的自动计算了！



01

—

实现相关性矩阵显著性检验计算函数框架


首先，我们要知道R语言其实提供了一个基本的相关系数显著性检验的函数cor.test()。但是此函数的使用只能支持两个变量之间的计算，如

cor.test(mt$mpg,mt$cyl)
  Pearson's product-moment correlation
data:  mt$mpg and mt$cylt = -8.9197, df = 30, p-value = 6.113e-10alternative hypothesis: true correlation is not equal to 095 percent confidence interval: -0.9257694 -0.7163171sample estimates:      cor -0.852162
‍


通过计算我们其实是可以得到两个向量之间相关系数的多种信息，主要包括统计量，自由度，p-value，置信区间以及相关系数的值。因此，我们的任务就是配对的计算所有的值，并且提取结果，保存结果，输出结果等等。


接下来我会给出代码的整体流程框架，一步一步的注释以及实现功能和输出：

mt = mtcarsmcor = cor(mt)> mcor            mpg        cyl       disp         hp        drat         wt        qsecmpg   1.0000000 -0.8521620 -0.8475514 -0.7761684  0.68117191 -0.8676594  0.41868403cyl  -0.8521620  1.0000000  0.9020329  0.8324475 -0.69993811  0.7824958 -0.59124207disp -0.8475514  0.9020329  1.0000000  0.7909486 -0.71021393  0.8879799 -0.43369788hp   -0.7761684  0.8324475  0.7909486  1.0000000 -0.44875912  0.6587479 -0.70822339drat  0.6811719 -0.6999381 -0.7102139 -0.4487591  1.00000000 -0.7124406  0.09120476wt   -0.8676594  0.7824958  0.8879799  0.6587479 -0.71244065  1.0000000 -0.17471588


逻辑框架：

# 提取特征数量 一般都是列nc = ncol(mt)# 外圈循环 用于抽取cor.test()的xfor (ri in 1:nc) {  # 内圈循环 获取cor.test()的y  for (ci in ri:nc) {    # 测试输出    print(paste(ri,'-',ci))  }  }
[1] "1 - 1"[1] "1 - 2"[1] "1 - 3"[1] "1 - 4"[1] "1 - 5"[1] "1 - 6"[1] "1 - 7"[1] "1 - 8"[1] "1 - 9"[1] "1 - 10"[1] "1 - 11"[1] "2 - 2"[1] "2 - 3"[1] "2 - 4"[1] "2 - 5"[1] "2 - 6"[1] "2 - 7"[1] "2 - 8"[1] "2 - 9"[1] "2 - 10"[1] "2 - 11"


下面我们的任务就是抽取数据并计算：

# 提取特征数量 一般都是列nc = ncol(mt)# 外圈循环 用于抽取cor.test()的xfor (ri in 1:nc) {  # 内圈循环 获取cor.test()的y  for (ci in ri:nc) {    # 捕获cor.test的结果    test_res = cor.test(mt[,ri],mt[,ci])    print(paste(test_res$p.value,test_res$conf.int))  }  }[1] "0 1" "0 1"[1] "6.11268714258096e-10 -0.925769361912065" "6.11268714258096e-10 -0.716317141481634"[1] "9.3803265373813e-10 -0.923359365005218" "9.3803265373813e-10 -0.708137602946971"[1] "1.78783525412106e-07 -0.885268613495523" "1.78783525412106e-07 -0.586099412685315"[1] "1.77623992875242e-05 0.436048383960751" "1.77623992875242e-05 0.83220104009279" [1] "1.29395870135052e-10 -0.933826413284994" "1.29395870135052e-10 -0.744087196460113"[1] "0.0170819884965196 0.0819548679090099" "0.0170819884965196 0.669618639209138" [1] "3.41593725441997e-05 0.410363014099788" "3.41593725441997e-05 0.822326236635468"


这里数据已经抽取出来了(注意test_res是list对象不是一个df对象)，我们现在要考虑存储的问题，构造三个矩阵来存储结果。

> pvalue_mat = lowCI_mat = upCI_mat = matrix(data = NA,nrow = nc,ncol = nc)> pvalue_mat      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [1,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA [2,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA [3,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA [4,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA [5,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA


通过使用矩阵创建函数来创建充满缺失值的矩阵，然后使用R的特定语法连等来一次性创建三个同样的矩阵（Py也有同样的语法）。


下面存值：

# pvalue lowCI upCIpvalue_mat = lowCI_mat = upCI_mat = matrix(data = NA,nrow = nc,ncol = nc)pvalue_mat
# 外圈循环 用于抽取cor.test()的xfor (ri in 1:nc) {  # 内圈循环 获取cor.test()的y  for (ci in ri:nc) {    # 捕获cor.test的结果    test_res = cor.test(mt[,ri],mt[,ci])    # 抽取结果并存储    pvalue_mat[ri,ci] = pvalue_mat[ci,ri] = test_res$p.value    lowCI_mat[ri,ci] = lowCI_mat[ci,ri] = test_res$conf.int[1]    upCI_mat[ri,ci] = upCI_mat[ci,ri] = test_res$conf.int[2]  }}> pvalue_mat              [,1]         [,2]         [,3]         [,4]         [,5]          [,6] [1,] 0.000000e+00 6.112687e-10 9.380327e-10 1.787835e-07 1.776240e-05  1.293959e-10 [2,] 6.112687e-10 0.000000e+00 1.802838e-12 3.477861e-09 8.244636e-06  1.217567e-07 [3,] 9.380327e-10 1.802838e-12 0.000000e+00 7.142679e-08 5.282022e-06  1.222320e-11 [4,] 1.787835e-07 3.477861e-09 7.142679e-08 0.000000e+00 9.988772e-03  4.145827e-05 [5,] 1.776240e-05 8.244636e-06 5.282022e-06 9.988772e-03 0.000000e+00  4.784260e-06 [6,] 1.293959e-10 1.217567e-07 1.222320e-11 4.145827e-05 4.784260e-06 2.271770e-236 [7,] 1.708199e-02 3.660533e-04 1.314404e-02 5.766253e-06 6.195826e-01  3.388683e-01 [8,] 3.415937e-05 1.843018e-08 5.235012e-06 2.940896e-06 1.167553e-02  9.798492e-04


到这里，我们的计算框架其实就已经写完了！




02

—

包装函数并支持结果自动保存为CSV


一般情况下，我们的代码的每一个功能模块都要包装为一个类或者一个函数，由于我们这里只是一个简单的计算P矩阵，所以包装为一个函数就足够了，下面开始包装为函数并定义返回值

CalPvalueMat <- function(mt){  # 提取特征数量 一般都是列  nc = ncol(mt)    # pvalue lowCI upCI  pvalue_mat = lowCI_mat = upCI_mat = matrix(data = NA,nrow = nc,ncol = nc)  pvalue_mat    # 外圈循环 用于抽取cor.test()的x  for (ri in 1:nc) {    # 内圈循环 获取cor.test()的y    for (ci in ri:nc) {      # 捕获cor.test的结果      test_res = cor.test(mt[,ri],mt[,ci])      # 抽取结果并存储      pvalue_mat[ri,ci] = pvalue_mat[ci,ri] = test_res$p.value      lowCI_mat[ri,ci] = lowCI_mat[ci,ri] = test_res$conf.int[1]      upCI_mat[ri,ci] = upCI_mat[ci,ri] = test_res$conf.int[2]    }  }  # 添加特征名  colnames(pvalue_mat) = rownames(pvalue_mat) = colnames(mt)  colnames(lowCI_mat) = rownames(lowCI_mat) = colnames(mt)  colnames(upCI_mat) = rownames(upCI_mat) = colnames(mt)  list(p = pvalue_mat,lCI = lowCI_mat,uCI = upCI_mat)}
pCI = CalPvalueMat(mt)


这里的返回值是一个list，使用$运算符取出结果就能使用我们的函数啦。这里要注意，R语言可以隐式返回，就是我们不需要显示的调用return，代码的最后一行变量默认作为返回值，我们只需要使用一个变量接受它就可以了。


下面用我们的结果随便画个corrplot吧

pCI = CalPvalueMat(mt)corrplot(mcor,         tl.col = 'black',         col = rev(COL2('RdBu', 200)),         method = 'square',          p.mat = pCI$p,         sig.level = .05)

图片


下面提供第三个功能，自动保存。其实我很多时候使用别人的函数或者自己写代码的时候都会主动添加这个功能，因为运行以后就给你存储起来不要太方便好嘛！！！尤其是大型计算，一次跑三四天的那种，忘了保存或者环境崩溃了就是噩梦，这比细胞养死了或者跑深度学习忘了存断点模型痛苦好几倍！！！

CalPvalueMat <- function(mt,fn = substitute(mt),is_save = T){    # 提取特征数量 一般都是列  nc = ncol(mt)    # pvalue lowCI upCI  pvalue_mat = lowCI_mat = upCI_mat = matrix(data = NA,nrow = nc,ncol = nc)  pvalue_mat    # 外圈循环 用于抽取cor.test()的x  for (ri in 1:nc) {    # 内圈循环 获取cor.test()的y    for (ci in ri:nc) {      # 捕获cor.test的结果      test_res = cor.test(mt[,ri],mt[,ci])      # 抽取结果并存储      pvalue_mat[ri,ci] = pvalue_mat[ci,ri] = test_res$p.value      lowCI_mat[ri,ci] = lowCI_mat[ci,ri] = test_res$conf.int[1]      upCI_mat[ri,ci] = upCI_mat[ci,ri] = test_res$conf.int[2]    }  }    colnames(pvalue_mat) = rownames(pvalue_mat) = colnames(mt)  colnames(lowCI_mat) = rownames(lowCI_mat) = colnames(mt)  colnames(upCI_mat) = rownames(upCI_mat) = colnames(mt)      if(is_save){    write.csv(pvalue_mat,paste(fn,"_pvalue.csv",sep = ""))    write.csv(lowCI_mat,paste(fn,"_lowCI.csv",sep = ""))    write.csv(upCI_mat,paste(fn,"_upCI.csv",sep = ""))}    list(p = pvalue_mat,lCI = lowCI_mat,uCI = upCI_mat)}
pCI = CalPvalueMat(mt)


这里面我们使用了通过变量symbol实现自动命名的技巧：

mt,fn = substitute(mt),is_save = T

这里的substitute函数可以获取变量的symbol名称，然后后面代码会根据is_save以及symbol自动命名。

图片

pCI = CalPvalueMat(mt_new)

图片


这个技巧的原理就在于R作为解释语言，R解释器会维护全局符号表，可以由此玩一些比较有趣的操作，不过用c/cpp实现此功能就有些痴心妄想了。



03

—

优化代码过程


终于进行到文章的最后部分了，很多人可能会问，代码已经写好了，还要做什么呢？我们虽然已经有了一组代码，它能用，但是它依旧还有很多小问题。比如：

代码能不能再快一些，循环次数是不是太多了？

相关系数检验还有其它的参数，如置信度设置，我们该如何传参来支持此功能呢？


我们首先解决第一个问题，我们知道任何特征对自己的相关系数都是1，p值是0，致信区间是1不变，我们应该跳过此位置并尽可能的初始化结果为合适结果。

CalPvalueMat <- function(mt,fn = substitute(mt),is_save = T){    # 提取特征数量 一般都是列  nc = ncol(mt)    # pvalue lowCI upCI  # 新的初始化方案 p矩阵直接初始化为0  pvalue_mat = matrix(data = 0,nrow = nc,ncol = nc)  # CI矩阵初始化为1  lowCI_mat = upCI_mat = matrix(data = 1,nrow = nc,ncol = nc)    # 优化循环次数  for (ri in 1:(nc-1)) {    for (ci in (ri+1):nc) {      # 捕获cor.test的结果      test_res = cor.test(mt[,ri],mt[,ci])      # 抽取结果并存储      pvalue_mat[ri,ci] = pvalue_mat[ci,ri] = test_res$p.value      lowCI_mat[ri,ci] = lowCI_mat[ci,ri] = test_res$conf.int[1]      upCI_mat[ri,ci] = upCI_mat[ci,ri] = test_res$conf.int[2]    }  }    colnames(pvalue_mat) = rownames(pvalue_mat) = colnames(mt)  colnames(lowCI_mat) = rownames(lowCI_mat) = colnames(mt)  colnames(upCI_mat) = rownames(upCI_mat) = colnames(mt)      if(is_save){    write.csv(pvalue_mat,paste(fn,"_pvalue.csv",sep = ""))    write.csv(lowCI_mat,paste(fn,"_lowCI.csv",sep = ""))    write.csv(upCI_mat,paste(fn,"_upCI.csv",sep = ""))}    list(p = pvalue_mat,lCI = lowCI_mat,uCI = upCI_mat)}

注意循环次数的优化！


接下来是额外参数的问题，R拥有一种神奇的语法...，这和Py中的*args和**kwargs有异曲同工之妙，不过R的更加简单直白。

CalPvalueMat <- function(mt,fn = substitute(mt),is_save = T,...){    # 提取特征数量 一般都是列  nc = ncol(mt)    # pvalue lowCI upCI  pvalue_mat = matrix(data = 0,nrow = nc,ncol = nc)  lowCI_mat = upCI_mat = matrix(data = 1,nrow = nc,ncol = nc)
  # 外圈循环 用于抽取cor.test()的x  for (ri in 1:(nc-1)) {    # 内圈循环 获取cor.test()的y    for (ci in (ri+1):nc) {      # 捕获cor.test的结果      test_res = cor.test(mt[,ri],mt[,ci],...)      # 抽取结果并存储      pvalue_mat[ri,ci] = pvalue_mat[ci,ri] = test_res$p.value      lowCI_mat[ri,ci] = lowCI_mat[ci,ri] = test_res$conf.int[1]      upCI_mat[ri,ci] = upCI_mat[ci,ri] = test_res$conf.int[2]    }  }    colnames(pvalue_mat) = rownames(pvalue_mat) = colnames(mt)  colnames(lowCI_mat) = rownames(lowCI_mat) = colnames(mt)  colnames(upCI_mat) = rownames(upCI_mat) = colnames(mt)      if(is_save){    write.csv(pvalue_mat,paste(fn,"_pvalue.csv",sep = ""))    write.csv(lowCI_mat,paste(fn,"_lowCI.csv",sep = ""))    write.csv(upCI_mat,paste(fn,"_upCI.csv",sep = ""))}    list(p = pvalue_mat,lCI = lowCI_mat,uCI = upCI_mat)}
pCI = CalPvalueMat(mt_new,is_save = F,conf.level = 0.95)


这样就简单有效的解决了问题。注意这里的...是可以输出的，将其转化为list对象就可以输出查看了。



这次的文章就先到这里了，下一期的话可能做一期涉及正则表达式的sed或者grep教程或者花式绘制一些常见的统计图。



最后，如果你觉得这篇文章让你学到了知识，就请你动动小手帮忙点一点赞和在看，赠人玫瑰，手有余香，你对我的认可就是我更新下去的最大动力！


-结束-



点击关注公众号，学习资料免费送，学生信，不迷路图片


喜欢此内容的人还喜欢

微信扫一扫
关注该公众号