「らせん対称をもつ粒子の３次元再構成法」の版間の差分

2017年1月7日 (土) 02:55時点における版

らせん対称構造のリファインメント

こちらはある学生が作成したらせん対称構造のリファインメントのためのmakefileです。 一般性が少し欠けているので、公開するには書き直した方がよいものですね。それでもmakefileの書き方などの例になるかと思い、解説してみます。

Makefile

下記に示したものは、アクトミオシンの構造解析の際に使っていたMakefileである。

 .SUFFIXES:  .roi .roi2 .mrc3d .ds6 \
			.pad2 .mrc2d .pad3d .corinfo \
			.3dlist_part .3dlist .3dorglist \
			.refined3d .refinedds6 

include ROI
include ROI2
include MRC3D
include PAD2
include MRC2D
include PAD3D
include CORINFO
include 3DLIST_PART
include 3DLIST
include 3DORGLIST
include REFINED3D
include REFINEDDS6

#Refinement param
SOURCE = F-actin#          #SOURCE: Reference for refinement
ROI_WIDTH = 128
 # Parameter for mrc3Dto2D
rXmin = 89
rXmax = 91
rXdel = 1
rYmin = 0#                 #rXmin,rXmax,rXdel: Used in mrc3Dto2D.
rYmax = 355#               #                   rXdel is delta of Angle
rYdel = 5#                 #                   rYmin,rYmax,rYdel are also same
#MODE  = YOYS#			   #Euler mode Y-X-Z 
#MODE  = ZOYS#			   #Euler mode Z-Y-X
MODE  = XEYS#			   #Euler mode X-Y-Z: 1st, Filament axis lap over Y-axis,
						   #				  2nd, Rotate filament about filament axis
 # Parameter for mrcImageAutoRotationCorrelation
PADWIDTH_FORDATA = 128#    #
PADHEIGHT_FORDATA = 128#   #
Xmin = 89
Xmax = 91
Xstep = 3 
Ymin = 0
Ymax = 355
Ystep = 72 
Zmin = 0
Zmax = 360
Zstep = 1

roi2:		$(ROI:.roi=.roi2)
pad2:		$(ROI2:.roi2=.pad2)
mrc2d:		$(MRC3D:.mrc3d=.mrc2d)
pad3d:		$(MRC2D:.mrc2d=.pad3d)
corinfo:	$(PAD2:.pad2=.corinfo)
3dlist_part:$(CORINFO:.corinfo=.3dlist_part)
3dlist:		$(3DLIST_PART:.3dlist_part=.3dlist)
3dorglist:	$(3DLIST:.3dlist=.3dorglist)
refined3d:	$(3DORGLIST:.3dorglist=.refined3d)
refinedds6:	$(REFINED3D:.refined3d=.refinedds6)

FRAGMENTATION::
	make roi2;	  make ROI2;

REFINEMENT::
	make MRC3D;	  make ROI;
	make pad2;	  make PAD2;
	make mrc2d;	  make MRC2D;
	make corinfo; make CORINFO;
	make 3dlist_part; make 3DLIST_PART;
	make 3dlist; make 3DLIST;
	make 3dorglist; make 3DORGLIST;
	make refined3d; make REFINED3D;
	make refinedds6; make REFINEDDS6;

ROI::
	echo "ROI=\\" > ROI
	ls -1 *.roi | sed s/roi/roi\\\\/ >> ROI
	echo "" >> ROI

ROI2::
	echo "ROI2=\\" > ROI2
	ls -1 *.roi2 | sed s/roi2/roi2\\\\/ >> ROI2
	echo "" >> ROI2

MRC3D::
	echo "MRC3D=\\" > MRC3D
	ls -1 *.mrc3d | sed s/mrc3d/mrc3d\\\\/ >> MRC3D
	echo "" >> MRC3D

PAD2::
	echo "PAD2=\\" > PAD2
	ls -1 *.pad2 | sed s/pad2/pad2\\\\/ >> PAD2
	echo "" >> PAD2

MRC2D::
	echo "MRC2D=\\" > MRC2D
	ls -1 *.mrc2d | sed s/mrc2d/mrc2d\\\\/ >> MRC2D
	echo "" >> MRC2D

CORINFO::
	echo "CORINFO=\\" > CORINFO
	ls -1 *.corinfo | sed s/corinfo/corinfo\\\\/ >> CORINFO
	echo "" >> CORINFO

