{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c41951cc-868a-4399-99c0-b0f674ec67c7",
   "metadata": {},
   "source": [
    "# Null models\n",
    "Null models are powerful tools for investigating community assembly mechanisms by comparing observed patterns to those expected under random processes.\n",
    "In qdiv, several widely used indices are implemented:\n",
    "- Raup–Crick index for assessing beta diversity deviations from randomness\n",
    "- Net Relatedness Index (NRI) and Nearest Taxon Index (NTI) for both alpha and beta diversity, which evaluate phylogenetic clustering or dispersion\n",
    "\n",
    "*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!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a842f143-d51f-4f73-abc9-dc4d7bcc3636",
   "metadata": {},
   "source": [
    "Let's load some example data and try:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "91ee1540-fd5d-47e4-8818-724e6c74b59b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<MicrobiomeData: 672 features, 16 samples, tax=yes, meta=yes, seq=yes, tree=yes>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import qdiv\n",
    "obj = qdiv.MicrobiomeData.load_example(\"Saheb-Alam_DADA2\") #First we load example data\n",
    "obj.rename_features(inplace=True, name_type=\"ASV\") #This is the name the features ASV1, ASV2\n",
    "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."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c2eca19-6164-4268-95fc-baf8b2a0bc85",
   "metadata": {},
   "source": [
    "First, we'll calculate the Raup-Crick index using the **model.rcq** function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a739944e-feb5-4965-8db6-c7366f54f193",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3fb503c2736b4eb1ac946d5b790ccc94",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Randomizing tables:   0%|          | 0/1000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2fdaa1f8c4c14ff183fe1076f95f6ef7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Comparing beta diversity:   0%|          | 0/1000 [00:00<?, ?iter/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "raup_crick_results = qdiv.model.rcq(obj, div_type=\"naive\", q=1, iterations=1000) #we specify diversity type (div_type) and diversity order (q)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c9f5b25-db45-4039-b6b1-74453124435d",
   "metadata": {},
   "source": [
    "This took a few minutes to run. The time it takes depends on the number of iterations. Let's look at the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "38818282-e17d-4a91-aee1-e53b47aee43b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_keys(['div_type', 'obs_d', 'p', 'null_mean', 'null_std', 'ses'])\n"
     ]
    }
   ],
   "source": [
    "print(raup_crick_results.keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d27bdb5-1f94-4699-984d-ab44b98b56d6",
   "metadata": {},
   "source": [
    "The result is a python dictionary. The first entry, div_type, is information about the diversity type that was calculate (naive in this case).  \n",
    "The rest of entries are pandas dataframe with the results. \n",
    "- **obs_d** is the observed dissimilarities between the samples\n",
    "- **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.\n",
    "- **null_mean** and **null_std** are the mean and standard deviation of the dissimilarity values for the randomizations.\n",
    "- **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.\n",
    "\n",
    "Let's look at a few of the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "4bd0e50f-9ba2-446a-9381-e41cf0814c42",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "--------\n",
      "Observed dissimilarity\n",
      "           S4        S5        S6        S7       S10\n",
      "S4   0.000000  0.058848  0.049525  0.074329  0.837658\n",
      "S5   0.058848  0.000000  0.057854  0.073645  0.778168\n",
      "S6   0.049525  0.057854  0.000000  0.067191  0.819105\n",
      "S7   0.074329  0.073645  0.067191  0.000000  0.853965\n",
      "S10  0.837658  0.778168  0.819105  0.853965  0.000000\n",
      "\n",
      "--------\n",
      "Mean dissimilarity of the randomized tables (i.e. the null expectation)\n",
      "           S4        S5        S6        S7       S10\n",
      "S4   0.000000  0.645231  0.734812  0.735250  0.623847\n",
      "S5   0.645231  0.000000  0.645889  0.655437  0.490451\n",
      "S6   0.734812  0.645889  0.000000  0.727974  0.624604\n",
      "S7   0.735250  0.655437  0.727974  0.000000  0.635891\n",
      "S10  0.623847  0.490451  0.624604  0.635891  0.000000\n",
      "\n",
      "--------\n",
      "The p index\n",
      "        S4     S5     S6     S7    S10\n",
      "S4   0.500  0.000  0.000  0.000  0.926\n",
      "S5   0.000  0.500  0.000  0.000  0.986\n",
      "S6   0.000  0.000  0.500  0.000  0.884\n",
      "S7   0.000  0.000  0.000  0.500  0.948\n",
      "S10  0.926  0.986  0.884  0.948  0.500\n",
      "\n",
      "--------\n",
      "The ses index\n",
      "           S4        S5        S6        S7       S10\n",
      "S4        NaN  3.413654  3.903284  3.717669 -1.327855\n",
      "S5   3.413654       NaN  3.550045  3.502462 -1.945836\n",
      "S6   3.903284  3.550045       NaN  3.642286 -1.169599\n",
      "S7   3.717669  3.502462  3.642286       NaN -1.351313\n",
      "S10 -1.327855 -1.945836 -1.169599 -1.351313       NaN\n"
     ]
    }
   ],
   "source": [
    "print(\"\\n--------\\nObserved dissimilarity\")\n",
    "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 \n",
    "\n",
    "print(\"\\n--------\\nMean dissimilarity of the randomized tables (i.e. the null expectation)\")\n",
    "print(raup_crick_results[\"null_mean\"].iloc[:5, :5]) \n",
    "\n",
    "print(\"\\n--------\\nThe p index\")\n",
    "print(raup_crick_results[\"p\"].iloc[:5, :5]) \n",
    "\n",
    "print(\"\\n--------\\nThe ses index\")\n",
    "print(raup_crick_results[\"ses\"].iloc[:5, :5]) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5910818b-5ca1-43f0-b505-0c20fec0b5e4",
   "metadata": {},
   "source": [
    "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\"."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "127565db-e30f-4bc4-acfa-80a5c70edeaf",
   "metadata": {},
   "source": [
    "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. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "3641b41e-86fc-4378-a316-fd64bc10c608",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computing pairwise tree distances...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "Leaves:   0%|                                                                                  | 0/672 [00:00<?, ?it/s]\u001b[A\n",
      "Leaves:  17%|████████████                                                          | 116/672 [00:00<00:00, 1106.65it/s]\u001b[A\n",
      "Leaves:  39%|███████████████████████████▍                                          | 263/672 [00:00<00:00, 1312.48it/s]\u001b[A\n",
      "Leaves: 100%|██████████████████████████████████████████████████████████████████████| 672/672 [00:00<00:00, 2091.95it/s]\u001b[A\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Tree distance matrix saved.\n",
      "Done!\n"
     ]
    }
   ],
   "source": [
    "dm = qdiv.sequences.tree_distance_matrix(obj, savename=\"Tree_distmat\") #The function will also save the distance matrix as a file"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9e70d4b-e146-48d1-bc49-d88df75fbcdc",
   "metadata": {},
   "source": [
    "Now, let's calculate NRI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "339e2cd1-05e0-4a23-8865-8888ad323abc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "02dee54746fa41a1b00438c2c0ba6213",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "NRIq null permutations:   0%|          | 0/1000 [00:00<?, ?iter/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nri_results = qdiv.model.nriq(obj, dm, q=1, iterations=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96875dad-2dc3-47de-9c40-10132e286457",
   "metadata": {},
   "source": [
    "The results are in a pandas dataframe. Let's have a look."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "94f21708-9c1e-4ca4-8525-fd95d07729dd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "         MPDq  null_mean  null_std      p       ses\n",
      "S4   0.198582   0.252812  0.063203  0.129  0.858044\n",
      "S5   0.269714   0.324084  0.074221  0.186  0.732540\n",
      "S6   0.219472   0.280447  0.068153  0.106  0.894674\n",
      "S7   0.201862   0.249928  0.062276  0.162  0.771805\n",
      "S10  0.407732   0.504229  0.071405  0.007  1.351419\n",
      "S11  0.383323   0.492184  0.077297  0.004  1.408337\n",
      "S12  0.440235   0.526294  0.051528  0.003  1.670146\n",
      "S13  0.418713   0.518590  0.062794  0.005  1.590538\n",
      "S20  0.235399   0.378035  0.109250  0.038  1.305596\n",
      "S21  0.230203   0.382119  0.108781  0.028  1.396537\n",
      "S22  0.203751   0.336986  0.105135  0.045  1.267286\n",
      "S23  0.208140   0.364889  0.107591  0.022  1.456892\n",
      "S26  0.489631   0.508709  0.068199  0.432  0.279745\n",
      "S27  0.453988   0.505221  0.072082  0.224  0.710756\n",
      "S28  0.437841   0.504393  0.074751  0.147  0.890305\n",
      "S29  0.402366   0.516238  0.065288  0.006  1.744144\n"
     ]
    }
   ],
   "source": [
    "print(nri_results)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64be62e7-b500-4c4f-8ab1-c8134d9b4dac",
   "metadata": {},
   "source": [
    "- First column: Shows the mean phylogenetic distance (MPD) calculated for the specified diversity order (q) in the function.\n",
    "- Next two columns: Show the mean and standard deviation of the MPD from the randomized communities (null model).\n",
    "- 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).\n",
    "- ses: Positive value indicate phylogenetic clustering. Negative value indicate phylogentic overdispersion.\n",
    "\n",
    "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."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}