Небольшой пример обработки микробных данных на примере анализа измкенчивости микробиома литторин (данные изменены)
Осторожно, встречаются неоптимизированные действия, неединообразный код и конфликты пакетов.
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 7920 taxa and 80 samples ]
## sample_data() Sample Data: [ 80 samples by 12 sample variables ]
## tax_table() Taxonomy Table: [ 7920 taxa by 8 taxonomic ranks ]
## [1] "L1_L_fab_gut2" "L1_L_nod_BF2" "L1_L_obt_gut1" "L1_L_obt_gut2"
## [5] "L1_L_obt_gut3" "L1_L_obt_gut5" "L1_L_obt_gut6" "L1_L_obt_gut7"
## [9] "L1_L_obt_gut8" "L1_L_ser_BF1" "L1_L_ser_BF2" "L1_L_stones_BF1"
## [13] "L1_L_ves_BF1" "L1_L_ves_BF2" "L1_L_ves_BF3" "L1_U_obt_gut1"
## [17] "L1_U_obt_gut2" "L1_U_obt_gut3" "L1_U_obt_gut7" "L11_L_fab_gut2"
## [21] "L11_L_fab_gut3" "L11_L_fab_gut4" "L11_L_fab_gut5" "L11_L_fab_gut6"
## [25] "L11_L_fab_gut7" "L11_L_fab_gut8" "L11_L_fab_tent1" "L11_L_fab_tent2"
## [29] "L11_L_fill_BF1" "L11_L_fill_BF2" "L11_L_fill_BF3" "L11_L_nod_BF2"
## [33] "L11_L_nod_BF3" "L11_L_obt_gut3" "L11_L_obt_gut4" "L11_L_obt_gut5"
## [37] "L11_L_obt_gut6" "L11_L_obt_gut7" "L11_L_obt_gut8" "L11_L_obt_tent1"
## [41] "L11_L_obt_tent2" "L11_L_stones_BF1" "L11_L_stones_BF2" "L11_L_stones_BF3"
## [45] "L11_U_fill_BF1" "L11_U_nod_BF2" "L11_U_obt_gut1" "L11_U_obt_gut2"
## [49] "L11_U_obt_gut3" "L11_U_obt_gut4" "L11_U_obt_gut5" "L11_U_obt_gut8"
## [53] "L11_U_obt_tent2" "L11_U_stones_BF1" "L11_U_stones_BF3" "L5_L_fab_gut1"
## [57] "L5_L_fab_gut3" "L5_L_fab_gut5" "L5_L_fab_gut6" "L5_L_fab_gut7"
## [61] "L5_L_fab_gut8" "L5_L_nod_BF2" "L5_L_obt_gut1" "L5_L_obt_gut3"
## [65] "L5_L_obt_gut4" "L5_L_obt_gut7" "L5_L_obt_gut8" "L5_L_stones_BF2"
## [69] "L5_L_ves_BF3" "L5_U_fill_BF1" "L5_U_fill_BF2" "L5_U_nod_BF1"
## [73] "L5_U_nod_BF3" "L5_U_obt_gut1" "L5_U_obt_gut2" "L5_U_obt_gut4"
## [77] "L5_U_obt_gut5" "L5_U_obt_gut6" "L5_U_obt_gut8" "L5_U_obt_tent2"
## [1] "Domain" "Phylum" "Class" "Order" "Family" "Genus" "Species"
## [8] "X"
## [1] "site" "species" "sp_type"
## [4] "subtype" "type" "sufficient"
## [7] "level" "sea" "season"
## [10] "season_time" "season_site" "season_site_level"
Фильтрация мягкая - удаляем ASV, встреченные 1 раз во всем наборе данных (singletones), ASV, присутствующие менее чем в 5% образцов. Порог обилия моджет быть иным, но даже если удалить фичи, встреченные менее, чем в 10%, это может помешать найти различия между образцами по редким микробам.
## L11_L_obt_gut5 L11_L_obt_gut7 L1_L_stones_BF1 L5_L_ves_BF3 L11_L_fill_BF2
## 0 0 2715 4108 4621
## L1_L_obt_gut7 L11_L_obt_gut6 L5_U_obt_gut2 L5_L_obt_gut7 L11_L_obt_gut3
## 5210 5366 5494 6242 6543
## L5_L_obt_gut8 L11_U_obt_gut1 L11_L_fab_tent2 L5_U_obt_gut5 L5_L_fab_gut6
## 7368 8125 8395 8413 8469
## L11_L_stones_BF2 L1_L_obt_gut5 L1_L_fab_gut2 L1_L_obt_gut2
## 34988 34010 29782 28611
## L1_U_obt_gut3 L1_U_obt_gut1 L1_L_obt_gut8 L1_U_obt_gut7
## 26706 25122 23868 23568
## L11_L_stones_BF3 L5_U_nod_BF3 L1_L_obt_gut1 L11_U_nod_BF2
## 22115 21893 21624 21035
## L1_L_obt_gut3 L11_U_stones_BF3 L11_U_stones_BF1
## 20357 19611 19070
первые два образца (с наиболее низкой глубиной прочтения, она НОЛЬ) далее в анализе не используются
кроме того, имеет смысл удалить лишние летние образцы (Лето 2021, неполные L.fabalis)
было принято волевое решение не удалять единственную пробу F.vesiculosus
## L5_L_ves_BF3 L11_L_fill_BF2 L11_L_obt_gut6 L5_U_obt_gut2 L5_L_obt_gut7
## 4108 4621 5366 5494 6242
## L11_L_obt_gut3 L5_L_obt_gut8 L11_U_obt_gut1 L11_L_fab_tent2 L5_U_obt_gut5
## 6543 7368 8125 8395 8413
## L5_L_fab_gut6 L5_L_obt_gut1 L11_L_nod_BF2 L5_L_obt_gut3 L11_L_obt_gut4
## 8469 8531 8536 8813 8861
Значения размеров библиотек для самых больших образцов если что такие:
## L5_U_obt_gut1 L11_L_stones_BF1 L11_U_obt_gut3 L11_U_obt_gut5
## 15521 15837 17379 17694
## L11_U_obt_gut8 L11_L_fill_BF3 L11_U_obt_gut4 L5_L_nod_BF2
## 17959 18104 18279 18738
## L11_L_fill_BF1 L11_U_stones_BF1 L11_U_stones_BF3 L11_U_nod_BF2
## 18896 19070 19611 21035
## L5_U_nod_BF3 L11_L_stones_BF3 L11_L_stones_BF2
## 21893 22115 34988
Мы в своем анализе делили на медианный размер библиотеки (здесь - 11539) для всех образцов, но можно хоть логарифмировать
Давайте сначала опишем простые вещи
альфа-ранообразие, полный датасет, объединено по типу образца
Наконец-то многомерный подход
Базово я привыкла делать nMDS по матрице различий Брея-Куртиса, но с некоторыми (напр, логарифмическими нормализациями он не дружит) Можно делать и Жаккара в качестве расстояний, и хоть главные компоненты в качестве ординации Можно делать хитмап или кластеризацию, главное понять, на какие расстояния относит композиция микробиома ваши образцы друг ото друга
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.200135
## Run 1 stress 0.2243491
## Run 2 stress 0.2147016
## Run 3 stress 0.2089865
## Run 4 stress 0.2105712
## Run 5 stress 0.205666
## Run 6 stress 0.2051715
## Run 7 stress 0.223335
## Run 8 stress 0.2075067
## Run 9 stress 0.2032072
## Run 10 stress 0.2146564
## Run 11 stress 0.2043048
## Run 12 stress 0.2136925
## Run 13 stress 0.2033042
## Run 14 stress 0.2151472
## Run 15 stress 0.2036829
## Run 16 stress 0.1987382
## ... New best solution
## ... Procrustes: rmse 0.02500797 max resid 0.165155
## Run 17 stress 0.2072265
## Run 18 stress 0.204302
## Run 19 stress 0.1969736
## ... New best solution
## ... Procrustes: rmse 0.03306505 max resid 0.2160054
## Run 20 stress 0.2051777
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 20: stress ratio > sratmax
## [1] 0.1969736
будем тестировать с помощью PERMANOVA, есть ли влияние вида для кишечного микробиома и есть ли влияние сезона
## [1] FALSE
##
## May-June_2022 November_2021
## 17 17
##
## L_fabalis L_obtusata
## 13 21
##
## Analysis of Variance, using Residual Randomization
## Permutation procedure: Randomization of null model residuals
## Number of permutations: 1000
## Estimation method: Ordinary Least Squares
## Sums of Squares and Cross-products: Type I
## Effect sizes (Z) based on F distributions
##
## Df SS MS Rsq F Z Pr(>F)
## species 1 0.3341 0.33409 0.09513 1.5769 1.3603 0.089 .
## Residuals 15 3.1779 0.21186 0.90487
## Total 16 3.5120
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call: lm.rrpp(f1 = guts_white_bray_snails ~ species, data = guts_white_df)