Calculation of repertoire metrics

VDJTOOLS的CalcDiversityStats函数提供了直接计算功能，详情参考免疫组库分析实战/mixcr+vdjtools+R实现，使用vdjtools进行免疫组库分析.
本文以TRUST4的结果为例（基本都大同小异），展示如何计算免疫组库的Richness/Diversity/Clonality。

trust4.report

• Diversity：香农-威纳指数（Shannon-WienerIndex ，H） 可以回答样本各物种均一性如何，一般选择以e为底的对数计算。
  在香农-威纳指数中，包含着两个成分：①种数（Richness）；②各种间个体分配的均匀性（Evenness）
  

  Shannon-WienerIndex

• Richness : 每个样本中的Number of Unique clone。

• Evenness：各种之间，个体分配越均匀，H值就越大。如果每一个体都属于不同的种，多样性指数就最大；如果每一个体都属于同一种，则其多样性指数就最小。可以通过估计群落的理论上的最大多样性指数（Hmax），然后以实际的多样性指数对Hmax的比率，从而获得均匀性指数，即为标准化的香农指数。
  Evenness=H/Hmax，当每一个Clone出现的概率相同时即可得到Hmax，此时，P=1/N，Hmax=-1/N*log(1/N)*N=-log(1/N)=log(N), Evenness=Diversity/log(Richness) 。

• Clonality: Clonality = 1 – Evenness = 1 - Diversity/log(Richness)。
  参考：Shannon–Wiener Index

代码实现（R）

Shannon_index <- function(list_p){
  sum=0
  for (p in list_p){
    sum <- p*log(p) + sum
  }
  H <- -sum
  return(H)
}
clonality <- function(diversity,richness){
  E<- diversity/log(richness)
  C <- 1-E
  return(C)
}

sample_id <- c(paste0("sample",1:10))
results <- data.frame(sample_id)

for (i in 1:length(sample_id)){
  df <- read.csv(file = paste0(sample_id[i],".txt"),
                    stringsAsFactors = F,sep = "\t")
  list_p <- df[,"frequency"]
  results[i,"richness"] <- nrow(df)
  results[i,"diversity"] <- Shannon_index(list_p)
  results[i,"clonality"] <- 1-clonality(Shannon_index(list_p),nrow(df))
}
write.table(results,filename,sep = "\t",row.names = F,quote = F)