目的

細菌や古細菌のゲノム上では、協働する遺伝子が近接して存在し、遺伝子クラスターを形成していることが多い（例：ラクトースオペロン）。

そのため、ある遺伝子の機能推定のために周辺遺伝子を調べることは有効であり、私の所属する研究室ではその作業はほぼルーティンワーク化している。
そこで、コードを書いてルーチンワークの自動化を行いたいと考えた。

また、少数のバクテリアのゲノムだけではその普遍性を担保できないため、周辺遺伝子を網羅的に解析したい。（遺伝子クラスターを形成していることが「多い」とはいえ、関係のなさそうな遺伝子がクラスター内にあったり、関係ありそうな遺伝子がなかったりする）

私の所属する研究室では、糖質関連酵素（CAZyme）の研究を行っているため、CAZyme 遺伝子クラスターを解析するためのコードを書いた。

スキーム

糖質関連酵素のアノテーションに特化したプログラムである dbCAN、CUPP と、タンパク質のアノテーションに広く用いられているHMMER によってアノテーションを行い、同じクラスターに存在しているたんぱく質を最後に集計する。

テストデータ

まず、この手法の有効性を調べるために、すでに遺伝子クラスターを形成していることが知られている、ガラクトマンナンを分解する酵素をターゲットとして解析を行うことにした。

ガラクトマンナンは、マンノースとガラクトースから成る多糖で、ガラクトマンナンの分解には、α-ガラクトシダーゼ（図中黄色い丸を切り離す酵素）と β-マンナナーゼ（図中緑色の丸を切り離す酵素）が関与していることが知られている¹。

α-ガラクトシダーゼも β-マンナナーゼも多くの種類が存在するが、ここではデータベースに登録されている遺伝子数が手頃だった糖質加水分解酵素ファミリー36（GH36）のα-ガラクトシダーゼをターゲットとする。

バクテリアゲノムの取得

GH36ファミリーは、同じファミリー内に、基質特異性（機能）が異なる複数の酵素が存在している。
ガラクトマンナンの分解に関わる α-ガラクトシダーゼのみを抽出するために、CUPP²プログラムによるサブファミリー分類を採用した。

CUPP の web サイトで、論文にある GH27 酵素が属するサブファミリーを調べ（GH36:9.1 だった）、そのサブファミリーに属する酵素の表からバクテリアの種名を抽出した。

NCBI_genome_download

遺伝子の取得には、NCBI Genome Downloading Scripts を利用した。（私は conda の方を使ってインストールした。）

注意点としては、最初に --dry-run オプションをつけて実行し、事前にダウンロードする
ファイル数を確認した方が良い（ファイル数に比例して実行に時間がかかる）。

ncbi-genome-download コマンド例


ncbi-genome-download --genera /user/genome/species.txt --assembly-level complete --formats fasta --dry-run bacteria

また、ダウンロードしたゲノムのファイル名（GCF_000009685.1 などの名前が付く）と種名の対応が分からなくならないように、dry-run コマンドの最後に > file_neme.txt をつけて標準出力の内容をファイルに記録しておいた。（今回の解析ではこのファイルも後で使った。）

抽出したバクテリアの種名をクエリとして、対象を bacteria の complete ゲノムに絞ってウンロードした。

ファイルの解凍

FASTA形式でダウンロードしたファイルを解凍するためにシェルスクリプトを作成した。

.zip ファイルが存在すれば重複してダウンロードしないので、元のファイルを残してディレクトリ中のすべてのファイルを解凍するようにした。

#!/bin/sh

dir_path="path/refseq/bacteria/*"#解凍するファイルがあるフォルダのパス
dirs=`find $dir_path -maxdepth 1 -type f -name *.gz`

for dir in $dirs;
do
    if [ ! -e ${dir%.*} ]; then
       echo $dir
       gzip -d -k $dir
    else
        print="$dir_path already extracted."
        echo $print
    fi
done

解凍するファイルがあるフォルダのパスを自分の環境に合わせて書き換えて、
bash unzip.sh で実行する。

ファイルの確認

解凍すると、FASTAファイルは ”〇〇_genomic.fna”　という名前になっていた。
フォルダ中の ”〇〇_genomic.fna”　の数を調べると、すべて１個だった。複数のコンティグに分かれているものも１つのファイルにまとめられているらしい。

翻訳とアノテーション

