Diversity

Diversity in ecological studies refers to different aspects of variation within and across communities:

Alpha diversity: The richness and evenness of taxa (or “features”) within a single sample or site. It measures how many different taxa are present and how evenly they are distributed.
Beta diversity: The degree of difference in community composition between two or more samples. It captures how similar or dissimilar communities are across environments.
Gamma diversity: The overall diversity across a collection of samples or an entire region, considering both the number of taxa and their distribution at a broader spatial scale.

Naive-, phylogenetic-, or functional diversity
One way to calculate diversity is to consider all taxa as separate entities and not care about their relationships. We call these diversity values naive.
A second way to calculate diversity takes the relationships in a phylogenetic tree into account. Two taxa that share many branches in the tree contribute less to diversity than two that are very far from each other in the tree. We call the diversity values that are calculated with a tree as input phylogenetic.
A third way to calculate diversity takes the pairwise distances between taxa into account. Two taxa that have a small pairwise distance contribute less to diversity than two that have a high pairwise distance. We call the diversity values that are calculated with a distance matrix as input functional.

Naive alpha diversity is simply the effective number of taxa, phylogenetic alpha diversity is the effective branch diversity in the phylogenetic tree associated with the taxa in the sample, and functional alpha diversity is the effective total distance between taxa in the sample. The word “effective” used above refers to the relative abundance weighting by the diversity order, q.

Hill numbers (or effective numbers)

Imagine a community with 10 species, all equally abundant. The most intuitive measure of diversity here is 10, and Hill numbers give exactly that result. Now consider a community with 2 dominant species and 8 rare ones. What should the diversity be?

If we ignore relative abundances, the answer is 10.
If we ignore rare species, the answer is 2.

Hill numbers solve this by introducing a single parameter, the diversity order (q), which controls how much weight is given to species’ relative abundances:

When q=0, all species count equally (richness).
When q=1, species are weighted by their proportional abundance (Shannon diversity).
When q=2, common species dominate the measure (Simpson diversity).

This flexibility makes Hill numbers a powerful and intuitive framework for comparing communities. And it doesn’t stop at 0, 1, or 2 — diversity can be explored across a continuum of q values, offering a nuanced view of community structure.

Alpha diversity

Let’s calculate alpha diversity using qdiv. First, we will load and rarefy some example data.

[1]:

import qdiv
obj = qdiv.MicrobiomeData.load_example("Saheb-Alam_DADA2")
rarefied_obj = obj.rarefy() #I rarefy the object before calculating alpha diversity

Next, we will calculate the naive alpha diversity values for the samples using the diversity.naive_alpha function.

[2]:

alpha = qdiv.diversity.naive_alpha(rarefied_obj, q=1)
print(alpha)

S4      4.905520
S5      8.761457
S6      5.746725
S7      4.922199
S10    46.192626
S11    40.920506
S12    89.238965
S13    54.704046
S20     7.215461
S21     6.875755
S22     4.723119
S23     6.007144
S26    47.436514
S27    48.542889
S28    39.325907
S29    59.500728
dtype: float64

This was the alpha diversity for diversity order q=1. Perhaps we want to calculate it for multiple q values and save it into a csv file together with the meta data.

[3]:

meta = rarefied_obj.meta.copy() #We make a copy of the metadata to avoid modifying the object
for q in [0, 1, 2]: #A for loop to calculate alpha diversity for several values of q
    alpha = qdiv.diversity.naive_alpha(rarefied_obj, q=q)
    meta["Div_q"+str(q)] = alpha #Add the calculated values to the meta dataframe
meta.to_csv("Metadata_with_alpha_diversity.csv") #This lines saves the dataframe as a CSV file.
print(meta) #This prints the file content to show it below

    location     feed mfc  Div_q0     Div_q1     Div_q2
S4     anode  acetate   B    74.0   4.905520   1.863492
S5     anode  acetate   B   136.0   8.761457   2.415367
S6     anode  acetate   B    78.0   5.746725   2.040309
S7     anode  acetate   B    74.0   4.922199   1.844776
S10  cathode  acetate   B   165.0  46.192626  13.296620
S11  cathode  acetate   B   196.0  40.920506  10.276571
S12  cathode  acetate   B   314.0  89.238965  31.613072
S13  cathode  acetate   B   188.0  54.704046  21.043678
S20    anode  glucose   D   153.0   7.215461   3.295074
S21    anode  glucose   D   110.0   6.875755   3.368090
S22    anode  glucose   D    57.0   4.723119   2.681520
S23    anode  glucose   D   110.0   6.007144   3.084631
S26  cathode  glucose   D   165.0  47.436514  16.851130
S27  cathode  glucose   D   218.0  48.542889  15.017341
S28  cathode  glucose   D   144.0  39.325907  13.886542
S29  cathode  glucose   D   241.0  59.500728  19.756538

