Hilldiv: Github使用说明及Hill numbers在多样性分析中应用代码

但是是基于R包的帮助文档。Github上目前已经有了更详尽的说明。

https://github.com/anttonalberdi/hilldiv/wiki

本文是Hilldiv的最新用法说明。

0. Hill numbers

Richness

代码语言：javascript复制

#Load data
data(bat.diet.otutable)
#Compute richness for each sample (all three options yield the same result)
index.div(bat.diet.otutable)
index.div(bat.diet.otutable, index="richness")
colSums(bat.diet.otutable != 0)

Evenness

代码语言：javascript复制

#Load data
data(bat.diet.otutable)
#Compute Shannon index for each sample 
index.div(bat.diet.otutable, index="shannon")
#Compute Simpson index for each sample
index.div(bat.diet.otutable, index="simpson")

Regularity,衡量OTU之间相似的程度。Faith’s PD考虑richness, Allen’s H同时考虑richness和evenness

代码语言：javascript复制

# Computing diversities that account for phylogeny-based regularity require an OTU phylogenetic tree.
data(bat.diet.tree)
#Compute Faith's PD for each sample 
index.div(bat.diet.otutable, tree=bat.diet.tree, index="faith")
#Compute Allen's H for each sample 
index.div(bat.diet.otutable, tree=bat.diet.tree, index="allen")
#Compute Rao's Q for each sample 
index.div(bat.diet.otutable, tree=bat.diet.tree, index="rao")

Phylogenetic indices and Hill numbers系统发育指数和Hill之间可以相互转化

代码语言：javascript复制

#Load data
data(bat.diet.otutable)
#Subset the first sample and add OTU names
bat.diet.sample <- bat.diet.otutable[,1]
names(bat.diet.sample) <- rownames(bat.diet.otutable)

#Diversity indices
richness <- index.div(bat.diet.sample, index ="richness")
shannon <- index.div(bat.diet.sample, index ="shannon")
simpson <- index.div(bat.diet.sample, index ="simpson")

faith <- index.div(bat.diet.sample, tree=bat.diet.tree, index="faith")
allen <- index.div(bat.diet.sample, tree=bat.diet.tree, index="allen")
rao <- index.div(bat.diet.sample, tree=bat.diet.tree, index="rao")

#Hill numbers
hill.div.q0 <- hill.div(bat.diet.sample, qvalue=0)
hill.div.q1 <- hill.div(bat.diet.sample, qvalue=1)
hill.div.q2 <- hill.div(bat.diet.sample, qvalue=2)
hill.phylodiv.q0 <- hill.div(bat.diet.sample, tree=bat.diet.tree, qvalue=0)
hill.phylodiv.q1 <- hill.div(bat.diet.sample, tree=bat.diet.tree, qvalue=1)
hill.phylodiv.q2 <- hill.div(bat.diet.sample, tree=bat.diet.tree, qvalue=2)

#Richness vs. neutral hill number of q=0
richness
hill.div.q0

#Shannon index vs. neutral hill number of q=1
shannon
exp(shannon)
hill.div.q1

#Simpson index vs. neutral hill number of q=2
simpson
1/(1-simpson)
hill.div.q2

#Compute phylogenetic tree depth to transform indices to Hill numbers
T <- tree.depth(bat.diet.tree,bat.diet.sample)

#Faith's PD vs. phylogenetic hill number of q=0
faith
faith/T
hill.phylodiv.q0

#Allen's H vs. phylogenetic hill number of q=1
allen
exp(allen/T)
hill.phylodiv.q1

#Rao's Q vs. phylogenetic hill number of q=2
rao
1/(1-rao/T)
hill.phylodiv.q2

1. Data files

读入数据

代码语言：javascript复制

#Load data
data(bat.diet.otutable)
data(bat.diet.tree)
data(bat.diet.hierarchy)

2. Data preprocessing

Filter samples by sequencing depth

代码语言：javascript复制

#Remove samples with less than 20000 sequence reads (the threshold employed here is arbitrary)
bat.diet.otutable.filt <- depth.filt(bat.diet.otutable, 20000)

Filter OTUs by sequence copy number

代码语言：javascript复制

#Remove all OTUs with less than 10 copies
copy.filt(bat.diet.otutable, 10)
#Remove all OTUs with less copies than 0.01% of the total number of reads of each sample
copy.filt(bat.diet.otutable, 0.0001)

Convert OTU table into relative abundances

代码语言：javascript复制

#Convert into relative
bat.diet.otutable.rel <- tss(bat.diet.otutable)

Convert abundance table into incidence table

代码语言：javascript复制

