微生物来源分析

[toc]

微生物来源分析

写在前面

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

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

准备

不单单是这一次，我在以后所有的分析都会将整个群落封装到phylsoeq，只是为了更好的更加灵活的对微生物群落数据进行分析，固然你们若是初次见面，可能须要安装依赖极多的phyloseq包。须要熟悉phylsoeq封装的结构和调用方法。

为了让你们更容易操做，我把数据保存为csv，方便还没有解除phylsoeq的小伙伴进行无压力测试。

微生物来源分析

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

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

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

# 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设置目标样品和来源样品的对应关系。是否不一样目标对应不一样来源样品，仍是不一样目标对应相同来源样品

##导入主函数
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)

出图，简单出一张饼图供你们参考

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个以上重复了，若是展现九个饼图，那就显得太夸张了，求均值，展现均值饼图

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的物种分类树