3DLIST_PART::
	echo "3DLIST_PART=\\" > 3DLIST_PART
	ls -1 *.3dlist_part | sed s/3dlist_part/3dlist_part\\\\/ >> 3DLIST_PART
	echo "" >> 3DLIST_PART

3DLIST::
	echo "3DLIST=\\" > 3DLIST
	ls -1 *.3dlist | sed s/3dlist/3dlist\\\\/ >> 3DLIST
	echo "" >> 3DLIST

3DORGLIST::
	echo "3DORGLIST=\\" > 3DORGLIST
	ls -1 *.3dorglist | sed s/3dorglist/3dorglist\\\\/ >> 3DORGLIST
	echo "" >> 3DORGLIST

REFINED3D::
	echo "REFINED3D=\\" > REFINED3D
	ls -1 *.refined3d | sed s/refined3d/refined3d\\\\/ >> REFINED3D
	echo "" >> REFINED3D

REFINEDDS6::
	echo "REFINEDDS6=\\" > REFINEDDS6
	ls -1 *.refinedds6 | sed s/refinedds6/refinedds6\\\\/ >> REFINEDDS6
	echo "" >> REFINEDDS6

#Triming of Noize
CENTERGET:
	mrcImageCenterGet -i $(SOURCE).mrc3d -o $(SOURCE).centre3d \
					  -Nx 60 -Ny 60 

#FRAGMENTATION
.roi.roi2:
	((i=1));\
	((Shift=64));\
	((ty=857));\
	((by=729));\
	while((by >= 0));\
	do \
		echo "i=$$i";\
		echo "ty=$$ty by=$$by";\
		mrcImageROI -i $*.roi -o $*-$$i.roi2\
					-r 0.0 $$by 121.0 $$by 121.0 $$ty 0.0 $$ty -m 2;\
		((i++));\
		((ty-=64));\
		((by-=64));\
	done

#REFINEMENT
.roi2.pad2:
	mrcImageCenterGet -i $*.roi2 -o $*.roiget
	mrcImagePad -i $*.roiget -o $*.pad2 \
				-W $(PADWIDTH_FORDATA) -H $(PADHEIGHT_FORDATA) \
				-m 2 -V 0

.mrc3d.mrc2d:
	mrc3Dto2D -i $(SOURCE).mrc3d -o $(SOURCE).mrc2d \
			  -EulerMode $(MODE) -InterpolationMode 0 \
			  -Rot1 $(rXmin) $(rXmax) $(rXdel) \
			  -Rot2 $(rYmin) $(rYmax) $(rYdel) \
			  -Rot3 0 0 5 -m 1

.pad2.corinfo:
	mrcImageAutoRotationCorrelation -i $*.pad2 -r $(SOURCE).mrc2d \
									-O $*.corinfo -fit $*.fit \
									-nRot1 $(Xstep) \
									-nRot2 $(Ystep) \
									-nRot3 $(Zstep) \
									-range -5 5 -n 10 -Iter 2 -m 2
#									-nRot1Area $(Xmin) $(Xmax) $(Xstep) \
									-nRot2Area $(Ymin) $(Ymax) $(Ystep) \
									-nRot3Area $(Zmin) $(Zmax) $(Zstep) \

.corinfo.3dlist_part:
	awk 'BEGIN {cor=0} \
			   {if(cor<$$18){cor=$$18; \
			   				 maxRot1=$$3; \
							 maxRot2=$$4; \
							 maxRot3=$$5;}} \
		 END {sub(".corinfo",".fit",FILENAME); \
		 	  printf("%s $(MODE) %f %f %f Cor %f\n", FILENAME\
			  									, maxRot1\
												, maxRot2\
												, maxRot3\
												, cor)}'\
		$*.corinfo > $*.3dlist_part

.3dlist_part.3dlist:
	cat $*.3dlist_part >> $(SOURCE).3dlist

.3dlist.3dorglist:
	awk '{print $$7 " " $$1 " " $$2 " " $$3 " " $$4 " " $$5 " " $$6}' \
		$(SOURCE).3dlist | \
	sort -r | \
	awk '{print $$2 " " $$3 " " $$4 " " $$5 " " $$6 " " $$7 " " $$1}' \
		> $(SOURCE).3dorglist

