I am working with RNA-seq from Patient Derived Xenografts (PDX), which means that I want to filter the mouse reads from the BAM files and keep only human reads. I have mapped a FastQ file (1) to the human genome and (2) to the mouse genome. To keep things as efficient as possible, the BAM files are not sorted. Instead, the reads are kept in the original FastQ order so that the human and mouse BAM files contain the same reads in the same order. Alternatively, the BAM files could be name-sorted.

The following process can be performed very quickly in Python using the pysam package:

1. Open the human and mouse BAM files for input and open one new BAM file for output.
2. Read one complete record (i.e., one pair of reads and all associated tags and flags) from each BAM file.
3. Compute a "goodness of mapping" metric from the cigar strings and tags.
4. If the human metric is better than the mouse metric, output the human record.
5. Repeat until end of file.

This process is fast because it passes through the files once only and because it requires no name-matching to match up reads between the human and mouse alignments.

Can the same process be implemented in R?

In principle, Rsamtools::filterBam is relevant, but it seems to want location-sorted BAM files and it seems as if I would have read in the two BAM files entirely and then use name-matching in order to setup the filter. That's slow.

I know how to read a BAM file one record at a time using BamFile() with yieldSize=1 followed by open() and scanBam(), but this produces a list and I can't see how to write a BAM file one record at a time. Also, scanBam seems to read only specified tags rather than importing an entire record.

The idea will be to use filterBam() to 'filter' the human BAM file, but to process the mouse BAM file in parallel. You'll want to do the I/O in chunks of many (100,000's?) reads so that you can vectorize the R calculations.

library(Rsamtools)
human = BamFile("human.bam", yieldSize = 100000)
mouse = BamFile("mouse.bam", yieldSize = 100000)
param = ScanBamParam(what = c("seq", "qual")) # minimum info for statistic

filter <- FilterRules(list(fun = function(human_reads) {
    ## create a parallel data frame from mouse
    aln <- scanBam(mouse, param = param)
    aln_df <- do.call(DataFrame, aln[[1]])
    ## calculate statistic and return logical(nrow(human_reads)) filter
    sample(c(TRUE, FALSE), nrow(human_reads), TRUE)
}))

destination <- tempfile(fileext=".bam")
open(mouse)
filterBam(human, destination, filter = filter, param = param)
close(mouse)

filterBam() copies the entire human record to destination, even though only some fields are input (for efficiency) into R; it would be more challenging if you were to want to write records from either the human or the mouse bam file based on some criterion...