背景
大背景是研究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)
四、结果文件
- 回归结果文件
lm_result:
行名为菌名,第一列为R2值,第二列为校正后R2,第三列为p值,第4列为显著性标记。
- 环境变量贡献度文件
env_importance:
行名为环境变量,列名为各个物种,NA表示不相关,列之和表示所有环境变量对相应物种回归的R2
输出文件中第一行行首没有制表符,建议自己加上;对于其中回归不清楚的可以查看背景里给出的链接,很详细。