多样性计算

Reps 包括gini, shannon.entropy, metric.entropy, simpson, d50, dxx
immunarch包括 chao1, hill, diversity, Gini-Simpson, inv.simp, gini, raref, d50, dxx.

GINI系数

基尼系数是本来是一个国际通用的经济学概念，用来衡量贫富差距。基尼系数介于0-1之间，基尼系数越大，表示不平等程度越高。

基尼系数最大为1，表示居民之间的收入分配绝对不平均，即100%的收入被一个单位的人全部占有了；
基尼系数最小为0，表示居民之间的收入分配绝对平均，即人与人之间收入完全平等，没有任容何差异。
假定一定数量的人口按收入由低到高顺序排队，分为人数相等的n组，从第1组到第i组人口累计收入占全部人口总收入的比重

image.png

参考immunarch

# Reps 中的代码
gini_coef <- function(.data, .do.norm = NA, .laplace = 0) {
  .data <- sort(check_distribution(.data, .do.norm, .laplace, .warn.sum = FALSE))
  n <- length(.data)
  res <- 1 / n * (n + 1 - 2 * sum((n + 1 - 1:n) * .data) / sum(.data))
  add_class(res, "immunr_gini")
}
# 等同于
gini.index  <-function(x){ 
  x <- sort(x)  
  G <- sum(x * 1L:length(x)) 
  G <- 2 * G/sum(x) - (length(x) + 1L)
  G/length(x)
}

D50

CDR3序列（元素）按照占比比例排序，从最高往最低累加，达到50%的总序列时候的CDR3序列种类比例，就是类似于饱和度曲线。
D50最大为0.5，意味着全部的CDR3序列占比一致，多样性好；
D50最小为0，意味着有且只有一种CDR3序列，多样性差。

# Reps 中的代码
dXX <- function(.data, .perc = 10) {
  if (has_class(.data, "list")) {
    return(t(sapply(.data, dXX, .perc = .perc)))
  }
  prop <- 0
  n <- 0
  col <- .data
  col <- sort(col, decreasing = TRUE)
  col.sum <- sum(col)
  while (prop < col.sum * (.perc / 100)) {
    n <- n + 1
    prop <- prop + col[n]
  }
  res <- c(Clones = n, Percentage = 100 * signif((prop / col.sum), 3), Dxx_DiversityIndex = n*100 / nrow(x))
  add_class(res, "immunr_dxx")
}

# 类似的
d50.index <-function(x,type='raw'){
  if(type=='raw'){
    myfreqs <- table(x)/length(x) 
    myvec <- sort(as.numeric(myfreqs),decreasing = T)
  }else{
    myvec=sort(x,decreasing = T)
  }
  len=length(myvec)
  state=cumsum(myvec)>sum(myvec)/2  
  (len -sum(state))/len
}

Shannon

参考 https://zhuanlan.zhihu.com/p/157346850

# Reps 中的代码
shannon.entropy <-function(x,type='raw'){
  if(type=='raw'){
    x <- x$copy
    myfreqs <- table(x)/length(x)
    myvec <- as.data.frame(myfreqs)[,2]
  }else{
    myvec=x
  }
  -sum(myvec * log2(myvec))
}

metric entropy

## modify shannon.entropy to metric entropy
metric.entropy <-function(x,type='raw'){
  if(type=='raw'){
    x <- x$copy
    myfreqs <- table(x)/length(x)
    myvec <- as.data.frame(myfreqs)[,2]
  }else{
    myvec=x
  }
  -sum(myvec * log(myvec,length(x)))
}

Simpson 指数

Vegan 包里的相关计算公式:

vegan.png

百度的说明:

Simpson baidu.png

# Reps 中的代码
get.simpson <-function(x,type='raw'){
  if(type=='raw'){
    x <- x$copy
    myfreqs <- table(x)/length(x)
    myvec <- as.data.frame(myfreqs)[,2]
  }else{
    myvec=x
  }
 # -sum(myvec * log(myvec,length(x)))
  1-sum( myvec ^2)
}

