読者です 読者をやめる 読者になる 読者になる

NGSデータ解析まとめ

サカナ研究者の手探りNGS解析(おもに進化生物学)

Perlスクリプトの記法を試してみる fasta_to_nex.pl

FASTA形式(.fa|.aln|.fasta)の配列ファイル(アラインメント済)をNEXUS形式 (.nex)に変換するスクリプトです。あとでMrBayesで解析を行うため、outfileにはmrbayes block(パラメータは固定)が付きます。引数はワイルドカードに対応していません。繰り返し処理を行う場合は、shell scriptを使用する必要があります。Infile, outfileの拡張子は任意(何でもいい)です。

#!/usr/bin/perl
# fasta_to_nex.pl
# 2015-03-20 Y Hashiguchi

# fasta format をnexus format (for MrBayes)に変更する
# interleaveしない
# sequenceの開始位置をそろえる
# ntax, ncharをカウント

# modelその他MrBayes blockを書き加える
# comment outはあり

## usage ##
# ./fasta_to_nex.pl [infile.fa|.aln|.fasta|.fas] > [outfile.nex]

#

use strict;
use warnings;


### Define variables ##

my $fname = $ARGV[0];

my @lines = ();		# infileの各行を入れる配列
my @seq_with_name = ();	#名前付きsequenceを入れる配列

my $ntax = 0;		# taxa数
my $nchar = 0;		# character数(gapやNも含むsequence長さ)
my $seq_data = 0;	# sequence with name

my %sequence;	#sequence nameおよびsequenceを入れるhash

##

### MAIN ###
# input file (FASTA) の読み込み @linesに入れる
open (FASTA, "$fname") || die "Cannot open file:$!";
while(<FASTA>) {
	chomp;
	push(@lines, $_);
}
close(FASTA);

# taxa数をcount (= ">"の数を数える)
foreach my $line (@lines) {
	$ntax++ if($line =~ /\>/);
}

$sequence{"Name"} = shift(@lines);	#最初の1配列の名前(1行目)は$sequence{"Name"}に入れる
$sequence{"Name"} =~ s/\>//;	# ">"を除く

while ($#lines > -1) {
	my $data = shift(@lines);
	if($data =~ /\>/) {
		$nchar = length($sequence{"Sequence"});	# test
		while(length($sequence{"Name"}) < 10) {		#長さが10文字まで // ここでは、nameの長さが10文字より少ないことを仮定している
			$sequence{"Name"} .= " ";	# 後ろに" "を加える
		}
#		print $ntax, "\t", $nchar, "\n";		# test
		$seq_data = $sequence{"Name"} . $sequence{"Sequence"};	#nameとsequence統合
		push(@seq_with_name, $seq_data);	#配列 @seq_with_nameに結果入れる
#		print $sequence{"Name"}, $sequence{"Sequence"}, "\n";		# test
		$sequence{"Sequence"} = "";		# sequence中身を削除
		$sequence{"Name"} = $data;	#sequence_name代入
		$sequence{"Name"} =~ s/\>//;	# ">"を除く
	} else {
		$sequence{"Sequence"} .= $data;
	}
}

# 最後の1配列について
while(length($sequence{"Name"}) < 10) {		#長さが10文字まで 
	$sequence{"Name"} .= " ";	# 後ろに" "を加える
}
$seq_data = $sequence{"Name"} . $sequence{"Sequence"};	#nameとsequence統合
push(@seq_with_name, $seq_data);	#配列 @seq_with_nameに結果入れる

#print $sequence{"Name"}, $sequence{"Sequence"}, "\n";	#最後のsequenceをprint

#print $ntax, "\t", $nchar, "\n";		# test
#foreach(@seq_with_name) {		# test
#	print $_, "\n";
#}


## print out as ".nex" format
print '#NEXUS', "\n\n";
print '[', "\n", 'Data from:', "\n", $ARGV[0], "\n", ']', "\n\n\n";
print 'begin data;', "\n    ", 'dimensions ntax=', $ntax,  ' nchar=', $nchar, ';', "\n";
print '    format datatype=dna interleave=no gap=-;', "\n";
print '    matrix', "\n";
foreach(@seq_with_name) {		# test
	print $_, "\n";
}
print '    ;', "\n";
print 'end;', "\n";


### add MrBayes block ###
print 'begin mrbayes;', "\n";
print '    lset nst=6 rates=invgamma;', "\n";
print '    set quitonerror=no;', "\n";
print '    mcmc ngen=200000 samplefreq=100 printfreq=100 diagnfreq=1000;', "\n";
print '    sumt burnin=50000;', "\n";
print '    quit;', "\n";