<html><head></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; "><div style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; ">Hi all,<div><br></div><div>I'm a little stumped as to how to successfully count the depths of all reads at a specific locus in a sam/bam file. I know I can do this with GATK DepthOfCoverage but I wanted to do some more customized things with my script yet I haven't figured out how to get the right base. I was a bit surprised there wasn't (or that it's not well documented) a method to get the individuals read's base at a specific position while getting the $refbase is quite easy. (I'm betting such a method exists and is just not documented well)</div><div><br></div><div>At any rate, gaps in the alignment are the cause for my problems, so if anyone knows a simpler way to do call the bases correctly, or a clever algorithm to deal with this issue, it would be much appreciated. Here's what I have for code and it works except in cases where there are multiple gaps in the reference sequence, e.g. the alignment below should be T-T here not C-C but is shifted due to the second gap. </div><div><br></div><div><img id="16a2bc7a-138b-4b1a-930f-a3511792bc6f" height="56" width="607" apple-width="yes" apple-height="yes" src="cid:4D6D308D-BABF-4405-953C-1709CE992486@knome.com"></div><div><br></div><div><div>#!/progs/bin/perl</div><div>use strict;</div><div>use warnings;</div><div>use Bio::DB::Sam;</div><div>use Bio::DB::Bam::AlignWrapper;</div><div>use Pod::Usage;</div><div>use Getopt::Long;</div><div>use Bio::DB::Bam::Pileup;</div><div>use Term::ANSIColor;</div></div><div><br></div><div><div>my $sam = Bio::DB::Sam->new(-bam =>$BAM,</div><div> -fasta=> $FASTA);</div><div>getBases($chr, $pos, $pos);</div><div><br></div><div><br></div><div>sub getBases {</div><div> my $print = 1;</div><div> my ($chr, $start_query, $end_query) = @_;</div><div> my @alignments = $sam->get_features_by_location(-seq_id => $chr,</div><div> -start => $start_query,</div><div> -end => $end_query);</div><div> my $refbase;</div><div> my ($a_count, $t_count, $g_count, $c_count, $n_count, $del_count, $ins_count) = (0, 0, 0, 0, 0, 0, 0);</div><div> for my $a (@alignments) {</div><div><br></div><div> my $start = $a->start;</div><div> my $end = $a->end;</div><div><br></div><div> my $query_start = $a->query->start; </div><div> my $query_end = $a->query->end;</div><div> my $ref_dna = $a->dna; # reference sequence bases</div><div> my ($ref, $matches, $query) = $a->padded_alignment;</div><div> my $offset = 0;</div><div> if ($ref =~ /^([-]+)[ATCG]+/){</div><div> $offset = length($1);</div><div> }</div><div> #print "$offset\n";</div><div> $refbase = $sam->segment($chr,$start_query,$start_query)->dna; </div><div> printAlignment($ref, $matches, $query, $start_query, $start, $offset);</div><div> my $base = substr($query, $start_query-$start+$offset, 1);</div><div> if (!$base){</div><div> next;</div><div> }</div><div> $a_count++ if ($base eq "A");</div><div> $t_count++ if ($base eq "T");</div><div> $c_count++ if ($base eq "C");</div><div> $g_count++ if ($base eq "G");</div><div> $n_count++ if ($base eq "N");</div><div> $del_count++ if ($base eq "-");</div><div> my @scores = $a->qscore; # per-base quality scores</div><div> my $match_qual= $a->qual; # quality of the match</div><div> }</div><div> my $total_depth = $a_count + $t_count + $c_count + $g_count + $n_count + $del_count;</div><div> if ($print == 1){</div><div> # print "$start_query\tref base: $ref_base\n"; </div><div> print "$chr:$start_query($refbase)\t";</div><div> print "A:$a_count\t";</div><div> print "T:$t_count\t";</div><div> print "C:$c_count\t";</div><div> print "G:$g_count\t";</div><div> print "N:$n_count\t";</div><div> print "D:$del_count\t";</div><div> print "Total:$total_depth\n";</div><div> }</div><div> return ($a_count, $t_count, $c_count, $g_count, $n_count, $del_count, $ins_count);</div><div>}</div><div>sub printAlignment{</div><div> my ($ref, $matches, $query, $start_query, $start, $offset) = @_;</div><div> print substr($ref, 0, $start_query-$start+$offset);</div><div> print (color("red"), substr($ref, $start_query-$start+$offset, 1), color("reset"));</div><div> print substr($ref, $start_query-$start+$offset+1),"\n";</div><div> print substr($matches, 0, $start_query-$start+$offset);</div><div> print (color("red"), substr($matches, $start_query-$start+$offset, 1), color("reset"));</div><div> print substr($matches, $start_query-$start+$offset+1),"\n";</div><div> print substr($query, 0, $start_query-$start+$offset);</div><div> print (color("red"), substr($query, $start_query-$start+$offset, 1), color("reset"));</div><div> print substr($query, $start_query-$start+$offset+1),"\n";</div><div>}</div></div><div><br></div></div></body></html>