文档章节

微生物来源分析

o
 osc_g8254g7s
发布于 2019/08/19 18:03
字数 1792
阅读 26
收藏 0

精选30+云产品,助力企业轻松上云!>>>

[toc]

微生物来源分析

写在前面

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

fig1

刘老师对这片文章进行了详细的解读: 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)

image

就正常样品而言,我们都会测定重复,这里基于多个样品的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)

image

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


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()

image

image

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

我们作为生物可能测定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()

image

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

image

历史目录

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的物种分类树
o
粉丝 0
博文 500
码字总数 0
作品 0
私信 提问
加载中
请先登录后再评论。
Nature综述:Rob Knight等大佬手把手教你开展菌群研究

自然微生物综述(2017 IF:31.851)于2018年5月23日在线发表了Rob Knight亲自撰写(一作兼通讯)的微生物组领域研究方法综述,不仅系统总结了过去,更为未来3-5年内本领域研究方法的选择,提供了...

osc_a6tqrvul
2018/11/12
4
0
浙大绘制首个地球微生物“社会关系”网络

来源:浙江大学 图片由课题组提供 单个微生物看不见、摸不着,但却无时不在、无处不在。但微生物的功能绝非“分解者”这么简单,影响到温室气体的、绿色生产的、人体健康的方方面面,其群落组...

人工智能学家
06/25
0
0
57种药物成分可被肠道菌改变,同类研究已在国内展开丨专家解读

     为什么在服用药物后,有些人疗效好,有些人疗效不佳、甚至会产生副作用?要回答这个问题,或许要首先考虑药物在人体内长时间停留的地方——肠道。   此前,已经有若干个研究发现了...

osc_g91p39eg
06/14
6
0
57种药物成分可被肠道菌改变,同类研究已在国内展开丨专家解读

     为什么在服用药物后,有些人疗效好,有些人疗效不佳、甚至会产生副作用?要回答这个问题,或许要首先考虑药物在人体内长时间停留的地方——肠道。   此前,已经有若干个研究发现了...

DeepTech深科技
06/13
0
0
想知道自己的自然寿命?未来AI或许可以帮到你

云栖号:https://yqh.aliyun.com 第一手的上云资讯,不同行业精选的上云企业案例库,基于众多成功案例萃取而成的最佳实践,助力您上云决策! 预估自己的自然寿命,或许不再是一个奇幻的命题,...

云栖号资讯小编
02/20
0
0

没有更多内容

加载失败,请刷新页面

加载更多

科技人文丨玻璃心:承受阈值与表达

大家好,我是SKODE。 有趣的灵魂,聊科技人文。 本系列博客地址:传送门 本文转载自B站:安慰记传送门 玻璃心是网络用语,意思是: 对负面事件的接受度很低 还有对别人可能给出的负面评价非常...

osc_u9mt0sus
57分钟前
20
0
迅睿CMS 游客不允许上传附件

游客不允许上传附件 迅睿CMS系统:https://www.xunruicms.com/ 本文档原文地址:https://www.xunruicms.com/doc/752.html...

迅睿CMS-PHP开源CMS程序
57分钟前
12
0
代理,注解,接口和实现类的小测验

* retention : 保留* policy : 策略 ps : 简单测试了一下手写代理,jdk动态代理,和cglib动态代理,根据其不同的结果分析其原理 一:测试目的 主要想看一下不同的代理模式对于目标类中方法上注...

岁一
58分钟前
12
0
V-Ray 5 For 3ds Max 正式发布:超越渲染 - 知乎

15个新功能,V-Ray5助你时间更节省,渲染更出色! 作者:ChaosGroup VRay 5 For 3ds Max 已正式发布! 2分钟视频,抢先预览新功能↓ 知乎视频 www.zhihu.com V-Ray 5 for 3ds Max 新增功能 ...

osc_o9u1um45
58分钟前
0
0
毕业的笑容和悲伤永远是校园的回忆

校园的风轻轻的拂过我的脸庞,风儿显得更加凉爽, 开满火红的凤凰树,染遍了校园的每个角落, 晚上那枝头蝉儿的竞相鸣奏,唱满了令人不舍的毕业歌, 它们彷彿告诉了我们要毕业了。 毕业典礼那...

瑾123
58分钟前
7
0

没有更多内容

加载失败,请刷新页面

加载更多

返回顶部
顶部