博文
技术贴 | 宏基因组Binning（七）进化分析，基因组圈图
|
本文由阿童木根据实践经验而整理，希望对大家有帮助。
原创微文，欢迎转发转载。
只做宏基因组太单调？为什么不试试宏基因组Binning呢？一次测序，“宏基因组”+“Binning”两种分析，微生太帮您一站式处理宏基因组难题。现在，微生太免费向所有人分享Binning的整套分析流程，包含：生信分析代码和R语言绘图代码。我们一共设计了7个课时，每周一次，课表（进度）如下。
图1
对Binning分析、R语言绘图感兴趣的朋友千万别错过。错过也没关系，每次课程不仅有回放，还有技术贴带您回顾课程内容。
下面我们一起来回顾第七节课的主体内容吧。
一、 分析内容与方法
1 分类学注释：
PhyloPhlAn是发表在Nature Communications上的一个用于分析微生物之间进化关系的软件。PhyloPhlAn能通过分析400种通用微生物蛋白展示微生物之间的进化关系。软件本身自带3000+种已知微生物的进化关系信息，因此PhyloPhlAn也常用于研究新发现的微生物与已知微生物进化关系。
2 进化分析：
使用PhyloPhlAn可以给物种或者Bin进行De novo进化关系分析，利用基因组数据分析某些物种或者Bin之间的进化关系。下图是Bin之间的进化关系树状图。树状图由R语言ggplot2和ggtree绘制。
3 Circos圈图：
Circos是由加拿大生物信息学科学家Martin Krzywinski利用Perl语言开发的一款可以用于描述关系型数据和多维数据的基因组圈图可视化软件。2009年Circos发表在Genome research。Circos不仅能将一个物种的整个基因组呈现在一张图片中，还可以给基因组添加丰富的注释信息，如功能注释信息，差异统计信息等等。
二、分类注释流程
1 phylophlan进化分析: insert
python phylophlan.py -i -t $project.insert --nproc $thread
输入文件：所有bin的蛋白序列文件，存放于同一个文件夹中
图2
输出文件：nwk bin进化树文件和物种注释信息
图3
2 合并insert分析物种注释结果
cat imputed_conf_* | sort -t $'\t' -k 2r | sed '1 iid\ttax' > bin_taxonomy.txt
输出文件：
图4
三、Bin进化关系分析流程
1 进化分析：de novo
python phylophlan.py -u $project.denovo --nproc $thread
输出文件：bin进化关系树文件
图6
2 制作map文件
sed '1d' ../$project.insert/bin_taxonomy.txt | sed 's/.p__/\t/g' | sed 's/.c__/\t/g' | awk '{print $1,$3}' | sed 's/ /\t/g' | sed '1 ibin\tPhylum' > map.txt
输出文件：bin门水平注释信息，用作绘图参数
图7
3 绘制BIN进化树
R CMD BATCH --args /home/cheng/huty/softwares/script_bin/tree_bins_circ.R
绘图代码：
图8
绘图结果：
图9
四、Circos绘制基因组圈图流程
文件准备主程序
bash ~/huty/softwares/script_circos/circos.sh
准备文件：取contig长度信息
grep '^k141' $i > Bin_all/Bin_circos/$base/bin.gene
# 取gene mapping信息
图10
grep '^##sequence' $i | awk '{print $2,$3,$4}' | sed 's/ /\t/g' | sed '1 icontig\tstart\tend' > Bin_all/Bin_circos/$base/bin.contig
# 取contig长度信息
图11
0. 文件：chr
R CMD BATCH --args /home/cheng/huty/softwares/script_circos/conf.chr.R
图12
1. 文件：+CDS -CDS tRNA rRNA tmRNA
awk '{if($3=="CDS" && $7=="+") print $1,$4,$5,$7}' $i/bin.gene > $i/bin.positive.cds
awk '{if($3=="CDS" && $7=="-") print $1,$4,$5,$7}' $i/bin.gene > $i/bin.negative.cds
图13
awk '{if($3=="tRNA") print $1,$4,$5,$3}' $i/bin.gene > $i/bin.tRNA
awk '{if($3=="rRNA") print $1,$4,$5,$3}' $i/bin.gene > $i/bin.rRNA
awk '{if($3=="tmRNA") print $1,$4,$5,$3}' $i/bin.gene > $i/bin.tmRNA
图14
2. 文件：GC GH CE AA PL CBM
awk '{print $1,$2}' $i | uniq | sed 's/[_|]/ /g' | awk '{print $1,$2,$4}' | sed 's/ /_/' | sed 's/[0-9]*$//' | sed '1 i id cazy' > Bin_all/Bin_circos/$bold/bin.cazy
# cazy注释结果
图15
sed 's/ID=//g' $i/bin.gene | sed 's/;/\t/g' | awk '{print $1,$4,$5,$9}' | sed '1 i contig start end id' > $i/bin.func
# gff文件
图16
R CMD BATCH --args /home/cheng/huty/softwares/script_circos/conf.merge.cazy_func.R # bin.cazy.conf
awk '{if($4 == "GH") print $0}' $i/bin.cazy.conf > $i/bin.cazy.GH
awk '{if($4 == "GT") print $0}' $i/bin.cazy.conf > $i/bin.cazy.GT
awk '{if($4 == "CE") print $0}' $i/bin.cazy.conf > $i/bin.cazy.CE
awk '{if($4 == "PL") print $0}' $i/bin.cazy.conf > $i/bin.cazy.PL
awk '{if($4 == "CBM") print $0}' $i/bin.cazy.conf > $i/bin.cazy.CBM
awk '{if($4 == "AA") print $0}' $i/bin.cazy.conf > $i/bin.cazy.AA
图17
3. GC含量计算
perl /home/cheng/huty/softwares/script_circos/contig_gc.pl $i > Bin_all/Bin_circos/$bold/bin.gc
# 计算gc含量
图18
sed -i '1 icontig\tgc\tcompare' Bin_all/Bin_circos/$bold/bin.gc
# 加header
R CMD BATCH --args /home/cheng/huty/softwares/script_circos/conf.merge.gc_contig.R # 加start end
图19
sed '1d' bin.gc.conf | awk '{if($5 == "small") print $0}' | awk '{print $1,$2,$3,$4}' > bin.gc.small
# 选择small
sed '1d' bin.gc.conf | awk '{if($5 == "big") print $0}' | awk '{print $1,$2,$3,$4}' > bin.gc.big
# 选择big
4. GC skew
sed -i '1 icontig\tvalue' bin.gc.skew
# 加header
R CMD BATCH --args /home/cheng/huty/softwares/script_circos/conf.merge.gc.skew_contig.R
# 加start end
图20
sed '1d' bin.gc.skew.conf | awk '{if($4 > 0) print $0}' > bin.gc.skew.positive
# 选择>0
sed '1d' bin.gc.skew.conf | awk '{if($4 < 0) print $0}' > bin.gc.skew.negative
# 选择<0
输出文件（准备文件）：
图21
5画circus
/home/cheng/huty/softwares/circos-0.69-9/bin/circos -conf /home/cheng/huty/softwares/script_circos/circos.conf --outputdir ../Figure -outputfile $base
6加图例
R CMD BATCH --args /home/cheng/huty/softwares/script_circos/circos.result.R
输出文件：
图22
图解：
图23
推荐一个Cicos绘图教程：
技术贴 | 手把手教你如何画Circos图https://mp.weixin.qq.com/s/j8KRzYE5FilzcMIdFlmaKQ
你可能还喜欢
2 技术贴 | 宏基因组Binning（二）质控、分箱、质检、可视化
3 技术贴 | 宏基因组Binning（三）丰度计算、差异分析
4 技术贴 | 宏基因组 Binning（四）COG EC RNA注释统计
5 技术贴 | 宏基因组Binning（五）KEGG GO注释统计
6 技术贴 | 宏基因组Binning（六）CAZyme注释统计
微生态科研学术群期待与您交流更多微生态科研问题
（联系微生态老师即可申请入群）
了解更多菌群知识，请关注“微生态”。
微信扫一扫
关注该公众号
https://blog.sciencenet.cn/blog-3474220-1285704.html
上一篇：科研 | Nature Communications：新生儿接触抗生素会通过扰乱肠道微生物的定植影响儿童六岁前的生长
下一篇：科研 | Microbiome：（国人佳作）Keystone分类群的特殊代谢功能维持了土壤微生物群落的稳定性
