2011年11月8日火曜日

PDFファイルの系統樹の順番に従ってアライメントをソートする

系統解析を行った後、結果として得られた系統樹の横に各OTUのアライメント(の一部)を表示して、系統樹の情報を加味しつつ配列同士の比較を行いたい(例えばインデルの有無など)。あるいはそれを誰かに示す図を作りたい。

あるいはとりあえずファーストトライで系統樹を推測し、次の解析を行うために配列データのリファインをやりたい。

こういうとき、扱う配列の数が多ければ多いほど、配列の並び替えや配列のリファインを行うために、いちいち配列データのファイルから目的の配列を抽出するのは面倒です。まあ配列のリファインは、一回一回配列名で検索かけてもさほど苦ではないかもしれませんが。


今回作成したスクリプトは、系統樹(ただしPhylogramあるいはCladogramに限る)で上から順に示されているOTUのオーダに従って、アライメントファイルのソートを行うものです。この処理を行えば、同じ順番でOTUが並んだ系統樹とアライメントを簡単に見比べることができ、上の作業のややこしさは軽減されると思います。

OTUの順序は、「PDF形式の系統樹」から取得します。これは、系統樹をrerootや昇順降順の並び替えによって編集しNEWICK形式のファイルを出力した場合には、「系統樹を可視化したときのOTUの順序」がそのままには反映されないからです。

PDF形式からのOTUの順序情報の取得は、PDFファイルからテキスト抽出を行うことで可能となります。これには「pdftotext」というプログラムを使います。

pdftotextはLinuxとMac、Windows全てで使用可能です。「Xpdf」が入っていれば、コマンドラインから使用することができます(Macの場合、Mac Portsでインストール可能)。Cygwinで使えるかは知りません。Cygwinでもpdftotextは使えます。

テキスト抽出は次のコマンド
$pdftotext -layout 系統樹ファイル.pdf 出力ファイル.txt
*Macではxpdf-pdftotext

配列のソートには、pdftotextの出力ファイルと、自前のアライメントファイルが必要です。*アライメントファイルはFASTA形式で、全てのOTUの名前は系統樹に示されているものと同じでなければいけません。

以下に記述するソースコードを適当なファイル(.pl)にコピペして、上記二つのファイルがあるディレクトリで次のコマンドを実行して下さい。

$perl ~~~.pl アライメントファイル名 pdftotext出力ファイル名

これで、PDFの系統樹のOTU順序(上から順)と同じ順番でソートされた配列ファイルが生成されます。

自前のデータでは、FigtreeでPDFファイルを作成し、これを元に並び替えを行った結果、きちんと配列がソートされることが確認できました。

何か要望があれば使って下さい(同じことができるプログラムが既にあったりして ^^;)。


ーーーーーー以下ソースコードーーーーーー

#!/usr/local/bin/perl
# BGM ~~ Extend Ash ~~
use strict;
use warnings;

my $t1 = $ARGV[0];
my $t2 = $ARGV[1];

my @taxa;
open(FILE,"$t1") or die "$!";
my @lines = <FILE>;
close (FILE);
foreach my $lines (@lines) {
    if ($lines =~ /\>/) {
        $lines =~ s/\>//;
        chomp $lines;
        push(@taxa, $lines);
    }
}

my @order;
open(FILE,"$t2") or die "$!";
my @pdf = <FILE>;
close (FILE);
foreach my $pdf (@pdf) {
    chomp $pdf;
    foreach my $taxa (@taxa) {
        if ($pdf =~ /$taxa/) {
            chomp $pdf;
            push(@order, $pdf);
        }
    }
}

open(FILE, "$t1") or die "$!";
my @sequence = <FILE>;
close(FILE);
my @names;
foreach my $sequence (@sequence) {
    if ($sequence =~ /\>/) {
        my $name = $sequence;
        $name =~ s/\>//;
        chomp $name;
        my $count = -1;
        foreach my $order (@order) {
            $count ++;
            if ($order =~ /$name/) {
                $name =~ s/\s+//;
                $names[$count] = $name;
            }                
        }
    }
}

my @result;
foreach my $names (@names) {
    open(FILE, "$t1") or die "$!";
    my @seqs = <FILE>;
    close(FILE);
    my $selection = 0;
    foreach my $seqs (@seqs) {
        if ($seqs =~ /$names/) {
            $selection = 1;
        } elsif ($seqs =~ /\>/) {
            $selection = 2;
        }
        if ($selection == 1) {
            push(@result,$seqs);
        }
    }
}

open(NEWFILE, "> $t1\_ordered.fas") or die "$!";
print NEWFILE @result;
close(NEWFILE);

0 件のコメント: