目录

• 微生物来源分析
  • 写在前面
  • 准备
  • 微生物来源分析
• rm(list = ls())
• gc()
  • 导入主函数
• 导入分组文件和OTU表格
• Load OTU table
  • 下面区分目标样品和来源样品。
• Extract the source environments and source/sink indices
• 对两组样品进行抽平
• Estimate source proportions for each sink
  • 就正常样品而言，我们都会测定重复，这里基于多个样品的sourceracker分析
  • 导入主函数
• 导入分组文件和OTU表格
• Load OTU table
  • 出图，简单出一张饼图供大家参考
• layouts = as.character(unique(metadata$SampleType))
  • 基于多个重复，我们合并饼图展示
  • 历史目录
    • R语言分析技术
    • 扩增子专题
    • 基于phyloseq的微生物群落分析
    • 代谢组专题
    • 当科研遇见python
    • 科学知识图谱
    • 杂谈

微生物来源分析

写在前面

最近由于老板有分析项目，我实在是进展缓慢，一直苦恼并艰难的探索和进展，所以很长时间没有和大家见面了，今天我为大家带来的source tracker分析，使用前一段时间刚出来的工具FEAST。

刘老师对这片文章进行了详细的解读：
Nature Methods：快速准确的微生物来源追溯工具FEAST。跟着刘老师的步伐，今天我对这个工具进行一个尝试。为什么作者不将这个工具封装到R包呢这样不就更容易了吗？可能好多小伙伴都没有从github上克隆过项目。

准备

不仅仅是这一次，我在之后全部的分析都会将整个群落封装到phylsoeq，只是为了更好的更加灵活的对微生物群落数据进行分析，当然大家如果初次见面，可能需要安装依赖极多的phyloseq包。需要熟悉phylsoeq封装的结构和调用方法。

为了让大家更容易操作，我把数据保存为csv，方便尚未解除phylsoeq的小伙伴进行无压力测试。

微生物来源分析

FEAST提供两种方式来做微生物来源分析。

1. 基于单个目标的来源。单个样品的来分析。
   2.基于多个目标和多个来源。多个样品进行来源分析。

首先我们来演示基于单个目标样品和来源样品的来源分析

```{r cars1212}

rm(list = ls())

gc()

path = "./phyloseq_7_source_FEAST"
dir.create(path)

导入主函数

source("./FEAST-master/FEAST_src//src.R")

ps = readRDS("./a3_DADA2_table/ps_OTU_.ps")

导入分组文件和OTU表格

metadata <- as.data.frame(sample_data(ps))
head(metadata)

write.csv(metadata,"metadata.csv",quote = F)

Load OTU table

vegan_otu <- function(physeq){
OTU <- otu_table(physeq)
if(taxa_are_rows(OTU)){
OTU <- t(OTU)
}
return(as(OTU,"matrix"))
}
otus <- as.data.frame(t(vegan_otu(ps)))
write.csv(otus,"otus.csv",quote = F)
otus <- t(as.matrix(otus))

下面区分目标样品和来源样品。

envs <- metadata$SampleType

metadata<- arrange(metadata, SampleType)
metadata$id = rep(1:6,4)
Ids <- na.omit(unique(metadata$id))
it = 1

train.ix <- which(metadata$SampleType%in%c("B","C","D")& metadata$id == Ids[it])
test.ix <- which(metadata$SampleType==‘A‘ & metadata$id == Ids[it])

Extract the source environments and source/sink indices

num_sources <- length(train.ix) #number of sources
COVERAGE = min(rowSums(otus[c(train.ix, test.ix),])) #Can be adjusted by the user

对两组样品进行抽平

sources <- as.matrix(rarefy(otus[train.ix,], COVERAGE))
sinks <- as.matrix(rarefy(t(as.matrix(otus[test.ix,])), COVERAGE))

dim(sinks)
print(paste("Number of OTUs in the sink sample = ",length(which(sinks > 0))))
print(paste("Seq depth in the sources and sink samples = ",COVERAGE))
print(paste("The sink is:", envs[test.ix]))

Estimate source proportions for each sink

EM_iterations = 1000 # number of EM iterations. default value

FEAST_output<-FEAST(source=sources, sinks = t(sinks), env = envs[train.ix], em_itr = EM_iterations, COVERAGE = COVERAGE)
Proportions_est <- FEAST_output$data_prop[,1]
names(Proportions_est) <- c(as.character(envs[train.ix]), "unknown")

print("Source mixing proportions")
Proportions_est
round(Proportions_est,3)

```

就正常样品而言，我们都会测定重复，这里基于多个样品的sourceracker分析

基于多个目标和来源的微生物来源分析:
different_sources_flags设置目标样品和来源样品的对应关系。是否不同目标对应不同来源样品，还是不同目标对应相同来源样品