The results seem to indicate higher diversity on cathode than anodes, and higher diversity in glucose-fed than acetate-fed microbial fuel cells. Higher diversity orders (q values) result in lower diversity. The change in diversity as q increases is an indication of the evenness. We can calculate evenness parameters using the function below.

[4]:

evenness = qdiv.diversity.evenness(rarefied_obj, index="pielou")
meta["Evenness"] = evenness
print(meta)

    location     feed mfc  Div_q0     Div_q1     Div_q2  Evenness
S4     anode  acetate   B    74.0   4.905520   1.863492  0.369502
S5     anode  acetate   B   136.0   8.761457   2.415367  0.441790
S6     anode  acetate   B    78.0   5.746725   2.040309  0.401365
S7     anode  acetate   B    74.0   4.922199   1.844776  0.370291
S10  cathode  acetate   B   165.0  46.192626  13.296620  0.750658
S11  cathode  acetate   B   196.0  40.920506  10.276571  0.703212
S12  cathode  acetate   B   314.0  89.238965  31.613072  0.781181
S13  cathode  acetate   B   188.0  54.704046  21.043678  0.764247
S20    anode  glucose   D   153.0   7.215461   3.295074  0.392854
S21    anode  glucose   D   110.0   6.875755   3.368090  0.410171
S22    anode  glucose   D    57.0   4.723119   2.681520  0.383985
S23    anode  glucose   D   110.0   6.007144   3.084631  0.381440
S26  cathode  glucose   D   165.0  47.436514  16.851130  0.755862
S27  cathode  glucose   D   218.0  48.542889  15.017341  0.721042
S28  cathode  glucose   D   144.0  39.325907  13.886542  0.738837
S29  cathode  glucose   D   241.0  59.500728  19.756538  0.744966

The cathode samples have, indeed, higher evenness than the anode samples.

Phylogenetic diversity

Since a phylogenetic tree is included in the example data object, we can also calculate phylogenetic diversity using the diversity.phyl_alpha function.

[5]:

for q in [0, 1, 2]: #A for loop to calculate alpha diversity for several values of q
    alpha = qdiv.diversity.phyl_alpha(rarefied_obj, q=q)
    meta["PD_q"+str(q)] = alpha #Add the calculated values to the meta dataframe
print(meta) #Show the results

    location     feed mfc  Div_q0     Div_q1     Div_q2  Evenness      PD_q0  \
S4     anode  acetate   B    74.0   4.905520   1.863492  0.369502  18.234736
S5     anode  acetate   B   136.0   8.761457   2.415367  0.441790  26.729741
S6     anode  acetate   B    78.0   5.746725   2.040309  0.401365  21.271491
S7     anode  acetate   B    74.0   4.922199   1.844776  0.370291  20.367709
S10  cathode  acetate   B   165.0  46.192626  13.296620  0.750658  30.481421
S11  cathode  acetate   B   196.0  40.920506  10.276571  0.703212  35.071286
S12  cathode  acetate   B   314.0  89.238965  31.613072  0.781181  46.723518
S13  cathode  acetate   B   188.0  54.704046  21.043678  0.764247  35.070279
S20    anode  glucose   D   153.0   7.215461   3.295074  0.392854  35.169738
S21    anode  glucose   D   110.0   6.875755   3.368090  0.410171  27.664994
S22    anode  glucose   D    57.0   4.723119   2.681520  0.383985  19.681168
S23    anode  glucose   D   110.0   6.007144   3.084631  0.381440  27.988865
S26  cathode  glucose   D   165.0  47.436514  16.851130  0.755862  36.015415
S27  cathode  glucose   D   218.0  48.542889  15.017341  0.721042  41.776060
S28  cathode  glucose   D   144.0  39.325907  13.886542  0.738837  28.993524
S29  cathode  glucose   D   241.0  59.500728  19.756538  0.744966  44.483627

        PD_q1     PD_q2