dbCAN

ローカル版の dbCAN （run_dbcan4）を GitHub からインストールした。
大量のゲノムを dbCAN に流し込みたいので、シェルスクリプトを作成した。（2~14 行目まではヘルプを表示するための部分。）ncbi-genome-download の標準出力のファイルを引数として、ファイル中にあるディレクトリを dbCAN のコマンドの引数とする。
ニッチ過ぎて需要がある気がしないが、一応コードを載せておく……。（16 行目と 22 行目のパスを自分の環境に合わせて書き換えるのと、実行前に conda activate でdbCAN の仮想環境を立ち上げておく必要がある。）

dbcan.sh

#!/bin/sh
#
# dbcan.sh
# ============
# 
# Usage:
# 
# 	dbcan.sh file.txt
# 

#Function help() shows help
if [ "$1" = "--help" ] || [ "$1" = "-h" ] || [ -z "$1" ]; then
cat $0 | grep '^\s*#' |grep -v '^\s*#!' | sed -e "s/^# //g" | grep -v '^\s*#'
fi

dir_path="/home/user/refseq/bacteria" #ncbi-genome-download の結果が出力されるディレクトリ
DATA=`tail -n +2 $1` #１行目を飛ばして２行目から読み取る
	
while read line
do
	dir=${line:0:15}
	outdir="/home/user/genecluster_ver2/dbcan/$dir" #結果を出力するディレクトリ＋/$dir
	echo $outdir
#アウトプットのフォルダがなければdbcanを実行
	if [ ! -d $outdir ]; then
		COMMAND="run_dbcan $dir_path/$dir/*genomic.fna prok -c cluster --out_dir $outdir"
    		echo $COMMAND
   		eval ${COMMAND}
   	else
        echo "$dir already exists"
    	fi

done << FILE
$DATA
FILE

私の環境では、ゲノム一つあたり大体 15 ~ 20 分位掛かった。

CUPP

ローカル版の CUPP もここからダウンロードできる。
自分の OS に合ったファイルをダウンロードして展開する。
実行する時は、CUPP ファイルがあるディレクトリに移動して ./cupp と打つか、そのディレクトリにパスを通す。

ニッチ過ぎて需要ry な、dbCANで翻訳されたタンパク質をCUPPに流し込むシェルスクリプト。これもncbi-genome-download のファイルを引数とする。（24、30 行目のパスを書き換える必要がある。）

cupp.sh

#!/bin/sh
#
# cupp.sh
# ============
# 
# Usage:
# 
# 	cupp.sh file.txt
# 
# file.txt	stdout of NCBI download dry-run.
# 
#Function help() shows help
if [ "$1" = "--help" ] || [ "$1" = "-h" ] || [ -z "$1" ]; then
cat $0 | grep '^\s*#' |grep -v '^\s*#!' | sed -e "s/^# //g" | grep -v '^\s*#'
fi

dir_path="/home/user/genecluster_ver2/dbcan" #dbCANの結果を出力したディレクトリ
DATA=`tail -n +2 $1`

eval cd /home/user/x86_64-unknown-linux-gnu #パスを通していないのでディレクトリを移動した
while read line
do
	dir=${line:0:15}
    outdir="/home/user/genecluster_ver2/cupp" #cuppの結果を出力するディレクトリ
	isoutput="$outdir/$dir/uniInput.fasta.tsv"
#アウトプットのファイルがなければ実行
	if [ ! -e $isoutput ]; then
	cp $dir_path/$dir/uniInput $dir_path/$dir/uniInput.fasta
	mkdir -p $outdir/$dir
	COMMAND="./cupp --query $dir_path/$dir/uniInput.fasta -c $outdir/$dir --library /home/user/lib-20230602T055732Z-001/lib/2023_CUPP_lib_CUPPlibrary.lib" #--libraryのパスを自分の環境に合わせて書き換える
    	echo $COMMAND
		eval ${COMMAND}
	else
		echo "$dir already exists"
	fi

done << FILE
$DATA
FILE

CUPPは計算がすごく早いが、ライブラリを読み込むのに時間が掛かる。毎回ライブラリを読み込むのではなく、FASTAファイルをまとめて処理できると高速化できる気がする。

HMMER

このサイトをほぼそのまま参考にさせて頂き、HMMERのインストールとPfam-Aのデータベースの準備をした。

