1简介

PhILR是“系统发育等距对数比变换”的缩写(Silverman et al. 2017)．此包提供用于分析的函数成分数据(例如，表示不同变量/部分比例的数据)。具体来说，该软件包允许分析成分数据，其中各部分可以通过系统发育树进行关联(在微生物区系调查数据中很常见)，并提供从系统发育树构建的等长对数比转换，并使用加权参考测量(Egozcue and Pawlowsky-Glahn 2016)．

2PhILR分析综述

PhILR的目标是将组合数据转换为具有系统发育/进化解释的正交无约束空间(实空间)，同时保留原始组合中包含的所有信息。与原始的合成空间不同，在变换后的真实空间中，可以应用标准的统计工具。对于由类群测量值组成的给定样本集，我们将数据转换为一个新的样本空间和称为“平衡”的标准正交坐标。每一个平衡都与系统发生树的一个内部节点相关联，分类单元为叶。余额表示从给定的内部节点下降的两组类群的几何平均丰度的对数比。有关此方法的详细信息，请参见西尔弗曼等人(2017)（链接）.

分析使用丰度表和系统发育树。这些数据可以作为单独的数据对象提供，也可以嵌入到标准R/Bioconductor数据容器中。philr R包支持两个可选的微生物组数据容器，TreeSE(Huang et al. 2021)而且phyloseq(McMurdie and Holmes 2013)．

3.数据集的加载和预处理

我们通过使用最初发布的Global Patterns数据集来演示PhILR分析卡波拉索等人(2011)．

让我们首先加载必要的库。

图书馆(philr);packageVersion(“philr”)

## [1] '1.25.0'

图书馆(猿);packageVersion(“猿”)

## [1] '5.6.2'

图书馆(ggplot2);packageVersion(“ggplot2”)

## [1] '3.3.6'

4数据准备:TreeSE

中最初概述的GlobalPatterns示例工作流(McMurdie and Holmes 2013)．

中检索示例数据TreeSummarizedExperiment（TreeSE)数据格式(Huang et al. 2021)，再举例子也为之phyloseq格式。属性提供了GlobalPatterns数据的TreeSE版本米娅包(Lahti et al. 2020)．

让我们加载数据。

图书馆(mia);packageVersion (mia)

## [1] '1.7.0'

图书馆(dplyr);packageVersion(“dplyr”)

## [1] '1.0.10'

data(GlobalPatterns, package = "mia")

##选择流行类群tse <- GlobalPatterns %>% subsetByPrevalentTaxa(detection = 3，患病率= 20/100,as_relative = FALSE) ##选择样本变量丰度变化显著的类群。Taxa <- apply(assays(tse)$counts, 1, function(x) sd(x)/mean(x) > 3.0) tse <- tse[变量。#倒树!否则将保留包含所有节点的原始树(包括那些从rowData中过滤出来的节点)。tip(phy = rowTree(tse)， tip = rowLinks(tse)$nodeNum) rowTree(tse) <- tree ##添加一个带有pseudocount assays(tse)$counts的新化验。移位<- assays(tse)$counts + 1

##类:treesummarizeexperiment ## dim: 1248 26 ##元数据(0):## assays(2): counts计数。## rownames(1248): 540305 108964…516119 145149 ##行数据名称(7):王国门…属种## colnames(26): CL3 CC1…## colData name (7): X.SampleID Primer…SampleType Description ## reducedDimNames(0): ## mainExpName: NULL ## altExpNames(0): ## rowLinks: a LinkDataFrame(1248行)## rowTree: 1 phylo tree(s)(1248个叶子)## colLinks: NULL ## colTree: NULL

图书馆(猿);packageVersion(“猿”)

## [1] '5.6.2'

Is . roots (tree) #树有根了吗?

##[1]真

is.binary(tree) #所有multichotomies都已解决?

##[1]真

tree <- makeNodeLabel(tree, method="number"， prefix='n') #将修改后的树添加回(' TreeSE ')数据对象rowTree(tse) <- tree

#从TreeSE对象tax中提取分类表<- rowData(tse)[，taxonomyRanks(tse)] #获取名称平衡名称。余额(树，税，'n1')

