【for循环对丰度表格中各个物种与所有环境变量进行多元线性回归，获得回归R2及各变量对R2的贡献度】

最新推荐文章于 2023-09-10 14:13:41 发布

jnope

最新推荐文章于 2023-09-10 14:13:41 发布

阅读量1.6k

点赞数 14

文章标签： r语言线性回归矩阵

本文链接：https://blog.csdn.net/jnope/article/details/124251130

版权

文章目录

背景
一、数据
二、原理
三、代码
四、结果文件

背景

大背景是研究16S扩增子测序的群落结构，已经获得不同物种水平（如门水平）在不同样品中的丰度表格和各个样本的环境因子（pH和各种离子浓度等），需要利用R分析环境因子与环境中物种相关性，参考教程（小白鱼的生统笔记，链接：

https://mp.weixin.qq.com/s?__biz=MzIxNzc1Mzk3NQ==&mid=2247486604&idx=1&sn=7b2bb39b8ee238b377760ebec1822885&chksm=97f5bc94a0823582790a238a5f13f113ffb01380204ca81b1c510aa652d451985eb1207bfaab&token=85792652&lang=zh_CN&scene=21#wechat_redirect

）。该过程涉及到利用多元线性回归分析——使用所有环节因子对每个物种的丰度进行拟合。需要对数据框进行循环操作，教程中采用手动一个个计算，本文对其进行循环操作。

一、数据

门水平OTU表格

前面步骤根据教程进行，将门水平物种丰度表格读入R后进行转至得到如下形式，行名为样品名（未全部显示），列名为物种名（门水平，未全部显示），由于已经求过相对丰度值，各行之和均为1。

phylum <- read.delim('zotu_phylum_per_nocoal.txt',sep = "\t",check.names = F,
                     row.names = 1,stringsAsFactors = FALSE)
phylum <- t(phylum)

在这里插入图片描述

环境因子表格

index <- index <- read.delim('index.txt',sep = "\t",check.names = F,
                             row.names = 1,stringsAsFactors = FALSE)

读入R后格式如下，行为样品名，列为各个理化指标值
在这里插入图片描述

二、原理

采用前向选择来减少环境变量；每个物种可能得出与0、1和多个环境变量相关，其中0个即无环境变量与其变化显著相关；1表示仅1个变量相关，此时求得的R2即为该环境变量重要性；多个环境变量相关时，可以计算每个变量的贡献度，其和为回归的R2。
代码使用for循环对丰度表格按列进行回归分析，并根据相关性变量数量分类后写入原本空的数据框内。

三、代码

#如果没有安装相应包，运行以下几句
# install.packages('relaimpo')
# library(BiocManager)
# BiocManager::install('packfor')
# install.packages('adespatial')
library(relaimpo)
library(adespatial)
#创建空的回归结果数据框
lm_result <- data.frame(matrix(,nrow=ncol(phylum),ncol=3))
rownames(lm_result) <- colnames(phylum)
colnames(lm_result) <- c('r.squared','adj.r.squared','p.value')
#创建空的环境因子重要性数据框
env_importance <- data.frame(matrix(,nrow=ncol(index),ncol=ncol(phylum)))
rownames(env_importance) <- colnames(index)
colnames(env_importance) <- colnames(phylum)
#循环计算环境因子与各个物种的拟合关系
for(i in colnames(phylum)){
  tryCatch({
	#前向选择
    forward.select <- forward.sel(Y=phylum[,i],X=index,
                                  R2more = 0.01,alpha = 0.05,nperm = 999)
    #无环境变量相关时跳过本次循环
    if(nrow(forward.select)==0){next}
    #至少一个变量相关时：
    env_sp <- cbind(index[forward.select$variables],phylum[,i])
    colnames(env_sp)[ncol(env_sp)] <- 'phylum'
    #进行线性回归
    fit_sample <- lm(phylum~.,data=env_sp)
    lm_stat <- summary(fit_sample)
    pF <- lm_stat$fstatistic
    #获取回归p值
    p_value <- pf(pF[1], pF[2], pF[3], lower.tail = FALSE)
    #结果写入回归表格中
    lm_result[i,] <- c(lm_stat$r.squared,lm_stat$adj.r.squared,p_value)
    #仅一个变量相关时
    if(nrow(forward.select)==1){
      #将单个变量的r2作为重要性
      env_importance[forward.select$variables,i] <- lm_stat$r.squared
      next
    }
    #多个变量相关时，获取每个变量贡献度并写入环境因子贡献度表格
    crf <- calc.relimp(fit_sample,rela=F)
    for(j in names(crf$lmg)){
      env_importance[j,i] <- crf$lmg[j]
    }
  }, 
  error=function(e){cat('error: ',i,conditionMessage(e),'\n')})
}
#在回归表格中添加显著性标记
lm_result[which(lm_result$p.value<0.001),'sig'] <- '***'
lm_result[which(lm_result$p.value<0.01 & lm_result$p.value>0.001),'sig'] <- '**'
lm_result[which(lm_result$p.value<0.05 & lm_result$p.value>0.01),'sig'] <- '*'
#输出文件
write.table(lm_result,'lm_resul.txt',sep='\t',
            row.names = T, col.names = T,quote = F)
write.table(env_importance,'env_importance.txt',sep='\t',
            row.names = T, col.names = T,quote = F)