#Incidence data for the overall system
to.incidence(bat.diet.otutable)
#Incidence data for each group
to.incidence(bat.diet.otutable,hierarchy=bat.diet.hierarchy)
#Output relative rather than absolute values 
to.incidence(bat.diet.otutable,hierarchy=bat.diet.hierarchy,relative=TRUE)
#Which is the same as
tss(to.incidence(bat.diet.otutable,hierarchy=bat.diet.hierarchy))

Match OTU table and tree

代码语言：javascript复制

data(bat.diet.otutable)
data(bat.diet.tree)
#Check if OTU table and OTU tree match
match.data(bat.diet.otutable,bat.diet.tree)

#Modify OTU table (retain>bat.diet.otutable2 <- bat.diet.otutable[c(1:300),]
#Check if OTU table and OTU tree match
match.data(bat.diet.otutable2,bat.diet.tree)
#Modify OTU tree (remove OTU1 and OTU2) to test the functionality of match.data()
bat.diet.tree2 <- drop.tip(bat.diet.tree,c("OTU1","OTU2"))
#Check if OTU table and OTU tree match
match.data(bat.diet.otutable,bat.diet.tree2)

#Check if OTU table and OTU tree match
match.data(bat.diet.otutable2,bat.diet.tree2)
#Filter OTU table
bat.diet.otutable3 <- match.data(bat.diet.otutable2,bat.diet.tree2,output="otutable")
#Filter OTU tree
bat.diet.tree3 <- match.data(bat.diet.otutable2,bat.diet.tree2,output="tree")
#Check if OTU table and OTU tree match
match.data(bat.diet.otutable3,bat.diet.tree3)

3.1 Diversity computation of a single system

Diversity computation of a single sample vs. multiple individual samples vs. overall system

代码语言：javascript复制

#Load data
data(bat.diet.otutable)

#Hill numbers of a single sample
#Select the first samples of the count table
bat.diet.sample <- bat.diet.otutable[,1]
#Compute Hill number of q=1 (Shannon diversity) using hill.div()
hill.div(bat.diet.sample,qvalue=1)

#Hill numbers of multiple samples
hill.div(bat.diet.otutable,qvalue=1)

#incidence数据
hill.div(tss(bat.diet.otutable),qvalue=1)

#Compute the Hill number of q=1 (Shannon diversity) of the overall system
gamma.div(bat.diet.otutable,qvalue=1)
#Compute the Hill number of q=0 (richness) of the overall system, which is exactly the number of OTU/ASVs in the system
gamma.div(bat.diet.otutable,qvalue=0)

Hill numbers vs diversity indices。Hill和多样性可以相互转化

代码语言：javascript复制

#Load data
data(bat.diet.otutable)

#NEUTRAL MEASURES

#Richness (identical for the index or the Hill number)
index.div(bat.diet.otutable,index="richness") 
index.div(bat.diet.otutable) # also works without using the index argument, as richness is the default option.
hill.div(bat.diet.otutable,qvalue=0)

#Neutral Hill number of q=1 (Shannon diversity) vs Shannon index
hill.div(bat.diet.otutable,qvalue=1)
index.div(bat.diet.otutable,index="shannon")
#To transform the Shannon index to Hill number of q=1 (Shannon diversity)
exp(index.div(bat.diet.otutable,index="shannon"))

#Neutral Hill number of q=2 (Simpson diversity) vs Simpson index
hill.div(bat.diet.otutable,qvalue=2)
index.div(bat.diet.otutable,index="simpson")
#To transform the Shannon index to Hill number of q=2 (Shannon diversity)
1/(1-index.div(bat.diet.otutable,index="simpson"))

#PHYLOGENETIC MEASURES
data(bat.diet.otutable)
data(bat.diet.tree)
T <- tree.depth(bat.diet.tree,bat.diet.otutable)

#Phylogenetic Hill number of q=0 (phylorichness) vs Faith's PD
index.div(bat.diet.otutable,tree=bat.diet.tree,index="faith") 
hill.div(bat.diet.otutable,qvalue=0,tree=bat.diet.tree)
#To transform the Faith's PD to phylogenetic Hill number of q=0 (phylorichness)
index.div(bat.diet.otutable,tree=bat.diet.tree,index="faith")/T

#Phylogenetic Hill number of q=1 vs Allens'H
index.div(bat.diet.otutable,tree=bat.diet.tree,index="allen") 
hill.div(bat.diet.otutable,qvalue=1,tree=bat.diet.tree)
#To transform the Allen's H to phylogenetic Hill number of q=1
exp(index.div(bat.diet.otutable,tree=bat.diet.tree,index="allen")/T)

#Phylogenetic Hill number of q=2 vs Rao's Q
index.div(bat.diet.otutable,tree=bat.diet.tree,index="rao") 
hill.div(bat.diet.otutable,qvalue=2,tree=bat.diet.tree)
#To transform the Allen's H to phylogenetic Hill number of q=1
1/(1-index.div(bat.diet.otutable,tree=bat.diet.tree,index="rao")/T)

