[Bioperl-l] Are Bio::DB::Bam::Alignment objects read-only ?

[Fields, Christopher J]

According to the docs the low-level API for Bio-Samtools, both read and write are allowed:

http://search.cpan.org/perldoc?Bio::DB::Sam#The_low-level_API

Using the low-level API for this purpose isn't documented as well, though (the high-level API is read only AFAICT).  

The error message is a standard one generated from the XS bindings where the passed argument passed isn't mapped correctly.  Looking through the Sam.xs file, qseq() is only prototyped as a reader; the only arg is a Bio::DB::Bam::Alignment (e.g. $self).  However, it appears there is a function specified for Bio::DB::Bam::Alignment names l_qseq() that might be the setter, wheras qseq() is maybe to be the getter (ignore the 'bama_' prefix):

....

int
bama_l_qseq(b,...)
    Bio::DB::Bam::Alignment b
PROTOTYPE: $;$
CODE:
    if (items > 1)
      b->core.l_qseq = SvIV(ST(1));
    RETVAL=b->core.l_qseq;
OUTPUT:
    RETVAL

SV*
bama_qseq(b)
Bio::DB::Bam::Alignment b
PROTOTYPE: $
PREINIT:
    char* seq;
    int   i;
CODE:
    seq = Newxz(seq,b->core.l_qseq+1,char);
    for (i=0;i<b->core.l_qseq;i++) {
      seq[i]=bam_nt16_rev_table[bam1_seqi(bam1_seq(b),i)];
    }
    RETVAL = newSVpv(seq,b->core.l_qseq);
    Safefree(seq);
OUTPUT:
    RETVAL


-chris

On Nov 23, 2011, at 10:02 AM, Cook, Malcolm wrote:

> Charles,
> 
> I suggest you reconsider your approach to rather, use `samtools view` to pipe your reads to stdout in sam format, then stream edit the barcode and pipe it back to samtools for conversion back to .bam file.
> 
> I know this is not what you're asking.  I'm pretty sure that direct answer to your question is, "yes - they are read-only".
> 
> ~Malcolm
> 
> 
>> -----Original Message-----
>> From: bioperl-l-bounces at lists.open-bio.org [mailto:bioperl-l-
>> bounces at lists.open-bio.org] On Behalf Of Charles Plessy
>> Sent: Wednesday, November 23, 2011 4:28 AM
>> To: bioperl-l at bioperl.org
>> Subject: [Bioperl-l] Are Bio::DB::Bam::Alignment objects read-only ?
>> 
>> Dear BioPerl developers,
>> 
>> I am trying to process some unaligned paired-end reads with Bio::DB::Sam.
>> For
>> each pair, I want to detect a sequence index and a unique molecular
>> identifier in
>> the linker, record them as auxiliary flags, and trim the linker from the read.
>> 
>> I collect the pairs through a features iterator, and can access all their data
>> through the high-level Bio::DB::Bam::Alignment API.  After modifying them
>> (linker trimming and adding flags), I want to write the resulting pairs as a
>> new unaligned BAM file.
>> 
>> I apologise if the solution is trivial, but my problem is that I do not manage to
>> modify the Bio::DB::Bam::Alignment objects.  Typically, attempts such as
>> “$pair[0]->qseq("GATACA")” give errors like
>> “Usage: Bio::DB::Bam::Alignment::qseq(b) at
>> /usr/lib/perl5/Bio/DB/Bam/AlignWrapper.pm line 80”.
>> 
>> Since I did not find explanations or portsions of source code indicating how to
>> modify Bio::DB::Bam::Alignment objects, I wonder if they are read-only…
>> 
>> Have a nice day,
>> 
>> --
>> Charles Plessy
>> Tsurumi, Kanagawa, Japan
>> _______________________________________________
>> Bioperl-l mailing list
>> Bioperl-l at lists.open-bio.org
>> http://lists.open-bio.org/mailman/listinfo/bioperl-l
> 
> _______________________________________________
> Bioperl-l mailing list
> Bioperl-l at lists.open-bio.org
> http://lists.open-bio.org/mailman/listinfo/bioperl-l