##[1]“王国古菌/王国细菌”

otu。Table <- t(as(assays(tse))$counts。移位，“矩阵”))树<- rowTree(tse) metadata <- colData(tse) tax <- rowData(tse)[，taxonomyRanks(tse)] otu。表[1:2,1:2]# OTU表

## 540305 108964 ## cl3 1 1 ## cc1 1 2

树#系统发育树

## ##系统发育树，有1248个尖端和1247个内部节点。## ##提示标签:## 540305,108964,175045,546313,54107,71074，…##节点标签:## n1, n2, n3, n4, n5, n6，…## ##扎根;包括分支长度。

head(metadata,2) #元数据

数据框架与2行7列## X.SampleID引物Final_Barcode条形码_truncated_plus_t# # <因子> <因子> <因子> <因子> ## CL3 CL3 ILBC_01 AACGCA TGCGTT ## CC1 CC1 ILBC_02 AACTCG CGAGTT Barcode_full_length样本类型描述## <因子> <因子> <因子> ## CL3 CTAGCGTGCGT土壤Calhoun南卡罗来纳州松树土，pH 4.9 ## CC1 CATCGACGAGT土壤，明尼苏达州雪松河，草地，pH 6.1

头(税，2)#分类表

##王国门类顺序科## <字符> <字符> <字符> <字符> <字符> ## 540305古菌Crenarchaeota Thaumarchaeota Cenarchaeales Cenarchaeaceae ## 108964古菌Crenarchaeota Thaumarchaeota Cenarchaeales Cenarchaeaceae ##属种## <字符> <字符> ## 540305 NA NA ## 108964亚基sopumilus pIVWA5

人类。样本<- factor(colData(tse)$SampleType %in% c(“粪便”，“Mock”，“皮肤”，“舌头”))

5使用PhILR转换数据

这个函数philr: philr ()为philr转换中的关键步骤实现用户友好的包装器。

1. 使用函数将系统发育树转换为其顺序二进制分区(SBP)表示philr: phylo2sbp ()
2. 计算分类单元(又名部分)的权重或使用用户指定的权重
3. 利用该函数从SBP和类群权重建立对比矩阵philr: buildilrBasep ()
4. 转换OTU表为相对丰度(使用philr: miniclo ())和使用权重的“shift”数据集(Egozcue and Pawlowsky-Glahn 2016)使用函数philr: shiftp ()．
5. 使用函数将数据转换为PhILR空间philr: ilrp ()
6. (可选)使用系统发育距离对得到的PhILR空间进行加权。这些权重可以由用户提供，也可以由函数计算philr: calculate.blw ()．

注意:应将预处理过的OTU表传递给函数philr: philr ()在对相对丰度进行封闭(归一化)之前，由于某些分类单元的预设权重使用原始计数数据来降低低丰度分类单元的权重。

这里我们将使用与主论文中相同的权重。

你可以跑philr用丰度表和系统发育树。

全科医生。Philr <- Philr (otu。表，树，part.weights=' enormous .x.gm。计数”,ilr.weights = ' blw.sqrt”)全科医生。philr[1:5,1:5]

## n1 n2 n3 n4 n5 ## CL3 -1.3638521 1.9756259 2.61196 -3.3174292 0.08335109 ## CC1 -0.9441168 3.9054807 2.9804522 -4.7771598 -0.05334306 ## SV1 5.8436901 5.9067782 6.7315081 -8.8020849 0.08335109 ## M31Fcsw -3.9010427 -0.1816618 -0.5432099 0.08335109 ## M11Fcsw -5.4554073 0.5398249 -0.5647474 0.5551616 -0.02389182

或者，您可以直接在TreeSE格式。