```{r more step}

导入主函数

source("./FEAST-master/FEAST_src//src.R")

ps = readRDS("./a3_DADA2_table/ps_OTU_.ps")

导入分组文件和OTU表格

metadata <- as.data.frame(sample_data(ps))
head(metadata)

Load OTU table

vegan_otu <- function(physeq){
OTU <- otu_table(physeq)
if(taxa_are_rows(OTU)){
OTU <- t(OTU)
}
return(as(OTU,"matrix"))
}
otus <- as.data.frame(t(vegan_otu(ps)))
otus <- t(as.matrix(otus))

head(metadata)

metadata<- arrange(metadata, SampleType)
metadata$id = rep(1:6,4)
EM_iterations = 1000 #default value
different_sources_flag = 1

envs <- metadata$SampleType
Ids <- na.omit(unique(metadata$id))
Proportions_est <- list()
it = 1

for(it in 1:length(Ids)){

# Extract the source environments and source/sink indices
if(different_sources_flag == 1){

train.ix <- which(metadata$SampleType%in%c("B","C","D")& metadata$id == Ids[it])
test.ix <- which(metadata$SampleType=='A'  & metadata$id == Ids[it])

}

else{

train.ix <- which(metadata$SampleType%in%c("B","C","D"))
test.ix <- which(metadata$SampleType=='A' & metadata$id == Ids[it])

}

num_sources <- length(train.ix)
COVERAGE = min(rowSums(otus[c(train.ix, test.ix),])) #Can be adjusted by the user

# Define sources and sinks

sources <- as.matrix(rarefy(otus[train.ix,], COVERAGE))
sinks <- as.matrix(rarefy(t(as.matrix(otus[test.ix,])), COVERAGE))

print(paste("Number of OTUs in the sink sample = ",length(which(sinks > 0))))
print(paste("Seq depth in the sources and sink samples = ",COVERAGE))
print(paste("The sink is:", envs[test.ix]))

# Estimate source proportions for each sink

FEAST_output<-FEAST(source=sources, sinks = t(sinks), env = envs[train.ix], em_itr = EM_iterations, COVERAGE = COVERAGE)
Proportions_est[[it]] <- FEAST_output$data_prop[,1]

names(Proportions_est[[it]]) <- c(as.character(envs[train.ix]), "unknown")

if(length(Proportions_est[[it]]) < num_sources +1){

tmp = Proportions_est[[it]]
Proportions_est[[it]][num_sources] = NA
Proportions_est[[it]][num_sources+1] = tmp[num_sources]

}

print("Source mixing proportions")
print(Proportions_est[[it]])

}

print(Proportions_est)

went = as.data.frame(Proportions_est)
colnames(went) = paste("repeat_",unique(metadata$id),sep = "")
head(went)

filename = paste(path,"/FEAST.csv",sep = "")
write.csv(went,filename,quote = F)

```

出图，简单出一张饼图供大家参考

```{r ming}

library(RColorBrewer)
library(dplyr)
library(graphics)

head(went)

plotname = paste(path,"/FEAST.pdf",sep = "")
pdf(file = plotname,width = 12,height = 12)
par(mfrow=c((length(unique(metadata$SampleType))%/%2 +2 ),2), mar=c(1,1,1,1))

layouts = as.character(unique(metadata$SampleType))

for (i in 1:length(colnames(went))) {

labs <- paste0(row.names(went)," \n(", round(went[,i]/sum(went[,i])*100,2), "%)")

pie(went[,i],labels=labs, init.angle=90,col = brewer.pal(nrow(went), "Reds"),
border="black",main =colnames(went)[i] )
}

dev.off()

```

基于多个重复，我们合并饼图展示

我们作为生物可能测定9个以上重复了，如果展示九个饼图，那就显得太夸张了，求均值，展示均值饼图

```{r cars123}
head(went)

asx = as.data.frame(rowMeans(went))

asx = as.matrix(asx)
asx_norm = t(t(asx)/colSums(asx)) #* 100 # normalization to total 100
head(asx_norm)

plotname = paste(path,"/FEAST_mean.pdf",sep = "")
pdf(file = plotname,width = 6,height = 6)
labs <- paste0(row.names(asx_norm)," \n(", round(asx_norm[,1]/sum(asx_norm[,1])*100,2), "%)")

pie(asx_norm[,1],labels=labs, init.angle=90,col = brewer.pal(nrow(went), "Reds"),
border="black",main = "mean of source tracker")
dev.off()

```

欢迎加入微生信生物讨论群，扫描下方二维码添加小编微信，小编带你入伙啦，大牛如云，让交流变得简单。

历史目录

R语言分析技术

1. 《扩增子16s核心OTU挑选-基于otu_table的UpSet和韦恩图》
2. 《分类堆叠柱状图顺序排列及其添加合适条块标签》
3. 《R语言绘制带有显著性字母标记的柱状图》
4. 《使用R实现批量方差分析(aov)和多重比较(LSD)》
5. Rstudio切换挂载R版本及本地安装一些包
6. 玩转R包
7. 用ggplot画你们心中的女神
8. ggplot版钢铁侠
9. 想用ggplot做一张致谢ppt，但是碰到了520，画风就变了

扩增子专题

1. 《16s分析之Qiime数据整理+基础多样性分析》
2. 《16s分析之差异展示（热图）》
3. 迅速提高序列拼接效率，得到后续分析友好型输入，依托qiime
4. https://mp.weixin.qq.com/s/6zuB9JKYvDtlomtAlxSmGw》
5. 16s分析之网络分析一(MENA)
6. 16s分析之进化树+差异分析（一）
7. 16s分析之进化树+差异分析（二）
8. Qiime2学习笔记之Qiime2网站示例学习笔记
9. PCA原理解读
10. PCA实战
11. 16s分析之LEfSe分析
12. 16s分析之FishTaco分析
13. PICRUSt功能预测又被爆出新的问题啦！

基于phyloseq的微生物群落分析

1. 开年工作第一天phyloseq介绍
2. phyloseq入门
3. R语言微生物群落数据分析:从原始数据到结果（No1）
4. R语言微生物群落数据分析:从原始数据到结果（No2））
5. phyloseq进化树可视化
6. 基于phyloseq的排序分析

代谢组专题

1. 非靶向代谢组学数据分析连载（第零篇引子）
2. 非靶向代谢组学分析连载（第一篇：缺失数据处理、归一化、标准化）

当科研遇见python

1.当科研遇见python

科学知识图谱

1.citespace 的基本术语与下载安装

杂谈

1. 我的生物信息之路
2. graphlan进一步学习笔记之进化树
3. 如约 为大家带来了graphlan的物种分类树

原文地址：https://www.cnblogs.com/wentaomicro/p/11378682.html