Jeans & Development コンピューターのことなどを綴ったメモ (旧：目から鱗 w/SQLite)

生物学の少し古い文献では、遺伝子のクローニングの後に制限酵素地図(restriction enzyme map)を論文に載せてはいるが、最終的に塩基配列(DNA sequence)まで発表していないケースがまれにある。一部でも塩基配列が載せられていれば、全ゲノム配列のどの部分に相当するかを調べるだけで、クローニングされたものがどの遺伝子でどんな塩基配列なのかが分かるのだが、制限酵素地図だけが載せられているようなケースでは、そういった調べ方はできない。

昨日、調べていた文献の中に、上記事項に相当するものが出てきた。数年前にも同じようなことがあって、あるプログラムを書いて、調べ上げたことがある。当時は、こんなことは一回こっきりだろうと思い、作ったプログラムをちゃんと保存していなかったので、一回使ったきり、無くしてしまった。今回は後でまた使えるように、ここにメモをしておく。

PHPを使って、作成した。テストランの実行結果は、次のとおり。ここでは、この論文で扱われている制限酵素地図を元に、該当のDNAフラグメントが大腸菌ゲノムのどの部分に相当するかを求めている。

Generating the list of genes...
Generating the genome sequence string...
Pick up candedates...
................................................................................
................................................................................
................................................................................
........................
264 candidates found
Check out each candidates
12888-16541 (3654 bp)
E---B--------C--------------------------P----------------------B
  EcoRI 0
  BglII 203
  ClaI 705
  PvuII 2293
  BamHI 3648
  14168-15298/gene="dnaJ"
  15445-16557/gene="insL-1"

648275-652175 (3901 bp)
E-B-------C--P------------------------------------------------B
  EcoRI 0
  BglII 122
  ClaI 583
  PvuII 766
  BamHI 3895
  648805-649713/gene="citE"
  649710-650006/gene="citD"
  650021-651079/gene="citC"
  651458-653116/gene="citA"

2590297-2593828 (3532 bp)
E---------B-------C----------------P--------------------------B
  EcoRI 0
  BglII 560
  ClaI 1006
  PvuII 1964
  BamHI 3526
  2590784-2590984/gene="ypfN"
  2591094-2591792/gene="ypfH"
  2591866-2593881/gene="ypfI"

3533569-3536763 (3195 bp)
E-B----------------------------C------P----------------------B
  EcoRI 0
  BglII 80
  ClaI 1621
  PvuII 1990
  BamHI 3189
  3533887-3534606/gene="ompR"
  3534834-3535310/gene="greB"
  3535407-3537728/gene="yhgF"

3881686-3885787 (4102 bp)
E---B--C-------------------P-------B
  EcoRI 0
  BglII 244
  ClaI 387
  PvuII 1703
  BamHI 2202
  3882359-3882499/gene="rpmH"
  3882516-3882875/gene="rnpA"
  3882839-3883096/gene="yidD"
  3883099-3884745/gene="yidC"
  3884851-3886215/gene="trmE"

9097-12893 (3797 bp)
B--------------------------------PC-B------------------------E
  BamHI 0
  PvuII 2064
  ClaI 2104
  BglII 2227
  EcoRI 3791
  9306-9893/gene="mog"
  9928-10494/gene="yaaH"
  10643-11356/gene="yaaW"
  10830-11315/gene="htgA"
  11382-11786/gene="yaaI"
  12163-14079/gene="dnaK"

874401-878990 (4590 bp)
B----P------------C-------------------------------------B-----E
  BamHI 0
  PvuII 372
  ClaI 1313
  BglII 4182
  EcoRI 4584
  874558-875886/gene="yliF"
  875933-877258/gene="yliG"
  877471-877854/gene="yliH"
  877965-879080/gene="yliI"

2277274-2280705 (3432 bp)
B---PC---------B----------------------------------------------E
  BamHI 0
  PvuII 205
  ClaI 211
  BglII 752
  EcoRI 3426
  2277810-2278505/gene="rsuA"
  2278654-2280414/gene="yejH"
  2280539-2280823/gene="rplY"

2529653-2533655 (4003 bp)
B------------------------P----C-------------------B-----------E
  BamHI 0
  PvuII 1603
  ClaI 1918
  BglII 3204
  EcoRI 3997
  2530431-2531402/gene="cysK"
  2531786-2532043/gene="ptsH"
  2532088-2533815/gene="ptsI"

3667606-3669993 (2388 bp)
B----------P---------------C---------------------B------------E
  BamHI 0
  PvuII 416
  ClaI 1052
  BglII 1895
  EcoRI 2382
  3667615-3669264/gene="treF"
  3669315-3669917/gene="yhjB"

いくつか候補が表示されているが、フラグメントの長さとそれぞれの制限酵素の位置から、ゲノム中の3533569-3536763の位置が該当するフラグメントであることが分かる。

PHPスクリプトのソースコードは、以下のとおり。使用する際は、「// Configurations」と「// Restrictiop map」の部分を書き換える。GENBANK_FILEで指定されたファイルを、PHPファイルを同じディレクトリに置くこと。

<?php

// Configurations
define('GENBANK_FILE','EcoliGenome.txt');
enzyme('EcoRI','gaattc');
enzyme('BglII','agatct');
enzyme('ClaI','atcgat');
enzyme('PvuII','cagctg');
enzyme('BamHI','ggatcc');
enzyme('BglI','gccnnnnnnggc');
enzyme('HincII','gtyrac');
enzyme('FokI','ggatg');

// Restrictiop map
$map=array();
$map[]='EcoRI';
$map[]='BglII';
$map[]='ClaI';
$map[]='PvuII';
$map[]='BamHI';
$length=3100;

