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。
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！
回头找到我们下载的 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 合并为大类，生成一个新的矩阵。
二话不说，上代码。
#!/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
接下来是要做柱状图还是 heatmap，就随便了。
脚本获取
关注公众号 “生信之巅”，聊天窗口回复 “9052” 获取下载链接。
 |  |
敬告：使用文中脚本请引用本文网址，请尊重本人的劳动成果，谢谢！
