GenBankフォーマットの処理

GenBank フォーマットでは、各エントリーは「LOCUS」の行から開始し、「//」で終了する。従って、GenBank フラットファイルを処理する場合は LOCUS 行に遭遇したときに新規エントリーの処理を開始するシグナルを設定して、「//」行に遭遇したときに終了シグナルを設定してデータの保存作業などを行う。

以下に Perl でこのようなファイルを処理する例を挙げているが、正確性や速さを求めるなら BioPerl などの SeqIO モジュールなどを利用した方がいい。

フィールド情報の取得

GenBank形式のファイルでは各行の最初の12文字にフィールド名が記述されている。以下の例では、LOCUSフィールドの内容と、ORGANISMフィールドの内容を取得する例である。

use strict;
use DBI;
use Data::Dumper;
use Compress::Zlib;

#複数のGenBankフォーマットファイルが置かれているディレクトリパス
my $datadir = "/data";
opendir (DATADIR , $datadir);
    my @gbff = readdir(DATADIR);
    foreach my $gbff (@gbff){
        print "処理開始:$gbff\n";
        &readGbffGz("$datadir/$gbff");
        print "処理終了:$gbff\n";
    }
closedir(DATADIR);


#gbffフォーマットファイルを処理するサブルーチン
#複数のエントリーが含まれても可
#「LOCUS」に遭遇すれば、新規エントリーの処理を開始する
#「//」に遭遇すれば、現在のエントリーの終了処理を行う
sub readGbffGz(){
	my $path = shift;    #フラットファイルへのパス
	my $i = 0;           #featureキーの順番を記録する
	my $j = 0;           #qualifierをの順番を記録する
	my $gb = {};         #データを保存するためのリファレンス
	my $tb = '';         #テーブル名を記録する(各行の最初の12文字)
	my $s  = '';         #1つのテーブルが複数行にわたって記載される場合、$tbの補助として働く
	my $buff = '';       #圧縮ファイルの行データを保存するための一時変数
	my $gz = gzopen($path, 'rb') or die "cannot unzip $path\n";

	#圧縮ファイルを一行ずつ読み込んで、その行のデータを$buffに保存する
	while( $gz->gzreadline($buff) > 0 ){
		chomp($buff);    #改行文字を切り取る
		#gbff形式は最初の12文字がテーブル名を示す
		$tb = substr($buff,0,12);  #テーブル名を記録する
		$tb =~ s/\s//g;            #空白を消す

		#エントリーのヘッダー行であれば、変数の初期化を行う
		if($tb eq 'LOCUS'){
			$gb = {};    #データを保存するリファレンスを初期化する
			my @value = split(/\s+/,substr($buff,12));    #空白で文字を分割する
			$gb->{leng} = $value[1];    #divisionを取得する
			$gb->{leng} = $value[5];
			next;
		}
    
		#生物の属名種名を記録する
		if($tb eq 'SOURCE'){
			unless(exists $gb->{genus}){    #属名種名は1行めに記載されているため、存在確認
				my @value = split(/\s+/,substr($buff,2));
				$gb->{genus} = $value[0];
				$gb->{species} = $value[1];
			}
			next;
		}

		#FEATURESテーブル情報を保存するための変数を初期化する
		if($tb eq 'FEATURES'){
			$i = 0;    #features key の順番
			$j = 0;    #そのキーに登録されているqualifierの順番
			$s = 'features';  #FEATURESテーブル処理開始シグナル
			next;
		}
        
		#FEATURESテーブルでは、各feature keyは各行の6文字目から記載されるため、最初の5文字は空白文字である
		#従って$s=='features'かつ行の先頭5文字が空白であればFEATURESテーブルにいると判断できる
		#$s=='features'かつ行の先頭5文字が空白でなければFEATURESテーブルは終了したと判断できる
		if( ($s eq 'features') && (substr($buff,0,5) eq '     ')){
		#行の5文字目~15文字目を記録する。$keyはfeature keyか空白かである
			my $key = substr($buff,5,15);
			$key =~ s/\s+//g;
			#keyが空白でないとき、すなわち新規にfeature keyが開始した
			#この行に含まれるlocation情報を0番目のqualifierとして処理する
			if($key ne ''){
				$j = 0;
				$gb->{ft}->[$i]->{key} = $key;
				#この行にlocation(0番目のqualifier)が存在するので、これを処理する
				$gb->{ft}->[$i]->{qual}->[$j] = substr($buff,21);

				#次(1番目)の新規qualifierが現れたときに使うために添字を1増やす
				$j++;
				#次に新規キーが現れたときに使うために添字を1増やす
				$i++;
			}
			#keyが空白であるとき、すなわちこの行のqualifierは現在のkeyのものである
			#上のif文の最後で$i++を処理したために、「現在のキー」の添字は[$i-1]となる
			else{
				#新規qualifierの開始
				if(substr($buff,21,1) eq '/'){
					$gb->{ft}->[$i-1]->{qual}->[$j] = substr($buff,22); #21文字目は「/」なので無視
					#次(2番目以降)の新規qualifierが現れたときに使うために添字を1増やす
					$j++;
				}
				#新規qualifierでない時、現在のqualifierとして処理する、if文では$j++としたために、添字は[$j-1]である
				else{
					$gb->{ft}->[$i-1]->{qual}->[$j-1] .= substr($buff,21);
				}
			}
			next;
		}
     
		#ORIGINテーブルの次の行から始まる塩基配列を保存するシグナルをONにする
		if($tb eq 'ORIGIN'){
			$gb->{seq} = "";
			$s = 'origin';
			next;
		}

		#配列が記載されている行の最初の10文字は、配列中における塩基の開始位置を数字で書かれている。
		#この行は数字、空白を含むが英文字と記号を含まれないので、塩基配列を抽出するシグナルとし利用する。
		if(($s eq 'origin') && ($tb =~ /[^a-zA-Z\/]/)){
			my $seq = substr($buff, 10);
			$seq =~ s/\s+//g;
			$gb->{seq} .= $seq;
			next;
		}
		#「//」に遭遇すれば終了処理を開始する
		#$gbに格納されたデータを選択・加工しMySQLに保存する
		if($tb eq '//' ){

			#ここで終了処理を行う
			#Data::Dumperを利用して、現在$gbに保存されているデータの構造を確認することができる。
			warn Dumper $gb;
			#位置情報を元に塩基配列を切り出すコードは次のURLを参照
			#http://bi.biopapyrus.net/perl/app/gbff-location.html
			#変数をすべて初期化する
			$gb = {};
			$i = 0;
			$j = 0;
			$s = '';
			$tb = '';
			next;
		}
	}

	$gz->gzclose();
}