.3dorglist.refined3d:
	mrc2Dto3D -I $(SOURCE).3dorglist -o $(SOURCE).refined3d \
			  -InterpolationMode 2 \
			  -CounterThreshold 1e-6 \
			  -DoubleCounter $(SOURCE).counter \
			  -WeightMode 2 \
			  -Double -m 1

.refined3d.refinedds6:
	mrc2map -i $(SOURCE).refined3d -o $(SOURCE).refinedds6

DispROIs::
	for i in $(ROI) ;\
	do \
		echo $$i ; \
		Display2 -i $$i -Zoom -0.3 -geometry 786x786+0+0 ; \
	done

Mafileの解説

これについて、簡単に解説する。

サフィックスの定義

 .SUFFIXES:  .roi .roi2 .mrc3d .ds6 \
			.pad2 .mrc2d .pad3d .corinfo \
			.3dlist_part .3dlist .3dorglist \
			.refined3d .refinedds6 

ここまでは、サフィックスを定義したものである。切り出した繊維構造のファイル名が.roiであることを前提にしている。 既に粒子化（フラグメンテーションが終わっている場合には、roi2となっていることが前提になる。その際のネーミングに注意が必要である。下記を参照してください。 フラグメンテーションを個々で行った方が検索が楽かもしれない。

ファイルターゲットの定義

include ROI
include ROI2
include MRC3D
include PAD2
include MRC2D
include PAD3D
include CORINFO
include 3DLIST_PART
include 3DLIST
include 3DORGLIST
include REFINED3D
include REFINEDDS6

ここまでは内部で作成されるリストファイルである。 例えば、

$ make ROI2

によりROI2に含まれるファイルリストが内部に含まれる。

設定するべきパラメータ

ここから下に設定するためのパラメータのリストで有り、この部分を別のファイルにしておくことも出来る。

#Refinement param
SOURCE = F-actin#          #SOURCE: Reference for refinement
ROI_WIDTH = 128
 # Parameter for mrc3Dto2D
rXmin = 89
rXmax = 91
rXdel = 1
rYmin = 0#                 #rXmin,rXmax,rXdel: Used in mrc3Dto2D.
rYmax = 355#               #                   rXdel is delta of Angle
rYdel = 5#                 #                   rYmin,rYmax,rYdel are also same
#MODE  = YOYS#			   #Euler mode Y-X-Z 
#MODE  = ZOYS#			   #Euler mode Z-Y-X
MODE  = XEYS#			   #Euler mode X-Y-Z: 1st, Filament axis lap over Y-axis,
						   #				  2nd, Rotate filament about filament axis
 # Parameter for mrcImageAutoRotationCorrelation
PADWIDTH_FORDATA = 128#    #
PADHEIGHT_FORDATA = 128#   #
Xmin = 89
Xmax = 91
Xstep = 3 
Ymin = 0
Ymax = 355
Ystep = 72 
Zmin = 0
Zmax = 360
Zstep = 1

この部分をファイル(Mafile.config)にしておき、

-include Makefile.config

とすることで、読み込むことができる。上記の記述の下でincludeすれば、デフォルト値を書き換えることができる。 -includeはファイルがなくてもエラーにならないので便利です。上述の各includeファイルもそのほうがよいかもしれません。

ターゲットの定義

roi2:		$(ROI:.roi=.roi2)
pad2:		$(ROI2:.roi2=.pad2)
mrc2d:		$(MRC3D:.mrc3d=.mrc2d)
pad3d:		$(MRC2D:.mrc2d=.pad3d)
corinfo:	$(PAD2:.pad2=.corinfo)
3dlist_part:$(CORINFO:.corinfo=.3dlist_part)
3dlist:		$(3DLIST_PART:.3dlist_part=.3dlist)
3dorglist:	$(3DLIST:.3dlist=.3dorglist)
refined3d:	$(3DORGLIST:.3dorglist=.refined3d)
refinedds6:	$(REFINED3D:.refined3d=.refinedds6)

この部分は、ファイルリストをターゲットにしている部分です。

 

FRAGMENTATION::
	make roi2;	  make ROI2;

下記のコマンドにより繊維構造から切り出される。

$ make FRAGMENTATION

