Calculation of repertoire metrics
VDJTOOLS的CalcDiversityStats
函数提供了直接计算功能,详情参考免疫组库分析实战/mixcr+vdjtools+R实现,使用vdjtools进行免疫组库分析.
本文以TRUST4
的结果为例(基本都大同小异),展示如何计算免疫组库的Richness
/Diversity
/Clonality
。
-
Diversity:香农-威纳指数(Shannon-WienerIndex ,H) 可以回答样本各物种均一性如何,一般选择以e为底的对数计算。
Shannon-WienerIndex
在香农-威纳指数中,包含着两个成分:①种数(Richness);②各种间个体分配的均匀性(Evenness)
-
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)
网友评论