前一阵我们翻译Rob Knight的综述，1.8万字，让你熟读2遍轻松握掌微生物组领域分析框架、把握未来分析趋势。目前在宏基因组平台累计1.9万人次，热心肠平台首发阅读8500+，科学网加精置顶阅读8700+，CSDN阅读1200+，四大平台阅读人数近4万次。

还没有仔细学习的，赶快读两遍吧！

• Nature综述：大佬手把手教你分析菌群数据——1.8万字

其中提到了一款追踪微生物来源的软件SourceTracker，被多位朋友留言，今天为大家分享一下软件的结果解读。同时详细说明输入文件的准备和代码的注释，方便大家上手应用于自己的数据。

SourceTracker有什么用？

用途是可以识别相关各组间来源的分析，如具体的问题：

• 婴儿的肠道菌群有哪些继承了母亲的肠道菌群、哪些来自阴道菌群、哪些来自皮肤

• 法医学的应用，尸体中的菌群与来源土壤的鉴定、腐败菌来自本身，还是周围环境

• 河流污染物的来源分析、周围工厂、农田、养殖厂对河流污染的贡献和来源追溯。

• 分析植物菌组形成过程：植物根际菌在土壤中来源和种子来源；叶际菌群的土壤来源比例等。

软件简介

Bayesian community-wide culture-independent microbial source tracking 于2011年发表于Nature Methods

由圣地亚哥大学的Scott教授及Rob Knight团队合作完成。第一作者Dan Knights。

Google统计目前引用299次。

该软件中目标样本为Sink，微生物污染源或来源的样品为Source；基于贝叶斯算法，探究目标样本（Sink）中微生物污染源或来源（Source）的分析。根据Source样本和Sink样本的群落结构分布，来预测Sink样本中来源于各Source样本的组成比例。

我们之前解读过Rob Knight的一篇Sciences文章中图2A就使用此软件分析确定尸体腐败过程中主要菌来自于土壤的结果。

• 16S+功能预测发Sciences

软件结果解读

SourceTracker分析图a，预测样本来源比例柱状图。一幅图代表一个预测样本，用不同颜色的柱子表示该样本中各来源的比例，Unknow代表未知来源分类，误差线代表100次Gibbs采样的标准差。

SourceTracker分析图b，预测样本来源比例面积图。一幅图代表一个预测样本，不同颜色代表不同来源的比例，每一列代表一次Gibbs采样结果，100次Gibbs采样结果按照相近的排列顺序进行展示。

SourceTracker分析图c，预测样本来源比例饼图。一个饼图代表一个预测样本，不同颜色扇形的比例代表该预测样本中各来源的比例。

原文解读

DOI: 10.1038/nmeth.1650

图1 SourceTracker和其他模型的比较。所示模型估计模拟样本中两个源环境的比例，Jensen-Shannon差异表示环境之间的重叠程度从0（完全相同，因此不可能消除歧义）到1（完全不重叠，因此容易区分）。绘制了估计比例的确定系数（r2）。每个点代表100个样本的三次试验的平均R2；误差条显示均值的标准误s.e.m.（n=3）。

Figure 1 | Comparison of SourceTracker and other models. Indicated models estimate the proportions of two source environments in simulated samples, as the degree of overlap between the environments was varied from a Jensen-Shannon divergence of 0 (completely identical and thus impossible to disambiguate) to 1 (completely non-overlapping and thus trivial to disambiguate). The coefficients of determination (R2) of the estimated proportions are plotted. Each point represents the mean R2 for three trials of 100 samples each; error bars show s.e.m. (n = 3).

Nature Method [1]原文：每个图形代表一个样本Sink，分别是Lab1的PCR水、NICU桌子、办公室电话；不同颜色表示不同样本来源Source，所占面积为在Sink样本中各来源的比例。

图2. SourceTracker对洗碗池样本子集的比例估计。（a–c）使用SourceTracker估计的三个洗碗池样本的源环境比例，每个源环境中包括45个训练样本：吉布斯采样中100次取样的平均比例（a），相同样本的数据，包括比例估计的标准变异S.D.（b）。100次吉布斯绘制的可视化；每个列显示一次采样的结果，列按保持相似的混合物进行排列在一起，使图形看起来看美观、更容易观察和比较规律（c）。

