Hello,

I'm using readGAlignmentsList from package "GenomicAlignments" to read a paired-end sequencing BAM file, which contains both paired mates and singleton reads (R1 or R2 unmapped).

According to the manual, "readGAlignmentsList pairs records into mates according to the pairing criteria described below. The 1st mate will always be 1st in the GAlignmentsList list elements that have mate_status set to "mated", and the 2nd mate will always be 2nd." This is easy to understand which one is mate1 or mate2. However, when we only have one mate (i.e. a singleton), readGAlignmentsList shows only one element for that read pair. In this case, I want to know how to tell which mate it is.

The only way I can think of is to use flag = ScanBamFlag(isFirstMateRead = TRUE) or flag = ScanBamFlag(isSecondMateRead = TRUE)in ScanBamParam. But this requires me to read BAM file twice, and moreover, I'll lose track of paired reads.

Does anybody have an idea? Thanks.

Maybe like this?

> param <- ScanBamParam(what = "flag", which = GRanges("chr1:1-1e7"))
> bf <- BamFile("aligned_nodups/A1-4Aligned.sortedByCoord.out.bam", asMates = TRUE)
> param <- ScanBamParam(what = "flag", which = GRanges("chr1:1-1e7"))
> z <- readGAlignmentsList(bf, param = param)
> z
GAlignmentsList object of length 58654:
[[1]] 
GAlignments object with 2 alignments and 1 metadata column:
      seqnames strand    cigar qwidth   start     end width njunc | flag
  [1]     chr1      -     101M    101 3118856 3118956   101     0 |   83
  [2]     chr1      + 53M4D48M    101 3118781 3118885   105     0 |  163

[[2]] 
GAlignments object with 2 alignments and 1 metadata column:
      seqnames strand cigar qwidth   start     end width njunc | flag
  [1]     chr1      +  101M    101 9509913 9510013   101     0 |   99
  [2]     chr1      -  101M    101 9510028 9510128   101     0 |  147

[[3]] 
GAlignments object with 2 alignments and 1 metadata column:
      seqnames strand cigar qwidth   start     end width njunc | flag
  [1]     chr1      -  101M    101 6287680 6287780   101     0 |   83
  [2]     chr1      +  101M    101 6287603 6287703   101     0 |  163

...
<58651 more elements>
-------
seqinfo: 66 sequences from an unspecified genome
> probablyStupidFun(z)
GAlignmentsList object of length 58654:
[[1]] 
GAlignments object with 2 alignments and 2 metadata columns:
      seqnames strand    cigar qwidth   start     end width njunc | flag
  [1]     chr1      -     101M    101 3118856 3118956   101     0 |   83
  [2]     chr1      + 53M4D48M    101 3118781 3118885   105     0 |  163
      whichmate
  [1]     first
  [2]    second

[[2]] 
GAlignments object with 2 alignments and 2 metadata columns:
      seqnames strand cigar qwidth   start     end width njunc | flag whichmate
  [1]     chr1      +  101M    101 9509913 9510013   101     0 |   99     first
  [2]     chr1      -  101M    101 9510028 9510128   101     0 |  147    second

[[3]] 
GAlignments object with 2 alignments and 2 metadata columns:
      seqnames strand cigar qwidth   start     end width njunc | flag whichmate
  [1]     chr1      -  101M    101 6287680 6287780   101     0 |   83     first
  [2]     chr1      +  101M    101 6287603 6287703   101     0 |  163    second

...
<58651 more elements>
-------
seqinfo: 66 sequences from an unspecified genome
> 


flagFun <- function(numbers){
    out <- matrix(0, length(numbers), 12)
    for(i in 1:12){
        out[,i] <- numbers %% 2
        numbers <- numbers %/%2
    }
    out
}

probablyStupidFun <- function(galst){
    tmp <- unlist(galst)
    flags <- flagFun(mcols(tmp)$flag)
    mcols(tmp)$whichmate <- NA
    mcols(tmp)$whichmate[as.logical(flags[,7])] <- "first"
    mcols(tmp)$whichmate[as.logical(flags[,8])] <- "second"
    relist(tmp, galst)
}