# 类似与
gini_simpson <- function(.data, .do.norm = NA, .laplace = 0) {
  .data <- check_distribution(.data, .do.norm, .laplace)
  res <- 1 - sum(.data^2)
  add_class(res, "immunr_ginisimp")
}

inverse_simpson 指数 将增加

Inv-Simpson

https://mothur.org/wiki/invsimpson/

Inv-Simpson.png

inverse_simpson <- function(.data, .do.norm = NA, .laplace = 0) {
  .data <- check_distribution(.data, .do.norm, .laplace)
  res <- 1 / sum(.data^2)
  add_class(res, "immunr_invsimp")
}

Diversity 将增加

参考immunarch

diversity_eco <- function(.data, .q = 5, .do.norm = NA, .laplace = 0) {
  .data <- check_distribution(.data, .do.norm = NA, .laplace = 0)
  if (.q == 0) {
    res <- length(.data)
  } else if (.q == 1) {
    res <- exp(-sum(.data * log(.data)))
  } else if (.q > 1) {
    res <- 1 / (sum(.data^.q)^(1 / (.q - 1)))
  } else {
    res <- NA
  }
}

hill 将增加

参考immunarch

hill_numbers <- function(.data, .max.q = 6, .min.q = 1,
                         .do.norm = NA, .laplace = 0) {
  .data <- check_distribution(.data, .do.norm = .do.norm, .laplace)
  if (.min.q < 0) {
    .min.q <- 0
  }
  res <- c()
  for (q in .min.q:.max.q) {
    res <- c(res, diversity_eco(.data, q))
  }
  names(res) <- paste0("q", .min.q:.max.q)
  add_class(res, "immunr_hill")
}

chao1 将增加

参考immunarch

chao1 <- function(.data) {
  counts <- table(.data)
  e <- NA
  v <- NA
  lo <- NA
  hi <- NA
  n <- sum(.data)
  D <- length(.data)
  f1 <- counts["1"]
  f2 <- counts["2"]
  # f1 == 0 && f2 == 0
  if (is.na(f1) && is.na(f2)) {
    e <- D
    i <- unique(.data)
    v <- sum(sapply(i, function(i) sum(.data == i) * (exp(-i) - exp(-2 * i)))) - (sum(sapply(i, function(i) i * exp(-i) * sum(.data == i))))^2 / n
    P <- sum(sapply(i, function(i) sum(.data == i) * exp(-i) / D))
    lo <- max(D, D / (1 - P) - qnorm(1 - .05 / 2) * sqrt(v) / (1 - P))
    hi <- D / (1 - P) + qnorm(1 - .05 / 2) * sqrt(v) / (1 - P)
  }
  # f1 != 0 && f2 == 0
  else if (is.na(f2)) {
    e <- D + f1 * (f1 - 1) / 2 * (n - 1) / n
    v <- (n - 1) / n * f1 * (f1 - 1) / 2 + ((n - 1) / n)^2 * f1 * (2 * f1 - 1)^2 / 4 - ((n - 1) / n)^2 * f1^4 / 4 / e
    t <- e - D
    K <- exp(qnorm(1 - .05 / 2) * sqrt(log(1 + v / t^2)))
    lo <- D + t / K
    hi <- D + t * K
  }
  # f1 != && f2 != 0
  else {
    const <- (n - 1) / n
    e <- D + f1^2 / (2 * f2) * const
    f12 <- f1 / f2
    v <- f2 * (const * f12^2 / 2 + const^2 * f12^3 + const^2 * f12^4 / 4)
    t <- e - D
    K <- exp(qnorm(.975) * sqrt(log(1 + v / t^2)))
    lo <- D + t / K
    hi <- D + t * K
  }
  res <- c(Estimator = e, SD = sqrt(v), Conf.95.lo = lo, Conf.95.hi = hi)
  add_class(res, "immunr_chao1")
}