Figure 2 | SourceTracker proportion estimates for a subset of sink samples. (a–c) Source environment proportions for three sink samples estimated using SourceTracker and 45 training samples from each source environment: mean proportions for 100 draws from Gibbs sampling (a), data for the same samples, including s.d. of the proportion estimates (b), and visualization of the 100 Gibbs draws; each column shows the mixture from one draw, with columns ordered to keep similar mixtures together (c).

图3 常见污染操作分类单元（OTU）的相对丰度。SourceTracker可以为水槽样本中的每个OTU观测序列分配不同的源环境。这十个OTU源对在水槽环境中具有最高的平均相对丰度，不包括未知源。图例给出了OTU的属级分类、OTU标识符和分配给这些观测的源环境。值得注意的是，被归类为肠杆菌的OTU，一种常见于肠道的谱系，在皮肤训练样本中比在肠道训练样本中更为普遍。

Figure 3 | Relative abundance of common contaminating operational taxonomic units (OTUs). SourceTracker may assign a different source environment to each observation (sequence) of an OTU in the sink samples. These ten OTU-source pairs had the highest average relative abundance across sink environments, excluding the unknown source. The legend gives the genus-level taxonomic classification of the OTU, the OTU identifier and the source environment assigned to these observations of the OTU. Note that the OTU classified as Enterobacter, a lineage commonly seen in the gut, was more prevalent in the skin training samples than the gut training samples.

文章实战解读

Rob Knight-2016-Sciences[2]文章中图(A) 动态贝叶斯推理网络: 尸体分解过程微生物分类群神经信息流动网络，土壤是其主要来源。小鼠尸体四种取样位置分别为头、躯干、腹部和土壤。颜色为3种环境，分别为沙漠、草地和森林，且均与土壤来源微生物非常显著相关。

更多内容，请阅读下文：

• 16S+功能预测发Sciences：尸体降解过程

[3]产道菌群移植对剖腹产婴儿缺失菌群的恢复：图中展示了婴儿三个部位皮肤Skin、口腔Oral、肛门Anal中，肠道菌群组成的来源，随时间的推移1-30天而发生的改变。

剖腹产的婴儿患免疫和代谢疾病的风险增高，被认为可能由于缺乏了与母亲生殖道分泌物（包括微生物）的接触；母亲生殖道的分泌物会覆盖顺产婴儿的全身，促进了婴儿口腔、肠道、皮肤菌群的定殖，以及对婴儿的保护作用；对剖腹产婴儿涂抹分泌物，随时间推移，其各部位菌群特征逐渐趋向于顺产婴儿。该方法可以部分的恢复剖腹产婴儿菌群，但对健康的长期影响有待观察，以及样本量也需要扩大。

软件安装

SourceTracker是一个R脚本， 最新版本地址：https://github.com/danknights/sourcetracker ，版本1.0，2016年9月18日更新

# 下载脚本和测试数据
git clone git@github.com:danknights/sourcetracker.git
cd sourcetracker/

example.r使用data目录中的OTU表和mapping实现文章中的分析实例

# 运行测试数据
Rscript example.r

运行中显示如下结果:

Rarefying training data at 1000Gut              Oral            Skin            Soil         Unknown
..........   1 of  125, depth=   12:  0.12 (0.09)       0.07 (0.03)     0.11 (0.06)     0.02 (0.04)     0.68 (0.09)
..........   2 of  125, depth=   20:  0.22 (0.03)       0.08 (0.03)     0.29 (0.05)     0.03 (0.04)     0.39 (0.05)
..........   3 of  125, depth=   10:  0.60 (0.00)       0.00 (0.00)     0.25 (0.13)     0.00 (0.00)     0.15 (0.13)
..........   4 of  125, depth=   10:  0.28 (0.06)       0.09 (0.03)     0.03 (0.05)     0.01 (0.03)     0.59 (0.07)
..........   5 of  125, depth=    2:  0.00 (0.00)       0.00 (0.00)     0.45 (0.16)     0.00 (0.00)     0.55 (0.16)

