vlambda博客
学习文章列表

SourceTracker

SourceTracker是一个基于贝叶斯方法的估计一个群落中来自另外一系列群落的物种的比例的工具。

从SourceTracker的名字就可以看出来，这个工具能够识别目的微生物群落中来自于其它来源微生物群落的比例。

实际应用场景：

• 室内空气中来自于人体不同部位的微生物比例；

• 污染源下游来自于点源排放的微生物比例；

• 新生婴儿或动物来源于母体的微生物比例；

• 肠道菌群中来源于食物的微生物比例；

• 通过食物链溯源微生物的传播过程；

• ……

结果解读

SourceTracker的结果包含3种不同的图形，用来解释目的样本中来源于不同Source的微生物比例。

• a图以饼图展示一组样本群落中不同来源微生物的比例。

• b图以条形展示一组样本群落中不同来源微生物的比例。

• c图通过面积图一组样本中单个样本不同来源微生物的比例。

分析实战

本文中的代码和作为来源的微生物群落信息来自于《Antimicrobial Chemicals Are Associated with Elevated Antibiotic Resistance Genes in the Indoor Dust Microbiome》，文章中所有的代码参见https://github.com/uo-green-lab/dust-2015。

数据准备

既然是溯源分析，那么自然就需要两个微生物群落信息，一个作为待分析的目的微生物群落，另一个作为微生物的来源信息。

用于溯源分析的目的样本的微生物群落数据使用属水平丰度表，在进行分析之前需要去除只在一个样本中出现的OTU。

⚠️一定要用抽平后的物种绝对丰度表，相对丰度不行。

用于作为来源的微生物群落数据同样使用属水平丰度表，这里使用上面提到文章中的示例文件作为微生物群落来源，文件中包括来自于人类不同组织以及室内空气和部分土壤的微生物群落数据。

接下来要制作一个mapping文件，这一步非常重要。

mapping文件的格式与QIIME要求的mapping文件基本一致，只是在最后需要添加一个SourceSink列。

该列规定了样本的类型，作为来源的样本在此处填写source，如果部分样本不希望包含在分析内，则填写NA，待分析的样本在此处填写sink，具体信息可以参考提供的示例mapping文件。

制作SourceTracker输入文件

有了sink、source和mapping文件之后，使用如下命令制作用于SourceTracker的输入文件。

## read in mapping file (present study + HMP)map.ger <- read.table("mapping_sourcetracker.txt",header = T, sep = '\t', row.names = 1, comment.char = '')## read in OTU data for GENUS (L6) level**## ** Note: this is filtered to OTUs that occur only in 2 or more samplesotu.ger.genus <- read.table("table_mc31752_sorted_L6.txt",header = T, sep = '\t', row.names = 1, skip = 1, comment.char = '')
otu.hmp.genus <- read.table('hmp_otu_table_L6_nospaces.txt',header = T, sep = '\t', row.names = 1, skip = 1, comment.char = '')

length(rownames(otu.ger.genus)[rownames(otu.ger.genus) %in% rownames(otu.hmp.genus)]) # 254length(rownames(otu.hmp.genus)[rownames(otu.hmp.genus) %in% rownames(otu.ger.genus)]) # 254## merge OTU data for present study and HMPotu.merge.genus <- merge(otu.ger.genus, otu.hmp.genus, by = 'row.names', all = TRUE)
row.names(otu.merge.genus) <- otu.merge.genus$Row.names## reduce OTU table to samples present in mapping tableotu.merge.genus <- otu.merge.genus[ ,(colnames(otu.merge.genus) %in% row.names(map.ger))]## write merged OTU table## *** NOTE *** you need to edit the output file to add a comment on the first line, and change ## the first cell of the table to '#OTU ID'. I also replaced NAs with zeros in the output table.write.table(otu.merge.genus, "sourcetracker_input.txt", sep = '\t', quote = FALSE)

得到的文件是将sink和source微生物数据合并在一起，得到一个包含所有信息的统一文件。

SourceTracker

在进行分析前，需要对上一步生成的sourcetracker输入文件进行一些修改。

• 在开始添加一行，填入# QIIME-formatted OTU table；

• 在样品名一行最开始添加一个单元格，填入#OTU ID。

接下来需要将sourcetracker_input.txt、mapping文件以及提供的SourceTracker.r原始函数文件放在同一个文件夹下。

PS：SourceTracker.r是封装在一起的SourceTracker运行及出图函数，文末提供下载方式。

首先载入mapping文件和用于分析的SourceTracker输入文件。

metadata <- read.table('mapping_sourcetracker.txt',sep='\t',h=T,row.names=1,check=F,comment='')

otus <- read.table('sourcetracker_input.txt',sep='\t', header=T,row.names=1,check=F,skip=1,comment='')
otus <- t(as.matrix(otus))
otus[is.na(otus)] <- 0

提取mapping文件和输入的分析文件中共有的样本信息。

common.sample.ids <- intersect(rownames(metadata), rownames(otus))
otus <- otus[common.sample.ids,]
metadata <- metadata[common.sample.ids,]#如果mapping文件和输入文件中没有共有样本，则会如下错误信息if(length(common.sample.ids) <= 1) {
    message <- paste(sprintf('Error: there are %d sample ids in common '),                    'between the metadata file and data table')    stop(message)
}

提取作为Source的样本分组以及样本对应的Source/Sink信息。

train.ix <- which(metadata$SourceSink=='source')
test.ix <- which(metadata$SourceSink=='sink')
envs <- metadata$Envif(is.element('Description',colnames(metadata))) desc <- metadata$Description

载入刚才提到的SourceTracker函数文件并相关设置和训练。

source('SourceTracker.r')
alpha1 <- alpha2 <- 0.001st <- sourcetracker(otus[train.ix,], envs[train.ix])

运行此步骤有可能会有warning信息，但不影响分析。

进行SourceTracker分析。

results <- predict(st,otus[test.ix,], alpha1=alpha1, alpha2=alpha2)
write.table(results$proportions,"sourcetracker_results.txt",sep = "\t",quote = FALSE)

运行过程会以“……”的方式表示运行的百分比，每个样本分别展示，运行完一个样本即可得到一个样本的结果，结果中不同来源对应的数值即为样本中来源于该Source的微生物比例。

最后对结果进行图像的绘制。

plot(results, type='pie', include.legend=TRUE, env.colors=c('#47697E','#5B7444','#CC6666','#79BEDB','#885588'))

示例的sink数据是我随便找的，与人类样本没有关系，所有结果中几乎所有的微生物来源都是Unknown。

⚠️如果需要条形图，将type设置为bar，面积图设置为dist。

扩展阅读

  
    
    
  

   
     
     
   
 
    
 
    
   
     
     
   
 
    
 
     
 
    
   
     
     
   
 
    
 
     
 
    
   
     
     
   
 
    
 
    
   
     
     
   
 
    
 
    
   
     
     
   
 
    
 
    
   
     
     
   
 
    
 
    
   
     
     
   
 
    
 
    
   
     
     
   
 
    
 
    
   
     
     
   
 
    
 
     
 
    
  
    
    
  
  
    
    
  

   
     
     
   
加群、交流和投稿