全科医生。Philr <- Philr (tse，丰度值= "计数。转移”,part.weights = ' enorm.x.gm。计数”,ilr.weights = ' blw.sqrt”)

或者，您可以在phyloseq格式。为了简单起见，让我们只转换TreeSE对象phyloseq对象来给出一个简单的例子。

gp. pseq <- makephyloseqfromtreesummarize实验(tse， _values="counts. moved ")Philr <- Philr (pseq, part.weights=' enormn .x.gm。计数”,ilr.weights = ' blw.sqrt”)

在运行philr转换后的数据以余额表示，由于每个余额都与树的单个内部节点相关联，因此我们使用分配给内部节点的相同名称表示余额(例如，n1）.

6PhILR空间的圣职仪式

PhILR空间中的欧氏距离可用于排序分析。让我们先计算距离，然后计算标准MDS排序。

#基于philr转换数据gp的样本之间的距离。Dist <- Dist (gp。philr，method="euclidean") # Calculate MDS for the distance matrix d <- as.data.frame(cmdscale(gp.dist)) colnames(d) <- paste0("PC", 1:2)

7用TreeSE进行可视化

接下来让我们想象一下圣职。本例使用标准工具进行排序和可视化，无论首选的数据容器是什么，都可以使用这些工具。注意phyloseq而且TreeSE框架可以提供对额外排序和可视化方法的访问。

#为可视化添加一些元数据d$SampleType <- factor(元数据$SampleType) #创建一个plot ggplot(data = d, aes(x=PC1, y=PC2, color=SampleType)) + geom_point() +实验室(title = "与phILR的欧几里得距离")

8确定区分人类/非人类的平衡

不仅仅是排序分析，PhILR还提供了一个可以使用标准多元工具的完整坐标系。在这里，我们将使用稀疏逻辑回归(从glmnet包装)，以确定少量的平衡，最好地区分人类和非人类样本。

现在我们将拟合一个稀疏逻辑回归模型(逻辑回归与\ (l1 \)处罚)

图书馆(glmnet);packageVersion(“glmnet”)

## [1] '4.1.4'

Glmmod <- glmnet(gp。philr,人类。样本，alpha=1, family="binomial")

的硬阈值\ (l1 \)处罚的\(\lambda = 0.2526\)我们选择哪一个来得到非零系数的个数\(\ \)约5(为了便于本教程中的可视化)。

上面。坐标<- as。矩阵(系数(glmmod, s=0.2526))坐标<- rownames(top. coordinates)[which(top. coordinates)]坐标!= 0)](上。坐标<- top. coordinates [2:length(top. coordinates)]) #删除截距作为坐标

##[1]“n16”“n106”“n122”“n188”“n730”

9名字余额

为了找到与这些平衡相对应的分类标签，我们可以使用这个函数philr: name.balance ()．这个函数使用一个简单的投票方案来分别命名给定平衡的两个后代分支。对于一个给定的分支，分类表是子集，只包含来自该分支的分类单元。从最好的分类等级(例如，物种)开始，检查子集分类表，看看是否有标签(例如，物种名称)代表≥该分类等级表项的阈值(默认95%)。如果没有找到一致的标识符，则在下一个最特定的分类秩(等等)检查表。

Tc.names <- sapply(top。坐标，函数(x)名称。Balance (tree, tax, x)) tc.names

# #但# #“Kingdom_Bacteria / Phylum_Firmicutes”# #不行了n106 # #“Order_Actinomycetales / Order_Actinomycetales”# # n122 # #“Kingdom_Bacteria / Phylum_Cyanobacteria”# # n188 # #“Genus_Campylobacter / Phylum_Proteobacteria”# # n730 # #“Order_Bacteroidales / Order_Bacteroidales”

我们还可以通过直接查看投票来获得更多关于命名的信息。

