How to safely and efficiently convert subset of bam to fastq?

Question

Question

How can I extract reads from a bam file (produced by bwa-mem) to fastq given a list of reference sequences to filter out?

Potential difficulties

• maintaining FR orientation of pair end reads (in bam all the sequences are reference sequences)
• keeping R1 and R2 reads
• keeping quality scores in the same encoding as original fastq (default illumina phred scores in my case)
• bam can be (ana usually is) sorted by coordinates

Almost there solutions

1. Sujai's perl solution in blobology does exact opposite - getting reads from list of references (so I could just reverse the list). The disadvantage is that the script outputs an interleaved fq file; requires unique names of mates, otherwise R1/R2 information is lost.

2. samtools + grep them all from fastq files

create a list of read names that do not map to filtered scaffolds. (cut will extract just read names, uniq will collapse pair end read names if they are the same). Then grep read names from fastq files and remove -- separator between hits

samtools view foo.bam | grep -vf list_of_scaffols_filter \
  | cut -f 1 | uniq > list_of_reads_to_keep

grep -A 3 -f list_of_reads_to_keep foo_R1.fq | grep -v "^--$" > foo_R1_filtered_bash.fq
grep -A 3 -f list_of_reads_to_keep foo_R2.fq | grep -v "^--$" > foo_R1_filtered_bash.fq

3. filter bam & picard tools

Or I could do just the filtering part and use Picard-tools (Picard.SamToFastq), but as usual I am avoiding java as much as I can. I guess

samtools view foo.bam | grep -vf list_of_scaffols_filter \
  | java -jar picard.jar SamToFastq INPUT=/dev/stdin \
  FASTQ=foo_R1_filtered_bash.fq SECOND_END_FASTQ=foo_R2_filtered_bash.fq

---

The first solution does not really work for me since I do not want to rename all the reads in bam file and I want to keep the R1/R2 information (since R1 and R2 have different error profiles). Both solutions 2 and 3 I find bit clumsy and I am not sure if they are general, I might get some unexpected behaviour if one of reads is mapped second is not... They both relay on the same filtering step.

I was wondering about some pysam solution. I guess it will be much slower, but at least it will be much clearer and perhaps more general. Something like in Convert Bam File To Fasta File - there is pysam solution for fasta (not fastq), almost there...

Story

I have very fragmented reference genome. Some of scaffolds them are too small to works with, and some of them are contaminants (identified using blobtools). I want to separate reads that are mapping to different groups to separate contaminants, short scaffolds and scaffold that will be used for downstream analysis. The reason is that if we remap all the reads to filtered reference (0.7 - 0.8 of original genome), the most of them (0.95 - 0.99) will still find a place where they map, therefore there is 0.2 - 0.3 of misplaced reads that will obviously have to affect downstream analysis, like variant calling.

This filtering idea is based logic that if the filtered duplicated genomic region will contain some small differences, they will attract their reads (and if I filter them I will improve variant calling) and if they will be exactly same, they will get reads assigned at random, so there is no harm in doing that.

Answer 1

I'm not aware of any pre-made program to do this, so I wrote one for you. This will take a BAM file with any ordering and produce properly ordered gzipped fastq files with the filtering as you requested. Internally, this iterates over all of the entries in the BAM file (ignoring secondary/supplemental entries and those where both mates map to your filter list), store the properly oriented sequence/quality/read name in a buffer, and then dumps that buffer entry to disk once the mate is found. This should be reasonably performant (hey, it's python, to don't expect too much), though if you happen to have indexed BAM files then one could think of ways to make this run faster.

Do check the output, since I've only run one test.

Answer 2

Ok, I wrote a bit brute pysam / BioPython parser that uses index of bam to get proper order of read pair for R1 / R2 files and bitwise flag. It should not be too difficult to add more sophisticated filtering rules now.

#!/usr/bin/env python3
# 1. arg - indexed bam file
# 2. arg - list of headers to filter
# 3. arg - name pattern for the output reads

import os
import sys
import pysam
from Bio import SeqIO, Seq, SeqRecord

samfile = pysam.AlignmentFile(sys.argv[1], "rb")
header_set = set(line.strip() for line in open(sys.argv[2]))
base = sys.argv[3]

out_R1 = base + 'R1_filtered.fq'
out_R2 = base + 'R2_filtered.fq'

with open(out_R1, mode='w') as R1, open(out_R2, mode='w') as R2:
    for entry in samfile:
        if entry.is_read1 and not entry.reference_name in header_set:
            # pait pair
            entry_R2 = samfile.mate(entry)
            # do some sequence gymnastics with R1
            seq_R1 = Seq.Seq(entry.seq)
            if entry.is_reverse :
                seq_R1 = seq_R1.reverse_complement()
            # do some sequence gymnastics with R2
            seq_R2 = Seq.Seq(entry_R2.seq)
            if entry_R2.is_reverse :
                seq_R2 = seq_R2.reverse_complement()
            R1.write('@' + entry.qname + '\n' + str(seq_R1) + '\n+\n' + entry.qqual + '\n')
            R2.write('@' + entry_R2.qname + '\n' + str(seq_R2) + '\n+\n' + entry_R2.qqual + '\n')


samfile.close()

There are some ugly bits, please feel free to make it nicer.