The ref2seq function should give a valid index that can be used to find the corresponding base in a BAM.sequence.
Current Behavior
ref2seq returns an invalid index position that is greater than the total length of the sample sequence.
Possible Solution / Implementation
My hunch is that ref2seq isn't smart enough to deal with the 'Hard Clip' (H) operation when counting anchors, as I've only encountered this bug on BAM alignments that start with H. I've looked at the source, but can't determine without attaching a debugger if that is indeed the case.
using BioAlignments
using FASTX
using XAM
# Get the reference sequence
reference_reader = open(FASTA.Reader, "reference.fasta")
reference_record = read(reference_reader)
reference_sequence = sequence(reference_record)
close(reference_reader)
# Read the problematic alignment
bam_reader = open(BAM.Reader, "sample.bam")
bam_record = read(bam_reader)
aligned_sequence = AlignedSequence(reference_sequence, BAM.alignment(bam_record))
# Let's look at postion 50 in the reference
reference_sequence[50]
# Is reference position 50 within the BAM sequence?
(50 in BAM.position(bam_record):BAM.rightposition(bam_record)) && println("it is")
# Find where reference position 50 is relative to the sequence
pos = ref2seq(aligned_sequence, 50)[1]
# Look at position 50 within the sequence
BAM.sequence(bam_record)[pos]
close(bam_reader)
Output
DNA_G
it is
ERROR: BoundsError: attempt to access BioSequences.LongDNASeq at index [284]
Stacktrace:
[1] checkbounds
@ ~/.julia/packages/BioSequences/k4j4J/src/biosequence/indexing.jl:117 [inlined]
[2] getindex(seq::BioSequences.LongDNASeq, i::Int64)
@ BioSequences ~/.julia/packages/BioSequences/k4j4J/src/biosequence/indexing.jl:142
[3] top-level scope
@ REPL[37]:1
Context
I'm trying to identify and compare variants at specific locations along the reference genome. Basically, I'd like to call a function with the reference genome position and the BAM alignment, and get the base that was called at that position. The ultimate goal is using this data in variant calling-type applications.
Your Environment
Package Version used: 2.0.0
Julia Version used: 1.6.1
Operating System and version (desktop or mobile): Fedora 34 (x86_64)
Expected Behavior
The
ref2seq
function should give a valid index that can be used to find the corresponding base in aBAM.sequence
.Current Behavior
ref2seq
returns an invalid index position that is greater than the total length of the sample sequence.Possible Solution / Implementation
My hunch is that
ref2seq
isn't smart enough to deal with the 'Hard Clip' (H) operation when counting anchors, as I've only encountered this bug on BAM alignments that start with H. I've looked at the source, but can't determine without attaching a debugger if that is indeed the case.Steps to Reproduce (for bugs)
Output
Context
I'm trying to identify and compare variants at specific locations along the reference genome. Basically, I'd like to call a function with the reference genome position and the BAM alignment, and get the base that was called at that position. The ultimate goal is using this data in variant calling-type applications.
Your Environment