技术贴 | 宏基因组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