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
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.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
- 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.