【原】微生物组大数据分析如何矫正不同研究误差？

写在前面

• github地址：https://github.com/wdl2459/ConQuR

• 英文教程地址：https://wdl2459./ConQuR/ConQuR.Vignette.html

简介

基于条件分位数回归（ConQuR）的微生物组数据批次效应去除，概述该软件包通过条件分位数回归方法从物种读数计数表中去除批次效应。同时考虑微生物组数据的分布属性，如零膨胀和过度离散。将生成一个去除批次效应后的物种读数计数表。

实战

数据介绍

我们将使用包中的样本数据来演示如何使用ConQuR。该数据集包含来自273个样本的100个微生物，这些微生物信息在3个批次中进行了测序。数据集包括批次ID（样本名字）和分组数据：sbp（血压，关键变量，连续性），性别（协变量1，二元），种族（协变量2，二元）和年龄（协变量3，连续性）。

批次差异使用的方法

• 逻辑回归是一种统计学习方法，用于建立分类模型，它将输入变量与二元输出变量之间的关系建模。逻辑回归可以用于解决二元分类问题，例如根据某些特征将个体分类为阳性或阴性。

• 标准分位数回归是一种回归分析方法，它将一个或多个自变量与因变量之间的关系建模，并且可以对不同的分位数进行建模。在标准分位数回归中，使用的是分位数回归系数，而不是最小二乘回归系数。这种方法常用于处理非正态分布的数据，特别是对于处理异常值具有较好的鲁棒性。标准分位数回归可以用于解决连续变量的回归问题，例如预测一个人的收入与其教育水平、工作经验等因素之间的关系。

• Logistic LASSO回归是一种统计学习方法，是基于逻辑回归模型的一种变量选择技术，它可以用于处理高维数据集，避免模型过拟合问题。在逻辑LASSO回归中，采用L1正则化约束来对模型进行稀疏化处理，即将某些系数置为零，从而实现特征选择。LASSO模型可以通过最小化目标函数来实现，该目标函数是对逻辑回归损失函数和L1正则项的加权和。与传统的逻辑回归模型相比，Logistic LASSO回归可以帮助减少不必要的变量，提高模型的预测精度和解释能力，并且在具有多重共线性的情况下可以更好地处理数据。它在分类问题中得到了广泛的应用，例如在医学研究中用于预测疾病风险、在金融领域中用于信用评估等。

• Quantile LASSO回归是一种回归分析方法，是基于标准分位数回归模型的一种变量选择技术。它通过在回归模型中添加L1正则化约束来实现特征选择和模型压缩，从而避免过拟合问题。相比于传统的最小二乘回归模型，Quantile LASSO回归可以更好地处理离群值和数据分布的偏斜问题。在Quantile LASSO回归中，采用的是分位数回归系数，而不是最小二乘回归系数，因此可以通过调整分位数来获得不同的回归结果。这种方法对于处理非正态分布的数据集具有较好的鲁棒性和稳健性，并且可以用于预测中位数、高分位数或低分位数等目标变量。Quantile LASSO回归在许多领域中得到了广泛的应用，例如在金融领域中用于股票价格预测、在医学研究中用于预测疾病风险等。

参数解析：

总的来说，推荐使用ConQuR的优化结果。我们有以下选项：

• batch_ref_pool：候选参考批次，例如c（“0”，“1”）。

• logistic_lasso_pool：是否使用逻辑LASSO回归，例如c（T，F）。

• quantile_type_pool：分位数回归类型，例如c（“standard”，“lasso”，“composite”）。

• simple_match_pool：是否使用简单的分位数匹配，例如c（T，F）。

• lambda_quantile_pool：分位数LASSO或复合分位数回归的惩罚参数的候选项，例如c（NA，“2p/n”，“2p/logn”）。注意：如果quantile_type_pool中包括“standard”，则始终包括NA。

• interplt_pool：是否在分段估计中使用插值，例如c（T，F）。

• frequencyL：需要优化的分类单位频率下限，例如0。

• frequencyU：需要优化的分类单位频率上限，例如1。

• cutoff：进行优化的分类单位频率截止值，例如0.1。

代码实战

# # 安装R包
# # library("doParallel")
# devtools::install_github("wdl2459/ConQuR")
# BiocManager::install("coda.base")

library(ConQuR)
library(doParallel)

内置数据介绍

