GenBankのlocation情報処理

GenBank/EMBL/DDBJ の FEATURES テーブルには配列の特徴などが記載されている。例えば配列の名前(CDS、mRNA、exonなど)と位置情報(12..399、join(23..100,320..559)など)などである。以下にあるPerlスクリプトは、配列の名前と位置情報を入力して、該当部分の塩基配列を出力するものである。

まず、位置情報(location)のフォーマットは次のようなものがある。

表示形式意味
467単一塩基だけからなる配列の場合は、その位置の番号だけが記載される。
340..565連続した配列の範囲を示す。340 番目の塩基から 565 番目の塩基までの配列。
<345..500345 番目の塩基から 500 番目の塩基までの配列。ただし、正確な開始位置は345よりも前(5'方向)にあり、正確には分かっていない。
<1..8881 よりも前にある塩基から始まり、888 番目塩基までの配列を示す。1 よりも前にある塩基は、現在のエントリーに記載されてない。
1..888>1 番目の塩基から始まり 888 番目よりも後の塩基で終わる。ただし、正確な位置は分かっていない。
102.110塩基配列の正確な位置情報が分かっていない。ただし、102 番目から 110 番目の間にある。
123^1242 塩基だけからなる塩基配列の場合は「^」によて記載される。
join(12..78,134..202)12 番目~78 番目の配列の後に、134 番目~202 番目の配列を結合させた配列。
complement(34..126)34 番目~125 番目の塩基配列の相補鎖。プログラムで処理するときは、まず 34..125 として処理し、続いて、塩基配列を逆転 (125..34) させて、A を T、C を G、G を C、T を A に転換する必要がある。
complement(join(269..451,550..790))269..451 と 550..790 の結合配列の相補鎖。
join(complement(118..168),complement(35..89))35..89 の相補鎖と、118..168 の相補鎖を結合させた配列。順序は 168->118+89->35となっている。
J00194.1:100..202accession 番号が J00194 に記載されている塩基配列中の 100 番目~202 番目までの配列。J00194.1 の小数点以降はバージョンを示す。
join(1..100,J00194.1:100..202)現在エントリーの 1..100 配列と、J00194 エントリーに含まれる 100..202 配列を結合した配列。

処理スクリプトの考え方としては次のようにする。ただし、id は Accession 番号で、J00194.1:100..202 などを含む位置情報を処理するために用いる。

  1. mainProcess(id, location, seq) を入力する。location は complement ならば complementProcess()、join ならば joinProcess() を呼び出す。それ以外ならば generalProcess() を呼び出す。

  2. joinProcess() に入力された位置情報に対して、外側の「join()」を削除する。すると中身は 1..10,21..30,complement(12..231) などのようにカンマ区切りになる。カンマを区切り文字として子クラスの location を作成して mainProcess() に代入する。

  3. complementProcess() も join() と同様に処理する。ただし、return $seg の直前に配列の相補鎖を求める。

  4. 上記を繰り返すことでいずれ、12..31、12^23、<34..110 などのような形式になる。これらは generalProcess() で処理する。

代入する塩基配列は何十万にも及ぶ可能性があるため、ここでは直接文字列を代入せずに、代わりに塩基配列のアドレスを代入する(リファレンス変数)。


my $segment = clipSequence($accession, $location, \$seq);
print $segment;



#塩基配列とlocation情報を入力して、locationに基づいて該当部分を切り出してかえす。
#location情報にはAccessionが記載されている場合があるので、Accessionの一致の確認を
#行うためにAccession情報を入力する。
sub clipSequence(){
	my ($id, $location, $seq) = @_;
	my $seg = '';
	$location =~ s/\s//g;

	if($location =~ /:/){
		#そのaccessionと現在エントリーのaccessionをチェックし、一致しなければ処理を行わない
		my @tmp = split(/\:/, $location);
		my $locate_id = $tmp[0];
		if($locate_id =~ /$id/){
			$seg = &clipSeq_mainProcess($id, $tmp[1], $seq);
		}else{
			$seg = '(('.$locate_id.'))';
		}
	}else{
		$seg = &clipSeq_mainProcess($id, $location, $seq);
	}

	return $seg;
}

#locationのタイプを識別して下位関数を呼び出す
#complete() , join() , その他 の3タイプ
sub clipSeq_mainProcess(){
	my ($id, $location, $seq) = @_;

	my $seg = '';
	my $s = substr($location, 0, 1);
	if($s eq 'c'){
		#complement()を含むlocationの処理
		$seg = &clipSeq_complementProcess($id, $location, $seq);
	}elsif($s eq 'j'){
		#join()を含むlocationの処理
		$seg = &clipSeq_joinProcess($id, $location,  $seq);
	}else{
		#complement()とjoin()を含まない処理
		#">", "<", "^", "..", "." などを含む
		$seg = &clipSeq_generalProcess($id, $location, $seq);
	}

	return $seg;
}


#complement(),join()を含まないlocationは次のような場合がある
# 21
# 21^22
# 21.80
# 21..80
# <21..80
# 21..80>
# NX00001.2:21..80
#location中のaccession情報が存在すればそれをチェックし、
sub clipSeq_generalProcess(){
	my ($id, $location, $seq) = @_;
	my $seg = '';
	my $pre = 0;

	#location2は「>」と「<」がついている位置を処理するときに用いる
	my $location2 = $location;

	$location =~ s/[\>|\<]//g;
	my @lim = split(/\.+/ , $location );
		
	#「^」を含むlocationの処理
	if($location =~ /\^/){
		$seg = substr($$seq, $lim[0]-1, 2);
	}
	#「.」を含むlocationの処理、ピリオドが1つの場合、不正確な範囲を示す
	elsif(($location =~ tr/\./\./) == 1){
		$seg = substr($$seq, $lim[0]-1, $lim[1]-$lim[0]+1);
		$seg = '['.$seg.']';
	}
	#「..」を含むlocationの処理、ピリオドが2つの場合、正確な範囲を示す
	elsif(($location =~ tr/\./\./) == 2){
		$seg = substr($$seq, $lim[0]-1, $lim[1]-$lim[0]+1);
	}
	#「22」数字のみのlocation処理
	else{
		$seg = substr($$seq, $lim[0]-1, 1);
	}

	#locationに「>」または「<」が含まれる場合は、そのコードを書き込む
	if($location2 =~ /\>/){$seg .= '>';}
	if($location2 =~ /\</){$seg = '<'.$seg;}

	return $seg;

}

#complement()を処理する関数
#complement()の文字列を削除して、中身を取り出して、mainProcessを利用してjoinタイプかgeneralタイプかを判断する
sub clipSeq_complementProcess(){
	my ($id, $location, $seq) = @_;
	my $seg = '';

	#complement()の場合、一番外側のcomplement( と )を消す
	$location =~ s/^complement\(//;
	$location =~ s/\)$//;

	#「21..50,60..80」などのように複数の場合があるのでカンマで区切られている
	my @location = split(/,/ , $location);
	foreach my $lim (@location){
		$seg .= &clipSeq_mainProcess($id, $lim, $seq);
	}

	#切り抜いた配列を逆転させ、相補鎖を求める
	$seg = reverse($seg);
	$seg =~ tr/ACGTacgt/TGCAtgca/;

	return $seg;
}

#join()を処理する関数
#join()の文字列を削除して、中身を取り出して、mainProcessを利用してcomplementタイプかgeneralタイプかを判断する
sub clipSeq_joinProcess(){
	my ($id, $location, $seq) = @_;
	my $seg = '';
	#join()の場合、一番外側のjoin( と )を消す
	$location =~ s/^join\(//;
	$location =~ s/\)$//;

	#「21..50,60..80」などのように複数の場合があるのでカンマで区切られている
	my @location = split(/,/ , $location);
	foreach my $lim (@location){
		$seg .= &clipSeq_mainProcess($id, $lim, $seq);
	}

	return $seg;
}