中文应该只有我写了这个。。
“ DivE”估计量是一种启发式方法，用于估计种群中的种类数量或物种数量(物种丰富度)。
文章
Quantification of HTLV-1 Clonality and TCR Diversity (plos.org)
免疫学和微生物多样性的估计对于我们理解感染和免疫反应是至关重要的。例如，t 细胞指令表的多样性是什么？这些问题通过高通量测序技术得到部分解决，从而能够在样本中识别免疫学和微生物“物种”。根据样本多样性估计种群多样性需要估计未知物种的数量。在这里，我们测试了五个广泛使用的非参数估计，并发展和验证了一种新的方法，DivE，估计物种丰富度和分布。我们使用了3个独立的数据集: (i)来自感染1型人类嗜T细胞病毒病毒的人群; (ii) t 细胞抗原受体克隆体型; 和(iii)来自婴儿粪便样本的微生物数据。当应用于具有稀疏曲线且不会出现高原现象的数据集时，现有的估计量随样本容量的增加而系统地增加。相比之下，DivE 一致而准确地估计所有数据集的多样性。我们确定了限制 DivE 应用的条件。我们还表明，DivE 可以用来准确估计潜在的人口频率分布。我们已经开发了一种新的方法，这种方法比微生物和免疫群体中常用的生物多样性估计方法要精确得多。
免疫学和微生物学数据在重要方面与生态学数据不同。首先，在许多免疫和微生物群体中，可以合理地假定”物种”在分类学上是相似的，个体的空间分布是均匀的，个体是随机、独立和概率相等的抽样。如果做出这些简化的假设，就可以对基于个体的稀疏曲线进行外推，这些曲线描述了预期的物种数量与抽样的个体数量之间的关系[11]-[14]。然而，在生态种群[14]-[18]中，上述假设常常被违反。在这些生态种群中，未被观察到的个体在颜色、物理大小、地理分布、迁移、生境的多样性以及与其他物种的关系等方面可能与观察到的个体有所不同，因此，即使随后进行了大量抽样，仍然未被观察到。其次，许多关于种群结构的常见假设对于免疫学和微生物种群来说是不适当的，例如所有物种的频率都是相等的，或者种群分布的功能形式是已知的。因此，我们考虑非参数估计。
单样本丰度与多样性指数 (rstudio-pubs-static.s3.amazonaws.com)
非参数估计的方法包括，如 Chao1和基于丰度的覆盖估计(ACE)。ACE 被认为是目前最好的方法，在微生物学和免疫学中得到了广泛的应用，例如用于评估人类胃肠道菌群、人类肠道宏基因组、小鼠 TCR 全谱、[33]、真菌[34]以及 htlv-1感染细胞克隆的数量[35]。虽然它们最初是用来估计下界的方法，但是 chao1估计量和经过修正的偏倚校正的 Chao1bc [36] ，已经被用来对 TCR 克隆型[37]、[38]、丙型肝炎病毒感染的 OTUs [1]、疟疾感染的寄生虫多样性[39]宏基因组大小[40]、基因治疗载体的整合位点数[41]、土壤多样性[42]以及 htlv-1感染的细胞克隆数[35]进行点估计。除了 ACE 和 Chao 估计，我们还考虑了另外两个非参数估计: Bootstrap [44]和 Good-Turing 估计[45]。
大多数多样性估计器的目的是估计两个相关群体中的一个群体的物种丰富度: 或者是从样本中提取的群体(例如，肠道中的微生物种类数量，取自肠道的样本) ，或者是稀疏曲线达到饱和的数值(例如，进一步取样不产生任何新物种时的物种数量)。这些有关利益群体的定义缺乏灵活性，而且可能不适合或者不能很好地界定当前的问题。事实上，如果某些物种由单个个体表示，稀疏曲线将不会饱和。对于许多微生物学和免疫学问题，最好能有一个估计量，使用者能够确定所关注的人群的大小。例如，我们可能希望了解血液和整个身体的 t 细胞指令系统的多样性。
https://cran.r-project.org/web/packages/DivE/DivE.pdf
R包描述
DivE 将许多数学模型拟合到基于个体的稀疏的多个嵌套子样本曲线。这些曲线描绘了预期的物种数量与个体数量的函数关系。（例如 T 细胞、病毒粒子、微生物）。每个模型都适合所有嵌套的子样本，产生多模型拟合。新标准用于对每个模型的拟合一致性进行评分从嵌套的子样本中导出完整的观察到的稀疏曲线，即仅从不完整的数据。性能最好的模型被外推到所需的总量规模，并且它们的评估聚合以估计总体中的种类。
这个包内容包括：
- 生成基于个体的稀疏（物种积累）数据的功能，并评估它们的曲率
- 将数学模型拟合到稀疏数据及其嵌套子样本的函数。这些函数广泛使用 R 包 FME (http://cran.r-project.org/web/packages/FME/index.html)
- 评估每个模型的新标准的功能。这些函数使用了 R 包rgeos (http://cran.r-project.org/web/packages/rgeos/index.html)
- 对竞争模型进行评分的功能
- 产生类数（多样性）最终估计值的函数
- 示例候选模型、拟合参数、参数范围和示例数据集
- 一个示例脚本。我们试图使代码对需要不同的用户具有灵活性细节和控制水平。
使用该包的最简单方法是 DiveMaster 函数。这个函数是 DivE 包提供的其他函数的封装器，能够一次性完成：创建子样本（函数 DivSubsamples）、拟合模型（函数 FitSingleMod）、评分模型（函数ScoreSingleMod）并产生最终的多样性/物种丰富度估计（函数 PopDiversity）等函数的全部功能。（相当于一个不需要动脑的pepline）
对每个模型拟合进行评分的新标准是：
Discrepancy差异——数据点和模型预测之间的平均百分比误差。
Accuracy准确度——全样本物种丰富度之间的百分比误差,来自给定子样本的样本物种丰富度。
Similarity 相似性——拟合子样本的曲线与拟合完整样本的曲线之间的面积，归一化为来自完整数据的曲线下面积，在区间 [0, Nobs] 上，其中 Nobs 是完整数据的大小。
Plausibility合理性 ——预测的物种数量必须单调增加或平稳增长，物种积累的预测速率必须减少或平稳（即对于 S(x) 和 x ≥1，其中 x 是个体数，S'(x) ≥0，且 S''(x) ≤0)
一步运行完毕
但是注意这里使用的模型数量，它只使用了两个，我们可以一次性使用58个模型，但是这样就会比较慢，至于一次怎么使用全部模型，去掉循环就行了
######### Test script for DivE ###########
# This script runs the DiveMaster function from the DivE package on the example dataset 'Bact1' (using only 2 models)
# Required package: DivE
###########################################
#### Step 1. Load package ###
###########################################
require(DivE)
###########################################
#### Step 2. Load data files ####
###########################################
data(Bact1) # Example sample dataset
data(ModelSet) # 58 models to fit
data(ParamRanges) # Parameter ranges for the 58 models
data(ParamSeeds) # 58 sets of candidate initial parameters
###########################################
#### Step 3. Truncate 58 models to 2 models (for quick demonstration purposes only) ####
###########################################
testmodels <- list()
testmeta <- list()
paramranges <- list()
Mods <- c(1,2)
for (i in 1:length(Mods))
{
testmodels <- c(testmodels , ModelSet[ Mods[i] ] )
testmeta [[i]] <- ParamSeeds [[ Mods[i] ]]
paramranges [[i]] <- ParamRanges [[ Mods[i] ]]
}
###########################################
#### Step 4 (simple version). Run the DiveMaster function ####
###########################################
################
# SIMPLEST METHOD, USING DiveMaster
################
# With two samples (Main sample + one subset)
result <- DiveMaster(models=testmodels, init.params=testmeta, param.ranges = paramranges,
main.samp=Bact1, subsizes=2, NResamples=50, nrf=10, fitloop=1, numit=50) # default parameters are: nrf=1, NResamples=1e3, numit=1e5. Values here chosen to speed example up
逐步运行，一步步分解运行它的命令
################
# LONGER METHOD, using component functions
################
#### Step 4 (component function version). ####
# Again, with two samples (Main sample + one subset)
###########################################
### 4.1 create rarefaction data (DivSubsamples object)
###########################################
Bact1length = sum(Bact1$Count)
dss_1 <- DivSubsamples(Bact1, nrf=2, minrarefac=1, maxrarefac=0.5*Bact1length , NResamples=30)
dss_2 <- DivSubsamples(Bact1, nrf=2, minrarefac=1, maxrarefac=Bact1length , NResamples=30)
dss <- list(dss_2, dss_1)
###########################################
### 4.2 fit models (create FitSingleMod object)
###########################################
fmm <- list() ## list of fitted model data
for (i in 1:length(Mods))
{
fsm.temp <- FitSingleMod(model.list = testmodels[i] , init.param = testmeta[[i]],
param.range = paramranges[[i]],
numit = 10^2,
varleft = 1e-8,
fitloops = 1,
minplaus = 10 ,
tot.pop = 100 * Bact1length,
nrf = 2,
minrarefac = 1,
NResamples = 30,
main.samp = Bact1,
# Option 1: provide previously calculated rareafaction data - keeps rarefaction data consistent across models
subsizes = c(Bact1length, 0.5*Bact1length),
data.default = FALSE,
dssamps = dss
## Option 2: calculate rarefaction data separately for each model - not recommended.
#subsizes = 2,
#data.default = TRUE
)
fmm[[fsm.temp$modelname]] <- fsm.temp
}
###########################################
### 4.3 score models (create list of class ScoreSingleMod)
###########################################
num.mod <- length(Mods)
ssm <- matrix(rep(NA, num.mod*4), nrow=num.mod, ncol=4) ## list of model scores
colnames(ssm) <- c("fit", "accuracy", "similarity", "plausibility")
mod.rownames <- matrix(rep(NA, num.mod), nrow = num.mod, ncol = 1)
mod.score <- matrix(rep(NA, num.mod), nrow = num.mod, ncol = 1)
TopX <- 5 ##
for (i in 1:length(Mods))
{
ssm.temp <- ScoreSingleMod(fsm = fmm[[i]], precision.lv = c(1e-04,
0.005, 0.005), plaus.pen = 500 )
mod.score.temp <- combine.criteria(ssm = ssm.temp, crit.wts = c(1, 1, 1, 1))
ssm[i, 1] <- ssm.temp$fit
ssm[i, 2] <- ssm.temp$accuracy
ssm[i, 3] <- ssm.temp$similarity
ssm[i, 4] <- ssm.temp$plausibility
mod.score[i, ] <- mod.score.temp
mod.rownames[i] <- fmm[[i]]$modelname
}
rownames(mod.score) <- mod.rownames
colnames(mod.score) <- "Combined score"
rownames(ssm) <- mod.rownames
###########################################
### 4.3 Produce Estimates
###########################################
m <- min(TopX, num.mod)
topX_scores <- sort(unique(mod.score))[1:m]
lenX <- length(topX_scores)
topX_index <- which(mod.score %in% topX_scores)
predX.vector <- c()
for (i in 1:lenX)
{
tmp <- ((fmm[[topX_index[i]]])$global)[1, ]
predX.vector <- c(predX.vector, tmp)
}
PointEstimate <- geo.mean(predX.vector)
UpperBound <- max(predX.vector)
LowerBound <- min(predX.vector)
###########################################
#### Step 5. View outputs ####
###########################################
result # Comparison of combined scores
summary(result) # Summary of combined score
result$ssm # Detailed comparison of scores
result$fmm$logistic # Fit details (model 1 - logistic model)
result$fmm$negexp # Fit details (model 2 - negative exponential model)
summary(result$fmm$logistic) # Fit summary (model 1)
result$fmm$logistic$param # Individual fit details ($param example)
plot(result$fmm$logistic) # Local plot of model 1 fit
plot(result$fmm$logistic, range="global") # Global plot of model 1 fit
plot(result$fmm$negexp) # Local plot of model 2 fit
plot(result$fmm$negexp, range="global") # Global plot of model 2 fit
PopDiversity(result, 10^6) ## calculate diversity at another population size.
###########################################
#### Step 6. Miscellaneous ####
###########################################
?DiveMaster # View main help file
?DivE # Package summary
# Component functions
?DivSubsamples
?FitSingleMod
?ScoreSingleMod
这是它写的示例，逐步运行的话，会有问题的,在倒数第二个循环那里会出现报错（暂时没有解决，估计我找的这个案例是之前的？以本人经验来看，一般是版本问题，估计cran上安装的是新一点的版本，不过我没有继续纠结这个，毕竟我只需要其中的最主要的函数而已）
Error in combine.criteria(ssm = ssm.temp, crit.wts = c(1, 1, 1, 1)) :
could not find function "combine.criteria"
DiveMaster(models, init.params, param.ranges, main.samp,
tot.pop=(100*(DivSampleNum(main.samp,2)[1])), numit=10^5,
varleft=1e-8, subsizes=6, dssamps=list(), nrf=1, minrarefac=1,
NResamples=1000, minplaus=10,
precision.lv=c(0.0001, 0.005, 0.005), plaus.pen=500,
crit.wts=c(1.0, 1.0, 1.0, 1.0), fitloops=2, numpred=5)
这个函数主要的参数
models
list of models; each model is written as a function: function(x, params) { with(as.list(params), <function of params>)}. Examples are given in the ModelSet data file as part of the DivE package.
init.params
list of matrices of initial seed model parameters. For each matrix, each row represents a given parameter set; each column represents a parameter value. Column names must match parameter names (params) in the corresponding model in the list models. Examples are given in the ParamSeeds data file as part of the DivE package.
param.ranges
list of matrices of lower and upper model parameters bounds. Used for the modFit function. The first and second row corresponds to the lower and upper bounds respectively; each column represents a parameter value. Column names must match parameter names (params) in the corresponding model in the list models. Examples are given in the param.ranges data file as part of the DivE package.
main.samp
the main sample, either as a 2-column data.frame (species ID, count of species), or a vector of species IDs.
tot.pop
total population (integer); default set to 100x the main.samp size.
numit
control argument passed to optimisation routine; the maximum number of iterations that modFit will perform. See modFit for details.
varleft
control argument passed to optimisation routine; see modFit for details.
subsizes
either number of nested subsamples (integer, must be 2 or greater), or a vector of nested sample lengths. If the former, then the vector of sample lengths will be created using the DivSampleNum function.
dssamps
list of user specified rarefaction data DivSubsamples objects. The length of each component vector of each object in the list must correspond to the vector of nested sample lengths (as defined by the user in subsizes).
nrf
difference between lengths of successive rarefaction datapoints.
minrarefac
minimum rarefaction x-axis value. This argument is not used if list of DivSubsamples object is specified in dssamps.
NResamples
number of resamples used to calculate the rarefaction data. This parameter is not used if list of DivSubsamples object is specified in dssamps.
minplaus
lower x-axis bound for plausibility check.
precision.lv
vector of precision level values for each criterion: 1. discrepancy – mean percentage error between rarefaction data points and model predicion, 2. Sample accuracy – percentage error between observed diversity of full rarefaction data and estimated diversity of full data from subsample, 3. local similarity. The scores for each criteria are defined as 1 + (multiples of bin sizes)
plaus.pen
penalty score for breaking the plausibility criterion: a model fit should be monotonically increasing and should have a slowing rate of species accumulation.
crit.wts
vector of weights of each of the four scoring criteria – fit, accuracy, similarity, plausibility. Default is c(1,1,1,1).
fitloops
number of fitting rounds performed for each model. In each round of fitting, the initial seed parameter values for each model will be the fitted parameters of the previous fitting run. This parameter has a significant impact on the computational time. The ‘sweet spot’ is 2.
numpred
number of topscoring models used for diversity prediction. Default is 5.
细节
这是 DivE estimator的主要功能。默认操作是四个步骤的组合。
- 从主样本生成嵌套样本长度列表。
- 对于每个嵌套子样本，生成稀疏数据向量及其相关的平均物种多样性。
- 为生成的数据拟合一组模型。
- 根据 DivE 多样性估计方法评估拟合并比较模型和拟合标准的分数。
可以使用 CombDM 函数将 DiveMaster 对象列表（每个对象表示适合不同模型集）组合成单个 DiveMaster 对象。当无法在单次运行中使用完整的 58 个模型集运行 DivE estimator时，这很有用。
可以使用 PopDiversity 函数估计给定种群的多样性，其中参数分别是 Divemaster 对象和种群大小。
