失眠网 > 技术贴 | 宏基因组Binning（三）丰度计算差异分析

技术贴 | 宏基因组Binning（三）丰度计算差异分析

时间：2023-06-03 09:14:23

点击蓝字↑↑↑“微生态”，轻松关注不迷路

本文由阿童木根据实践经验而整理，希望对大家有帮助。

原创微文，欢迎转发转载。

只做宏基因组太单调？为什么不试试宏基因组Binning呢？一次测序，“宏基因组”+“Binning”两种分析，微生太帮您一站式处理宏基因组难题。现在，微生太免费向所有人分享Binning的整套分析流程，包含：生信分析代码和R语言绘图代码。我们一共设计了7个课时，每周一次，课表（进度）如下。

对Binning分析、R语言绘图感兴趣的朋友千万别错过。错过也没关系，每次课程不仅有回放，还有技术贴带您回顾课程内容。

下面我们一起来回顾第三节课的主体内容吧。

一、课程内容

图2

说明：黑虚线框住的是上节课的内容，红虚线框住的是本节课将讲解的内容。

二、输入文件

图3

说明：1号文件是质控后的序列文件；2号文件是clean reads组装得到的contig序列文件；3号文件是contig分箱得到的初始bins再经质控得到的clean bins；4号文件是样品分组信息表。

三、分析方法

首先，用metawrap中的quant_bins模块（salmon算法），估计每个样品中每个contig的丰度，计算Bin平均丰度，获得样品-Bin丰度表。Bin丰度是“每百万序列基因组拷贝数”，是已经标准化的数据，类似RNAseq分析中的TPM（每百万转录本）。

接着，用R语言pheatmap绘图函数绘制样品-Bin丰度热图。热图可以以色块颜色深浅的方式表达Bin丰度的大小，还可以进行Bin-Bin聚类和样品-样品聚类，相似的丰度模式会被聚到一起。然后，统计每个Bin的丰度总数，用ggplot geom_bar绘制柱形图并排序，展示整批数据中所有Bin的丰度情况。如果客户提供的样品分组>=2且每组样品在三个以上，那么还可以用lefse进行组间差异分析寻找与分组有关的Bin。

四、分析流程

#1.计算bin丰度metawrapquant_bins-t$thread-oBin_all/Bin_quant/-bBin_all/Bin_pick/-aBin_all/Assembly_out/final.contigs.faBin_all/Clean_>2.绘制bin丰度热图RCMDBATCH--args/home/cheng/huty/softwares/script_bin/quant.pheatmap.Rconvert-density400-quality200Bin_all/Bin_quant/bin_abundance_pheatmap.pdfBin_all/Bin_quant/bin_abundance_pheatmap.png#3.绘制bin丰度柱形图RCMDBATCH--args/home/cheng/huty/softwares/script_bin/quant.bar.Rconvert-density400-quality200Bin_all/Bin_quant/bin_abundance_bar.pdfBin_all/Bin_quant/bin_abundance_bar.png#4.丰度lefse分析python3/home/cheng/huty/softwares/script_bin/quant_lefse.py$metadata$all_group/home/cheng/huty/databases/group_color.list

4.1 计算bin丰度

图4

说明：使用metawrap quant_bins计算Bin丰度，得到结果丰度表如上。

4.2 绘制bin丰度热图

R语言代码：

#!/usr/bin/envRscriptlibrary(pheatmap)>,header=T,row.names=1,sep="\t")pheatmap(data[1:20,],filename="bin_abundance_pheatmap.pdf",scale="column",cellwidth=20,cellheight=20)pheatmap(data[1:20,],filename="bin_abundance_pheatmap.png",scale="column",cellwidth=20,cellheight=20)

绘图结果：

图5

4.3 绘制bin丰度柱形图

R语言代码：

#!/usr/bin/envRscriptlibrary(ggplot2)>,header=T,sep="\t")data$sum=apply(data[,2:length(data[1,])],1,FUN=sum)#求和，加到列尾data=data[,c("Genomic.bins","sum")]#取sum信息data=data[order(data$sum,decreasing=T),]#倒序write.table(data,file="bin_abundance_sum.txt",quote=F,sep="\t",row.names=F)data$color=paste("color",1:length(data[,1]),sep="_")result=ggplot(data,aes(x=data[,1],y=data[,2],fill=color))+geom_bar(stat="identity",width=0.5)+labs(x="",y="Binabundance")+scale_y_continuous(expand=c(0,0))+scale_x_discrete(limits=factor(data[,1]))+theme(legend.position="none",axis.text.x=element_text(angle=45,hjust=1))+theme(panel.grid=element_blank(),panel.background=element_rect(color="black",fill="transparent"))wide=round(length(data[,1])/40)*7ggsave(result,filename="bin_abundance_bar.pdf",width=wide)

绘图结果：

图6

4.4 丰度LEfSe差异分析

绘图方法：

1 无需转发，学术课程免费领取 | 使用LEFSE寻找在组间有差异的特征微生物

2 LEfSE网页版：http://huttenhower.sph.harvard.edu/galaxy/

3 LEfSE命令行：/biobakery/biobakery/wiki/lefse

4 微生太云平台（即将发布）

绘图结果：

图7

说明：丰度表结合分组信息后可进行Lefse bIn丰度差异分析寻找组将差异Bins。分析结果可进行LDA柱形图绘制，展示差异bins的名称和LDA score值。

你可能还喜欢