1

私はかなり perl に不慣れで、BioPerl に非常に慣れていないので、これが些細な質問のように思えたら、申し訳ありません。Bio::AlignIO と Bio::SimpleAlign を使用して、参照配列 (この場合はヒト mtDNA 参照配列) に対する目的の配列のペアワイズ アラインメントを生成します。アラインメントを fasta 形式で問題なく生成できますが、ウェブツール Haplogrep ( http://haplogrep.uibk.ac .at/index.html)。これをループとして実行することを想像できます。これにより、参照配列の各塩基を調べ、参照内の位置がクエリ配列と一致しない場合は、クエリ配列に存在する位置とヌクレオチドを記録します. ただし、アライメントの一部として生成された Bio::AlignIO オブジェクトからこれを行う方法がわかりません。これが私がこれまでに持っているコードです:

#!/usr/bin/perl

use strict;
use warnings;
use Bio::AlignIO;
use Bio::SimpleAlign;

my $ref = "rCRS.fas";

my $comp = shift @ARGV or die "please provide comparison sequence: $!\n";

my $comp_id;

if ($comp =~ /(\S*)\.fas*/) {
    $comp_id = $1;
}

my $CAT = "cat $ref ".$comp." > ".$comp.".tmp";
try_cmd($CAT);

my $ALN = "muscle -in ".$comp.".tmp -out ".$comp.".aln";
try_cmd($ALN);

my $str = Bio::AlignIO->new(
     -file => $comp.".aln",
    -format => "fasta",
    );

my $aln = $str->next_aln();

####Subroutine####

sub try_cmd {
    my $cmd = shift;
    my $status;
    $status = system($cmd);
    if ( $status ) {
    print STDERR ( "problem running $cmd\nReturned $status\n" );
    }
    else {
    print STDERR ( "Successfully ran $cmd\n" );
    }
}

簡単な例を挙げると、次のようなアラインメントを生成したとします。

>Ref_seq   AATTTGGGCTACT 
>Query_seq AAATTCGGCTACA 

次に、次のような出力を生成します。

3A 6C 13A

その方法を理解するのを手伝ってくれる人はいますか?これを行うためのより良い/簡単な方法に関する一般的なコメントも大歓迎です。

4

0 に答える 0