CAZy碳水化合物活性酶預(yù)測

CAZy數(shù)據(jù)庫簡介

CAZy 全稱為Carbohydrate-Active enZYmes Database，碳水化合物酶相關(guān)的專業(yè)數(shù)據(jù)庫，內(nèi)容包括能催化碳水化合物降解、修飾、以及生物合成的相關(guān)酶系家族。其包含五個主要分類：糖苷水解酶（Glycoside Hydrolases, GHs）、糖基轉(zhuǎn)移酶（GlycosylTransferases, GTs）、多糖裂解酶（Polysaccharide Lyases, PLs）、糖酯酶（Carbohydrate Esterases, CEs）和氧化還原酶（Auxiliary Activities, AAs）。此外，還包含與碳水化合物結(jié)合結(jié)構(gòu)域（Carbohydrate-Binding Modules， CBMs）。五大分類和一個結(jié)構(gòu)域下，都分別建立了多個Family。

  • GHs：糖苷鍵的水解和/或重排

  • GTs：糖苷鍵的形成

  • PLs：糖苷鍵的非水解裂解

  • CEs：水解碳水化合物的酯類

  • AAs：與 CAZymes 協(xié)同作用的氧化還原酶

  • CBMs：與碳水化合物結(jié)合

CAZy數(shù)據(jù)庫的準(zhǔn)備

在進(jìn)行預(yù)測之前需要準(zhǔn)備數(shù)據(jù)庫，CAZy貌似沒有提供FASTA格式的序列數(shù)據(jù)庫，而僅提供了序列的Assenssion number，需要我們自己從NCBI數(shù)據(jù)庫中下載序列。下載方法參照我之前的文章《根據(jù)assession number批量從NCB下載數(shù)據(jù)》，在文章中提供了下載CAZy序列的方法和腳本，此處不再贅述。

在上一篇文章結(jié)尾獲得的“All.sequences.fas”文件包含了所有的CAZy數(shù)據(jù)庫序列，在正式預(yù)測之前需要完成數(shù)據(jù)庫的格式化。后面我們將通過Diamond軟件從基因組中預(yù)測CAZy蛋白，因此采用Diamond格式化數(shù)據(jù)庫。

  • 序列預(yù)處理

    不知道什么原因，下載的序列存在兩個問題，其一，下一條序列的ID連接著上一條序列的末尾，沒有斷行；其二，序列中存在著一段網(wǎng)頁代碼。因此，需要分兩步進(jìn)行修正。

    • 解決斷行問題

      撰寫腳本“add_linebreak.pl”，內(nèi)容如下：

      #!/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

    • 刪除無關(guān)內(nèi)容

      用EmEditor軟件打開CAZy.fas，Ctrl+F調(diào)出查找功能，搜索“www.” 可以看到如下內(nèi)容，手動將其刪除，并保存文件。

  • 構(gòu)建Diamond數(shù)據(jù)庫

diamond?makedb?--in?CAZy.fas?-d?CAZy

開始序列比對

當(dāng)然，我們選擇用Perl進(jìn)行批量比對

#!/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表示線程數(shù)，在筆記本上用6個即可
?system("diamond?blastp?-d?CAZy?-q?$_?-e?1e-5?-f?6?-o?$out?-k?1?--sensitive?-p?30?--query-cover?50");
}

將上述代碼復(fù)制到文件中，命名為“run_diamond_CAZy.pl”，將其和序列文件放在同一目錄下，并在終端中輸入如下命令，完成分析，得到“*.CAZy.diamond”：

perl?run_diamond_CAZy.pl

比對結(jié)果過濾

在比對過程中我們控制了evalue和query coverage，但是沒有控制identity。但是很多時候，需要設(shè)定一個identity的閾值，低于閾值的比對將會被刪除，該步驟可以將比對結(jié)果拷貝到Excel中根據(jù)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;
}

將上述代碼復(fù)制到文件中，命名為“filter_cazy_diamond.pl”，將其和上一步產(chǎn)生的文件放在同一目錄下，并在終端中輸入如下命令，完成過濾，保留identity >= 40% 的行，得到“*.CAZy.diamond.filtered”。

perl?filter_cazy_diamond.pl

你以為完了？還得mapping！

得到的結(jié)果如下圖所示，第二列的Hits是NCBI的Assession number，我們根本只知道這是什么CAZy家族，因此需要mapping！

回頭找到我們下載的cazy_data.txt，里面保存的是CAZy家族與Assession number的對應(yīng)關(guān)系。比較閑的兄弟可以用查找-復(fù)制-粘貼的方法將“*.CAZy.diamond.filtered”中的Assession number替換為CAZy家族。我為比較忙的兄弟準(zhǔn)備了下面的代碼，批處理。不過我輸出的是一個矩陣。

#!/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;

將上述代碼復(fù)制到文件中，命名為“assession2cazy.pl“，將其和”cazy_data.txt“，及上一步產(chǎn)生的文件“*.CAZy.diamond.filtered”放在同一目錄下，并在終端中輸入如下命令：

perl?assession2cazy.pl

得到一個矩陣“CAZy.Matrix.txt”，內(nèi)容如下，行為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;

將上述代碼復(fù)制到文件中，命名為“cazyfamily2categories.pl”，將其和上一步產(chǎn)生的文件“CAZy.Matrix.txt”放在同一目錄下，并在終端中輸入如下命令，得到文件“CAZy.Category.Matrix.txt”。

perl?cazyfamily2categories.pl

接下來是要做柱狀圖還是heatmap，就隨便了。

腳本獲取

“生信之巔” GZH。

敬告：使用文中腳本請引用本文網(wǎng)址，請尊重本人的勞動成果，謝謝！Notice: When you use the scripts in this article, please cite the link of this webpage. Thank you!

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容