# 运行一下教程（来自官方英文教程）#----------
# https://wdl2459./ConQuR/ConQuR.Vignette.html

data(Sample_Data)

taxa = Sample_Data[, 1:100]
taxa[146:150, 1:5]

#--批次因素
batchid = Sample_Data[, 'batchid']
summary(batchid)

# 提取分组信息
covar = Sample_Data[, c('sex', 'race', 'age', 'sbp')]
summary(covar)

colnames(taxa)

选择逻辑回归和分位数回归矫正误差：

options(warn=-1)
taxa_corrected1 = ConQuR(tax_tab=taxa, batchid=batchid, covariates=covar, batch_ref="0")
taxa_corrected1[146:150, 1:5]
colnames(taxa_corrected1)

加上lasso变量选择：

options(warn=-1)
taxa_corrected2 = ConQuR(tax_tab=taxa, 
                         batchid=batchid, 
                         covariates=covar, 
                         batch_ref="0",
                         logistic_lasso=TRUE,
                         quantile_type="lasso", 
                         interplt=TRUE)

taxa_corrected2[146:150, 1:5]

通过排序分析展示矫正效果

par(mfrow=c(2, 3))

Plot_PCoA(TAX=taxa, factor=batchid, main="Before Correction, Bray-Curtis")
Plot_PCoA(TAX=taxa_corrected1, factor=batchid, main="ConQuR (Default), Bray-Curtis")
Plot_PCoA(TAX=taxa_corrected2, factor=batchid, main="ConQuR (Penalized), Bray-Curtis")

Plot_PCoA(TAX=taxa, factor=batchid, dissimilarity="Aitch", main="Before Correction, Aitchison")
Plot_PCoA(TAX=taxa_corrected1, factor=batchid, dissimilarity="Aitch", main="ConQuR (Default), Aitchison")
Plot_PCoA(TAX=taxa_corrected2, factor=batchid, dissimilarity="Aitch", main="ConQuR (Penalized), Aitchison")

通过置换检验查看：

#--评估批次效应造成的误差是否降低，而由关注分组误差导致误差变化不大。# 具体细节：使用bray距离基于原始count和相对丰度转化的三个置换检验R方；#  第一行是批次误差的R方法，第二行是关注的处理
PERMANOVA_R2(TAX=taxa, batchid=batchid, covariates=covar, key_index=4)

PERMANOVA_R2(TAX=taxa_corrected1, batchid=batchid, covariates=covar, key_index=4)

PERMANOVA_R2(TAX=taxa_corrected2, batchid=batchid, covariates=covar, key_index=4)

通过随机森林回归预测查看：

#--矫正后表格是否可以更好预测关键变量
sbp = covar[, 'sbp']
taxa_result = list(taxa, taxa_corrected1, taxa_corrected2)

pred_rmse = matrix(ncol=3, nrow=5)
colnames(pred_rmse) = c("Original", "ConQuR (Default)", "ConQuR (Penalized)")

for (ii in 1:3){
  pred_rmse[, ii] = RF_Pred_Regression(TAX=taxa_result[[ii]], variable=sbp)$rmse_across_fold
}

par(mfrow=c(1,1))
boxplot(pred_rmse, main="RMSE of Predicting SBP")

推荐函数-组合一下方法都进行尝试选择最优方案

#---推荐函数-组合一下方法都进行尝试选择最优方案
result_tuned = Tune_ConQuR(tax_tab=taxa, 
                           batchid=batchid, 
                           covariates=covar,
                           batch_ref_pool=c("0", "1"),
                           logistic_lasso_pool=F, 
                           quantile_type_pool=c("standard", "lasso"),
                           simple_match_pool=F,
                           lambda_quantile_pool=c(NA, "2p/n"),
                           interplt_pool=F,
                           frequencyL=0,
                           frequencyU=1)

# 多种方法组合置换检验结果
result_tuned$method_final

#--选择最优标准化方案
taxa_optimal = result_tuned$tax_final

par(mfrow=c(1, 2))
Plot_PCoA(TAX=taxa_optimal, factor=batchid, main="Fine-Tuned ConQuR, Bray-Curtis")
Plot_PCoA(TAX=taxa_optimal, factor=batchid, dissimilarity="Aitch",
          main="Fine-Tuned ConQuR, Aitchison")
根际互作生物学研究室 简介