the problem is as follows:

I read in BWA mapped BAM with readGAlignmentPairs().
findOverlaps() does not produce a hit if only the soft-clipped part of the read pair is overlapping the specified GRanges.
Is there a (easy) way to make soft-clipped parts of the read count in overlap calculation?

Thank you,
Lovro

PS: if the interval in-between read pairs overlaps a certain range, does findOverlaps() produce a hit or not ? I would like it to.

Thank you both for help!

Let me explain the purpose of this:

It is not to calculate coverage, but to select reads, together with their sequences, that will be relevant for further analysis.

I have BAMs of probands and I want to genotype them for insertions of certain DNA sequences with breakpoints at certain known loci.
Only the read pairs spanning the breakpoints are of interest since they can prove or disprove the presence of insertion.
I can do this with presence of multiple reads that start clipping at the same loci or with read pairs in which one pair maps to the reference and the other to the insertion sequence.

I have the location of the breakpoints and want reads that span them with at least 35bp. Thus I have created a list of breakpoints:

lsBP <- split(BPgr,seq(1,98))    #BPgr is just a granges object with all breakpoints (width = 1) but I want the results separately for each of them, so I put them in a list.

I also made lists of loci 35 upstream/downstream of these loci.

 lsHP <- split(HPgr,seq(1,98))
 lsRP <- split(RPgr,seq(1,98))

In the last of the following lines I choose only the reads that overlap with all of those. #BAMs is a list of readGAlignmentPairs(BAM)
 
for(i in 1:length(BAMs)) {
 
 
  OverlapsBP <- lapply(lsBP, FUN = function(x){findOverlaps(BAMs[[i]],x,ignore.strand=T)})
  OverlapsHP <- lapply(lsHP, FUN = function(x){findOverlaps(BAMs[[i]],x,ignore.strand=T)})
  OverlapsRP <- lapply(lsRP, FUN = function(x){findOverlaps(BAMs[[i]],x,ignore.strand=T)})
 
  ntxsBP <- lapply(OverlapsBP,countQueryHits)
  ntxsHP <- lapply(OverlapsHP,countQueryHits)
  ntxsRP <- lapply(OverlapsRP,countQueryHits)
 
  readsINT <- mapply(function(BP,HP,RP){BAMs[[i]][which(BP!=0 & HP!=0 & RP!=0),]},ntxsBP,ntxsHP,ntxsRP)
 
...continues...

The code works perfectly except for getting me the reads that I'm actually interested in.
The parts that would map onto the insertion are clipped (insertion is not in the hg19 reference) and return HP==0.
You can now see why I also don't care if the NNNNN inbetween read pairs or the actual reads overlap with those 3 points.
To me it only matters if the whole insert somehow spans the putative insertion breakpoint with at least 35bp overhang.

It is important to know, that I then use the actual sequences from those reads, to call BLAT on them externally.
So the solution would be to use "granges(BAMs)" to get indexes and then use those indexes on BAMs to get the actual sequences, as you suggested:
  "Do this by calling granges() on your GAlignmentPairs object."

Reading your edit, I'm not sure will this work ? Does "Oops, not so fast! " concern this solution as well? Does granges() depend on rglist method?  

If it does, is there a way to setMethod in a way to do the job for me?
Maybe to use
ops <- c("M", "=", "X", "I", "D", "S")  # "S" added!
ans <- cigarRangesAlongQuerySpace(x@cigar, pos=x@start, ops=ops, reduce.ranges=TRUE) #cigarRangesAlongQuerySpace instead of cigarRangesAlongReferenceSpace

?

Thank you !
Lovro

Hi,

If we agree that it wouldn't make sense to have the interval in-between read pairs account in the overlap calculation without also having the so called "skipped regions" (i.e. N in the cigars), then you can completely ignore the fine range structure of each pair and replace it by a single range. Do this by calling granges() on your GAlignmentPairs object.

As for accounting soft-clipped regions in the overlap calculation, I agree with Michael that it's a questionable thing to do. [Edit: The following solution is incorrect. See my comment below this answer for a follow-up.] A quick and dirty and hacky solution would be for you to redefine the rglist method for GAlignments objects (which is called behind the scene by grglist() and findOverlaps()) just before calling findOverlaps() or grglist() on your GAlignmentPairs, and to restore the original method right after the call returns. Redefine the rglist method with:

setMethod("rglist", "GAlignments",
    function(x, use.names=TRUE, use.mcols=FALSE,
                order.as.in.query=FALSE, drop.D.ranges=FALSE)
    {
        if (!isTRUEorFALSE(use.names))
            stop("'use.names' must be TRUE or FALSE")
        if (!isTRUEorFALSE(use.mcols))
            stop("'use.mcols' must be TRUE or FALSE")
        if (!isTRUEorFALSE(order.as.in.query))
            stop("'order.as.in.query' must be TRUE or FALSE")
        ops <- c("M", "=", "X", "I", "D", "S")  # "S" added!
        ans <- cigarRangesAlongReferenceSpace(x@cigar, pos=x@start,
                                              ops=ops, reduce.ranges=TRUE)
        if (order.as.in.query)
            ans <- revElements(ans, strand(x) == "-")
        if (use.names)
            names(ans) <- names(x)
        if (use.mcols)
            mcols(ans) <- mcols(x)
        ans
    }
)

Untested!!

H.