REFINEMENT::
	make MRC3D;	  make ROI;
	make pad2;	  make PAD2;
	make mrc2d;	  make MRC2D;
	make corinfo; make CORINFO;
	make 3dlist_part; make 3DLIST_PART;
	make 3dlist; make 3DLIST;
	make 3dorglist; make 3DORGLIST;
	make refined3d; make REFINED3D;
	make refinedds6; make REFINEDDS6;

こちらはらせん対称性構造のリファインメントのための内部でのターゲット生成の順序になります。

ターゲットファイルリストの作成

ROI::
	echo "ROI=\\" > ROI
	ls -1 *.roi | sed s/roi/roi\\\\/ >> ROI
	echo "" >> ROI

ROI2::
	echo "ROI2=\\" > ROI2
	ls -1 *.roi2 | sed s/roi2/roi2\\\\/ >> ROI2
	echo "" >> ROI2

MRC3D::
	echo "MRC3D=\\" > MRC3D
	ls -1 *.mrc3d | sed s/mrc3d/mrc3d\\\\/ >> MRC3D
	echo "" >> MRC3D

PAD2::
	echo "PAD2=\\" > PAD2
	ls -1 *.pad2 | sed s/pad2/pad2\\\\/ >> PAD2
	echo "" >> PAD2

MRC2D::
	echo "MRC2D=\\" > MRC2D
	ls -1 *.mrc2d | sed s/mrc2d/mrc2d\\\\/ >> MRC2D
	echo "" >> MRC2D

CORINFO::
	echo "CORINFO=\\" > CORINFO
	ls -1 *.corinfo | sed s/corinfo/corinfo\\\\/ >> CORINFO
	echo "" >> CORINFO

3DLIST_PART::
	echo "3DLIST_PART=\\" > 3DLIST_PART
	ls -1 *.3dlist_part | sed s/3dlist_part/3dlist_part\\\\/ >> 3DLIST_PART
	echo "" >> 3DLIST_PART

3DLIST::
	echo "3DLIST=\\" > 3DLIST
	ls -1 *.3dlist | sed s/3dlist/3dlist\\\\/ >> 3DLIST
	echo "" >> 3DLIST

3DORGLIST::
	echo "3DORGLIST=\\" > 3DORGLIST
	ls -1 *.3dorglist | sed s/3dorglist/3dorglist\\\\/ >> 3DORGLIST
	echo "" >> 3DORGLIST

REFINED3D::
	echo "REFINED3D=\\" > REFINED3D
	ls -1 *.refined3d | sed s/refined3d/refined3d\\\\/ >> REFINED3D
	echo "" >> REFINED3D

REFINEDDS6::
	echo "REFINEDDS6=\\" > REFINEDDS6
	ls -1 *.refinedds6 | sed s/refinedds6/refinedds6\\\\/ >> REFINEDDS6
	echo "" >> REFINEDDS6

ここまでがリストをつくっているところです。基本的には特定のサフィックスをもつファイルリストをmakefile用の変数へと変更しています。

３次元画像の取り扱い

#Triming of Noize
CENTERGET:
	mrcImageCenterGet -i $(SOURCE).mrc3d -o $(SOURCE).centre3d \
					  -Nx 60 -Ny 60 

この部分は、三次元画像の周辺を切り取るための仕掛です。-Nx, -Nyが大きさを示していますが、この部分も変数にした方がよいですね。

断片化

#FRAGMENTATION
.roi.roi2:
	((i=1));\
	((Shift=64));\
	((ty=857));\
	((by=729));\
	while((by >= 0));\
	do \
		echo "i=$$i";\
		echo "ty=$$ty by=$$by";\
		mrcImageROI -i $*.roi -o $*-$$i.roi2\
					-r 0.0 $$by 121.0 $$by 121.0 $$ty 0.0 $$ty -m 2;\
		((i++));\
		((ty-=64));\
		((by-=64));\
	done

　この部分で断片化しています。ここも内部に直接数が書かれていますが、修正した方が便利ですね。 ty, byは、切り出すボックスのy方向の位置を示しています。ループ毎に値が変わります。 121.0は、ボックスのｘ方向の大きさ（フィラメントの太さ方向）

これらはもとの画像に依存する部分になります。

64がty, byのシフト量（Shift変数は定義されているけれど使われていない） tyとbyが最初の切り出し位置を示しています。これも画像に応じて帰るべきですね。

例えば、つぎのような遣り方がありますね。

#FRAGMENTATION

#こちらが、どれだけずらすかのpixel単位
YShift=64    
      
