アメリエフの技術ブログ

Amelieff Staff Blog

Bio-SamTools

Bio-SamToolsを使ってBAMから情報を抜き出す方法をご紹介します。

例えば、あるBAMファイルのchr21:19,660,000-19,660,600の領域に
マップされたリードの「ID」と「アライメント位置」と
「ペアワイズアライメント結果」を表示したい場合には、
次のようなPerlスクリプトを実行します。

#!/usr/bin/perl
use strict;
use warnings;
use Bio::DB::Sam;

my $bam = BAMファイル;
my $obj = Bio::DB::Sam->new(-bam => $bam);

my @ali = $obj->get_features_by_location(
__-seq_id => 'chr21', -start => 19660000, -end => 19660600);

foreach my $a(@ali){
__my $id = $a->query->name;
__my ($chr, $sta, $end) = ($a->seq_id, $a->start - 1, $a->end);
__my ($ref, $mat, $que) = $a->padded_alignment;

__print "ID: $id¥n";
__print "Position: $chr:$sta-$end¥n";
__print join("¥n", $ref, $mat, $que), "¥n¥n";
}


このような結果が返ってきます。



アライメントの詳細を見たい時に便利ですね。
この他にも、CIGARやクオリティスコアや「ペアリードがマップ
できているか」など、さまざまな情報が取得できます。