// Reach to the DNA sequence
echo "Generating the list of genes...\n";
$genes=array();
$gene=$pos=false;
$fh=fopen(GENBANK_FILE,'r');
while(!feof($fh)){
	$line=fgets($fh);
	if (preg_match('/^ORIGIN/',$line)) break;
	if ($pos) {
		$genes[$pos]=$gene.trim($line);
		$gene=$pos=false;
		continue;
	}
	if (preg_match('/^[\s]+gene[\s]+[^0-9]*([0-9]+)[^0-9]+([0-9]+)/',$line,$m)) {
		$pos=(int)$m[1];
		$gene="$m[1]-$m[2]";
	}
}

// Get the genome sequence
echo "Generating the genome sequence string...\n";
$genome='';
while(!feof($fh)){
	$genome.=preg_replace('/[^gatc]+/','',fgets($fh));
}

// List up fragments with the ends of first and last restriction enzyme sites.
echo "Pick up candedates...\n";
$found=array();
$search='/'.enzyme($map[0]).'[gatc]{'.(int)($length/2).','.(int)($length*3/2).'}'.enzyme($map[count($map)-1]).'/';
preg_match_all($search,$genome,$m);
foreach($m[0] as $sequence){
	$pos=strpos($genome,$sequence);
	echo '.';
	$found[$pos]=$sequence;
}
$search='/'.enzyme($map[count($map)-1]).'[gatc]{'.(int)($length/2).','.(int)($length*3/2).'}'.enzyme($map[0]).'/';
preg_match_all($search,$genome,$m);
foreach($m[0] as $sequence){
	$pos=strpos($genome,$sequence);
	echo '.';
	$found[$pos]=$sequence;
}
echo "\n".count($found)." candidates found\n";

// Check each candidates.
echo "Check out each candidates\n";
foreach($found as $pos=>$sequence){
	$size=strlen($sequence);
	$enzpos=array();
	foreach($map as $enzyme){
		if (!preg_match('/^([gatc]*?)'.enzyme($enzyme).'/',$sequence,$m)) break;
		$enzpos[]=strlen($m[1]);
	}
	if (count($enzpos)<count($map)) continue;
	if ($enzpos[0]==0) {
		for($i=count($enzpos)-2;0<=$i;$i--){
			if ($enzpos[$i]<$enzpos[$i+1]) continue;
			$i=1;
			break;
		}
		if (0<$i) continue;
		echo $pos.'-'.($pos+$size-1).' ('.$size." bp)\n";
		$text='';
		for($i=0;$i<count($map);$i++){
			if (0<$i) echo str_repeat('-',(int)(60*($enzpos[$i]-$enzpos[$i-1])/$size));
			echo substr($map[$i],0,1);
			$text.='  '.$map[$i].' '.$enzpos[$i]."\n";
		}
		echo "\n$text";
		foreach($genes as $genepos=>$gene){
			if ($genepos<$pos) continue;
			if ($pos+$size-1<$genepos) break;
			echo "  $gene\n";
		}
		echo "\n";
	} else {
		for($i=count($enzpos)-2;0<=$i;$i--){
			if ($enzpos[$i]>$enzpos[$i+1]) continue;
			$i=1;
			break;
		}
		if (0<$i) continue;
		echo $pos.'-'.($pos+$size-1).' ('.$size." bp)\n";
		$text='';
		for($i=count($map)-1;0<=$i;$i--){
			if ($i<count($map)-1) echo str_repeat('-',(int)(60*($enzpos[$i]-$enzpos[$i+1])/$size));
			echo substr($map[$i],0,1);
			$text.='  '.$map[$i].' '.$enzpos[$i]."\n";
		}
		echo "\n$text";
		foreach($genes as $genepos=>$gene){
			if ($genepos<$pos) continue;
			if ($pos+$size-1<$genepos) break;
			echo "  $gene\n";
		}
		echo "\n";
	}
}


// Functions follow
function enzyme($type,$sequence=false){
	static $enzymes=array();
	$type=strtolower($type);
	if (isset($enzymes[$type])) return $enzymes[$type];
	if ($sequence) $enzymes[$type]=recog_seq($sequence);
}
function recog_seq($sequence){
	static $change;
	if (!isset($change)) {
		$change=array(
			'r'=>'[ag]',
			'y'=>'[ct]',
			'm'=>'[ac]',
			'k'=>'[gt]',
			's'=>'[cg]',
			'w'=>'[at]',
			'h'=>'[act]',
			'b'=>'[cgt]',
			'v'=>'[acg]',
			'd'=>'[agt]',
			'n'=>'[acgt]');
	}
	$sequence=strtolower($sequence);
	if (preg_match('/[^a-z]/',$sequence)) exit('ERROR'.__LINE__);
	$comp=complementary($sequence);
	if ($comp==$sequence) return strtr($sequence,$change);
	else return '(?:'.strtr($sequence,$change).'|'.strtr($comp,$change).')';
}
function complementary($sequence){
	static $change;
	if (!isset($change)) {
		$change=array(
			'g'=>'c',
			'a'=>'t',
			't'=>'a',
			'c'=>'g',
			'r'=>'y',
			'y'=>'r',
			'm'=>'k',
			'k'=>'m',
			's'=>'s',
			'w'=>'w',
			'h'=>'d',
			'b'=>'v',
			'v'=>'b',
			'd'=>'h',
			'n'=>'n');
	}
	$comp='';
	for($i=strlen($sequence)-1;0<=$i;$i--){
		$comp.=$change[substr($sequence,$i,1)];
	}
	return $comp;
}