Finding a map between observed bases and positions in the reference genome

[scikal]

Hi,

Thank you for building and maintaining this great package --- I wouldn't start analyzing genomic data via Julia without it!

I have a steps that I used to perform using pysam, but could not figure how they can be done via XAM:

1. For each record/read, get a list of start and end positions of aligned gapless blocks.
2. For each record/read, get a list of reference positions that this read aligns to. The length of the list should be the same length as the read. This should allow mapping each observed base to a position in the reference genome, in the presence of insertions and deletions.

Would you be able to guide me how this can be done?

Sincerely,
Daniel

[kescobo]

As a first pass, I'd look here and set up the for record in reader... bit. As for what happens inside the loop, the easiest (though not necessarily most performant) thing would probably be to just make a couple of vectors and push! to them. Eg

startstops = []
refpos = []
for record  in reader
    # whatever to get start / end positions...
    push!(startstops, (starts, stops))
    push!(refpos, (position(record), rightposition(record)))
end

[scikal]

Thank you for the prompt response.

I'm trying to deal with the case of a read that contains short deletions, so the read length is smaller than the rightposition(record)-position(record).

[scikal]

I was wondering if there is a more easy solution than parsing the cigar.

[kescobo]

Ah, that I don't know. Will have to wait for someone more familiar with the format to weigh in, sorry.