Skip Menu |
 

This queue is for tickets about the Bio-SamTools CPAN distribution.

Report information
The Basics
Id: 66387
Status: new
Priority: 0/
Queue: Bio-SamTools

People
Owner: Nobody in particular
Requestors: Yi-Shiou.Chen [...] accelrys.com
Cc:
AdminCc:

Bug Information
Severity: (no value)
Broken in: (no value)
Fixed in: (no value)



MIME-Version: 1.0
X-Bigfish: VPS-4(zz4015Lzz1202hzz8275bh8275dhz2dh2a8h668h66h)
X-Spam-Flag: NO
Acceptlanguage: en-US
X-Virus-Checked: Checked by ClamAV on 16.mx.develooper.com
Content-Type: multipart/alternative; boundary="_000_C7666A11F4F3084195AEDDD0173198A147B8607DEXCH1COLOaccelr_"
X-Virus-Scanned: Debian amavisd-new at bestpractical.com
X-Spam-Score: -6.649
Received: from localhost (localhost [127.0.0.1]) by hipster.bestpractical.com (Postfix) with ESMTP id C60B624173C for <cpan-bug+Bio-SamTools [...] hipster.bestpractical.com>; Fri, 4 Mar 2011 18:29:26 -0500 (EST)
Received: from hipster.bestpractical.com ([127.0.0.1]) by localhost (hipster.bestpractical.com [127.0.0.1]) (amavisd-new, port 10024) with ESMTP id gw8Y3rVvQAqC for <cpan-bug+Bio-SamTools [...] hipster.bestpractical.com>; Fri, 4 Mar 2011 18:29:24 -0500 (EST)
Received: from la.mx.develooper.com (x1.develooper.com [207.171.7.70]) by hipster.bestpractical.com (Postfix) with SMTP id 36C36241729 for <bug-Bio-SamTools [...] rt.cpan.org>; Fri, 4 Mar 2011 18:29:23 -0500 (EST)
Received: (qmail 26884 invoked by uid 103); 4 Mar 2011 23:29:23 -0000
Received: from x16.dev (10.0.100.26) by x1.dev with QMQP; 4 Mar 2011 23:29:23 -0000
Received: from db3ehsobe006.messaging.microsoft.com (HELO DB3EHSOBE006.bigfish.com) (213.199.154.144) by 16.mx.develooper.com (qpsmtpd/0.80) with ESMTP; Fri, 04 Mar 2011 15:29:20 -0800
Received: from mail18-db3-R.bigfish.com (10.3.81.253) by DB3EHSOBE006.bigfish.com (10.3.84.26) with Microsoft SMTP Server id 14.1.225.8; Fri, 4 Mar 2011 23:29:16 +0000
Received: from mail18-db3 (localhost.localdomain [127.0.0.1]) by mail18-db3-R.bigfish.com (Postfix) with ESMTP id 2A91BDB06C7 for <bug-Bio-SamTools [...] rt.cpan.org>; Fri, 4 Mar 2011 23:29:16 +0000 (UTC)
Received: from mail18-db3 (localhost.localdomain [127.0.0.1]) by mail18-db3 (MessageSwitch) id 1299281355188121_1935; Fri, 4 Mar 2011 23:29:15 +0000 (UTC)
Received: from DB3EHSMHS005.bigfish.com (unknown [10.3.81.253]) by mail18-db3.bigfish.com (Postfix) with ESMTP id 29D34D28052 for <bug-Bio-SamTools [...] rt.cpan.org>; Fri, 4 Mar 2011 23:29:15 +0000 (UTC)
Received: from EXCH1-COLO.accelrys.net (216.240.168.124) by DB3EHSMHS005.bigfish.com (10.3.87.105) with Microsoft SMTP Server (TLS) id 14.1.225.8; Fri, 4 Mar 2011 23:29:14 +0000
Received: from EXCH1-COLO.accelrys.net ([fe80::353e:a7f7:a8e4:527f]) by EXCH1-COLO.accelrys.net ([fe80::353e:a7f7:a8e4:527f%13]) with mapi; Fri, 4 Mar 2011 15:29:02 -0800
Delivered-To: cpan-bug+Bio-SamTools [...] hipster.bestpractical.com
Subject: Need to process CIGAR N in dna function
Thread-Index: Acvaw+zDgs6MykF0QlqeUo3iPR9kkQ==
X-Spam-Check-BY: 16.mx.develooper.com
Date: Fri, 4 Mar 2011 15:28:53 -0800
X-Originatororg: accelrys.com
X-Spam-Level:
X-CR-Hashedpuzzle: Ai2+ BBbX BuMr Cd1t Co6k C3zb EM1i G8Jw IzMG JVe2 Jt5x LF2i LGrD LsEK L5At MV9b;1;YgB1AGcALQBiAGkAbwAtAHMAYQBtAHQAbwBvAGwAcwBAAHIAdAAuAGMAcABhAG4ALgBvAHIAZwA=;Sosha1_v1;7;{66F01E10-718C-4917-91D1-67282AF52CAD};eQBpAC0AcwBoAGkAbwB1AC4AYwBoAGUAbgBAAGEAYwBjAGUAbAByAHkAcwAuAGMAbwBtAA==;Fri, 04 Mar 2011 23:28:53 GMT;TgBlAGUAZAAgAHQAbwAgAHAAcgBvAGMAZQBzAHMAIABDAEkARwBBAFIAIABOACAAaQBuACAAZABuAGEAIABmAHUAbgBjAHQAaQBvAG4A
To: "bug-Bio-SamTools [...] rt.cpan.org" <bug-Bio-SamTools [...] rt.cpan.org>
X-Spamscore: -4
From Yi-Shiou.Chen [...] accelrys.com Fri Mar 4 18: 29:26 2011
X-Forefront-Antispam-Report: KIP:(null);UIP:(null);IPVD:NLI;H:EXCH1-COLO.accelrys.net;RD:lwdc.ar06.fa2-63.host27.24465.americanis.net;EFVD:NLI
X-Spam-Status: No, score=-6.649 tagged_above=-99.9 required=10 tests=[AWL=0.250, BAYES_00=-1.9, HTML_MESSAGE=0.001, RCVD_IN_DNSWL_HI=-5] autolearn=ham
Content-Language: en-US
Message-ID: <C7666A11F4F3084195AEDDD0173198A147B8607D [...] EXCH1-COLO.accelrys.net>
X-MS-Tnef-Correlator:
X-CR-Puzzleid: {66F01E10-718C-4917-91D1-67282AF52CAD}
Return-Path: <Yi-Shiou.Chen [...] accelrys.com>
X-Original-To: cpan-bug+Bio-SamTools [...] hipster.bestpractical.com
X-RT-Mail-Extension: bio-samtools
Thread-Topic: Need to process CIGAR N in dna function
X-MS-Has-Attach:
Accept-Language: en-US
X-Spam-TCS-SCL: 5:0
From: Yi-Shiou Chen <Yi-Shiou.Chen [...] accelrys.com>
Content-Length: 0
content-type: text/plain; charset="utf-8"
Content-Transfer-Encoding: quoted-printable
X-RT-Original-Encoding: us-ascii
Content-Length: 3445
Download (untitled) / with headers
text/plain 3.3k
According to samtools's calmd function, the MD tag of BAM doesn't carry base information for skipping segments (i.e. N in CIGAR). For example, a read with 40M100N35N in CIGAR can has 75 in MD, without any information about the 100N. The dna function of Bio::DB::Bam::AlignWrapper doesn't consider this case and returns wrong reference sequences for reads with N in CIGAR. Propose to fix the function as follows. Thanks, Yi-Shiou Chen Accelrys sub dna { my $self = shift; my $sam = $self->{sam}; if (my $md = $self->get_tag_values('MD')) { # try to use MD string my $qseq = $self->qseq; #preprocess qseq using cigar array my $cigar = $self->cigar_array; my $seq = ''; my @Npos = (); my @Ncount = (); for my $op (@$cigar) { my ($operation,$count) = @$op; if ($operation eq 'M') { $seq .= substr($qseq,0,$count,''); # include these residues } elsif ($operation eq 'N') { # record "N" positions push @Npos, length($seq); # 0-based push @Ncount, $count; my $ns = $self->start + length($seq); # insert bases in "N" segment my $nseq; if(defined $self->{sam}) { $nseq = $self->{sam}->seq($self->seq_id, $ns, $ns + $count - 1); } else{ $nseq = 'N' x $count; } $seq .= $nseq; } elsif ($operation eq 'S' or $operation eq 'I') { substr($qseq,0,$count,''); # skip soft clipped and inserted residues } } my $start = 0; my $result; while ($md =~ /(\d+)|\^([gatcn]+)|([gatcn]+)/ig) { if (defined $1) { my $c = $1; while($c > 0) { if(@Npos == 0) { $result .= substr($seq,$start,$c); $start += $c; last; } # may cross N segments, need to copy piece by piece if(($start + $c) > $Npos[0]) { my $toAdd = $Npos[0] - $start; $result .= substr($seq, $start, $toAdd).substr($seq, $Npos[0], $Ncount[0]); $start += $toAdd + $Ncount[0]; $c -= $toAdd; shift @Npos; shift @Ncount; } else { $result .= substr($seq, $start, $c); $start += $c; last; } } } elsif ($2) { $result .= $2; } elsif ($3) { $result .= $3; $start += length $3; # insert N segments while(@Npos > 0 and $start >= $Npos[0]) { my $offset = $start - $Npos[0]; $result = substr($result, 0, length($result) - $offset).substr($seq, $Npos[0], $Ncount[0]).substr($result, length($result) - $offset); $start += $Ncount[0]; shift @Npos; shift @Ncount; } } } return $result; } else { return $self->{sam}->seq($self->seq_id,$self->start,$self->end); } }
content-type: text/html; charset="utf-8"
Content-Transfer-Encoding: quoted-printable
X-RT-Original-Encoding: us-ascii
Content-Length: 18550
Download (untitled) / with headers
text/html 18.1k

Message body is not shown because it is too large.



This service is sponsored and maintained by Best Practical Solutions and runs on Perl.org infrastructure.

Please report any issues with rt.cpan.org to rt-cpan-admin@bestpractical.com.