Null models
Null models are powerful tools for investigating community assembly mechanisms by comparing observed patterns to those expected under random processes. In qdiv, several widely used indices are implemented:
Raup–Crick index for assessing beta diversity deviations from randomness
Net Relatedness Index (NRI) and Nearest Taxon Index (NTI) for both alpha and beta diversity, which evaluate phylogenetic clustering or dispersion
What’s special with qdiv? These indices can be calculated for any diversity order (q), allowing you to incorporate abundance sensitivity into null model analyses!
Let’s load some example data and try:
[1]:
import qdiv
obj = qdiv.MicrobiomeData.load_example("Saheb-Alam_DADA2") #First we load example data
obj.rename_features(inplace=True, name_type="ASV") #This is the name the features ASV1, ASV2
obj.tax_prefix(add=True, inplace=True) #This is to add prefix to the taxonomic classified, i.e., d__ for domain, p__ for phylum, etc.
[1]:
<MicrobiomeData: 672 features, 16 samples, tax=yes, meta=yes, seq=yes, tree=yes>
First, we’ll calculate the Raup-Crick index using the model.rcq function.
[3]:
raup_crick_results = qdiv.model.rcq(obj, div_type="naive", q=1, iterations=1000) #we specify diversity type (div_type) and diversity order (q)
This took a few minutes to run. The time it takes depends on the number of iterations. Let’s look at the results.
[5]:
print(raup_crick_results.keys())
dict_keys(['div_type', 'obs_d', 'p', 'null_mean', 'null_std', 'ses'])
obs_d is the observed dissimilarities between the samples
p is the “Raup-Crick” index. A p close to 0 means that observed dissimilarity is lower than the null expectation with randomized table while a p close to 1 means that it is significantly higher than the null expectation.
null_mean and null_std are the mean and standard deviation of the dissimilarity values for the randomizations.
ses is the Raup-Crick index expressed as standardized effect size \(\frac{(d_{null,mean}-d_{obs})}{d_{null,std}}\). A positive value indicates that the observed dissimilarity is lower than the null expectation while a negative value means that the community are most dissimilar than the null expectation.
Let’s look at a few of the results:
[9]:
print("\n--------\nObserved dissimilarity")
print(raup_crick_results["obs_d"].iloc[:5, :5]) #.iloc[:5, :5] means that we only print the first five rows and columns of the dataframe
print("\n--------\nMean dissimilarity of the randomized tables (i.e. the null expectation)")
print(raup_crick_results["null_mean"].iloc[:5, :5])
print("\n--------\nThe p index")
print(raup_crick_results["p"].iloc[:5, :5])
print("\n--------\nThe ses index")
print(raup_crick_results["ses"].iloc[:5, :5])
--------
Observed dissimilarity
S4 S5 S6 S7 S10
S4 0.000000 0.058848 0.049525 0.074329 0.837658
S5 0.058848 0.000000 0.057854 0.073645 0.778168
S6 0.049525 0.057854 0.000000 0.067191 0.819105
S7 0.074329 0.073645 0.067191 0.000000 0.853965
S10 0.837658 0.778168 0.819105 0.853965 0.000000
--------
Mean dissimilarity of the randomized tables (i.e. the null expectation)
S4 S5 S6 S7 S10
S4 0.000000 0.645231 0.734812 0.735250 0.623847
S5 0.645231 0.000000 0.645889 0.655437 0.490451
S6 0.734812 0.645889 0.000000 0.727974 0.624604
S7 0.735250 0.655437 0.727974 0.000000 0.635891
S10 0.623847 0.490451 0.624604 0.635891 0.000000
--------
The p index
S4 S5 S6 S7 S10
S4 0.500 0.000 0.000 0.000 0.926
S5 0.000 0.500 0.000 0.000 0.986
S6 0.000 0.000 0.500 0.000 0.884
S7 0.000 0.000 0.000 0.500 0.948
S10 0.926 0.986 0.884 0.948 0.500
--------
The ses index
S4 S5 S6 S7 S10
S4 NaN 3.413654 3.903284 3.717669 -1.327855
S5 3.413654 NaN 3.550045 3.502462 -1.945836
S6 3.903284 3.550045 NaN 3.642286 -1.169599
S7 3.717669 3.502462 3.642286 NaN -1.351313
S10 -1.327855 -1.945836 -1.169599 -1.351313 NaN
The results are as we could expect for this example data. Samples S4 and S5 have much lower observed dissilarity than in the randomized null case. Consequenctly, p is 0 and ses is high. Samples S4 and S10, on the other hand, have higher dissimilarity than the null case, resulting in high p and negative ses. Often ses values higher than 2 or lower than -2 would be considered “statistically significant”.
Let’s also try one of the phylogenetic indices. The net relatedness index is calculated using the model.nriq function and the nearest taxon index is calculated using model.ntiq. Both function require that we convert our phylogenetic tree data to a distance matrix with information about the phylogenetic distance between each pair of features in our data. This can be done using the sequences.tree_distance_matrix function.
[10]:
dm = qdiv.sequences.tree_distance_matrix(obj, savename="Tree_distmat") #The function will also save the distance matrix as a file
Computing pairwise tree distances...
Leaves: 0%| | 0/672 [00:00<?, ?it/s]
Leaves: 17%|████████████ | 116/672 [00:00<00:00, 1106.65it/s]
Leaves: 39%|███████████████████████████▍ | 263/672 [00:00<00:00, 1312.48it/s]
Leaves: 100%|██████████████████████████████████████████████████████████████████████| 672/672 [00:00<00:00, 2091.95it/s]
Tree distance matrix saved.
Done!
Now, let’s calculate NRI.
[14]:
nri_results = qdiv.model.nriq(obj, dm, q=1, iterations=1000)
The results are in a pandas dataframe. Let’s have a look.
[16]:
print(nri_results)
MPDq null_mean null_std p ses
S4 0.198582 0.252812 0.063203 0.129 0.858044
S5 0.269714 0.324084 0.074221 0.186 0.732540
S6 0.219472 0.280447 0.068153 0.106 0.894674
S7 0.201862 0.249928 0.062276 0.162 0.771805
S10 0.407732 0.504229 0.071405 0.007 1.351419
S11 0.383323 0.492184 0.077297 0.004 1.408337
S12 0.440235 0.526294 0.051528 0.003 1.670146
S13 0.418713 0.518590 0.062794 0.005 1.590538
S20 0.235399 0.378035 0.109250 0.038 1.305596
S21 0.230203 0.382119 0.108781 0.028 1.396537
S22 0.203751 0.336986 0.105135 0.045 1.267286
S23 0.208140 0.364889 0.107591 0.022 1.456892
S26 0.489631 0.508709 0.068199 0.432 0.279745
S27 0.453988 0.505221 0.072082 0.224 0.710756
S28 0.437841 0.504393 0.074751 0.147 0.890305
S29 0.402366 0.516238 0.065288 0.006 1.744144
First column: Shows the mean phylogenetic distance (MPD) calculated for the specified diversity order (q) in the function.
Next two columns: Show the mean and standard deviation of the MPD from the randomized communities (null model).
p-value: Close to 0 means that MPD is lower in the real community than in the randomized (i.e. clustering). Close 1 means that MPD in the real community is higher than expected (i.e. overdispersion).
ses: Positive value indicate phylogenetic clustering. Negative value indicate phylogentic overdispersion.
In this case, our results indicate phylogenetic clustering across most samples, as shown by positive ses values and low p-values. However, the clustering is not statistically significant under common thresholds (e.g., |ses| ≥ 2 or p < 0.05 for strong evidence). Most ses values fall between 0.7 and 1.7, suggesting a tendency toward clustering rather than strong phylogenetic structure.