3.2 Diversity computation and comparison of multiple systems

Computing and contrasting (sub)system mean diversities manually

代码语言：javascript复制

#Species 1 (MSc) can be compared to species 2 (Mmy) and species 3 (Rme) by subsetting the corresponding columns and averaging resulting Hill numbers.
MSc <- hill.div(bat.diet.otutable[,c(1:5)],qvalue=1)
MMy <- hill.div(bat.diet.otutable[,c(6:10)],qvalue=1)
RMe <- hill.div(bat.diet.otutable[,c(11:15)],qvalue=1)

##统计检验
#Two contrasting groups
t.test(MSc,MMy)
wilcox.test(MSc,MMy)

#Three contrasting groups (requires generation of a table containing values and grouping variable)
all.divs <- c(MSc,MMy,RMe)
all.divs <- as.data.frame(cbind(all.divs,c(rep("MSc",5),rep("MMy",5),rep("RMe",5))))
all.divs[,1] <- as.numeric(as.character(all.divs[,1]))
colnames(all.divs) <- c("diversity","species")

summary(aov(diversity ~ species, data = all.divs))
kruskal.test(diversity ~ species, data = all.divs)

Computing and contrasting (sub)system mean diversities using div.test()

代码语言：javascript复制

#Test whether the diversity values differ across all contrasting groups
div.test(bat.diet.otutable,qvalue=1,hierarchy=bat.diet.hierarchy)
#Test whether the diversity values differ across the first two contrasting groups (first 10 lines in hierarchy table)
div.test(bat.diet.otutable,qvalue=1,hierarchy=bat.diet.hierarchy[c(1:10),])

Post hoc analyses

代码语言：javascript复制

#To run internal post hoc analysis
div.test(bat.diet.otutable,qvalue=1,hierarchy=bat.diet.hierarchy,posthoc=TRUE)
#To run external post hoc analysis (e.g. Dunn test with Bonferroni correction using package FSA)
library(FSA)
div.test.result <- div.test(bat.diet.otutable,qvalue=1,hierarchy=bat.diet.hierarchy)
FSA::dunnTest(Value ~ Group, data=div.test.result$data, method="bonferroni")

Visualising results using div.test.plot()

代码语言：javascript复制

#Run and save div.test() results
divtest.noposthoc<- div.test(bat.diet.otutable,qvalue=1,hierarchy=bat.diet.hierarchy)
divtest.posthoc<- div.test(bat.diet.otutable,qvalue=1,hierarchy=bat.diet.hierarchy,posthoc=TRUE)

#Plot without and with posthoc data
div.test.plot(divtest.noposthoc)
div.test.plot(divtest.noposthoc, posthoc=TRUE) #returns an error because the object does not contain posthoc data
div.test.plot(divtest.posthoc, posthoc=TRUE)
#With custom threshold
div.test.plot(divtest.posthoc, posthoc=TRUE, threshold=0.5)
#With custom colors
div.test.plot(divtest.posthoc, colour=c("#35a849","#9d1923","#f7ab1b","#ed7125","#cc4323","#b6d134","#fcee21","#085ba7"), posthoc=TRUE, threshold=0.5)

3.3 Diversity profiles 多样性和Hill阶数的关系

代码语言：javascript复制

#Load data
data(bat.diet.otutable)
#Diversity profile of a single sample (default profile)
div.profile(bat.diet.otutable[,1])
#Diversity profile of a single sample (custom profile from 0 to 3)
div.profile(bat.diet.otutable[,1],qvalues=c(seq(from = 0, to = 3, by = (0.1))))

#Diversity profile of multiple samples
div.profile(bat.diet.otutable)
#Diversity profile of alpha diversity of multiple samples (entire count table). 
div.profile(bat.diet.otutable,level="alpha")
#Diversity profile of gamma diversity of multiple samples (entire count table)
div.profile(bat.diet.otutable,level="gamma")

#Diversity profiles of groups determined by hierarchy tables
data(bat.diet.hierarchy)
#Alpha diversity
div.profile(bat.diet.otutable,hierarchy=bat.diet.hierarchy,level="alpha")
#Gamma diversity
div.profile(bat.diet.otutable,hierarchy=bat.diet.hierarchy,level="gamma")

4. Diversity partitioning 多样性分解

2个样本

代码语言：javascript复制

#ModelA
S1 <- c(OTU1=1,OTU2=1,OTU3=1,OTU4=1,OTU5=0,OTU6=0,OTU7=0,OTU8=0)
S2 <- c(OTU1=0,OTU2=0,OTU3=0,OTU4=0,OTU5=1,OTU6=1,OTU7=1,OTU8=1)
modelA <- cbind(S1,S2)
div.part(modelA,qvalue=0)

