Hey Monks! I'm attempting to use BioPerl to process/filter a Bam alignment file from Bowtie, iterating through the file alignment by alignment, and deciding whether or not to write a given alignment based on the identity or characteristics of the target. I've traditionally used the high-level API to access the elements of an alignment and "do things" with them, generally as follows:
use Bio::DB::Sam;
$fh = Bio::DB::Sam->new(-fasta=> "fasta file", -bam =>"bam input file
+");
my $header = $fh->header();
my $bam_out = Bio::DB::Bam->open("bam out file", "w+");
my @alignments = $fh->features();
for my $align (@alignments) {
my $query = $align->qname;
my $target = $align->seq_id;
my $is_unmapped = $align->unmapped;
unless ($in_unmapped) {
$bam_out->write1($align);
}
}
However, this usage results in an error: "Bio::DB::Bam::write1: align is not of type Bio::DB::Bam::Alignment." Further investigation and reading RE: BioPerl suggests that the high-level API returns object in the AlignWrapper class, which are apparently not writeable with my current approach. Is there a way to extract the Alignment object from the AlignWrapper or a way to extract it using the BioPerl methods? Ultimately, I could try to work with the lower-level API, but I've been pretty happy with my interactions with the high-level API to date.
Any help/direction is greatly appreciated!! ~dark_matter