加载中...
CAZy碳水化合物活性酶预测
发表于:2021-09-30 | 分类: 生物信息
字数统计: 1.8k | 阅读时长: 8分钟 | 阅读量:

CAZy 数据库简介

CAZy 全称为 Carbohydrate-Active enZYmes Database,碳水化合物酶相关的专业数据库,内容包括能催化碳水化合物降解、修饰、以及生物合成的相关酶系家族。其包含五个主要分类:糖苷水解酶(Glycoside Hydrolases, GHs)、糖基转移酶(GlycosylTransferases, GTs)、多糖裂解酶(Polysaccharide Lyases, PLs)、糖酯酶(Carbohydrate Esterases, CEs)和氧化还原酶(Auxiliary Activities, AAs)。此外,还包含与碳水化合物结合结构域(Carbohydrate-Binding Modules, CBMs)。五大分类和一个结构域下,都分别建立了多个 Family。

  • GHs:糖苷键的水解和 / 或重排

  • GTs:糖苷键的形成

  • PLs:糖苷键的非水解裂解

  • CEs:水解碳水化合物的酯类

  • AAs:与 CAZymes 协同作用的氧化还原酶

  • CBMs:与碳水化合物结合

CAZy 数据库的准备

在进行预测之前需要准备数据库,CAZy 貌似没有提供 FASTA 格式的序列数据库,而仅提供了序列的 Assenssion number,需要我们自己从 NCBI 数据库中下载序列。下载方法参照我之前的文章《根据 assession number 批量从 NCB 下载数据》,在文章中提供了下载 CAZy 序列的方法和脚本,此处不再赘述。

在上一篇文章结尾获得的 “All.sequences.fas” 文件包含了所有的 CAZy 数据库序列,在正式预测之前需要完成数据库的格式化。后面我们将通过 Diamond 软件从基因组中预测 CAZy 蛋白,因此采用 Diamond 格式化数据库。

  • 序列预处理

    不知道什么原因,下载的序列存在两个问题,其一,下一条序列的 ID 连接着上一条序列的末尾,没有断行;其二,序列中存在着一段网页代码。因此,需要分两步进行修正。

    • 解决断行问题

      撰写脚本 “add_linebreak.pl”,内容如下:

      #!/usr/bin/perl
      use strict;
      use warnings;
      # Author: Liu hualin
      # Date: Sep 30, 2021
      
      local $/=">";
      open IN, "All.sequences.fas" || die;
      open OUT, ">CAZy.fas" || die;
      <IN>;
      while (<IN>) {
      	chomp;
      	print OUT ">$_\n";
      }
      close IN;
      close OUT;

      将脚本与 "All.sequences.fas" 放在同一目录下,在终端或者命令行中运行如下命令,得到 “CAZy.fas”。

      perl add_linebreak.pl
    • 删除无关内容

      用 EmEditor 软件打开 CAZy.fas,Ctrl+F 调出查找功能,搜索 “www.” 可以看到如下内容,手动将其删除,并保存文件。

      数据库中需要手动删除的网页信息

  • 构建 Diamond 数据库

diamond makedb --in CAZy.fas -d CAZy

开始序列比对

当然,我们选择用 Perl 进行批量比对

#!/usr/bin/perl
use strict;
use warnings;
# Author: Liu hualin
# Date: Sep 30, 2021

my @faa = glob("*.faa");# 读取所有后缀为“.faa”的文件,可以自己更改
foreach  (@faa) {
	$_=~/(.+).faa/;
	my $out = $1 . ".CAZy.diamond";
	# -p表示线程数,在笔记本上用6个即可
	system("diamond blastp -d CAZy -q $_ -e 1e-5 -f 6 -o $out -k 1 --sensitive -p 30 --query-cover 50");
}

将上述代码复制到文件中,命名为 “run_diamond_CAZy.pl”,将其和序列文件放在同一目录下,并在终端中输入如下命令,完成分析,得到 “*.CAZy.diamond”:

perl run_diamond_CAZy.pl

比对结果过滤

在比对过程中我们控制了 evalue 和 query coverage,但是没有控制 identity。但是很多时候,需要设定一个 identity 的阈值,低于阈值的比对将会被删除,该步骤可以将比对结果拷贝到 Excel 中根据 identity 排序,手动删除阈值以下的行,然而我选择用 Perl 批处理。

#!/usr/bin/perl
use strict;
use warnings;
# Author: Liu hualin
# Date: Sep 30, 2021

my @cazy = glob("*.CAZy.diamond");
foreach  (@cazy) {
	$_=~/(.+).CAZy.diamond/;
	my $out = $1 . ".CAZy.diamond.filtered";
	open IN, "$_" || die;
	open OUT, ">$out" || die;
	while (<IN>) {
		chomp;
		$_=~s/[\r\n]+//g;
		my @lines = split /\t/;
		if ($lines[2] >= 40) {
			print OUT $_ . "\n";
		}
	}
	close IN;
	close OUT;
}