投票<- name。Balance (tree, tax， 'n730'， return)投票= c('向上'，'向下'))投票[[c('向上'。投票'，'家庭')]]#家庭层面的分子

Porphyromonadaceae ## 1

投票[[c('。投票'，'家庭')]]#家庭层面的分母

##拟杆菌科卟啉单胞菌科普氏菌科Rikenellaceae ## 12 9 10 5

10可视化的结果

图书馆(ggtree);packageVersion(“ggtree”)

## [1] '3.7.0'

图书馆(dplyr);packageVersion(“dplyr”)

## [1] '1.0.10'

上面我们找到了前5个坐标(平衡)，区分样本是来自人类还是非人类来源。现在使用ggtree(Yu et al. 2016)包中，我们可以使用geom_balance对象。为了使用这些函数，我们需要知道树中这些余额的实际节点号(而不仅仅是我们给出的名称)。为了在节点号和名称之间进行转换，我们添加了函数philr: name.to.nn ()而且philr: nn.to.name ()．此外，重要的是，我们知道余额的哪一部分是对数比的分子(+)，哪一部分是分母(-)。为了便于跟踪，我们创建了函数philr: annotate_balance ()这让我们可以很容易地给这两支演化支贴上标签。

tc。nn<- name.to.nn(tree, top.coords) tc.colors <- c('#a6cee3', '#1f78b4', '#b2df8a', '#33a02c', '#fb9a99') p <- ggtree(tree, layout='fan') + geom_balance(node=tc.nn[1], fill=tc.colors[1], alpha=0.6) + geom_balance(node=tc.nn[2], fill=tc.colors[2], alpha=0.6) + geom_balance(node=tc.nn[3], fill=tc.colors[3], alpha=0.6) + geom_balance(node=tc.nn[4], fill=tc.colors[4], alpha=0.6) + geom_balance(node=tc.nn[5], fill=tc.colors[5], alpha=0.6) p <- annotate_balance(tree, 'n16', p=p, labels = c('n16+', 'n16-'), offset.text=0.15, bar=FALSE) annotate_balance(tree, 'n730', p=p, labels = c('n730+', 'n730-'), offset.text=0.15, bar=FALSE)

我们还可以查看人类/非人类来源的这5个平衡的分布。为了画出ggplot2我们首先需要将PhILR转换后的数据转换为长格式。我们包含了一个函数philr: convert_to_long ()为了这个目的。

gp.philr.long <- convert_to_long(gp. philr.long)philr,人类。samples) %>% filter(coord %in% top.coords)

##警告:' gather_() '在tidyr 1.2.0中已弃用。##ℹ请使用“gather()”代替。##ℹ已弃用的功能可能用于philr包。##请在报告问题。

ggplot (gp.philr。Long, aes(x=labels, y=value)) + geom_boxplot(fill='lightgrey') + facet_grid(。~coord, scales='free_x') + labs(x =' Human'， y =' Balance Value') + theme_bw()

11使用天平减维

让我们看看余额n16和余额n730(我们在上面的树中注释的那些)。

图书馆(tidyr);packageVersion(“tidyr”)

## [1] '1.2.1'

gp.philr.long %>% dplyr::rename(Human=labels) %>% dplyr::filter(coord %in% c('n16'， 'n730')) %>% tidyr::spread(coord, value) %>% ggplot(aes(x=n16, y=n730, color=Human)) + geom_point(size=4) + labs(x = tc.names['n16']， y= tc.names['n730']) + theme_bw()

12包版本

sessionInfo ()

## R正在开发中(不稳定)(2022-10-25 r83175) ##平台:x86_64-pc-linux-gnu(64位)##运行在Ubuntu 22.04.1 LTS ## ##矩阵产品:默认## BLAS: /home/biocbuild/bbs-3.17-bioc/R/lib/libRblas。so ## LAPACK: /usr/lib/x86_64-linux-gnu/ LAPACK /liblapack.so.3.10.0 ## ## locale: ## [1] LC_CTYPE=en_US。UTF-8 LC_NUMERIC= c# # [3] LC_TIME=en_GB LC_COLLATE= c# # [5] LC_MONETARY=en_US。utf - 8 LC_MESSAGES = en_US。UTF-8 ## [7] LC_PAPER=en_US。UTF-8 LC_NAME= c# # [9] LC_ADDRESS=C lc_phone = c# # [11] LC_MEASUREMENT=en_US。UTF-8 LC_IDENTIFICATION=C ## ##附加的基本包:## [1]stats4 stats graphics grDevices utils datasets methods ##[8]基础## ##其他附加包:# # # # [1] tidyr_1.2.1 ggtree_3.7.0 [3] glmnet_4.1-4 Matrix_1.5-1 # # [5] dplyr_1.0.10 mia_1.7.0 # # [7] MultiAssayExperiment_1.25.0 TreeSummarizedExperiment_2.7.0 # # [9] Biostrings_2.67.0 XVector_0.39.0 # # [11] SingleCellExperiment_1.21.0 SummarizedExperiment_1.29.0 # # [13] Biobase_2.59.0 GenomicRanges_1.51.0 # # [15] GenomeInfoDb_1.35.0 IRanges_2.33.0 # # [17] S4Vectors_0.37.0 BiocGenerics_0.45.0 # # [19] MatrixGenerics_1.11.0 matrixStats_0.62.0 # # [21] ggplot2_3.3.6 ape_5.6-2 # # [23] philr_1.25.0BiocStyle_2.27.0 ## ##通过命名空间加载(并且没有附加):# # # # [1] shape_1.4.6 jsonlite_1.8.3 [3] magrittr_2.0.3 magick_2.7.3 # # [5] ggbeeswarm_0.6.0 farver_2.1.1 # # [7] rmarkdown_2.17 zlibbioc_1.45.0 # # [9] vctrs_0.5.0 multtest_2.55.0 # # [11] memoise_2.0.1 DelayedMatrixStats_1.21.0 # # [13] rcurl_1.98 - 1.9 htmltools_0.5.3 # # [15] BiocNeighbors_1.17.0 Rhdf5lib_1.21.0 # # [17] rhdf5_2.43.0 gridGraphics_0.5-1 # # [19] sass_0.4.2 bslib_0.4.0 # # [21] plyr_1.8.7 DECIPHER_2.27.0 # # [23] cachem_1.0.6 igraph_1.3.5 # # [25] lifecycle_1.0.3 iterators_1.0.14 # # [27]## [31] GenomeInfoDbData_1.2.9 digest_0.6.30 ## [33] aplot_0.1.8 colorspace_2.18 ## [39] vegan_2.6-4 beachmat_2.15.0 ## [41] labeling_0.4.2 fansi_1.0.3 ## [43] mgcv_1. 1.8-41 compiler_4.3.0 ## [45] bit64_4.0.5 withr_2.5.0 ## [47] BiocParallel_1.33.0 viridis_0.6.2 ## [49] DBI_1.1.3 highr_0.9 ## [53] biomformat_1.27.0 ## [53] aplot_0.1.8 colorspace_2.18 ## [35]permute_0.9-7 ## [55] tools_4.3.0 vipor_0.4.5 ## [57] beeswarm_0.4.0 glue_1.6.2 ## [59] quadprog_1.5-8 nlme_1 .1-160 ## [61] rhdf5filters_1.11.0 grid_4.3.0 ## [63] cluster_2.1.4 reshape2_1.4.4 ## [65] ade4_1.7-19 generics_0.1.3 ## [67] gtable_0.3.1 data.table_1.14.4 ## [69] BiocSingular_1.15.0 ScaledMatrix_1.7.0 ## [71] utf8_1.2.2 ggrepel_0.9.1 ## [73] foreach_1.5.2 pillar_1.8.1 ## [75] string_1 .4.1 yulb .utils_0.0.5 ## [77] splines_4.3.0 treeio_1.23.0 ## [79] lattice_0.20-45幸存者3.4-0 ##[81] bit_4.0.4 tidyselect_1.2.0 ## [83] dirichlet多omial_1.41.0 scuttle_1.9.0 ## [85] knitr_1.40 gridExtra_2.3 ## [87] bookdown_0.29 phyloseq_1.43.0 ## [91] xfun_0.34 stringi_1.7.8 ## [93] yaml_2.3.6 evaluate_0.17 ## [95] codetools_0.2-18 tibble_3.1.8 ## [97] BiocManager_1.30.19 ggplotify_0.1.0 ## [99] cli_3.4.1 munsell_0.5.0 ## [101] jquerylib_0.1.4 Rcpp_1.0.9 ## [105] assertthat_0.2.1 blob_1.2.3 ## [107] sparseMatrixStats_1.11.0 bitops_1.0-7 ## [109] phangorn_2.10.0 decontam_1.19.0 ## [111] viridisLite_0.4.1 tidytree_0.4.1 ## [113] scales_1.2.1 purrr_0.3.5 ## [115] crayon_1.5.2 rlang_1.0.6 ## [117] fastmatch_1.1-3