Rarefaction
mothur calculates rarefaction curve data using the rarefaction.single command for the number of observed OTUs as a function of distance between sequences and the number of sequences sampled. By theory, the rarefaction curve should match the following expression:
\[S_n = S_t - \left (\frac{\sum_{i=1}^{S_t}{N - N_i \choose n} }{ {N \choose n}} \right )\]where,
\(S_n\) = Average number of OTUs observed after drawing n individuals.
\(S_{t}\) = Total number of OTUs in sample of N total individuals.
In the amazonian dataset the total number of sequences was 98, and there were 84 total OTUs with 75 singletons, 6 doubletons, 1 tripleton, and 2 quadrupletons when defining OTUs at a 0.03 distance cutoff. If we were interested in the predicted observed richness after sampling 50 individuals from this distribution we would evaluate the formula:
\[S_{50} = 84 - 75\frac{\binom{98 - 1}{50}}{\binom{98}{50}} - 6\frac{\binom{98 - 2}{50}}{\binom{98}{50}} - 1\frac{\binom{98 - 3}{50}}{\binom{98}{50}}- 2\frac{\binom{98 - 4}{50}}{\binom{98}{50}}\] \[S_{50} = 84 - 75\left(0.4898 \right) - 6 \left(0.2373 \right) - 1 \left(0.1137 \right) - 2 \left(0.0539 \right) = 45.620\]Although this closed form solution exists, it is actually very slow compared to averaging across a large number of randomized collector’s curves. The following table compares the output based on the formula and using 1,000 permutations using mothur:
n Theory mothur Diff. %Error n Theory mothur Diff. %Error
1 1.000 1.000 0.000 0.000 50 45.620 45.634 -0.014 0.031
2 1.996 1.995 0.001 0.049 51 46.461 46.477 -0.016 0.035
3 2.987 2.986 0.001 0.024 52 47.299 47.317 -0.017 0.037
4 3.974 3.971 0.003 0.069 53 48.136 48.154 -0.019 0.038
5 4.956 4.954 0.002 0.045 54 48.970 48.995 -0.025 0.051
6 5.935 5.934 0.001 0.010 55 49.802 49.822 -0.019 0.038
7 6.909 6.909 0.000 0.002 56 50.633 50.647 -0.014 0.028
8 7.880 7.881 -0.001 0.014 57 51.461 51.478 -0.016 0.032
9 8.846 8.847 -0.001 0.013 58 52.288 52.305 -0.018 0.034
10 9.808 9.809 -0.001 0.011 59 53.112 53.126 -0.014 0.027
11 10.767 10.768 -0.001 0.012 60 53.935 53.946 -0.011 0.020
12 11.721 11.721 0.000 0.001 61 54.755 54.767 -0.012 0.022
13 12.672 12.675 -0.003 0.028 62 55.574 55.586 -0.012 0.021
14 13.619 13.623 -0.004 0.031 63 56.391 56.401 -0.009 0.017
15 14.562 14.566 -0.004 0.027 64 57.206 57.213 -0.007 0.012
16 15.502 15.510 -0.008 0.051 65 58.020 58.026 -0.006 0.010
17 16.438 16.446 -0.008 0.046 66 58.832 58.840 -0.008 0.014
18 17.371 17.378 -0.007 0.040 67 59.642 59.647 -0.005 0.009
19 18.300 18.305 -0.005 0.027 68 60.450 60.448 0.002 0.003
20 19.225 19.233 -0.007 0.037 69 61.256 61.260 -0.004 0.006
21 20.148 20.158 -0.010 0.049 70 62.061 62.057 0.004 0.007
22 21.066 21.079 -0.013 0.061 71 62.865 62.861 0.004 0.006
23 21.982 21.995 -0.012 0.057 72 63.666 63.667 0.000 0.000
24 22.894 22.904 -0.010 0.042 73 64.467 64.469 -0.003 0.004
25 23.804 23.814 -0.010 0.042 74 65.265 65.273 -0.008 0.013
26 24.710 24.719 -0.010 0.039 75 66.062 66.070 -0.008 0.012
27 25.613 25.623 -0.011 0.041 76 66.857 66.862 -0.005 0.007
28 26.512 26.521 -0.008 0.032 77 67.651 67.654 -0.002 0.004
30 28.303 28.313 -0.010 0.034 79 69.235 69.241 -0.007 0.010
31 29.194 29.204 -0.010 0.033 80 70.024 70.035 -0.011 0.015
32 30.082 30.093 -0.010 0.034 81 70.812 70.824 -0.011 0.016
33 30.967 30.984 -0.016 0.052 82 71.599 71.612 -0.013 0.018
34 31.850 31.863 -0.013 0.041 83 72.384 72.391 -0.007 0.010
35 32.729 32.745 -0.015 0.047 84 73.168 73.175 -0.007 0.009
36 33.606 33.618 -0.012 0.035 85 73.950 73.958 -0.008 0.010
37 34.481 34.498 -0.017 0.050 86 74.731 74.746 -0.015 0.020
38 35.352 35.371 -0.019 0.052 87 75.511 75.530 -0.019 0.026
39 36.221 36.241 -0.020 0.055 88 76.289 76.302 -0.013 0.017
40 37.088 37.109 -0.021 0.056 89 77.066 77.081 -0.015 0.020
41 37.952 37.973 -0.022 0.057 90 77.842 77.856 -0.014 0.018
42 38.813 38.835 -0.021 0.055 91 78.616 78.619 -0.003 0.004
43 39.672 39.687 -0.014 0.036 92 79.389 79.391 -0.002 0.003
44 40.529 40.542 -0.013 0.032 93 80.161 80.161 0.000 0.000
45 41.383 41.391 -0.008 0.020 94 80.931 80.930 0.002 0.002
46 42.235 42.245 -0.010 0.023 95 81.700 81.700 0.000 0.000
47 43.085 43.093 -0.009 0.020 96 82.468 82.473 -0.005 0.006
48 43.932 43.945 -0.013 0.029 97 83.235 83.237 -0.003 0.003
49 44.777 44.792 -0.015 0.034 98 84.000 84.000 0.000 0.000