#ModelB
S1 <- c(OTU1=1,OTU2=1,OTU3=1,OTU4=1,OTU5=0,OTU6=0,OTU7=0,OTU8=0)
S2 <- c(OTU1=0,OTU2=0,OTU3=1,OTU4=0,OTU5=0,OTU6=1,OTU7=1,OTU8=1)
modelB <- cbind(S1,S2)
div.part(modelB,qvalue=0)

#ModelC
S1 <- c(OTU1=1,OTU2=1,OTU3=1,OTU4=1,OTU5=0,OTU6=0,OTU7=0,OTU8=0)
S2 <- c(OTU1=0,OTU2=1,OTU3=0,OTU4=1,OTU5=0,OTU6=1,OTU7=0,OTU8=1)
modelC <- cbind(S1,S2)
div.part(modelC,qvalue=0)

#ModelD
S1 <- c(OTU1=1,OTU2=1,OTU3=1,OTU4=1,OTU5=0,OTU6=0,OTU7=0,OTU8=0)
S2 <- c(OTU1=1,OTU2=1,OTU3=1,OTU4=1,OTU5=0,OTU6=0,OTU7=0,OTU8=0)
modelD <- cbind(S1,S2)
div.part(modelD,qvalue=0)

3个样本

代码语言：javascript复制

#ModelA
S1 <- c(OTU1=1,OTU2=1,OTU3=1,OTU4=1,OTU5=0,OTU6=0,OTU7=0,OTU8=0,OTU9=0,OTU10=0,OTU11=0,OTU12=0)
S2 <- c(OTU1=0,OTU2=0,OTU3=0,OTU4=0,OTU5=1,OTU6=1,OTU7=1,OTU8=1,OTU9=0,OTU10=0,OTU11=0,OTU12=0)
S3 <- c(OTU1=0,OTU2=0,OTU3=0,OTU4=0,OTU5=0,OTU6=0,OTU7=0,OTU8=0,OTU9=1,OTU10=1,OTU11=1,OTU12=1)
modelA <- cbind(S1,S2,S3)
div.part(modelA,qvalue=0)

#ModelB
S1 <- c(OTU1=1,OTU2=1,OTU3=1,OTU4=1,OTU5=0,OTU6=0,OTU7=0,OTU8=0,OTU9=0,OTU10=0,OTU11=0,OTU12=0)
S2 <- c(OTU1=0,OTU2=0,OTU3=0,OTU4=1,OTU5=0,OTU6=1,OTU7=1,OTU8=1,OTU9=0,OTU10=0,OTU11=0,OTU12=0)
S3 <- c(OTU1=0,OTU2=0,OTU3=0,OTU4=0,OTU5=0,OTU6=0,OTU7=1,OTU8=0,OTU9=1,OTU10=0,OTU11=1,OTU12=1)
modelB <- cbind(S1,S2,S3)
div.part(modelB,qvalue=0)

#ModelC
S1 <- c(OTU1=1,OTU2=1,OTU3=1,OTU4=1,OTU5=0,OTU6=0,OTU7=0,OTU8=0,OTU9=0,OTU10=0,OTU11=0,OTU12=0)
S2 <- c(OTU1=0,OTU2=1,OTU3=0,OTU4=1,OTU5=0,OTU6=1,OTU7=0,OTU8=1,OTU9=0,OTU10=0,OTU11=0,OTU12=0)
S3 <- c(OTU1=0,OTU2=0,OTU3=1,OTU4=1,OTU5=0,OTU6=0,OTU7=0,OTU8=0,OTU9=0,OTU10=0,OTU11=1,OTU12=1)
modelC <- cbind(S1,S2,S3)
div.part(modelC,qvalue=0)

#ModelD
S1 <- c(OTU1=1,OTU2=1,OTU3=1,OTU4=1,OTU5=0,OTU6=0,OTU7=0,OTU8=0,OTU9=0,OTU10=0,OTU11=0,OTU12=0)
S2 <- c(OTU1=1,OTU2=1,OTU3=1,OTU4=1,OTU5=0,OTU6=0,OTU7=0,OTU8=0,OTU9=0,OTU10=0,OTU11=0,OTU12=0)
S3 <- c(OTU1=1,OTU2=1,OTU3=1,OTU4=1,OTU5=0,OTU6=0,OTU7=0,OTU8=0,OTU9=0,OTU10=0,OTU11=0,OTU12=0)
modelD <- cbind(S1,S2,S3)
div.part(modelD,qvalue=0)

END

data mean numbers table using

0 人点赞