标签：uml   fse   adr   钢铁侠   扫描   otp   设置   bbb   nio   

微生物来源分析

写在前面

最近由于老板有分析项目，我实在是进展缓慢，一直苦恼并艰难的探索和进展，所以很长时间没有和大家见面了，今天我为大家带来的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的物种分类树

微生物来源分析

标签：uml   fse   adr   钢铁侠   扫描   otp   设置   bbb   nio   

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