jh13-wtsi
commented
10 years ago Here is an explicit example using a BAM file of PacBio reads aligned with GMAP.
$ samtools view 'ftp://ngs.sanger.ac.uk/scratch/project/jmg/MCF7_Boley/IsoSeq_MCF7.gmapdenovo.sorted.bam' chr4:0-100000 | grep 'm130604_203349_42161_c100519042550000001823081209281390_s1_p0/159470/ccs'
m130604_203349_42161_c100519042550000001823081209281390_s1_p0/159470/ccs 16 chr4 86163 0 9M3I4M2I2M2I12M1I15M1I2M1I9M3I15M4I5M3I14M10I7M1I20M2I10M2I4M3I3M1I17M3I4M1I4M1I4M1I14M3I15M1I18M2I4M2I22M1I5M3I11M1I8M3I8M1I10M1D10M1I15M1I16M2I14M1I18M1I40M1I19M2I17M1I2M1I6M2I8M1I8M1I7M1I6M3I5M1I7M1I5M1I7M1I2M2I18M1I16M1I9M1I23M1D14M1I5M2I19M1I26M1I13M4I3M348H * 0 0 AAATGTGAAAGGGAATTGGTTGGGCAAAGCCTTTAACAAGGTCCACAACAGCTGGAATGAACACTGCCAAAAAATTCATACTACTTGGAGGTCAGAAACCTACCAGACATTGCTGGTTATAAGGCAATGTGGCAAAGCTTTAGGAATGGGTCCACAAGGGCCTTGGGAATTTAACATAAGAATATTCGCAGTACTTGGAGGAGACAACCCTACAAATATGTAAAAGAATGGCGAAAAGGCCTTTAGACAGTCCAGGGCAGCCCTGGAATGAACATAAAAATATTCATTACTGGCTGCGAAAAACCCATACACATGGGTGGAAAAATGGTGGCAAAGCTTTAACCAATGCCTCAAGTCTTATTATTACACAGGAGCATTCATTTTCTGAACAAAAACTTTTACAAATGTGAAGAATGGTGGCAAAAGCTTTACTTTGTCCTCATCCCTTAATAAACATTAAGAGAATTCATACTGGAGGGAGAAACCCTACACATGCTGAAAGAATCGTTGGCAAAGGCTTTTTATTAGGTCCTTCACACCTGCTTGCTTAAGCATAAAGAGAAATTCATGACTTTGGAGGAGAACCGTACACGGTGCGAAGAATGTGGCAAAAGCTTTTAAACCAATCGTCAACTCTTATATTCACAAGAGAATCCATTTCTGGCGGCAAAAACCTTACAAATGTTGAAGAATGGGGCAAAGCCTTTACACGGGTCCACAACACTTGGAGCA * MD:Z:54A1G23CTA2A2G3A14C2TA2T17G15A23G9TG1C25T45T17^T81GC8G65G36A22AGA4C38C14^A20A26T30A1 NH:i:9 HI:i:1 NM:i:134 SM:i:40 XQ:i:0 X2:i:0 XO:Z:UT XT:Z:NN-NN,0.00,0.00,347..349 PG:Z:M
m130604_203349_42161_c100519042550000001823081209281390_s1_p0/159470/ccs 0 chr4 86644 0 73S31M2D13M2I20M1I11M1D16M1D14M161636N17M1I6M5I17060N61M3I13M2I4M2I9M2I13M2I8M2I17M732H * 0 0 AATCAGATAGACTTTTATTTATAAGGCTCTTTGCCAACACACACCATTGCTTAACTACAATATAATAACAAATTTCATACTGGAGAGAAACCCTACACGTGCGAGATGTGGCAAAGCCTTTTAAAACAATCCTCAACTCTTTATATTACACAGAGAATCCATTCTGGCAAAAACCTTACAAAATTGTGAAAGAAGTGGCCAAAGCCCTGATGCTGTAGGGTTTCTCCTCAGTATGAATATTCTTATGTTCATCAGGCTTGGCGGACCCCATATGAAAGCTTTTGCCACCCATTTCTTCCACATTTATGTAGGGTTTCTCCGCCCAGTATCGGAATTTTCTTATGTTCG * MD:Z:31^AG0A15T2C24^A16^A16A6GA1T10ATG13TC11C9T2T1A1GG1T5TTGTC5GC11G1T13C5T12C10A NH:i:9 HI:i:1 NM:i:58 SM:i:40 XQ:i:0 X2:i:0 XO:Z:UT XS:A:- XT:Z:NN-NN,0.00,0.00,347..349 PG:Z:MNote the trailing "732H" in the CIGAR of the second read.
Now let's use perl to examine the same reads.
#!/usr/bin/env perl
use strict;
use warnings;
my $bam = 'ftp://ngs.sanger.ac.uk/scratch/project/jmg/MCF7_Boley/IsoSeq_MCF7.gmapdenovo.sorted.bam';
my $qname = 'm130604_203349_42161_c100519042550000001823081209281390_s1_p0/159470/ccs';
use Bio::DB::Sam;
printf "%4d -> %4d\n", $_->query->start, $_->query->end
for ( grep { $_->qname eq $qname; }
Bio::DB::Sam
->new( '-bam' => $bam )
->get_features_by_location(
'-seq_id' => 'chr4',
'-start' => 0,
'-end' => 100000,
));==>
1 -> 384
74 -> -384The second read has a negative query end due to the large hard clip at the end of the CIGAR!
lstein
keiranmraine
It is possible to create an alignment from valid BAM data such that $align->query->end < 0 . The problem is that $query->end calls the samtools function bam_cigar2qlen() to get the total length and subtracts all clipping (hard and soft) at the end of the CIGAR. This is incorrect because bam_cigar2qlen() only accounts for soft clipping, so if there is hard clipping at the end as well then the return value is too small. If there is a lot of hard clipping then the result can be negative.