S4   1.945577  1.264527
S5   2.338791  1.368980
S6   2.082046  1.295584
S7   2.021452  1.271060
S10  3.117366  1.581757
S11  2.996559  1.543676
S12  3.614873  1.661804
S13  3.360075  1.624280
S20  2.133923  1.365343
S21  2.009238  1.350702
S22  1.840235  1.304523
S23  1.850619  1.308257
S26  4.067021  1.800328
S27  3.805963  1.741100
S28  3.465882  1.691958
S29  3.512326  1.641231

Functional diversity

To calculate functional diversity, we need a distance matrix with pairwise functional distances between features in the dataset. The distance matrix should be in the form of a pandas dataframe with feature names as row- and column indices. For the purpose of demonstrating the use of the function, we will calculate such as distance matrix based on the phylogenetic distances between features. We use the sequences.tree_distance_matrix function to get a distance matrix saved in the variable distmat.

[6]:

distmat = qdiv.sequences.tree_distance_matrix(rarefied_obj, save=False) #By setting save=False, no file is saved

Leaves: 100%|██████████████████████████████████████████████████████████████████████| 671/671 [00:00<00:00, 1435.52it/s]

Next, we call the diversity.func_alpha function to calculate functional alpha diversity.

[7]:

alpha = qdiv.diversity.func_alpha(rarefied_obj, distmat, q=1)
print(alpha)

S4      107.693085
S5      268.177859
S6      127.536179
S7      121.187460
S10    1624.876208
S11    1482.267854
S12    4958.065562
S13    2025.582666
S20      56.226573
S21      42.605367
S22      19.818575
S23      30.097781
S26    1471.048419
S27    1509.827585
S28     911.095470
S29    2085.256385
dtype: float64

Beta diversity

Relationship between alpha, beta, and gamma diversity
One common way to define beta diversity is as a ratio between gamma and alpha: \(D_\beta=\frac{D_\gamma}{D_\alpha}\)
Gamma is the diversity of the pooled samples and alpha is the mean diversity of the individual samples.
Beta diversity takes a value between 1 and N (the number of samples being compared). For functional beta diversity it goes between 1 and NxN.

Beta diversity can be converted to a dissimilarity index, which takes a value between 0 (if the compared samples are identical) and 1 (if the compared samples are completely different). Dissimilarity can be calculated with a local or regional viewpoint. The local viewpoint quantifies the effective proportion of species in a sample that is not shared across all samples. The regional viewpoint quantifies the effective proportion of species in the pooled samples that is not shared across all samples. Just like alpha diversity, beta diversity can be calculated for any diveristy order, q, which means we can decide the emphasis we want to put on the relative abundance of the species.

Just like the alpha diversity indices, beta diversity and dissimilarity can be calculated using naive-, phylogenetic-, and functional indices.

Let’s calculate beta diversity for the example data. First, we’ll calculate pairwise dissimilarities between all samples in our data using the diversity.naive_beta function.

[6]:

import qdiv
obj = qdiv.MicrobiomeData.load_example("Saheb-Alam_DADA2")
rarefied_obj = obj.rarefy()

dissimilarity = qdiv.diversity.naive_beta(rarefied_obj, q=1)
print(dissimilarity.iloc[:5, :5]) #Here we print the results for five samples

           S4        S5        S6        S7       S10
S4   0.000000  0.058442  0.049793  0.075557  0.839512
S5   0.058442  0.000000  0.057203  0.073835  0.783948
S6   0.049793  0.057203  0.000000  0.067338  0.821311
S7   0.075557  0.073835  0.067338  0.000000  0.853798
S10  0.839512  0.783948  0.821311  0.853798  0.000000

The dissimilarity matrix shows the pairwise dissimilarity between the samples. It can, e.g., be used later for plotting an ordination.

We can also check the beta diversity across sample groups with multiple samples. Let’s group the samples based on feed (one recieved acetate and one glucose). We’ll use the diversity.naive_multi_beta function.

[7]:

multi_beta = qdiv.diversity.naive_multi_beta(rarefied_obj, by="feed", q=1)
print(multi_beta)

           N      beta  local_dis  regional_dis
acetate  8.0  1.967056   0.325346      0.325346
glucose  8.0  1.666429   0.245587      0.245587

Each sample group contains 8 samples (N=8). The beta diversity value indicates the effective number of distinct communities represented by the data. In other words, it reflects how much differentiation exists among the samples within a group. To provide additional perspective, beta diversity is also expressed as dissimilarity indices, calculated using either:

local view, which considers differences relative to individual samples, or
regional view, which considers differences relative to the entire pooled set of samples.

The analysis can also be done with diversity.phyl_multi_beta or diversity.func_multi_beta.