How to tell a GAlignmentsList singleton is mate1 or mate2
1
1
Entering edit mode
foehn ▴ 100
@foehn-16281
Last seen 2.6 years ago
United States

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.

Rsamtools GenomicAlignments GAlignmentsList • 963 views
ADD COMMENT
2
Entering edit mode
@james-w-macdonald-5106
Last seen 1 hour ago
United States

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)
}

ADD COMMENT
0
Entering edit mode

Could use the flag utilities like bamFlagAsBitMatrix() and bamFlagTest().

ADD REPLY
0
Entering edit mode

Nice! thanks to both.

ADD REPLY

Login before adding your answer.

Traffic: 761 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