EMBLフォーマットの処理

EMBLフォーマットは GenBank フォーマットとほぼ同じである。EMBL フォーマットの場合はすべての行に大文字2文字からなる識別子が先頭にかかれている。次のスクリプトは、gzip によって圧縮された EMBL フォーマットのファイルを読み込んで必要なデータを抽出するサンプルである。

なお、正確性を求めるなら BioPerl の SeqIO を利用したほうが良い。


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

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



# 一つのgenbankを読み込んで、データを抽出する関数
# 複数のエントリーが含まれても可。「//」で終了処理を行っている
sub readEmblGz(){
	my $path = shift;	#フラットファイルへのパス
	my $i = 0;			#featureキーの順番を記録する
	my $j = 0;			#qualifierをの順番を記録する
	my $gb = {};		#データを保存するためのリファレンス
	my $tb = '';		#テーブル名を記録する(各行の最初の2文字)
	my $buff = '';
	my $gz = gzopen($path, 'rb') or die "cannot unzip $path\n";

	while( $gz->gzreadline($buff) > 0 ){
		chomp($buff);	#改行文字を切り取る
		#embl形式は最初の2文字がテーブル名を示す
		$tb = substr($buff,0,2);	#テーブル名を記録する

		#エントリーのヘッダー行であれば、変数の初期化を行う
		if($tb eq 'ID'){
			$gb = {};	#データを保存するリファレンスを初期化する

			my @value = split(/;/,substr($buff,2));	#セミコロンで文字を分割する
			$gb->{division} = $value[2];	#divisionを取得する
			next;
		}
	
		#生物の属名種名を記録する
		#OSテーブルは2行以上の場合もあるため、2行め以降の情報を無視する
		if($tb eq 'OS'){
			unless(exists $gb->{genus}){	#属名種名は1行めに記載されているため、存在確認
				my @value = split(/\s+/,substr($buff,2));
				$gb->{genus} = $value[1];
				$gb->{species} = $value[2];
			}
			next;
		}


		#FEATURESテーブル情報を保存するための変数を初期化する
		if($tb eq 'FH'){
			$i = 0;	#features key の順番
			$j = 0;	#そのキーに登録されているqualifierの順番
			next;
		}
		
		#FEATUERSテーブルの処理を開始する
		if($tb eq 'FT'){
			#FTテーブルにおいて新規にキーが現れば最初の20文字に書かれる
			my $key = substr($buff, 2, 18);
			$key =~ s/\s//g;
	
			#新規にキーが開始するとき、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++;
			}
			#新規にキーが開始でない時、qualifier属性は現在のキーのものとして処理する
			#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;
		}

		#SQテーブルの次の行から始まる塩基配列を保存するシグナルをONにする
		if($tb eq 'SQ'){
			$gb->{seq} = "";
			next;
		}

		#配列が記載されている行の先頭2文字は空白文字となっている
		if($tb eq '  '){
			my $seq = substr($buff, 0, 70);
			$seq =~ s/\s//g;
			$gb->{seq} .= $seq;
			next;
		}

		#「//」に遭遇すれば終了処理を開始する
		#$gbに格納されたデータを選択・加工しMySQLに保存する
		if($tb eq '//'){
			#
			#データの処理
			#
			warn Dumper $gb;
			#位置情報を元に塩基配列を切り出すコードは次のURLを参照
			#https://bi.biopapyrus.net/perl/app/gbff-location.html

			#変数をすべて初期化する
			$gb = {};
			$i = 0;
			$j = 0;
			$tb = '';
		}

	}

	$gz->gzclose();
}