これもシェルスクリプトを書いてdbCANで翻訳した配列をHMMERに流し込んだ。
書いたというか、ほぼ上のコードを流用した。

hmmscan.sh

#!/bin/sh
#
# hmmscan.sh
# ============
# 
# Usage:
# 
# 	hmmscan.sh file.txt
# 
# file.txt	stdout of NCBI download dry-run.
# 
#Function help() shows help
if [ "$1" = "--help" ] || [ "$1" = "-h" ] || [ -z "$1" ]; then
cat $0 | grep '^\s*#' |grep -v '^\s*#!' | sed -e "s/^# //g" | grep -v '^\s*#'
fi

dir_path="/home/user/genecluster_ver2/dbcan" #dbCANの結果を出力したディレクトリ
DATA=`tail -n +2 $1`

eval cd hmmer #hmmer のディレクトリに移動
while read line
do
	dir=${line:0:15}
    outdir="/home/user/genecluster_ver2/hmmer" #hmmerの結果を出力するディレクトリ
	isoutput="$outdir/$dir/hmm_out.txt"
#フォルダがなければ実行
	if [ ! -e $isoutput ]; then
	mkdir -p $outdir/$dir
	COMMAND="hmmscan --tblout $outdir/$dir/hmm_out.txt --cpu 6 -E 1e-5 Pfam-A.hmm $dir_path/$dir/uniInput"
    		echo $COMMAND
		eval ${COMMAND}
	else
		echo "$dir already exists"
	fi

done << FILE
$DATA
FILE

結果の集計と考察

Python を駆使して、データを集計した。
その過程で、もともとのデータがシークエンスがよく読まれている生物のゲノムの割合が高く、偏っていると感じた。（例えば、大腸菌が例外的な組み合わせの遺伝子クラスターを持っていた場合、ゲノムデータが重複しているので集計するとマジョリティに見えてしまう。）
そこで、全く同じ遺伝子構成のクラスターは一つとカウントすることにした。しかし、この方法だと種を超えて広く同じ遺伝子構成が保存されている場合に、逆に少なくカウントされすぎてしまうので、これをどう扱うかは今後の課題としたい。

とりあえず、重複をなくした集計結果でGH36:9.1 と同じクラスターに存在している遺伝子を見てみる。

dbCAN のファミリー分類の結果

一番左の棒がGH36。予想に反し、β-マンナナーゼ、α-ガラクトシダーゼの報告はない GH43、GH3 ファミリーが上位に来ていた。
β-マンナナーゼが属するGH130ファミリーも上位に来ていたが、GH36:9.1は必ずしも論文にあったGH130のβ-マンナナーゼと同じ遺伝子クラスターに存在する訳ではなかった。
その理由としては、他のファミリーのマンナナーゼと協働している、ガラクトマンナン以外の多糖に含まれているα-ガラクトース部分を加水分解する働きもある、などが考えられる。

CUPP によるサブファミリー分類の結果

CUPPの分類ではGH43がさらに細分化されていたが、上位に位置するファミリーは大体同じだった。

HMMERによるその他の酵素のアノテーション

HMMERの結果からは、糖のトランスポーターがGH36:9.1と多く同じ遺伝子クラスターにあることが分かった。トランスポーターは、糖質関連酵素よりも高い頻度で同じクラスターにあった。

糖関連酵素も、dbCAN、CUPP と同じような割合で検出されている。

HMMERの結果では、機能未知ドメインを持つタンパク質も複数ランクインしていたが、これらの機能はガラクトマンナンの分解に関わるものと推測できるかもしれない。

https://doi.org/10.1002/1873-3468.12250 ↩
タンパク質中の短い配列の保存性によって酵素をクラスタリングする手法である、「k-mer 法」で糖質関連酵素の分類を行っているデータベース。
Barrett, K., Hunt, C. J., Lange, L., & Meyer, A. S. (2020). Conserved unique peptide patterns (CUPP) online platform: peptide-based functional annotation of carbohydrate active enzymes. Nucleic Acids Research, 48(W1), W110–W115. https://doi.org/10.1093/nar/gkaa375
Barrett, K., & Lange, L. (2019). Peptide-based functional annotation of carbohydrate-active enzymes by conserved unique peptide patterns (CUPP). Biotechnology for Biofuels, 12(1). https://doi.org/10.1186/s13068-019-1436-5 ↩