先将训练样本抽平至1000条。再进行125次的重采样，来源包括肠、口腔、皮肤、土壤和未知五大类来源。结果为以上讲解的几种图形，建议在Rstudio中打开R脚本绘图，交互探索结果。

输入文件准备

按照data/metadata.txt标准的实验设计和data/otus.txt标准的OTU表，在example.r中修改对应的分组，即可分析自己的数据了，非常easy。最好按照示例的文件填写内容，减少代码修改直接运行，以下有实验设计信息说明：

实验设计

示例文件：data/metadata.txt

主要包括的列有：样品名、描述、环境Env、来源SourceSink、研究、细节

你的实验设计必须有前4列

1. #SampleID列：样品编号，以英文字母开头，最好只包括字母和数字，必须与OTU表行名一致

2. Description列：样品名，即样本名称，可以包括空格，非字符，展示为图为标签方便阅读理解

3. Env列：为样品来源注释列，包括上面输出的Gut、Oral、Skin、Soil，主要为采样来源

4. SourceSink列，要计算的样品标Sink，而来源数据标Source

#SampleID       Description     Env     SourceSink      Study   Details
Run20100430_H2O-1       PCR water 1     Lab 1   sink    Lab 1   PCR_water_sample_1_2010_04_30_run
Run20100430_H2O-2       PCR water 2     Lab 1   sink    Lab 1   PCR_water_sample_2_2010_04_30_run
Run20100430_H2O-3       PCR water 3     Lab 1   sink    Lab 1   PCR_water_sample_3_2010_04_30_run
BB2    Spodosol 1    Soil    source    88_Soils    NA
IT2    Spodosol 2    Soil    source    88_Soils    NA
CL3    Ultisol     Soil    source    88_Soils    NA

OTU表

示例文件：data/otus.txt

QIIME导出biom后的经典表格格式：行为OTU编号，列为样本名，矩阵对应为原始测序数量

# QIIME-formatted OTU table
#OTU ID Run20100430_ESC_C-1ss   Run20100430_ESC_C-2ss   Run20100430_ESC_C-3ss   Run20100430_ESC_C-4ss   Run20100430_ESC_C-5ss
1       0       0       0       0       0       0       0       0       0       0       0       0       0       0       0
6       0       0       0       0       0       0       0       0       0       0       0       0       0       0       0
7       0       0       0       0       0       0       0       0       0       0       0       0       0       0       0
8       0       0       0       0       0       0       0       0       0       0       0       0       0       0       0
12      0       0       0       0       0       0       0       0       0       0       0       0       0       0       0
16      0       0       0       0       0       0       0       0       0       0       0       0       0       0       0
17      0       0       0       0       0       0       0       0       0       0       0       0       0       0       0

代码中文解读