将上述代码复制到文件中,命名为 “filter_cazy_diamond.pl”,将其和上一步产生的文件放在同一目录下,并在终端中输入如下命令,完成过滤,保留 identity >= 40% 的行,得到 “*.CAZy.diamond.filtered”。

perl filter_cazy_diamond.pl

你以为完了?还得 mapping!

得到的结果如下图所示,第二列的 Hits 是 NCBI 的 Assession number,我们根本只知道这是什么 CAZy 家族,因此需要 mapping!

Diamond比对结果

回头找到我们下载的 cazy_data.txt,里面保存的是 CAZy 家族与 Assession number 的对应关系。比较闲的兄弟可以用查找 - 复制 - 粘贴的方法将 “*.CAZy.diamond.filtered” 中的 Assession number 替换为 CAZy 家族。我为比较忙的兄弟准备了下面的代码,批处理。不过我输出的是一个矩阵。

#!/usr/bin/perl
use strict;
use warnings;
# Author: Liu hualin
# Date: Sep 30, 2021

my %cazy;

open IN, "cazy_data.txt" || die;
while (<IN>) {
	chomp;
	my @lines = split /\t/;
	$cazy{$lines[2]} = $lines[0];
}
close IN;

my %hash;
my %samples;
my %ids;
my @filtered = glob("*.CAZy.diamond.filtered");
foreach  (@filtered) {
	$_=~/(.+).CAZy.diamond.filtered/;
	my $sample = $1;
	$samples{$1}++;
	open IN, "$_" || die;
	while (<IN>) {
		chomp;
		my @line = split /\t/;
		if (exists $cazy{$line[1]}) {
			$ids{$cazy{$line[1]}}++;
			$hash{$sample}{$cazy{$line[1]}}++;
		}
	}
	close IN;
}

open OUT, ">CAZy.Matrix.txt" || die;

my @samples = sort keys %samples;
my @ids = sort keys %ids;

print OUT "\t" . join("\t", @samples) . "\n";
for (my $i=0; $i<@ids ;$i++) {
	print OUT $ids[$i];
	for (my $j=0; $j<@samples ;$j++) {
		if (exists $hash{$samples[$j]}{$ids[$i]}) {
			print OUT "\t$hash{$samples[$j]}{$ids[$i]}";
		}else {
			print OUT "\t0";
		}
	}
	print OUT "\n";
}
close OUT;

将上述代码复制到文件中,命名为 “assession2cazy.pl“,将其和”cazy_data.txt“,及上一步产生的文件 “*.CAZy.diamond.filtered” 放在同一目录下,并在终端中输入如下命令:

perl assession2cazy.pl

得到一个矩阵 “CAZy.Matrix.txt”,内容如下,行为 CAZy 家族,列为基因组 / 样本名。拿到本文件后,可以做热图看 CAZy 家族在各样本中的分布情况,然而这个热图将会比鞋帮子脸还要长,可读性不高,因此我选择将这些 family 合并为大类,生成一个新的矩阵。

Mapping后的矩阵

二话不说,上代码。

#!/usr/bin/perl
use strict;
use warnings;
# Author: Liu hualin
# Date: Sep 30, 2021

my %category;
my %hash;
my @samples;
my $count = 0;
open IN, "CAZy.Matrix.txt" || die;
while (<IN>) {
	$count++;
	chomp;
	if ($count == 1) {
		@samples = split;
	}else {
		my @lines = split;
		$lines[0]=~/(.+?)\d+/;
		my $cate = $1;
		$category{$cate}++;
		for (my $i=0; $i<@samples; $i++) {
			my $j = $i + 1;
			$hash{$samples[$i]}{$cate} += $lines[$j];
		}
	}
}
close IN;


open OUT, ">CAZy.Category.Matrix.txt" || die;

my @category = sort keys %category;

print OUT "\t" . join("\t", @samples) . "\n";
for (my $i=0; $i<@category ;$i++) {
	print OUT $category[$i];
	for (my $j=0; $j<@samples ;$j++) {
		if (exists $hash{$samples[$j]}{$category[$i]}) {
			print OUT "\t$hash{$samples[$j]}{$category[$i]}";
		}else {
			print OUT "\t0";
		}
	}
	print OUT "\n";
}
close OUT;

将上述代码复制到文件中,命名为 “cazyfamily2categories.pl”,将其和上一步产生的文件 “CAZy.Matrix.txt” 放在同一目录下,并在终端中输入如下命令,得到文件 “CAZy.Category.Matrix.txt”。

perl cazyfamily2categories.pl

CAZy.Category.Matrix.txt内容概览

接下来是要做柱状图还是 heatmap,就随便了。

脚本获取

关注公众号 “生信之巅”,聊天窗口回复 “9052” 获取下载链接。

生信之巅微信公众号生信之巅小程序码

敬告:使用文中脚本请引用本文网址,请尊重本人的劳动成果,谢谢!Notice: When you use the scripts in this article, please cite the link of this webpage. Thank you!

上一篇:
Hexo渲染异常
下一篇:
构建样本vs基因矩阵
本文目录
本文目录