Perl : 线粒体tRNA二级结构的美化（mitos2）

Perl : 线粒体tRNA二级结构的美化（mitos2）

mitos2是一款动物线粒体注释的软件，同时也可以获得tRNA的二级结构，且获得的二级结构一般都是三叶草形。如果我们有tRNA的序列，也可以使用一些软件将其二级结构画出来，但是如果参数调整不好，画出来的不一定是三叶草形，这点没有mitos2生成的tRNA二级结构友好。
虽然mitos2的tRNA二级结构比较规范，但是它展示的图片信息并不美观。原图如下：

它给出的结果包含了svg的格式，svg是一种矢量图，我们可以修改svg文件的内容来改变图片的形状，修改后的的图片如下：

修改的策略是：

1. 去掉相邻碱基之间的连线
2. 减短配对碱基之间的连线（也可以删除）
3. 增加底色
4. 增加tRNA名称

好的，接下来我们使用perl脚本来完成以上操作：

#!/usr/bin/env perl -w 
use strict;
use File::Basename qw/basename dirname/;

#用于修改mitos2生成的tRNA二级结构svg图片
#从命令行读入svg文件，输出到标准输出，需要将输出内容重定向到文件
#perl $0 sample-trnname-start-end.svg > sample-polish-trnname-start-end.svg
#contact:xul@genepioneer.cn
#2020-06-24

open IN,$ARGV[0] or &usage;	#打开指定的文件，第一个参数

my ($name) = basename($ARGV[0]) =~ /^.*?-(.*?)-/;	
print STDERR "$name\n";	#注意不能把这个输出重定向到文件中，不然svg格式错误
#mitos2给的svg命名符合：sample-trnname-start-end.svg格式
#这里匹配到trna的名称

my $flag = 0;
while(<IN>){
	chomp;
	if(/<g\s*transform/){	#此处是修改图片内容的大小，使其变得小一些，便于以后图片的组合
		my ($size1) = /scale\((\S*?),/g;
		if($size1 > 1.5){
			s/scale\(\S*?\)/scale(1.5,1.5)/;
			print "$_\n";
		}
		next;
	}
	
	if(/<line/){	#这里是配对碱基的连线，使用get_mid求三等分点，然后取中间的一小段线段
					#配对的线如果直接删除也是可以的
		my ($x1,$y1,$x2,$y2);
		($x1) = /x1=.*?"\s*(\S+?)\s*"/;
		($y1) = /y1=.*?"\s*(\S+?)\s*"/;
		($x2) = /x2=.*?"\s*(\S+?)\s*"/;
		($y2) = /y2=.*?"\s*(\S+?)\s*"/;
		my ($new_x1,$new_x2,$new_y1,$new_y2) = get_mid($x1,$y1,$x2,$y2);
		s/$x1/$new_x1/;
		s/$y1/$new_x2/;
		s/$x2/$new_y1/;
		s/$y2/$new_y2/;
		print"$_\n";
	}else{
		if(/polyline/){	#这里是相邻碱基的连线
			$flag = 1;
			s/black/#FFFACD/;	#将连线给成背景色
			s/1.5/20/g;	#加粗
			print "$_\n";
			while($flag){
				my $tmp = <IN>;	#其他不用改动的地方直接输出
				print "$tmp";
				$flag = 0 if($tmp =~ />/);
			}	
			print'      <g style="font-family: SansSerif" id="name">',"\n",'      <text x="1" y="-4">',$name,"</text>\n      </g>\n";	#此处在左上角增加了tRNA的名称
		}else{
			print"$_\n";	##其他不用改动的地方直接输出
		}
	}
}

sub get_mid{
	my ($x1,$y1,$x2,$y2) = @_;
	my $mid_x = 0.5 * ($x1 + $x2);	#中点
	my $mid_y = 0.5 * ($y1 + $y2);

	my $new_x11 = 0.5 * ($x1 + $mid_x);	#二次中点
	my $new_x21 = 0.5 * ($x2 + $mid_x);	
	my $new_y11 = 0.5 * ($y1 + $mid_y);
	my $new_y21 = 0.5 * ($y2 + $mid_y);

	my $new_x1 = 0.5 * ($new_x11 + $mid_x);	#三次求中点
	my $new_x2 = 0.5 * ($new_x21 + $mid_x);
	my $new_y1 = 0.5 * ($new_y11 + $mid_y);
	my $new_y2 = 0.5 * ($new_y21 + $mid_y);

	return($new_x1,$new_y1,$new_x2,$new_y2);
}

sub usage{
	print "\nusage:\n\tperl $0 infile > outfile\n\t\tinfile is mitos2 trna secondary structure\n\n";
	exit;
}

单张的图片修改很容易，对于多张图片，我们可以使用一个for循环来实现此操作。

for i in *.svg ;do perl test.pl  $i > ${i%%-*}-polish-${i#*-};done

这样我们可以快速得到美化后的tRNA二级结构，但是图片是分开的，我们一般使用的都是组合后的图片，可以使用ps等工具将其组合。这里我使用的是convert命令组合。

a=(`ls *.svg`)	#修改好的svg图片，需要自己看下多少个
convert +append  ${a[@]:0:6} tmp1.png	#每行六个横向拼接，0表示从第一张图片开始，6表示一行放六张
convert +append  ${a[@]:6:6} tmp2.png	#从第七张图片开始，放六张
convert +append  ${a[@]:12:6} tmp3.png
...	#注意最后一组可能剩的少于6张
convert -append tmp1.png tmp2.png tmp3.png ... final.png	#将横向拼接好的纵向拼接成完整的

这样，就获得我们最终的图了：

好的，本期分享结束，对细胞器基因组研究感兴趣的可加群：636121790

Authors: Lei xu
contact : xul@genepioneer.cn