# This runs SourceTracker on the original "contamination" data set
# (data included in 'data' folder)# load sample metadata
metadata <- read.table('data/metadata.txt',sep='\t',h=T,row.names=1,check=F,comment='')# load OTU table
# This 'read.table' command is designed for a
# QIIME-formatted OTU table.
# namely, the first line begins with a '#' sign
# and actually _is_ a comment; the second line
# begins with a '#' sign but is actually the header
# 读取行、列名，跳过一行，注释为空即读#号行
otus <- read.table('data/otus.txt',sep='\t', header=T,row.names=1,check=F,skip=1,comment='')
# 读入数据框变矩阵，且转置为样本为行的旧格式
otus <- t(as.matrix(otus))# extract only those samples in common between the two tables
common.sample.ids <- intersect(rownames(metadata), rownames(otus))
otus <- otus[common.sample.ids,]
metadata <- metadata[common.sample.ids,]
# double-check that the mapping file and otu table
# had overlapping samples
# 判断是否存在共有样品，否则退出
if(length(common.sample.ids) <= 3) {message <- paste(sprintf('Error: there are %d sample ids in common ',length(common.sample.ids)),'between the metadata file and data table')stop(message)
}# extract the source environments and source/sink indices
# 筛选哪些是来源或目标真假T/F，which转化为位置编号
# 共筛选训练集180个，测试集125个
train.ix <- which(metadata$SourceSink=='source')
test.ix <- which(metadata$SourceSink=='sink')
# 测试集太多，只保留6个样品做演示
test.ix = head(test.ix)
envs <- metadata$Env
# 判断是否存在Description列，存在列保存于desc
if(is.element('Description',colnames(metadata))) desc <- metadata$Description# load SourceTracker package
# 加载软件包
source('src/SourceTracker.r')# tune the alpha values using cross-validation (this is slow!)
# 使用交叉验证调整alpha值，非常耗时
# tune.results <- tune.st(otus[train.ix,], envs[train.ix])
# alpha1 <- tune.results$best.alpha1
# alpha2 <- tune.results$best.alpha2
# note: to skip tuning, run this instead:
# 跳过优化alpha值步骤，直接设置为0.001继续计算
alpha1 <- alpha2 <- 0.001# train SourceTracker object on training data
# 基于训练集和对应描述获得预测模型
st <- sourcetracker(otus[train.ix,], envs[train.ix])# Estimate source proportions in test data
# 估计测试集中来源比例
results <- predict(st,otus[test.ix,], alpha1=alpha1, alpha2=alpha2)# Estimate leave-one-out source proportions in training data
# 在训练集中留一法（一种交叉验证方法）估计来源比例，计算次数等于训练集样本数量，极耗时
# results.train <- predict(st, alpha1=alpha1, alpha2=alpha2)# plot results
# 结果绘图, 将环境和描述列合并作为标签展示
labels <- sprintf('%s %s', envs,desc)
# 绘制饼形图比例
plot(results, labels[test.ix], type='pie')# other plotting functions
plot(results, labels[test.ix], type='bar')
plot(results, labels[test.ix], type='dist')
# plot(results.train, labels[train.ix], type='pie')
# plot(results.train, labels[train.ix], type='bar')
# plot(results.train, labels[train.ix], type='dist')# plot results with legend
# 添加图例，并人工指定颜色
plot(results, labels[test.ix], type='pie', include.legend=TRUE, env.colors=c('#47697E','#5B7444','#CC6666','#79BEDB','#885588'))
plot(results, labels[test.ix], type='pie', include.legend=TRUE, env.colors=rainbow(5))

最后结果绘图，可选饼形图(pie)、柱状图(bar)和堆叠图(dist)，如上面示例所示。

Reference

1. Knights D, Kuczynski J, Charlson ES, et al. Bayesian community-wide culture-independent microbial source tracking [J]. Nature Methods, 2011, 8(9): 761. DOI: 10.1038/nmeth.1650

2. Metcalf, J. L., et al. (2016). “Microbial community assembly and metabolic function during mammalian corpse decomposition.” Science 351(6269): 158-162.

3. Dominguez-Bello MG, De JKM, Nan S, et al. Partial restoration of the microbiota of cesarean-born infants via vaginal microbial transfer [J]. Nature Medicine, 2016, 22(3): 250-253.

4. https://www.nature.com/articles/nmeth.1650

5. https://www.ncbi.nlm.nih.gov/pubmed/21765408

6. 锐翌16S分析升级之③ SourceTracker — 寻找微生物的来源 https://mp.weixin.qq.com/s/eAD42C8ZZAcHBXWmr6HnUw

猜你喜欢

10000+：菌群分析 宝宝与猫狗 梅毒狂想曲 提DNA发Nature Cell专刊 肠道指挥大脑

系列教程：微生物组入门 Biostar 微生物组  宏基因组

专业技能：学术图表 高分文章 生信宝典 不可或缺的人

一文读懂：宏基因组 寄生虫益处 进化树

必备技能：提问 搜索  Endnote

文献阅读 热心肠 SemanticScholar Geenmedical

扩增子分析：图表解读 分析流程 统计绘图

16S功能预测   PICRUSt  FAPROTAX  Bugbase Tax4Fun

在线工具：16S预测培养基 生信绘图

科研经验：云笔记  云协作 公众号

编程模板: Shell  R Perl

生物科普:  肠道细菌 人体上的生命 生命大跃进  细胞暗战 人体奥秘

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外5000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。PI请明示身份，另有海内外微生物相关PI群供大佬合作交流。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。

学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”

点击阅读原文，跳转最新文章目录阅读