CAZy碳水化合物活性酶预测

作者: 了尘兰若 | 来源:发表于2022-04-21 20:15 被阅读0次

    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!

    回头找到我们下载的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,就随便了。

    脚本获取

    “生信之巅” GZH。

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

    相关文章

      网友评论

        本文标题:CAZy碳水化合物活性酶预测

        本文链接:https://www.haomeiwen.com/subject/vwdwertx.html