#こちらが、繊維方向の長さをしめしたもの
YLength=128

.roi.roi2:
	((i=1));\
	((ty=`mrcInfo -i $*.roi | awk '/^N / {printf("%d",$$4);}'`));\
	((by=$$ty-$(YLength)));\
	while((by >= 0));\
	do \
		echo "i=$$i";\
		echo "ty=$$ty by=$$by";\
		mrcImageROI -i $*.roi -o $*-$$i.roi2\
					-r 0.0 $$by 121.0 $$by 121.0 $$ty 0.0 $$ty -m 2;\
		((i++));\
		((ty-=$(YShift)));\
		((by-=$(YShift)));\
	done

リファインメント

パッド

#REFINEMENT
.roi2.pad2:
	mrcImageCenterGet -i $*.roi2 -o $*.roiget
	mrcImagePad -i $*.roiget -o $*.pad2 \
				-W $(PADWIDTH_FORDATA) -H $(PADHEIGHT_FORDATA) \
				-m 2 -V 0

ここで最終的に繊維（単粒子）の前処理を行います。

モデルの作成

SOURCEは予め設定されています。

.mrc3d.mrc2d:
	mrc3Dto2D -i $(SOURCE).mrc3d -o $(SOURCE).mrc2d \
			  -EulerMode $(MODE) -InterpolationMode 0 \
			  -Rot1 $(rXmin) $(rXmax) $(rXdel) \
			  -Rot2 $(rYmin) $(rYmax) $(rYdel) \
			  -Rot3 0 0 5 -m 1

角度の決定

.pad2.corinfo:
	mrcImageAutoRotationCorrelation -i $*.pad2 -r $(SOURCE).mrc2d \
									-O $*.corinfo -fit $*.fit \
									-nRot1 $(Xstep) \
									-nRot2 $(Ystep) \
									-nRot3 $(Zstep) \
									-range -5 5 -n 10 -Iter 2 -m 2
#									-nRot1Area $(Xmin) $(Xmax) $(Xstep) \
									-nRot2Area $(Ymin) $(Ymax) $(Ystep) \
									-nRot3Area $(Zmin) $(Zmax) $(Zstep) \

探索範囲に関してはすべて予め定数化されています。 -rangeは、面内回転分を示しています。最初が１度刻み、その後、０．２度、０．０２度と変化していきます。 この部分に時間がかかります。Relion等はフーリエ空間で常に比較しています。この辺りは今後プログラムの修正が必要かも知れません。

相関値による選択

.corinfo.3dlist_part:
	awk 'BEGIN {cor=0} \
			   {if(cor<$$18){cor=$$18; \
			   				 maxRot1=$$3; \
							 maxRot2=$$4; \
							 maxRot3=$$5;}} \
		 END {sub(".corinfo",".fit",FILENAME); \
		 	  printf("%s $(MODE) %f %f %f Cor %f\n", FILENAME\
			  									, maxRot1\
												, maxRot2\
												, maxRot3\
												, cor)}'\
		$*.corinfo > $*.3dlist_part

.3dlist_part.3dlist:
	cat $*.3dlist_part >> $(SOURCE).3dlist

.3dlist.3dorglist:
	awk '{print $$7 " " $$1 " " $$2 " " $$3 " " $$4 " " $$5 " " $$6}' \
		$(SOURCE).3dlist | \
	sort -r | \
	awk '{print $$2 " " $$3 " " $$4 " " $$5 " " $$6 " " $$7 " " $$1}' \
		> $(SOURCE).3dorglist

　この部分で、尤も高い相関値をもつ角度を見出しています。この辺りを尤度などを利用していくこともできます。

三次元再構成　

.3dorglist.refined3d:
	mrc2Dto3D -I $(SOURCE).3dorglist -o $(SOURCE).refined3d \
			  -InterpolationMode 2 \
			  -CounterThreshold 1e-6 \
			  -DoubleCounter $(SOURCE).counter \
			  -WeightMode 2 \
			  -Double -m 1

.refined3d.refinedds6:
	mrc2map -i $(SOURCE).refined3d -o $(SOURCE).refinedds6

DispROIs::
	for i in $(ROI) ;\
	do \
		echo $$i ; \
		Display2 -i $$i -Zoom -0.3 -geometry 786x786+0+0 ; \
	done

　この部分で三次元再構成を行っています。