mapToAlignments_issue
0
0
Entering edit mode
@davidporubsky-15691
Last seen 8 days ago
United States

Dear Bioc-team,

I'm trying to map genome coordinates back to local (contig/read) coordinates, however, I have an issue to do this properly in cases when alignment between contig and the genome (reference) is in minus (reverse) orientation. Is there a way to use function GenomicAlignments::mapToAlignments in strand aware fashion such that reverse alignments will be mapped in minus orientation?

Is there an example code on how to do this properly?

Thank you in advance,

David

GenomicAlignments • 121 views
ADD COMMENT
1
Entering edit mode

Do you have example code showing what you are trying to do and how it doesn't work the way you expect? Ideally a self-contained example so others can reproduce themselves?

ADD REPLY
0
Entering edit mode

Example alignment between query (Q) and target (T) sequences in PAF format.

paf.aln <- dplyr::tibble( q.name='Q', q.len=100, q.start=1, q.end=100, strand='*', t.name='T', t.len=1000, t.start=700, t.end=800, n.match=100, aln.len=100, mapq=60, cg='100=')

Target region to be lifted to the query

lift.gr <- as('T:701-710', 'GRanges')

PLUS

Setting GAlignments object in positive orientation (using PAF alignment)

plus.alignment <- GenomicAlignments::GAlignments(seqnames = paf.aln$t.name, pos = as.integer(paf.aln$t.start + 1), cigar = paf.aln$cg, strand=GenomicRanges::strand('+'), names = paf.aln$q.name)

Lifting target region to query coordinates

plus.gr <- GenomicAlignments::mapToAlignments(x = lift.gr, alignments = plus.alignment)

MINUS

Setting GAlignments object in negative orientation (using PAF alignment)

minus.alignment <- GenomicAlignments::GAlignments(seqnames = paf.aln$t.name, pos = as.integer(paf.aln$t.start + 1), cigar = paf.aln$cg, strand=GenomicRanges::strand('-'), names = paf.aln$q.name)

Lifting target region to query coordinates

minus.gr <- GenomicAlignments::mapToAlignments(x = lift.gr, alignments = minus.alignment)

Here both ranges lifted to query coordinates are the same position 1-10 of the query. I would expect in case of reverse alignment lifted ranges to correspond to position 91-100 of the query sequence. Am I doing something wrong here?

ADD REPLY
0
Entering edit mode

Any help here from the developers? I don't think this is a correct behavior from the 'mapToAlignments' function. I think it is important to make this mapping in strand aware fashion to obtain correct coordinates.

ADD REPLY

Login before adding your answer.

Traffic: 217 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6