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

[Fields, Christopher J]

Ah, okay, makes sense.  I thought it was oddly named. :)

Chris

Sent from my iPad

On Nov 23, 2011, at 4:05 PM, "Lincoln Stein" <lincoln.stein at gmail.com<mailto:lincoln.stein at gmail.com>> wrote:

Unfortunately l_qseq read/writes the length of the query sequence, not the sequence itself.

Lincoln

On Thu, Nov 24, 2011 at 3:26 AM, Fields, Christopher J <cjfields at illinois.edu<mailto:cjfields at illinois.edu>> wrote:
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> [mailto:bioperl-l-<mailto:bioperl-l->
>> bounces at lists.open-bio.org<mailto: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<mailto: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<mailto: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<mailto: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<mailto:Bioperl-l at lists.open-bio.org>
http://lists.open-bio.org/mailman/listinfo/bioperl-l



--
Lincoln D. Stein
Director, Informatics and Biocomputing Platform
Ontario Institute for Cancer Research
101 College St., Suite 800
Toronto, ON, Canada M5G0A3
416 673-8514
Assistant: Renata Musa <Renata.Musa at oicr.on.ca<mailto:Renata.Musa at oicr.on.ca>>