Question: Should GenomicAlignments::mapToAlignments() terminate when no alignments are found for one read?
gravatar for Shian Su
12 days ago by
Shian Su20
Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia
Shian Su20 wrote:

When using mapToAlignments() with multiple GAlignments, it appears to terminate if none of the query GRanges map to a GAlignment.

alignments <- GAlignments(
    c("chr1", "chr1"),
    c(20L, 10L),
    c("11M", "11M"),
    strand(c("*", "*")),
    names=c("read_A", "read_B")
x <- GRanges("chr1", IRanges(c(12, 12), width=c(6, 20)))
mapToAlignments(x, alignments) ## Empty output

However if we reverse the order

mapToAlignments(x, alignments[c(2, 1)])
# GAlignments object with 1 alignment and 0 metadata columns:
#         seqnames strand       cigar    qwidth     start       end     width     njunc
#             <Rle>  <Rle> <character> <integer> <integer> <integer> <integer> <integer>
# read_B     chr1      *         11M        11        10        20        11         0
# -------
# seqinfo: 1 sequence from an unspecified genome; no seqlengths

The behaviour appears to be that if a GAlignment cannot be mapped to by any of the GRanges, then function terminates and returns the results obtained so far.

Is this the intended behaviour of the function? It makes it very difficult to use as a batch request, because there's no efficient way to figure out ahead of time which reads will cause the function to prematurely terminate. The fact it terminates silently is also dangerous, because if the offending read occurred half way down a large query then it would be very difficult for a user to realise what happened.

EDIT: I found another strange termination issue shown here: Both reads 2 and 3 have can be mapped to from the query, but no results are returned for them unless read 1 is excluded, not sure what is the termination criterion here.

genomicalignments • 78 views
ADD COMMENTlink modified 7 days ago by Valerie Obenchain6.7k • written 12 days ago by Shian Su20

Hi Shian, Thanks for the report and clear example. I'll look into it. Valerie

ADD REPLYlink written 9 days ago by Valerie Obenchain6.7k
Answer: Should GenomicAlignments::mapToAlignments() terminate when no alignments are fou
gravatar for Valerie Obenchain
7 days ago by
United States
Valerie Obenchain6.7k wrote:

This is fixed in GenomicAlignments 1.22.1 and 1.23.1.


ADD COMMENTlink written 7 days ago by Valerie Obenchain6.7k
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 451 users visited in the last hour